Search code examples
pythonscipyfftfrequency

Python Frequency filtering with seemingly wrong frequencies

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.

Solution

  • 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])