2016-08-10 3 views
1

fft - 出力の逆FFTを手動で計算しようとしています。私は最初にfftを使ってデータセットのFFTを計算する次のスクリプトを使用しています。私は手動で逆FFTを見つけようとしますが、それはifftの結果と似ていません。手動で逆FFTを計算する

あなたは自分のエラーを発見できますか?私は単に

data = [ 
    -0.0005 
    -0.0004 
    -0.0003 
    -0.0002 
    -0.0001 
    -0.0000 
    0.0001 
    0.0001 
    0.0001 
    0.0002 
    0.0002 
    0.0002 
    0.0002 
    0.0002 
    0.0002 
    0.0002 
    0.0002 
    0.0002 
    0.0003 
    0.0004 
    0.0005 
    0.0006 
    0.0007 
    0.0009 
    0.0010 
    0.0011 
    0.0011 
    0.0012 
    0.0011 
    0.0011 
    0.0011 
    0.0010 
    0.0011]; 


delta = 0.0125; 
fs = 1/delta; 
x  = (0:1:length(data)-1)/fs; 

X=fft(data); 

%find fft 
N=length(data); 
ws = 2*pi/N; 
wnorm = -pi:ws:pi; 
wnorm = wnorm(1:length(x)); 
w = wnorm*fs; 
figure(2) 
plot(w/(2*pi),abs(fftshift(X))) 


%find inverse fft manually 
for m=1:length(X) 
    for k=1:length(data) 
     X_real(m) = X(k)*exp(i*k*ws*(m-1)); 
    end 
end 

figure(3) 
plot(1:length(data), abs(X_real), 1:length(data), ifft(X)) 
+1

まず、あなたの 'exp'の部分が実際に正しいかどうか分からないうちに、' m'ごとにX_real(k)の 'sum'が足りなくなりました。 – GameOfThrows

+1

@GameOfThrowsは、X_real(m)= X_real(m)+ exp(i * k * ws *(m-1)); (i * k * ws *(m-1)); '(' X_real'をあらかじめゼロの配列に初期化する必要があるかもしれません) –

答えて

2

https://en.wikipedia.org/wiki/Fast_Fourier_transform#Definition_and_speedは、以下のようなループのためにあなたを変更してください、FFTの標準逆式がここで提示し使用しています。

for m=1:length(X) 
    for k=1:length(data) 
     temp(k) = X(k)*exp(i*(m-1)*ws*(k-1)); 
    end 
    X_real(m)=(1/N)*sum(temp); 
end 

figure(3) 
plot(1:length(data), real(X_real)) 

ifftの方程式は、hereで見つけることができます。

あなたは2つのことを逃しました。

1つは正規化、もう1つは加算です。

+1

3番目のものは 'k-1'(' k'だけでなく、彼らがオフ・バイ・ワンのエラーについて何を言っているか知っています)です。 –

+0

@AhmedFasih私は '1..length(data)'から 'k 'を集計していますが、' -0.5fs'から '0.5fs'まで走る' w /(2 * pi) 'を使って' FFT'Iプロットをプロットしています。 。これらの2つの周波数の間の対応は何ですか?私は、プロットの '<0.25fs'に対応する逆行列にコンポーネントを含めるだけで、' k'はそれに対応しますか? – BillyJean

+0

@BillyJean私はこのような周波数ベクトルを作る: 'N = length(X_real); * 1秒あたりのサイクル数*(つまり、あなたの入力信号が毎秒 'fs'サンプルでサンプリングされたと仮定して、Hz)の単位を持つf =(0:N-1)/ N * fs'です。 'plot(f、abs(X_real))'が正しくなり、 'xlim(0、0.25 * fs)'を実行して制限を設定することができます。 –