scipy Fast fourier transform doesn't recognize the signal

418 Views Asked by At

i'm trying to get the frequency of a signal via fourier transform but it's not able to recognize it (sets the peak to f=0). Maybe something is wrong in my code (FULL reprudible code at the end of the page):

PF = fft.fft(Y[0,:])/Npoints #/Npoints to get the true amplitudes
ZF = fft.fft(Y[1,:])/Npoints
freq = fft.fftfreq(Npoints,deltaT)
PF = fft.fftshift(PF) #change of ordering so that the frequencies are increasing
ZF = fft.fftshift(ZF)
freq = fft.fftshift(freq)
plt.plot(freq, np.abs(PF))
plt.show()
plt.plot(T,Y[0,:])
plt.show()

where Npoints is the number of intervals (points) and deltaT is the time spacing of the intervals. You can see that the peak is at f=0 Fourier transform PF over freq

I show also a plot of Y[0,:] (my signal) over time where it's clear that the signal has a characteristic frequency the signal to be processed

FULL REPRUDICIBLE CODE

import numpy as np 
import matplotlib.pyplot as plt
#numerical integration
from scipy.integrate import solve_ivp
import scipy.fft as fft

r=0.5
g=0.4
e=0.6
H=0.6
m=0.15
#define a vector of K between 0 and 4 with 50 componets
K=np.arange(0.1,4,0.4)

tsteps=np.arange(7200,10000,5)
Npoints=len(tsteps)
deltaT=2800/Npoints #sample spacing
for k in K :
    i=0
    def RmAmodel(t,y):
        return [r*y[0]*(1-y[0]/k)-g*y[0]/(y[0]+H)*y[1], e*g*y[0]/(y[1]+H)*y[1]-m*y[1]]

    sol = solve_ivp(RmAmodel, [0,10000], [3,3], t_eval=tsteps) #t_eval specify the points where the solution is desired
    T=sol.t
    Y=sol.y
    vk=[]
    for i in range(Npoints):
        vk.append(k)
    XYZ=[vk,Y[0,:],Y[1,:]]
    #check periodicity over P and Z with fourier transform

    
    
#try Fourier analysis just for the last value of K        
    PF = fft.fft(Y[0,:])/Npoints #/Npoints to get the true amplitudes
    ZF = fft.fft(Y[1,:])/Npoints
    freq = fft.fftfreq(Npoints,deltaT)
    PF = fft.fftshift(PF) #change of ordering so that the frequencies are increasing
    ZF = fft.fftshift(ZF)
    freq = fft.fftshift(freq)
    plt.plot(T,Y[0,:])
    plt.show()
    plt.plot(freq, np.abs(PF))
    plt.show()
1

There are 1 best solutions below

1
Abdur Rakib On

I can't pinpoint where the problem is. It looks like there is some problem in the fft code. Anyway, I have little time so I will just put a sample code I made before. You can use it as reference or copy-paste it. It should work.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq

fs = 1000   #sampling frequency
T = 1/fs    #sampling period
N = int((1 / T) + 1) #number of sample points for 1 second
t = np.linspace(0, 1, N) #time array
pi = np.pi

sig1 = 1 * np.sin(2*pi*10*t)
sig2 = 2 * np.sin(2*pi*30*t)
sig3 = 3 * np.sin(2*pi*50*t)
#generate signal
signal = sig1 + sig2 + sig3
#plot signal
plt.plot(t, signal)
plt.show()

signal_fft = fft(signal)    #getting fft
f2 = np.abs(signal_fft / N) #full spectrum
f1 = f2[:N//2]              #half spectrum
f1[1:] = 2*f1[1:]           #actual amplitude

freq = fs * np.linspace(0,N/2,int(N/2)) / N     #frequency array
#plot fft result
plt.plot(freq, f1)
plt.xlim(0,100)
plt.show()