CODEWITHSUNDEEP
pythonfft
TheStrangeQuark
TheStrangeQuark

Reputation: 2405

blocked FFT of signal with frequency noise

I have a signal with a phase that is around 40 Hz, with small amounts of noise. I am trying to analyze it to find the frequency in 1 second blocks, I am using python to do a simulation of this. I found that the blocked FFT analysis was not accurate, so I looked at specific blocks to see why this was happening and found that at later times, the FFT was not making sense and having a spike around the 40 Hz where it should. I am using Python for this simulation and here is the code:

First create the signal:

# Time
time = 500
# sample spacing
T = 1.0 / 5000.0
sample_freq = 1.0/T
#number of points
N = time / T
x = np.linspace(0.0, N*T, N)
noise = np.random.normal(scale = 0.01, size = len(x))
noisy_freq = 40 + noise
y = 50.0 * np.sin(noisy_freq * 2.0*np.pi*x) + 1.0*np.sin(80.0 * 2.0*np.pi*x)

Then I looked at the FFT for the first second and receive:

first second

and zoomed in to around 40 Hz

zoom of first second

Then I looked at the 10th second block:

10th second

and zoomed to around 40 Hz:

10th second zoom

It is visible the signal is degrading already and this horizontal line is beginning to develop which I'm not sure what is is from.

Then I looked at the 100th second and this is what I found:

100th second

and zoomed to around 40 Hz:

100th second zoom

and the FFT response is showing barely any certainty of a spike around 40 Hz. I cannot figure out why the signal degrades at later times. I have tried using windowing functions, but this does not help.

Here is the code I use to create the FFT:

sample_freq = 1/T
time_step = 10
step = int(time_step * sample_freq)

x = x[100/T:101/T]
y = y[100/T:101/T]

flat = flat_top_windowing(len(y))
y = y*flat

yf = np.abs(np.fft.fft(y))
x_n = x.size

xf = np.fft.fftfreq(x_n,1/sample_freq)

plt.close()
plt.plot(xf, 2.0/N * yf[0:N/2])
plt.grid()
plt.show()

Upvotes: 2

Views: 314

Answers (1)

SleuthEye
SleuthEye

Reputation: 14577

The problem is with the generated signal:

np.sin(noisy_freq * 2.0*np.pi*x)

As the time advances (values in variable x), the multiplication noisy_freq*x means that noisy_freq variations have an ever increasing impact on the sin phase. Given increasing phase variations, the actual instantaneous frequency varies that much more (to the extent that it eventually looks like it randomly jumps all over the whole spectrum).

To generate a frequency modulated signal, you should instead integrate the frequency contributions with:

dphi = 2.0*np.pi*noisy_freq[:-1]*T; # per sample frequency contributions
dphi = np.insert(dphi, 0, 0);       # set the initial phase to 0
phi  = np.cumsum(dphi);             # integrate phase

y = 50.0 * np.sin(phi) + 1.0*np.sin(80.0 * 2.0*np.pi*x)

Upvotes: 2

Related Questions