Matlab模擬傅立葉變換
傅立葉變換是我們最早開始接觸的時頻域變換方法，雖然經常使用，知道怎麼用紙筆計算，但是還從來沒有在電腦中模擬過，正好現在開始學習數字訊號處理，藉著這個機會再學習如何在電腦上模擬傅立葉變換。
以下大部分內容來自 Digital Signal Processing Using Matlab 和 數字訊號處理教程 程佩青
此次選擇的軟體平臺為Matlab。
由於Matlab無法處理無限長序列，所以需要處理的訊號必須是有限長的。
連續時間傅立葉變換
傅立葉變換的公式為：
\[ X_a(j\Omega)=\int x_a(t)e^{-j\Omega t}dt \]
為了在計算機中模擬傅立葉變換，我們將積分變為求和的方式，上下限也從正無窮到負無窮變為一段長度M，dt需要儘可能小
\[ X_a(j\Omega) = \sum_m x_a(m\Delta t)e^{-j\Omega m\Delta t}\Delta t \]
在Matlab中，函式的自變數因變數的集合都是使用矩陣來儲存的，從矩陣的角度來看傅立葉變換的公式如下：
\[ [X_a(0)\ X_a(1)\ X_a(2)\ ..] = [x_a(0)\ x_a(1)\ x_a(2)\ ..] \left[ \begin{matrix} e^{-j\omega_0 t_0} & e^{-j\omega_1 t_0} & \cdots & e^{-j\omega_K t_0} \\ e^{-j\omega_0 t_1} & e^{-j\omega_1 t_1} & \cdots & e^{-j\omega_K t_1} \\ \vdots & \vdots & \ddots & \vdots \\ e^{-j\omega_0 t_N} & e^{-j\omega_1 t_N} & \cdots & e^{-j\omega_K t_N} \\ \end{matrix} \right] \]
角頻率向量定義為 \(\omega=[\omega_0\ \omega_1\ ...\ \omega_K]\)
時間向量定義為 \(t=[t_0 :\Delta t: t_N]\)
因此矩陣指數可寫為 \(-j*t'*\omega\)
整個傅立葉變換可寫為
Xa = xa * exp(-1j*t'*W) * Dt;
具體實現
其實下面這個例子是 Digital Signal Processing Using Matlab 中的，來自P64頁，不過想到都看到這裡了還要讀者翻書不太好，就一起放上來了。
定義 \(x_a(t) = e^{-1000|t|}\)
先進行數學上的分析，
\[ \begin {aligned} X_a(j \Omega) &= \int^\infty_{-\infty}x_a(t)e^{-j\Omega t}dt \\ &= \int^0_{-\infty}e^{1000t}e^{-j\Omega t}dt + \int^\infty_0 e^{-1000t}e^{-j\Omega t}dt \\ &= \frac{0.002}{1+(\frac{\Omega}{1000})^2} \end {aligned} \]
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定理，將使用兩種不同的取樣頻率 1. 使用Fs=5000sam/sec取樣來獲得x1(n) 2. 使用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取樣得到的序列的頻域波形和前面的頻域波形不同，這就是混疊導致的，而且過低的取樣率採集的訊號的變換的不可逆的。
執行效果如下：
