不使用ifft重建时间序列数据的FFT结果

Question

我分析了下面的sunspots.dat数据，使用了fft，这在这个领域是个经典的例子。我从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

mtrw的回答非常有帮助，让我解决了和提问者一样的问题，不过我在理解嵌套循环的时候差点儿头炸了。

这是最后一部分，不过用了numpy的广播功能（不确定在提问时这个功能是否已经存在），而不是调用f函数：

xs = np.arange(N)
omega = 2*np.pi/N
phase = omega * xs[:,None] * xs[None,:]
reconstruct = Y[None,:] * (np.cos(phase) + 1j*np.sin(phase))
reconstruct = (reconstruct).sum(axis=1).real / N

# same output
plt.plot(reconstruct)
plt.plot(wolfer)

Answer 2

当你调用 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()