1
$\begingroup$

I have been studying the following code but still unable to explain its output. I really hope someone could shed some light.

Let me begin with the baseline. Suppose I have a 20 seconds sinusoid with frequency components of [1,2,3,4....500Hz]. The signal has Fs=1024Hz, so it's FFT output looks as follows. The FFT amplitude is also scaled to preserve Parseval's Theorem. The frequency resolution is Fs/nfft = 1.0 Hz per bin in this case, as the result the FFT output smeared itself and make the output looks like a flat line, which is expected.

# Python code

import numpy as np
import matplotlib.pyplot as plt

dur = 20
fs = 1024
t1 = np.arange(0, dur, 1/fs)    
x1 = np.zeros_like(t1)
for f in np.arange(1, 500):
    x1 += np.sin(2*np.pi*f*t1)


plt.figure(figsize=(8,4))
plt.title('Signal')
plt.plot(t1[:100], x1[:100])
plt.xlabel('sec')
plt.grid()
plt.show()


nfft = 1024
df = fs / nfft
X1 = np.fft.fft(x1, nfft)[:nfft//2]

# Scale FFT amplitude to preserve Parseval's Theorem
X1 = X1 * np.sqrt(2 / nfft)
f1 = np.arange(0, nfft//2) * df

plt.figure(figsize=(8,4))
plt.title(f'Baseline FFT, Freq Resolution={df}Hz')
plt.plot(f1, np.abs(X1))
plt.grid()
plt.xlim([0, nfft/2])
plt.xlabel('Hz')
plt.show()

p1 = np.sum(x1[:nfft]**2)
p2 = np.sum(np.abs(X1)**2)
print('Parseval check:', np.allclose(p1, p2))

enter image description here

Here comes the tricky part. The signal is resampled in the following way, x2 = x1[::10] (downsampling by a factor of 10?). If that's really downsampling then new effective Fs should become 1024/10 = 102.4Hz, then the Nyquist frequency becomes $102.4/2 = 51.2Hz$, so the FFT output should only contain frequency components up to 51Hz. At least that's from the theoretical perspective. However the following output suggests otherwise. I can see the frequency resolution has improved significantly from 1.0Hz to 0.1Hz, but then the output only contain frequency components up to 12Hz instead of 51Hz. Why is that happen? What do I miss? Can anyone shed some light? Thanks

step = 10
x2 = x1[::step]

plt.figure(figsize=(8,4))
plt.title('Downsampled Signal')
plt.plot(t1[:100], x2[:100])
plt.xlabel('sec')
plt.grid()
plt.show()


nfft = 1024
fs2 = fs / step
df2 = fs2 / nfft
X2 = np.fft.fft(x2, nfft)[:nfft//2]

# To preserve Parseval's Theorem
X2 = X2 * np.sqrt(2/nfft)
f2 = np.arange(0, nfft//2) * df2

plt.figure(figsize=(8,4))
plt.title(f'Downsampled FFT, Freq Resolution={df2}Hz')
plt.plot(f2, np.abs(X2))
plt.grid()
plt.xlim([0, 50])
plt.show()

p1 = np.sum(x2[:nfft]**2)
p2 = np.sum(np.abs(X2)**2)
print('Parseval check:', np.allclose(p1, p2))

enter image description here

$\endgroup$
  • $\begingroup$ Does the Parseval theorem applies in your resampled signal? $\endgroup$ – Filipe Pinto May 28 at 15:33
  • $\begingroup$ I believe it is a kind of aliasing problem when you do the downsampling. I've tried your code and If you change the variable step to 9 or 11 then the output is more compatible with the theory. $\endgroup$ – Filipe Pinto May 28 at 15:38
  • $\begingroup$ Hi @FilipePinto, thanks for the reminder. I have rescaled the FFT amplitude to preserve Parseval Theorem. I hope it is done correctly. $\endgroup$ – Scoodood C May 28 at 21:21
5
$\begingroup$

I think there are a couple things going on here.

First, the flatness of your baseline plot is a little deceiving because the tones that you have added (1 - 500 Hz) are all almost exact multiples of the fundamental frequency of the fft (fs/nfft = 1.024 Hz/bin).

If you increased nfft by a factor of 10, for example, the finer frequency resolution of the zero-padded FFT would show something similar to your New FFT Plot. Here it is with nfft = 10240:

enter image description here

Because the tones are centered well on the FFT bins, there is no spectral leakage. If you were to change the tones so that they were for example at multiples of 0.75 Hz, you'd see something like this:

enter image description here

It looks like a mess, but what's really happening is that each tones is getting smeared across several bins and the tones are constructively and destructively interfering. Here's a single tone at 0.75 Hz with nfft=10240:

enter image description here

This spectral leakage occurs because the FFT assumes that its input is repeated forever, so if you have a sinusoid whose period is a multiple of the fundamental frequency of the FFT everything looks as you'd expect, but if not then there will be a discontinuity between the first sample and the last sample. This effect can be mitigated using window functions: https://en.wikipedia.org/wiki/Window_function.

However, the main issue that you're seeing is the result of aliasing. If you want to avoid aliasing, you need to low-pass filter your signal before downsampling. Otherwise, any frequency content in the original signal that is above the Nyquist rate of your down-sampled signal will alias. In your case, the sampling rate after down-sampling is 1024/10 = 102.4 Hz, so all of the tones you have above 51.2 Hz in your original signal will alias to frequencies below 51.2 Hz and interfere with the true signals at those frequencies.

Changing your down-sampling factor to 8 makes the effect easier to demonstrate because the aliased tones will land right on top of other true tones. If your original signal only has tones below the 64 Hz, all tones are present in the down-sampled signal:

for f in np.arange(1, 64):
...
step = 8
x2 = x1[::step]   # down-sampling by a factor of 10??
fs2 = fs/step

enter image description here

Now let's add a tone at 65 Hz:

for f in np.arange(1, 66):

There is no longer any tone at 63 Hz or 65 Hz. The 65 Hz tone has aliased and destructively interfered with the true tone at 63 Hz to cancel it out:

enter image description here

Now, let's change the phase of the original tone by 180 degrees:

for f in np.arange(1, 66):
x1 += np.sin(2*np.pi*f*t1 + np.pi/2)

enter image description here

There's still no tone at 65 Hz, but the one at 65 Hz is now constructively interfering, doubling the magnitude in the 63 Hz bin.

So if you extend this to 500 Hz, you can start to see what's going on. The fix is to to apply a low-pass filter before you down-sample. The combined steps of low-pass filtering and then down-sampling is called decimation. Here's the scipy function to do it: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.decimate.html.

$\endgroup$
  • 1
    $\begingroup$ Hi @rkor, thanks for your excellent explanation. I now understand the result of my FFT output was due to aliasing. I follow you advice and apply decimation to the signal and now all frequency components reveal themselves. I had to use signal.decimate(...ftype='fir', n=1500,...) in order to remove aliasing in this example. If the signal has unknown frequency component it would be a challenge to find the right setting for decimation. $\endgroup$ – Scoodood C May 28 at 21:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.