快速傅立葉變換（MATLAB實現）
轉自：http://blog.sina.com.cn/s/blog_a1d5b9ba0102vxj2.html
一、快速傅立葉介紹
傅立葉原理表明：任何連續測量的時序或訊號，都可以表示為不同頻率的餘弦（或正弦）波訊號的無限疊加。FFT是離散傅立葉變換的快速演算法，可以將一個訊號變換到頻域。那其在實際應用中，有哪些用途呢？
1.有些訊號在時域上是很難看出什麼特徵的，但是如果變換到頻域之後，就很容易看出特徵（頻率，幅值，初相位）；
2.FFT可以將一個訊號的頻譜提取出來，進行頻譜分析，為後續濾波準備；
3.通過對一個系統的輸入訊號和輸出訊號進行快速傅立葉變換後，兩者進行對比，對系統可以有一個初步認識。
假設取樣頻率Fs，訊號頻率F，訊號長度L，取樣點數N。那麼FFT之後結果就是一個為N點的複數。每一個點就對應著一個頻率點。這個點的模值，就是該頻率值下的幅度特性。
具體跟原始訊號的幅度有什麼關係呢？
1. 假設原始訊號的峰值為A，那麼FFT的結果的每個點（除了第一個點直流分量之外）的模值就是A的N/2倍，而第一個點就是直流分量（即0Hz），它的模值是直流分量的N倍；
2. 每個點的相位呢，就是在該頻率下的訊號的相位。第一個點表示直流分量，它的相位是該頻率的初相位，matlab以cos為底的，若訊號時正弦形式sin(t)，則變成cos(t-pi/2)即可。
取樣頻率Fs,被N-1個點平均分成N等份，每個點的頻率依次增加。為了方便進行FFT運算，通常N取大於訊號長度L的2的整數次方。
例如某點n所表示的頻率為：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到頻率為為Fs/N。如果取樣頻率Fs為1024Hz，取樣點數為1024點，則可以分辨到1Hz。
1024Hz的取樣率取樣1024點，剛好是1秒，也就是說，取樣1秒時間的訊號並做FFT，則結果可以分析到1Hz。如果取樣2秒時間的訊號，則N為2048，並做FFT，則結果可以分析到0.5Hz。
如果要提高頻率分辨力，則必須增加取樣點數，也即取樣時間。頻率解析度和取樣時間是倒數關係。
假設FFT之後某點n用複數a+bi表示，該複數的模就是An=sqrt（a*a+b*b），相位就是Pn=atan2(b,a)。根據以上的結果，就可以計算出n點（n≠1，且n<=N/2）對應的訊號的表示式為：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)；對於n=1點的訊號，是直流分量，幅度即為A1/N。
由於FFT結果的對稱性，通常我們只使用前半部分的結果，即小於取樣頻率一半的結果。
二、例子
假設我們有一個訊號，它含有5V的直流分量，頻率為50Hz、相位為-30度、幅度為7V的交流訊號以及一個頻率為90Hz、相位為90度、幅度為3V的交流訊號。數學表示式為：
x = 5 + 7*cos(2*pi*15*t - 30*pi/180) + 3*cos(2*pi*40*t - 90*pi/180)。
我們以128Hz的取樣率對這個訊號進行取樣，總共取樣256點。按照我們上面的分析，Fn=(n-1)*Fs/N，我們可以知道，每兩個點之間的間距就是0.5Hz。我們的訊號有3個頻率：0Hz、15Hz、40Hz
出於程式設計方便，因為直流分量的幅值A1/N，其他點幅值為An/(N/2)，故直流分量最後要除以2才是對的。
一般FFT所用資料點數N與原含有訊號資料點數L相同，這樣的頻譜圖具有較高的質量，可減小因補零或截斷而產生的影響。
三、Matlab程式碼
clear
Fs = 128; % 取樣頻率
T = 1/Fs; % 取樣時間
L = 256; % 訊號長度
t = (0:L-1)*T; % 時間
x = 5 + 7*cos(2*pi*15*t - 30*pi/180) + 3*cos(2*pi*40*t - 90*pi/180); %cos為底原始訊號
y = x + randn(size(t)); %新增噪聲
figure;
plot(t,y)
title('加噪聲的訊號')
xlabel('時間(s)')
>> N = 2^nextpow2(L); %取樣點數，取樣點數越大，分辨的頻率越精確，N>=L，超出的部分訊號補為0
Y = fft(y,N)/N*2; %除以N乘以2才是真實幅值，N越大，幅值精度越高
f = Fs/N*(0:1:N-1); %頻率
A = abs(Y); %幅值
P = angle(Y); %相值
figure;
subplot(211);plot(f(1:N/2),A(1:N/2)); %函式fft返回值的資料結構具有對稱性,因此我們只取前一半
title('幅值頻譜')
xlabel('頻率(Hz)')
ylabel('幅值')
subplot(212);plot(f(1:N/2),P(1:N/2));
title('相位譜頻')
xlabel('頻率(Hz)')
ylabel('相位')
