2

The script below filters frequencies by cutting of all frequencies larger than 6. However instead of using the seemingly correct function rfftfreq, fftfreq is being used.

To my understandnig rfftfreq should be used together with rfft. Why does this code work although it uses fftfreq with rfft?

import numpy as np
from scipy.fftpack import rfft, irfft, fftfreq

time   = np.linspace(0,10,2000)
signal = np.cos(5*np.pi*time) + np.cos(7*np.pi*time)

W = fftfreq(signal.size, d=time[1]-time[0])
f_signal = rfft(signal)

# If our original signal time was in seconds, this is now in Hz    
cut_f_signal = f_signal.copy()
cut_f_signal[(W<6)] = 0

cut_signal = irfft(cut_f_signal)

Background: rfft gives an array sorting the fourier modes with real and imaginery in separate entries. Such as [R(0), R(1), I(1), ... R(N/2),I(N/2)] for R(n) and I(n) being the real and imaginery part of the fourier mode respectively. (assuming an even entry array)

Therefore, rfftfreq yields an array of frequencies corresponding to this array such as (assuming an even entry array and sampling spacing of 1):

[0, 1/n, 1/n, 2/n, 2/n ... n/(2n), n/(2n)]

However this code works with fftfreq where the output of the function is 

[0, 1/n, ... n/(2n), -n/(2n), ..., -1/n]

Clearly, fftfreq should lead to wrong results when used together with rfft because the frequencies and FFT bins do not match.

MB-F
  • 22,770
  • 4
  • 61
  • 116
dann
  • 23
  • 4
  • does this work? Does it do what it is supposed to? What is the question? – wwii May 06 '18 at 16:57
  • It works, but it shouldnt imo...Therefore question: why does it work? – dann May 06 '18 at 17:01
  • If it works and you think it should not, maybe you have a misunderstanding. Unfortunately this isn't a discussion forum or tutorial service. Please take the time to read [ask] and the other links found on that page. – wwii May 06 '18 at 17:16
  • 1
    sorry, but whats wrong with a simple question starting with why!? Im not asking for a discussion neither for a tutorial. Just a straight answer why it works with fftfreq but not with rfftfreq. – dann May 06 '18 at 17:19
  • @dann nothing wrong in principle. People get frustrated from many bad questions, and an opening statement like *I have found the following script...* is due to trigger bad feelings ;) – MB-F May 07 '18 at 09:15
  • I think the question was stated clearly enough, with additional information demonstrating knowledge of the topic. I have applied some editing to make it more accessible. – MB-F May 07 '18 at 09:31
  • The question could be further improved by stating how you determine that the code works with one function and not with the other, and what results you would expect instead. – MB-F May 07 '18 at 09:32

1 Answers1

1

You mis-specified the frequencies in the original signal.

A sine wave is parameterized according to this equation (from Wikipedia): enter image description here

The factor 2 is missing in the definition of signal = np.cos(5*np.pi*time) + np.cos(7*np.pi*time). Therefore, the actual frequencies are 

5*pi*t = 2*pi*t * f
f = (5*pi*t) / (2*pi*t) = 5 / 2

7*pi*t = 2*pi*t * f
f = (7*pi*t) / (2*pi*t) = 7 / 2

In words, the two frequencies are half of what you think they are. Ironically, that is why it seems to work with fftfreq instead of rfftfreq. The former covers twice the frequency range (positive and negative frequencies), and therefore compensates for the missing factor 2.

This is the correct code:

signal = np.cos(5 * 2*np.pi * time) + np.cos(7 * 2*np.pi * time)
W = rfftfreq(signal.size, d=time[1]-time[0])
MB-F
  • 22,770
  • 4
  • 61
  • 116