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

我使用 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 生成的图（我在两侧使用相同数量的系数）。可能出什么问题了？

Data:

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

当你打电话时fft(wolfer)，您告诉变换假设基本周期等于数据的长度。要重建数据，您必须使用相同基本周期的基函数 =2*pi/N。同样的道理，你的时间指数xs必须在原始信号的时间样本范围内变化。

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

这是固定代码。注意我改名了i to ctr以避免混淆sqrt(-1), and n to 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()