（此帖引至網絡資源，僅供參考學習）
第一：頻譜
一.調用方法
X=FFT(x)；
X=FFT(x，N)；
x=IFFT(X);
x=IFFT(X,N)
用MATLAB進行譜分析時注意：
（1）函數FFT返回值的數據結構具有對稱性。
例：
N=8;
n=0:N-1;
xn=[4 3 2 6 7 8 9 0];
Xk=fft(xn)
→
Xk =
39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i
Xk與xn的維數相同，共有8個元素。Xk的第一個數對應於直流分量，即頻率值為0。
（2）做FFT分析時，幅值大小與FFT選擇的點數有關，但不影響分析結果。在IFFT時已經做了處理。要得到真實的振幅值的大小，只要將得到的變換后結果乘以2除以N即可。
二.FFT應用舉例
例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采樣頻率fs=100Hz，分別繪制N=128、1024點幅頻圖。
clf;
fs=100;N=128; %采樣頻率和數據點數
n=0:N-1;t=n/fs; %時間序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信號
y=fft(x,N); %對信號進行快速Fourier變換
mag=abs(y); %求得Fourier變換后的振幅
f=n*fs/N; %頻率序列
subplot(2,2,1),plot(f,mag); %繪出隨頻率變化的振幅
xlabel('頻率/Hz');
ylabel('振幅');title('N=128');grid on;
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %繪出Nyquist頻率之前隨頻率變化的振幅
xlabel('頻率/Hz');
ylabel('振幅');title('N=128');grid on;
%對信號采樣數據為1024點的處理
fs=100;N=1024;n=0:N-1;t=n/fs;
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信號
y=fft(x,N); %對信號進行快速Fourier變換
mag=abs(y); %求取Fourier變換的振幅
f=n*fs/N;
subplot(2,2,3),plot(f,mag); %繪出隨頻率變化的振幅
xlabel('頻率/Hz');
ylabel('振幅');title('N=1024');grid on;
subplot(2,2,4)
plot(f(1:N/2),mag(1:N/2)); %繪出Nyquist頻率之前隨頻率變化的振幅
xlabel('頻率/Hz');
ylabel('振幅');title('N=1024');grid on;
fs=100Hz，Nyquist頻率為fs/2=50Hz。整個頻譜圖是以Nyquist頻率為對稱軸的。並且可以明顯識別出信號中含有兩種頻率成分：15Hz和40Hz。由此可以知道FFT變換數據的對稱性。因此用FFT對信號做譜分析，只需考察0~Nyquist頻率范圍內的福頻特性。若沒有給出采樣頻率和采樣間隔，則分析通常對歸一化頻率0~1進行。另外，振幅的大小與所用采樣點數有關，采用128點和1024點的相同頻率的振幅是有不同的表現值，但在同一幅圖中，40Hz與15Hz振動幅值之比均為4：1，與真實振幅0.5：2是一致的。為了與真實振幅對應，需要將變換后結果乘以2除以N。
例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,繪制：
（1）數據個數N=32，FFT所用的采樣點數NFFT=32；
（2）N=32，NFFT=128；
（3）N=136，NFFT=128；
（4）N=136，NFFT=512。
clf;fs=100; %采樣頻率
Ndata=32; %數據長度
N=32; %FFT的數據長度
n=0:Ndata-1;t=n/fs; %數據對應的時間序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %時間域信號
y=fft(x,N); %信號的Fourier變換
mag=abs(y); %求取振幅
f=(0:N-1)*fs/N; %真實頻率
subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率之前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=32');grid on;
Ndata=32; %數據個數
N=128; %FFT采用的數據長度
n=0:Ndata-1;t=n/fs; %時間序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真實頻率
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率之前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=128');grid on;
Ndata=136; %數據個數
N=128; %FFT采用的數據個數
n=0:Ndata-1;t=n/fs; %時間序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真實頻率
subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率之前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=128');grid on;
Ndata=136; %數據個數
N=512; %FFT所用的數據個數
n=0:Ndata-1;t=n/fs; %時間序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真實頻率
subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率之前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;
結論：
（1）當數據個數和FFT采用的數據個數均為32時，頻率分辨率較低，但沒有由於添零而導致的其他頻率成分。
（2）由於在時間域內信號加零，致使振幅譜中出現很多其他成分，這是加零造成的。其振幅由於加了多個零而明顯減小。
（3）FFT程序將數據截斷，這時分辨率較高。
（4）也是在數據的末尾補零，但由於含有信號的數據個數足夠多，FFT振幅譜也基本不受影響。
對信號進行頻譜分析時，數據樣本應有足夠的長度，一般FFT程序中所用數據點數與原含有信號數據點數相同，這樣的頻譜圖具有較高的質量，可減小因補零或截斷而產生的影響。
例3：x=cos(2*pi*0.24*n)+cos(2*pi*0.26*n)
（1）數據點過少，幾乎無法看出有關信號頻譜的詳細信息；
（2）中間的圖是將x(n)補90個零，幅度頻譜的數據相當密，稱為高密度頻譜圖。但從圖中很難看出信號的頻譜成分。
（3）信號的有效數據很長，可以清楚地看出信號的頻率成分，一個是0.24Hz，一個是0.26Hz，稱為高分辨率頻譜。
可見，采樣數據過少，運用FFT變換不能分辨出其中的頻率成分。添加零后可增加頻譜中的數據個數，譜的密度增高了，但仍不能分辨其中的頻率成分，即譜的分辨率沒有提高。只有數據點數足夠多時才能分辨其中的頻率成分。
第二： 相譜
（相位譜和頻率普是回事兒，想着把頻譜中的幅值部分換成相角就可以了）
由於沒有找到具體的理論，就舉幾個例子說明一下。
比如要求y=sin(2*pi*60*t) 的相位譜，
程序如下：
fs=200;N=1024;n=0:N-1;t=n/fs;y=sin(2*pi*60*t);
Y=fft(y,N);
A=abs(Y);f=n*fs/N;
ph=2*angle(Y(1:N/2));
ph=ph*180/pi;
plot(f(1:N/2),ph(1:N/2));
xlabel('頻率/hz'),ylabel('相角'),title('相位譜');
grid on;
期中的 ph=2*angle(Y(1:N/2));ph=ph*180/pi;是利用angle函數求出每個點的角度，並由弧度轉化成角度！
angle函數解釋：
Phase angle
Syntax
P = angle(Z)
Description
P = angle(Z) returns the phase angles, in radians, for each element of complex array Z. The angles lie between ±π.
For complex Z, the magnitude R and phase angle theta are given by
R = abs(Z)
theta = angle(Z)
and the statement
Z = R.*exp(i*theta)
converts back to the original complex Z.
Examples
Z = [ 1 - 1i 2 + 1i 3 - 1i 4 + 1i
1 + 2i 2 - 2i 3 + 2i 4 - 2i
1 - 3i 2 + 3i 3 - 3i 4 + 3i
P = angle(Z)
P =
-0.7854 0.4636 -0.3218 0.2450
1.1071 -0.7854 0.5880 -0.4636
-1.2490 0.9828 -0.7854 0.6435
1.3258 -1.1071 0.9273 -0.7854
Algorithms
The angle function can be expressed as angle(z) = imag(log(z)) = atan2(imag(z),real(z)).
第三：功率譜
matlab實現經典功率譜估計
fft做出來是頻譜，psd做出來是功率譜；功率譜丟失了頻譜的相位信息；頻譜不同的信號其功率譜是可能相同的；功率譜是幅度取模后平方，結果是個實數
matlab中自功率譜密度直接用psd函數就可以求，按照matlab的說法，psd能實現Welch法估計，即相當於用改進的平均周期圖法來求取隨機信號的功率譜密度估計。psd求出的結果應該更光滑吧。
1、直接法：
直接法又稱周期圖法，它是把隨機序列x(n)的N個觀測數據視為一能量有限的序列，直接計算x(n)的離散傅立葉變換，得X(k)，然后再取其幅值的平方，並除以N，作為序列x(n)真實功率譜的估計。
Matlab代碼示例：
clear;
Fs=1000; %采樣頻率
n=0:1/Fs:1;
%產生含有噪聲的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));
2、間接法：
間接法先由序列x(n)估計出自相關函數R(n)，然后對R(n)進行傅立葉變換，便得到x(n)的功率譜估計。
Matlab代碼示例：
clear;
Fs=1000; %采樣頻率
n=0:1/Fs:1;
%產生含有噪聲的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
cxn=xcorr(xn,'unbiased'); %計算序列的自相關函數
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);
3、改進的直接法：對於直接法的功率譜估計，當數據長度N太大時，譜曲線起伏加劇，若N太小，譜的分辨率又不好，因此需要改進。
3.1、Bartlett法
Bartlett平均周期圖的方法是將N點的有限長序列x(n)分段求周期圖再平均。
Matlab代碼示例：
clear；
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(length(n)); %矩形窗
noverlap=0; %數據無重疊
p=0.9; %置信概率
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
figure(1)
plot(k,plot_Pxx);
pause;
figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
3.2、Welch法
Welch法對Bartlett法進行了兩方面的修正，一是選擇適當的窗函數w(n)，並再周期圖計算前直接加進去，加窗的優點是無論什么樣的窗函數均可使譜估計非負。二是在分段時，可使各段之間有重疊，這樣會使方差減小。
Matlab代碼示例：
clear;
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %數據無重疊
range='half'; %頻率間隔為[0 Fs/2]，只計算一半的頻率
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
figure(1)
plot(f,plot_Pxx);
pause;
figure(2)
plot(f,plot_Pxx1);
pause;
figure(3)
plot(f,plot_Pxx2);
第四： 相關性分析
1. 首先說說自相關和互相關的概念。
t=[0:dt:100];
x=cos(t);
[a,b]=xcorr(x,'unbiased');
plot(b*dt,a)
上面代碼是求自相關函數並作圖，對於互相關函數，稍微修改一下就可以了，即把[a,b]=xcorr(x,'unbiased');改為[a,b]=xcorr(x,y,'unbiased');便可。
2. 實現過程：
在Matalb中，求解xcorr的過程事實上是利用Fourier變換中的卷積定理進行的，即R(u)=ifft(fft(f)×fft(g))，其中×表示乘法，注：此公式僅表示形式計算，並非實際計算所用的公式。當然也可以直接采用卷積進行計算，但是結果會與xcorr的不同。事實上，兩者既然有定理保證，那么結果一定是相同的，只是沒有用對公式而已。下面是檢驗兩者結果相同的代碼：
dt=.1;
t=[0:dt:100];
x=3*sin(t);
y=cos(3*t);
subplot(3,1,1);
plot(t,x);
subplot(3,1,2);
plot(t,y);
[a,b]=xcorr(x,y);
subplot(3,1,3);
plot(b*dt,a);
yy=cos(3*fliplr(t)); % or use: yy=fliplr(y);
z=conv(x,yy);
pause;
subplot(3,1,3);
plot(b*dt,z,'r');
即在xcorr中不使用scaling。
（1）相關程度與相關函數的取值有什么聯系？
對於相關系數的大小所表示的意義目前在統計學界尚不一致，但通常按下是這樣認為的：
相關系數 相關程度
0.00-±0.30 微相關
±0.30-±0.50 實相關
±0.50-±0.80 顯著相關
±0.80-±1.00 高度相關
分別用這兩個函數對同一個序列計算，為什么結果不太一樣？因為xcorr是沒有將均值減掉做的相關，autocorr則是減掉了均值的。而且，用離散信號做自相關時，信號截取長度（采樣點N）不一樣，自相關函數就不一樣。
（3）xcorr是計算互相關函數，帶有一個option的參數:
a=xcorr(x,y,'option')
option=baised時，是計算互相關函數的有偏估計；
option=unbaised時，是計算互相關函數的無偏估計；
option=coeff時，是計算歸一化的互相關函數，即為互相關系數，在-1至1之間；
option=none，是缺省的情況。
所以想要計算互相關系數，可用'coeff'參數。
*************************************************************************
用這個xcorr函數作離散互相關運算時要注意，當x, y是不等長向量時，短的向量會自動填0與長的對齊，運算結果是行向量還是列向量就與x一樣。
互相關運算計算的是x,y兩組隨機數據的相關程度，使用參數coeff時，結果就是互相關系數，在-1至1之間，否則結果不一定在這范圍，有可能很大也有可能很小，這視乎x, y數據的大小，所以一般要計算兩組數據的相關程度，一般選擇coeff參數，對結果進行歸一化。
所謂歸一化簡單理解就是將數據系列縮放到-1到1范圍，正式的就是一種簡化計算的方式，即將有量綱的表達式，經過變換，化為無量綱的表達式，成為純量。變換式為X=(X實測--Xmin)/(Xmax-Xmin)。
一般來說選擇歸一化進行互相關運算后，得到結果絕對值越大，兩組數據相關程度就越高。