scipy 快速傅立叶变换无法识别信号

Question

我正在尝试通过傅立叶变换获取信号的频率，但它无法识别它（将峰值设置为 f=0）。也许我的代码有问题（页面末尾的完整代码）：

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()

其中 Npoints 是间隔（点）的数量，而 deltaT 是间隔的时间间隔。可以看到峰值在 f=0

我还展示了 Y[0,:]（我的信号）随时间变化的图，很明显信号具有特征频率

完整的可重复代码

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()

Answer 1

我无法确定问题出在哪里。 fft 代码中似乎存在一些问题。反正我时间不多，就放一个我之前做的示例代码。您可以将其用作参考或复制粘贴。它应该可以工作。

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()