ifftを使用せずにFFT結果を使って時系列データを再作成する

Question

私はこの領域の古典的な例であるfftを使ってsunspots.datデータ（以下）を解析しました。私はfftの結果を実際の部品とimaginery部品で得ました。その後、フーリエ変換の式に従ってデータを再作成するために、これらの係数（最初の20）を使用しようとしました。 B_NするA_Nとimagineryに対応し、実際の部品を考えて、私はいくつかの理由

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によって生成されたものを（私は両側の係数の同じ番号を使用）ミラー化されません。何が間違っていますか？

データ：あなたはfft(wolfer)を呼ばれると

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

Answer 1

、あなたはデータの長さに等しい基本周期を想定する変換を語りました。データを再構成するには、基本周期= 2*pi/Nの基本関数を使用する必要があります。同じトークンで、あなたの時間インデックスxsは元の信号の時間サンプルの範囲を超えなければなりません。

もう一つの間違いは、完全な複素乗算に忘れることでした。これはY[omega]*exp(1j*n*omega/N)と考えるのが簡単です。

ここに固定コードがあります。注意sqrt(-1)、n〜Nと混同しないように、iをctrに変更しました。サンプルの小文字と大文字のサンプルの長さの通常の信号処理規則に従ってください。また、整数除算についての混乱を避けるために__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*pi/N) + 1j*np.sin(x*ctr*2*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()