Why resampling change this FFT output?

Question

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

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

Answer 1

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:

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:

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:

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

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:

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)

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.