在不使用ifft的情况下使用FFT结果重新创建时间序列数据 [英] Recreating time series data using FFT results without using ifft

问题描述

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

I analyzed the sunspots.dat data (below) using fft which is a classic example in this area. I obtained results from fft in real and imaginery parts. Then I tried to use these coefficients (first 20) to recreate the data following the formula for Fourier transform. Thinking real parts correspond to a_n and imaginery to b_n, I have

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生成的图(我在两侧使用相同数量的系数).有什么问题吗?

For some reason however, my plot does not mirror the one generated by ifft (I use the same number of coefficients on both sides). What could be wrong ?

数据:

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

推荐答案

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

When you called fft(wolfer), you told the transform to assume a fundamental period equal to the length of the data. To reconstruct the data, you have to use basis functions of the same fundamental period = 2*pi/N. By the same token, your time index xs has to range over the time samples of the original signal.

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

Another mistake was in forgetting to do to the full complex multiplication. It's easier to think of this as Y[omega]*exp(1j*n*omega/N).

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

Here's the fixed code. Note I renamed i to ctr to avoid confusion with sqrt(-1), and n to N to follow the usual signal processing convention of using the lower case for a sample, and the upper case for total sample length. I also imported __future__ division to avoid confusion about integer division.

忘了早点添加:请注意，SciPy的fft在累加后不会除以N.在使用Y[n]之前，我没有对此进行划分；如果您想找回相同的数字，而不只是看到相同的形状，就应该这样做.

forgot to add earlier: Note that SciPy's fft doesn't divide by N after accumulating. I didn't divide this out before using Y[n]; you should if you want to get back the same numbers, rather than just seeing the same shape.

最后，请注意，我正在对整个频率系数范围进行求和.当我绘制np.abs(Y)时，它看起来在较高频率中有明显的值，至少直到样本70左右.我认为，通过对整个范围求和，查看正确的结果，然后回算系数并查看会发生什么，将更容易理解结果.

And finally, note that I am summing over the full range of frequency coefficients. When I plotted np.abs(Y), it looked like there were significant values in the upper frequencies, at least until sample 70 or so. I figured it would be easier to understand the result by summing over the full range, seeing the correct result, then paring back coefficients and seeing what happens.

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

这篇关于在不使用ifft的情况下使用FFT结果重新创建时间序列数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！