0

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()
Pancio
  • 23
  • 5
  • 1
    suggest you post a minimal reproducible example with data – Mitch Wheat Nov 24 '21 at 09:53
  • 2
    The `fft` result is complex numbers. You need to take their absolute values before plotting. Also, the code you provided is not enough. What is that `T`? – Abdur Rakib Nov 24 '21 at 09:56
  • @AbdurRakib sorry i forgot a bit of code, now it is edited. I plotted the abs value of the fft. T is the array with the time steps associated to the signal, used in the last plot – Pancio Nov 24 '21 at 10:23
  • @MitchWheat I reported at the bottom the full code so it can be runned with data – Pancio Nov 24 '21 at 10:30
  • 1
    The amplitude of the component at f=0 is the *mean* of the signal (i.e. the "DC component"). Look at the plot of your data; the mean is approximately 1.342. Now look at the size of the spike in your plot of the FFT result. The amplitudes of the sinusoidal components of your signal are all much smaller than 1.342, so the DC component dominates the plot. – Warren Weckesser Nov 24 '21 at 10:45
  • See, for example, https://stackoverflow.com/questions/65024650/fft-scipy-calculating-frequency and https://stackoverflow.com/questions/20634833/scipy-numpy-fft-on-data-from-file. Some searching would probably turn up more related questions. – Warren Weckesser Nov 24 '21 at 10:54

1 Answers1

0

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()
Abdur Rakib
  • 336
  • 1
  • 12