Matlab模擬傅立葉變換
傅立葉變換是我們最早開始接觸的時頻域變換方法，雖然經常使用，知道怎麼用紙筆計算，但是還從來沒有在電腦中模擬過，正好現在開始學習數字訊號處理，藉著這個機會再學習如何在電腦上模擬傅立葉變換。
以下大部分內容來自Digital Signal Processing Using Matlab和數字訊號處理教程 程佩青
此次選擇的軟體平臺為Matlab。
由於Matlab無法處理無限長序列，所以需要處理的訊號必須是有限長的。
連續時間傅立葉變換
傅立葉變換的公式為：
為了在計算機中模擬傅立葉變換，我們將積分變為求和的方式，上下限也從正無窮到負無窮變為一段長度M，dt需要儘可能小
在Matlab中，函式的自變數因變數的集合都是使用矩陣來儲存的，從矩陣的角度來看傅立葉變換的公式如下：
角頻率向量定義為
時間向量定義為
因此矩陣指數可寫為
整個傅立葉變換可寫為
Xa = xa * exp(-1j*t'*W) * Dt;
具體實現
其實下面這個例子是Digital Signal Processing Using Matlab中的，來自P64頁，不過想到都看到這裡了還要讀者翻書不太好，就一起放上來了。
定義
先進行數學上的分析，
MATLAB實現如下：
% Analog Signal Dt = 0.00005; t = -0.005:Dt:0.005; xa = exp(-1000*abs(t)); % Continuous-time Fourier Transform Wmax = 2*pi*2000; K = 500; k = 0:1:K; W = k*Wmax/K; Xa = xa * exp(-1j*t'*W) * Dt; Xa = abs(Xa); W = [-fliplr(W), W(2:501)]; Xa = [fliplr(Xa), Xa(2:501)]; subplot(2,1,1); plot(t*1000,xa); xlabel('t in msec.'); ylabel('xa(t)'); title('Analog Signal'); subplot(2,1,2); plot(W/(2*pi*1000),Xa*1000); xlabel('Frequency in KHz'); ylabel('Xa(jW)*1000'); title('Continuous-time Fourier Transform');
執行效果如下：
如果想確認變換的正確性，可以在執行完上面這個指令碼後，在命令列輸入
plot(W/(2*pi*1000),(0.002./(1+(W./1000).^2))*1000);
xlabel('Frequency in KHz'); ylabel('Xa(jW)*1000');
執行效果如下：
這時會發現，根據上面推導的變換公式直接plot出的圖形和變換後得到的圖形是一樣的，這樣可以確定變換的正確性。
存在問題
目前存在的問題是，對於複函式的變換結果不正確。我想了很多天都找不出問題所在，只能暫時放棄，等以後有機會再研究。
離散時間傅立葉變換
下面是對上一個例子中的模擬輸入訊號做離散化，然後再進行離散傅立葉變換。
為了體現Nyquist定理，將使用兩種不同的取樣頻率
-
使用Fs=5000sam/sec取樣來獲得x1(n)
-
使用Fs=1000sam/sec取樣來獲得x2(n)
% Analog Signal
Dt = 0.00005;
t = -0.005:Dt:0.005;
xa = exp(-1000*abs(t));
% Discrete-time Signal
Ts = 0.0002;
n = -25:1:25;
x = exp(-1000*abs(n*Ts));
% Discrete-time Fourier transform
K = 500;
k = 0:1:K;
w = pi*k/K;
X = x*exp(-j*n'*w); X = real(X);
w = [-fliplr(w), w(2:K+1)];
X = [fliplr(X), X(2:K+1)];
subplot(2,1,1);plot(t*1000,xa);
xlabel('t in msec.');
ylabel('x1(n)');
title('Discrete Signal');hold on;
stem(n*Ts*1000,real(x));gtext('Ts=0.2 msec');hold off;
subplot(2,1,2);plot(w/pi,X);
xlabel('Frequency in pi units');ylabel('X1(w)');
title('Discrete-time Fourier Transform');
Fs=5000sam/sec
xa(t)的頻率為2KHz，因此它的Nyquist頻率為4KHz，而它的取樣頻率為5KHz，所以是滿足Nyquist取樣定律的，此時不會發生混疊。
執行效果如下：
Fs=1000sam/sec
這裡使用的取樣頻率為1KHz，不滿足Nyquist條件，因此會發生混疊。觀察一下就會發生，1KHz取樣得到的序列的頻域波形和前面的頻域波形不同，這就是混疊導致的，而且過低的取樣率採集的訊號的變換的不可逆的。
執行效果如下：