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

Question

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

7

我使用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

- BBSysDyn

2

只是出于好奇，您在处理频谱方面做了什么？如果您正在尝试确定各个组件的相对光谱幅度，您可能需要使用数据窗口（en.wikipedia.org/wiki/Window_function）。例如，如果您绘制 np.abs(fft(wolfer*hanning(len(wolfer))))，则n=30附近的峰值显示比没有窗口时更多的结构。 - mtrw

2个回答

0

mtrw的回答非常有帮助，帮助我回答了与OP相同的问题，但是我几乎要被理解嵌套循环的过程搞炸了。

这里是最后一部分，但使用了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)

- SuperCodeBrah

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- mtrw · Accepted Answer

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