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