使用 FFT 结果重新创建时间序列数据而不使用 ifft

Question

我使用 fft 分析了 sunspots.dat 数据（如下），这是该领域的经典示例。我从 fft 获得了真实和想象部分的结果。然后我尝试使用这些系数（前 20 个）按照傅里叶变换的公式重新创建数据。认为真实部分对应a_n，想象部分对应b_n，我有

import numpy as np
from scipy import *
from matplotlib import pyplot as gplt
from scipy import fftpack

def f(Y,x):
    total = 0
    for i in range(20):
        total += Y.real[i]*np.cos(i*x) + Y.imag[i]*np.sin(i*x)
    return total

tempdata = np.loadtxt("sunspots.dat")

year=tempdata[:,0]
wolfer=tempdata[:,1]

Y=fft(wolfer)
n=len(Y)
print n

xs = linspace(0, 2*pi,1000)
gplt.plot(xs, [f(Y, x) for x in xs], '.')
gplt.show()       

然而，出于某种原因，我的情节并未反映 ifft 生成的情节（我在两侧使用相同数量的系数）。有什么问题？

数据：

http://linuxgazette.net/115/misc/andreasen/sunspots.dat

Answer 1

当您调用fft(wolfer) 时，您告诉转换假设一个基本周期等于数据的长度。要重建数据，您必须使用相同基本周期 = 2*pi/N 的基函数。同理，您的时间索引xs 必须在原始信号的时间样本范围内。

另一个错误是忘记做完全复数乘法。将其视为Y[omega]*exp(1j*n*omega/N) 更容易。

这是固定代码。注意我将i 重命名为ctr 以避免与sqrt(-1) 混淆，并将n 重命名为N 以遵循通常的信号处理约定，即对样本使用小写字母，对总样本长度使用大写字母.我还导入了 __future__ division 以避免混淆整数除法。

之前忘记添加了： 请注意，SciPy 的 fft 在累加后不会除以 N。在使用Y[n] 之前我没有把它分开；如果您想获得相同的数字，而不是仅仅看到相同的形状，您应该这样做。

最后，请注意，我是在整个频率系数范围内求和。当我绘制np.abs(Y) 时，看起来高频中存在显着值，至少在样本 70 左右之前是这样。我认为通过在整个范围内求和，查看正确的结果，然后缩减系数并查看会发生什么，会更容易理解结果。

from __future__ import division
import numpy as np
from scipy import *
from matplotlib import pyplot as gplt
from scipy import fftpack

def f(Y,x, N):
    total = 0
    for ctr in range(len(Y)):
        total += Y[ctr] * (np.cos(x*ctr*2*np.pi/N) + 1j*np.sin(x*ctr*2*np.pi/N))
    return real(total)

tempdata = np.loadtxt("sunspots.dat")

year=tempdata[:,0]
wolfer=tempdata[:,1]

Y=fft(wolfer)
N=len(Y)
print(N)

xs = range(N)
gplt.plot(xs, [f(Y, x, N) for x in xs])
gplt.show()