I have a 10 second signal of a pure 50 Hz sine wave sampled at 40.96 kSa/s (409,600 samples). I apply a Hanning window and the magnitude plot of my (numpy) FFT shows (ignoring negative image components in this post) 3 bins at 49.9, 50.0 and 50.1 Hz. Behaving as expected:
FFT magnitude plot of input (zoomed y-axis)
Now when I resample this signal (taking every 50th value) to get 8192 samples (effective sampling rate of 819.2 Sa/s), the magnitude plot of this FFT suddenly shows spectral leakage at components with a distance of approx. n*5 Hz from 50 Hz (n being an integer), see attached image.
FFT magnitude plot of resampled input (zoomed y-axis)
I observed that they decrease in magnitude, when:
I increase the resampling frequency.
Increase the signal duration or do zero-padding (of course they are just distributed to more bins).
Apply a Flattop window instead.
My questions are:
Where do they come from?
Why are they completely gone, when I use 8193 samples for my FFT? Is this approach used in practice?
Is there a way to prevent this leakage, but still have 2^n (n being an integer) samples in my FFT? This is a requirement for my C++ implementation.
Find a MVE here for reproduction:
import numpy as np
import plotly.express as px
# ___ Input ___
m = 1 # signal amplitude (-)
H = 1 # harmonic order (-)
phi = 0 # angle shift (°)
frequency_H1 = 50 # frequency of fundamental H1 (Hz)
T_max = 10 # signal duration (s)
sampling_rate = 40960 # sampling rate of raw input (Sa/s)
fft_size_res = int(8192) # fft size of resampled signal (-)
fft_y_axis_limit = 0.01 # y-axis limit of fft magnitude plot (-)
# ___ Derived Values ___
omega = 2*np.pi*frequency_H1
t = np.arange(0, T_max, 1/sampling_rate)
# 1. raw signal x(t) (10 s at 40960 Hz) _______________________________________________________________________
# 1.1 signal creation _________________________________________________________________________________________
x = m * np.cos(H * omega * t + phi)
#px.line(x=t, y=x, labels={"x": "Time (s)", "y": "Amplitude"}, title="x").show()
# 1.2 Windowing _______________________________________________________________________________________________
x_windowed = x*np.hanning(len(x))
#px.line(x=t, y=x_windowed, labels={"x": "Time (s)", "y": "Amplitude"}, title="x_windowed").show()
# 1.3 FFT _____________________________________________________________________________________________________
x_windowed_fft = np.fft.fft(x_windowed)
# 1.3.1 FFT Magnitude _________________________________________________________________________________________
x_windowed_fft_abs = np.abs(x_windowed_fft)
x_windowed_fft_abs_normalized = x_windowed_fft_abs/len(x)*2 # factor of "2" is introduced due to the hanning
# window, which decreases the signal magnitude by a factor of approx. 2
x_windowed_fft_frequencies = np.fft.fftfreq(len(x), d=1/sampling_rate)
(px.scatter(x=x_windowed_fft_frequencies, y=x_windowed_fft_abs_normalized,
labels={"x": "Frequencies (Hz)", "y": "Amplitude"}, title="x_windowed_fft_abs_normalized").
update_xaxes(range=[30, 70]).update_yaxes(range=[0, fft_y_axis_limit]).show())
# 2. resampled signal x_res(t) (10 s at 819.2 Hz) ______________________________________________________________
# 2.1 resampling / signal creation _____________________________________________________________________________
resample_indices = np.linspace(0, len(x)-1, fft_size_res, dtype=int)
x_res = x[resample_indices]
t_res = t[resample_indices]
sampling_rate_resampling = sampling_rate/(len(x)/len(x_res)) # effective sample rate after resampling
#px.line(x=t_res, y=x_res, labels={"x": "Time (s)", "y": "Amplitude"}, title="x_res").show()
# 2.2 Windowing _______________________________________________________________________________________________
x_res_windowed = x_res*np.hanning(len(x_res))
#px.line(x=t_res, y=x_res_windowed, labels={"x": "Time (s)", "y": "Amplitude"}, title="x_res_windowed").show()
# 2.3 FFT _____________________________________________________________________________________________________
x_res_windowed_fft = np.fft.fft(x_res_windowed)
# 2.3.1 FFT Magnitude _________________________________________________________________________________________
x_res_windowed_fft_abs = np.abs(x_res_windowed_fft)
x_res_windowed_fft_abs_normalized = x_res_windowed_fft_abs/len(x_res)*2 # factor of "2" is introduced due to the hanning
# window, which decreases the signal magnitude by a factor of approx. 2
x_res_windowed_fft_frequencies = np.fft.fftfreq(len(x_res), d=1/sampling_rate_resampling)
(px.scatter(x=x_res_windowed_fft_frequencies, y=x_res_windowed_fft_abs_normalized,
labels={"x": "Frequencies (Hz)", "y": "Amplitude"}, title="x_res_windowed_fft_abs_normalized").
update_xaxes(range=[30, 70]).update_yaxes(range=[0, fft_y_axis_limit]).show())