第1節.Matlab中的fft

05-12

傅里葉變化（FFT）是光學數值模擬（特別是光學傳播中）的重要工具。

在matlab中實現FFT是依靠fft函數。那麼我們在實現FFT時，一定要知道matlab中的fft做了什麼，以及與理論上傅里葉變換的對應關係。

在 matlab Documentation中，可以看到

Y = fft(X) computes the discrete Fourier transform (DFT) of X using a fast Fourier transform (FFT) algorithm.

基於Discrete Fourier transform實現的是

$Y(k)=sum_{j=1}^{n}{X(j)W_n^{(j-1)(k-1)}}$

where $W_n=e^{(-2pi i)/n}$

而FFT對應的積分描述是

$Y(omega)=int_{-infty}^{+infty}x(t)e^{-iomega t}dt$

我們來看一下：

整數 $k$ 對應角頻率 $omega$ ，整數 $j$ 對應時間 $t$ . 要知道，matlab中的下標是從1開始的，所以 $j=1$ 對應的時刻是0時刻。同時，我們可以看到積分的上下限是 $-infty$ 到 $+infty$ 。而Discrete Fourier transform（DFT）確只用了定義域非 $-infty$ 到 $+infty$ 的數據。

其實這裡DFT用了周期延拓（函數的延拓就是把一個區間上的函數拓展到整個區間.方法是利用周期函數的性質，其中原區間的長度為一個周期）。

我們可以這樣子理解，圖一的是我們的輸入信號。

圖1.信號

而在DFT看來，這裡信號是長這樣子（如圖2）。

圖2.周期延拓後的信號

也就是，把這個信號複製無窮個，並讓他們首尾相接了。因為這個首尾相接，其實會導致一些問題，我們在編程中需要避免的。（主要是首尾相接可能會照成突變，而面對這種情況可以採用切趾函數來實現，或者說加窗口，就不展開了。）

離散化後看是這樣子的（如圖3）。

圖3.周期延拓後的離散信號

所以時域信號 $x(t)$ 的採樣數據 $X(j)$ 中， $t$ 和 $j$ 滿足 $t/T=(j-1)/n$ 。

我們回去看DFT的定義，以及傅里葉變換的定義，可以得到 $omega t=2pi(j-1)(k-1)/n$ 。

所以可以得到 $omega=2pi(k-1)/T$ 。

下面講講負時間和負角頻率。

負時間( $t<0$ )其實是比較好從圖3看出來的，因為是周期延拓，我們可以用 $j=1$ 之前的數據來當作是 $t<0$ 的數據，我一般喜歡用 $[n/2,n]$ 這個區間的數據當作是負時間。

我們從前面好像比較難看到負角頻率的跡象，其實負頻率是由於採樣導致的。我們可以用一個梳狀函數去描述採樣數據， $sum_{1}^{n}{X(j)delta (t-(j-1)Delta t))}=x(t)comb(t/Delta t)$ 。

圖4.梳狀函數去採樣數據

由於梳狀函數的存在，他的譜分布是這樣子的。其實把頻譜分布在頻率域周期延拓了。

圖4.採樣數據的頻譜

所以，而 $2pi/Delta t$ 正是 $k=n+1$ 時的頻率 $omega=2pi(n+1-1)/T=2pi/Delta t$ 。負角頻率( $omega<0$ )圖4看出來的，因為是周期延拓，我們可以用 $k=1$ 之前的數據來當作是 $omega<0$ 的數據，我一般喜歡用 $[n/2,n]$ 這個區間的數據當作是負角頻率。

其實我平時是比較喜歡用頻率而不是角頻率。所以，按定義 $f=omega/(2pi)$ 就可以了。

所以，就有了下面這個函數：

function [t,ft]=domain(window,sample);%構建變數空間% [t,ft]=domain(window,sample);% 輸入% window - 空間/時間窗口% sample - 採樣點數目%% 輸出% t - 空間/時間變數% ft - 空間頻率/時間頻率變數if(sample <= 1) t = 0; ft = 0;else if mod(sample,2)==1 error([sample should be even number ]) else t = linspace(-window/2,window/2,sample+1); t(end) = []; t = fftshift(t); ft =[0sample/2-1) -sample/2:-1]/window; endend

採樣點數我強制是2的倍數，也是避免一些麻煩的東西。

用這個函數生成變數空間。

[t,ft] = domain(1,512);subplot(2,1,1)plot(t)title(t)subplot(2,1,2)plot(ft)title(ft)

圖5.變數數據

生成的數據可能不太直觀，但是其實這樣子做使得後期模擬不容易出錯。我們只要在繪圖時統一加上fftshift就可以了，雖然在後期處理比較麻煩，但是，我們後面模擬可以看到，這樣子計算處理的結果跟理論值計算是一樣。但是，如果不這麼做，為了跟理論計算一致，就需要不斷在公式中加入fftshift來實現。而這樣子定義變數，我們在模擬時，使用公式就不需要加入fftshift，只需要編寫公式就可以，這個其實可以帶來巨大的方便，以及不容易出錯。

[t,ft] = domain(1,512);subplot(2,1,1)plot(fftshift(t))title(t)subplot(2,1,2)plot(fftshift(ft))title(ft)

圖6.ffthsift後變數數據

最後用一個計算結束第1節。生成若干頻率的信號疊加，計算頻率。根據我們的相位設置，得到頻率2Hz的是cos函數，5Hz的是sin函數，7Hz的是-cos函數。計算結果的實部和虛部都是跟信號分析之類的書上面的結果是一樣的。

[t,ft] = domain(1,512);% 生成若干頻率的信號疊加deg = pi/180;freq = [2,5,7];fai = [0 90 180]*deg;X= zeros(size(t));for index = 1:length(freq) freq_temp = freq(index); fai_temp = fai(index); X = X+cos(2*pi*freq_temp*t+fai_temp);end%傅里葉變換Y= fft(X);subplot(2,1,1)plot(fftshift(t),fftshift(X));title(X(t))subplot(2,1,2)bar(fftshift(ft),fftshift(real(Y)),red);hold onbar(fftshift(ft),fftshift(imag(Y)),blue);hold offxlim([-10,10])title(Y(f))