2013-10-10 249 views
2

我有一個奇怪的問題與離散fft。我知道高斯函數exp(-x^2/2)的傅里葉變換又是相同的高斯函數exp(-k^2/2)。我試圖用MatLab和FFTW中的一些簡單代碼來測試,但我得到了奇怪的結果。FFTW和fft與MatLab

首先,結果的虛部不是零(在MatLab中),因爲它應該是。

其次,實部的絕對值是高斯曲線,但沒有絕對值的一半的模式具有負係數。更確切地說,每個第二模式都有一個係數,它是應該是負數的係數。第三,得到的高斯曲線的峯值(在取實部的絕對值之後)不是一而是高得多。其高度與x軸上的點數成正比。但是,比例因子不是1,但接近1/20。

任何人都可以解釋我做錯了什麼?

這裏是MATLAB代碼,我用:

function [nooutput,M] = fourier_test 

    Nx = 512;  % number of points in x direction 

    Lx = 50;  % width of the window containing the Gauss curve 

    x = linspace(-Lx/2,Lx/2,Nx);  % creating an equidistant grid on the x-axis 

    input_1d = exp(-x.^2/2);     % Gauss function as an input 
    input_1d_hat = fft(input_1d);   % computing the discrete FFT 
    input_1d_hat = fftshift(input_1d_hat); % ordering the modes such that the peak is centred 

    plot(real(input_1d_hat), '-') 
    hold on 
    plot(imag(input_1d_hat), 'r-') 
+0

您可能將* continuous * FT與* discrete * FT混合在一起。高斯的連續FT是高斯函數,但我不確定DFT(FFT)是如此。 –

+2

當然,離散性使得事物不精確,但離散FT被創建爲連續FT的近似值。當點數和間隔足夠大時,這裏就是這種情況,近似值應該是好的。會有一些偏差,但他們絕對不應該那麼大,如此係統。另外,當我增加Nx和Lx時,結果不會改變,這意味着收斂是可以實現的。 –

+1

我認爲在你的DFT情況下可能會出現一個隱含的時間偏移,這就意味着頻域的相位旋轉 - 我認爲幅度是正確的,相位是一條直線(如果你是模2π不解開它,即鋸齒)。 –

回答

4

答案基本上是保羅 - [R表明在他的第二個意見,你介紹的相移(線性依賴於頻率),因爲中心由input_1d_hat描述的高斯實際上在k>0,其中k+1input_1d_hat的索引。相反,如果你Center的資料(這樣input_1d_hat(1)對應中心)如下你在頻域中的相位校正高斯:

Nx = 512;  % number of points in x direction 
    Lx = 50;  % width of the window containing the Gauss curve 

    x = linspace(-Lx/2,Lx/2,Nx);  % creating an equidistant grid on the x-axis 

    %%%%%%%%%%%%%%%% 
    x=fftshift(x); % <-- center 
    %%%%%%%%%%%%%%%% 

    input_1d = exp(-x.^2/2);     % Gauss function as an input 
    input_1d_hat = fft(input_1d);   % computing the discrete FFT 
    input_1d_hat = fftshift(input_1d_hat); % ordering the modes such that the peak is centered 

    plot(real(input_1d_hat), '-') 
    hold on 
    plot(imag(input_1d_hat), 'r-') 

從DFT的定義,如果高斯不集中等那最大值發生在k=0,你會看到相位扭曲。 fftshift的效果是執行數據集左側和右側的循環移位或交換,這相當於將峯的中心移至k=0

至於幅度縮放,這是在Matlab中實現的DFT的定義的問題。從對FFT的文檔:

For length N input vector x, the DFT is a length N vector X, 
with elements 
       N 
    X(k) =  sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N. 
       n=1 
The inverse DFT (computed by IFFT) is given by 
       N 
    x(n) = (1/N) sum X(k)*exp(j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N. 
       k=1 

注意的是,在向前邁進了一步的總和由N.標準化因此如果增加點Nx的數量的總和,同時保持寬度高斯函數常數的Lx,您將按比例增加X(k)

至於信號泄漏到假想頻率維度,也就是說是由於DFT,這導致截斷和其它效果的離散形式,如再由Paul R.注意到如果減少Lx同時保持恆定Nx,你應該看到在虛擬維數相對於的信號量減少到實際維度(比較光譜,同時保持實際維度中的峯值強度相等)。

你會發現類似問題的其他答案herehere

+0

謝謝你的提示。但是,我仍然不明白爲什麼我應該用fftshift來更改x陣列。你能詳細解釋一下嗎?我認爲這是通過在輸出上做fftshift來提供的。爲什麼我的版本中的高斯有效地以x> 0爲中心? –

+0

你對我有其他兩個問題的原因有什麼想法嗎?我實施你的解決方案。它修復了輸出兩部分的振盪,但虛部仍然不爲零。而且,當我改變點數Nx時,它保持不變,而實部的峯值與Nx成比例變化。任何想法爲什麼這是這樣的? –

+0

任何想法如何解決與FFTW使用C++時的問題?在那裏,我得到了同樣令人討厭的震盪。 你能給出關於這個問題的數學工具參考嗎?我沒有看到關於fft例程解釋頁面上的評論。 –