CODEWITHSUNDEEP
pythonpython-3.xscipysignal-processingfft
CyberMathIdiot
CyberMathIdiot

Reputation: 303

Can't find the frequency of the signal using Fourier transform

I integrate a system of 400 differential equations using odeint (the first 200 equations are the equation of x component of 200 neurons, and the other 200 of y component). So the main body of my code that does the integration is this

t_final = 100.0
dt = 0.01
t = np.arange(0, t_final, dt)
sol = odeint(full_derivative, z0, t)
x10 = sol[:,9]
y10 = sol[:,209]

It doesn't matter which is the model that I use (I don't want to make things more complicated), but the integration part is correct. In x10 there is the signal of x-component for the 10th oscillator of my system, which looks like that

enter image description here

It is obvious that this is a periodic signal with a specific period and frequency. So I want to do a Fourier transform to find this frequency. I use this code to do the transform

from scipy import fftpack

f_s = len(t)//2

X = fftpack.fft(x10)
freqs = fftpack.fftfreq(len(x10)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
#ax.set_ylim(-5, 110)

and the result that I take is this (which is not very beautiful because it shows that the frequency is approximately zero).

enter image description here

What can I do to fix the bug in my code?

p.s. Maybe in this example it is relatively obvious which is the frequency of the system, but if I change the parameters of my problem I can end up in more complicated solutions. This is the reason why I want to do a fourier transform.

Upvotes: 0

Views: 158

Answers (1)

Bob
Bob

Reputation: 14654

The plot makes sense to me if the time unit in the first plot is second, because then you should have an important frequency component close to 0.1Hz.

I see in the first part you are using dt = 0.01 and I understand this is the sampling interval. In second you are using f_s = len(t) // 2 that should be the one 1.0/dt this will actually make the frequency you will find even smaller since now f_s will be 100 instead of 5000, but the frequency you are searching is still ~0.2% of the sampling frequency, so you have will have to zoom in to the region of interest. Other thing to pay attention is that if the signal has non-zero mean there will be a peak corresponding to frequency zero.

Upvotes: 0

Related Questions