Python：为什么使用FFT方法去除频率似乎没有消除数据中的周期性？

Question

我有各种气象数据的时间序列，当目测时，它们都显示出明显的日周期和季节周期。我试图用FFT方法从时间序列中去除这些季节和日周期。我已经通过线性插值从数据集中删除了nans。每个变量以1分钟的统一间隔测量，并存储在一个数组（y）中，该数组的长度为1012320x1（703天的数据值）

到目前为止，我有代码：

#Calculate the anomalies from the mean of time series y
y = y - y.mean()

#Calculate FFT 
fftVar = np.fft.fft(y)

#n = 1012320
n = len(y)

#Just take positive values
n2 = np.divide(n,2)
rng = range(1,int(n2))

#Define sampling rate and normalization factor 
#1440 = number of minutes in a day (to get result in cycles/per day) 
Fs = np.divide(1,1440) 
norm = n*Fs

#Calculate frequency
frq = np.divide(rng,norm)

#Calculate power 
n3 = np.floor(n2
rng2 = range(1,int(n3))
pwr = np.abs(fftVar[rng2])**2 

#Remove frequencies not needed 
frqThreshold = 0.9 
for j in range(len(frq)):
    if frq[j] >= frqThreshold: 
        fftVar[j] = 0.0

        #Found this section online, not sure what this line does
        #When uncommented, I get bizzare result (see 3rd example below)
        #fftVar[int(n2) + j] = 0.0 

postFFT = np.fft.ifft(fftVar)

绘制功率与频率的关系图给出了以下内容（例如一个变量）：

Before diurnal/seasonal cycles were removed

After diurnal/seasonal cycles were removed

所以我有理由相信这一部分是有效的（我知道这是一个非常粗糙的方法，在这个阶段是可以接受的）。如果我做得对的话，y轴应该是功率，x轴应该是频率，白天周期的峰值是1个周期/天，然后另一个峰值显示了更长的周期。在

这个问题（我认为）发生在我处理去掉频率后产生的时间序列中。以下是我删除频率前后的一个变量（与上述相同的示例变量）的曲线图：

Example time series of variable before

Same variable and time after the FFT analysis

Same variable and time after the FFT analysis, with specified line of code uncommented (see code block)

请注意，这两个图的x轴都是时间（长度为1天），而y轴是原始值（在本例中是Pa），应该是第二个曲线图中没有日/季节循环的平均值的异常。在

我是一个信号处理的初学者，所以我甚至不能开始理解第三个绘图显示的是什么，或者这是“更正确”还是“不太正确”。在

我用以下简单代码绘制了图：

^{pr2}$

然而，频率似乎没有被移除。在代码中，我认为只要频率高于某个阈值，就将fftVar的结果设置为0，然后执行逆变换就可以做到这一点，但我错了（或者至少编码错误）。在

我是python的第一次用户，毫无疑问，有很多方法可以改进上面的代码（如果有人有时间，我也很感谢这里的任何提示），但我最感兴趣的是为什么这段代码没有删除所需的频率（或者可能只是一个绘图问题）？在

谢谢你的帮助！在

Answer 1

我在理解你的代码时遇到了一些问题，但我认为问题是你没有把负频率下的傅立叶变换也设为零。这里有一些编造的数据代码，其中包括季节性、双季节性和日变化，从中过滤出日变化。在

import numpy as np
import matplotlib.pyplot as plt
plt.ion()


t = np.linspace(0,10*365,1e5)

y = np.sin(2*np.pi*t) + np.sin(2*np.pi*t/365) + np.sin(2*np.pi*t/365/2)
y = y - y.mean()

#Calculate FFT
freqs = np.fft.fftshift(np.fft.fftfreq(len(y),t[1]-t[0]))
fftVar = np.fft.fftshift(np.fft.fft(y))

fft_filtered = fftVar.copy()
# set spectrum above frequency treshold to zero. It is important to also set the negative frequencies to zero, 
# which probably caused your problem
fft_filtered[np.abs(freqs)>1e-1] = 0 


postFFT = np.fft.ifft(np.fft.ifftshift(fft_filtered))
plt.figure(1,figsize=(9,7))
ax = plt.subplot(211)

ax.plot(t,y)
ax.plot(t,postFFT)
ax.set_xlim(t.min(),t.max())
ax.set_xlabel('Time [days]')

ax1 = plt.subplot(212)
ax1.loglog(freqs,np.abs(fftVar)**2)
ax1.plot(freqs,np.abs(fft_filtered)**2)
ax1.set_xlim(1e-3,freqs.max())
ax1.set_ylim(1e-6,1e10)
ax1.set_xlabel('Frequency [1/days]')

结果如下图所示，显示了过滤和未过滤的时间序列和频谱。厚厚的蓝色斑点是昼间振荡，由于时间跨度太大（10年），这一点在图中没有得到解决。

Answer 2

您的数据以1分钟的间隔获取，因此Fs=1/60而不是1/1440。在

我不清楚什么是完全不适合你，但试试这个修改后的版本。在

from scipy.fftpack import fftfreq

# Sampling rate
Fs = 1.0/60.0

#Calculate the anomalies from the mean of time series y
y = y - y.mean()

#Calculate FFT 
fftVar = np.fft.fft(y)

#calculate the frequencies.
#Here 1/Fs is the sampling interval, which is 1 minute (60 seconds)
f = fftfreq(len(y), 1/Fs)

#Remove frequencies not needed 
frqThreshold = 0.9 * Fs/2

fftVar[(f >= frqThreshold) | (f <= -frqThreshold)] = 0.0

postFFT = np.real(np.fft.ifft(fftVar))

另外，要进行数字滤波，您可能需要考虑使用butter和{}from{}。在