1、第一：頻譜一.調用方法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 

2、   5.0000              4.7782 + 7.7071i         0 + 5.0000i -10.7782 - 6.2929iXk與xn的維數相同，共有8個元素。Xk的第一個數對應于直流分量，即頻率值為0。（2）做FFT分析時，幅值大小與FFT選擇的點數有關，但不影響分析結果。在IFFT時已經做了處理。要得到真實的振幅值的大小，只要將得到的變換后結果乘以2除以N即可。二.FFT應用舉例例1：x=0

3、.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;&

4、#160;    %頻率序列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;%對信號

5、采樣數據為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

6、,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對信號做譜分析，只需考察0Nyquist頻

7、率范圍內的福頻特性。若沒有給出采樣頻率和采樣間隔，則分析通常對歸一化頻率01進行。另外，振幅的大小與所用采樣點數有關，采用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=

8、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頻率之

9、前的振幅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

10、(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)*f

11、s/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*

12、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）由于在時間域內信號加零，致使振

13、幅譜中出現很多其他成分，這是加零造成的。其振幅由于加了多個零而明顯減小。（3）FFT程序將數據截斷，這時分辨率較高。（4）也是在數據的末尾補零，但由于含有信號的數據個數足夠多，FFT振幅譜也基本不受影響。      對信號進行頻譜分析時，數據樣本應有足夠的長度，一般FFT程序中所用數據點數與原含有信號數據點數相同，這樣的頻譜圖具有較高的質量，可減小因補零或截斷而產生的影響。例3：x=cos(2*pi*0.24*n)+cos(2*pi*0.26*n)（1）數據點過少，幾乎無法看出有關信號頻譜的詳細信息；（2）中間的圖是將x(n)補90個零，幅度頻譜的數據相

14、當密，稱為高密度頻譜圖。但從圖中很難看出信號的頻譜成分。（3）信號的有效數據很長，可以清楚地看出信號的頻率成分，一個是0.24Hz，一個是0.26Hz，稱為高分辨率頻譜。         可見，采樣數據過少，運用FFT變換不能分辨出其中的頻率成分。添加零后可增加頻譜中的數據個數，譜的密度增高了，但仍不能分辨其中的頻率成分，即譜的分辨率沒有提高。只有數據點數足夠多時才能分辨其中的頻率成分。第二： 相譜（相位譜和頻率普是回事兒，想著把頻譜中的幅值部分換成相角就可以了）  由于沒有找到具體的理論，就舉幾個例子說明一下。&

15、#160;   比如要求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/p

16、i;是利用angle函數求出每個點的角度，并由弧度轉化成角度！angle函數解釋：Phase angle SyntaxP = angle(Z)DescriptionP = 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 byR = abs(Z)theta = angle(Z)and

17、the statementZ = R.*exp(i*theta)converts back to the original complex Z.ExamplesZ = 1 - 1i   2 + 1i   3 - 1i   4 + 1i      1 + 2i   2 - 2i   3 + 2i   4 - 2i      1 - 3i   2 + 3i   

18、3 - 3i   4 + 3iP = 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

19、   -1.1071    0.9273   -0.7854AlgorithmsThe angle function can be expressed as angle(z) = imag(log(z) = atan2(imag(z),real(z). 第三：功率譜matlab實現經典功率譜估計fft做出來是頻譜，psd做出來是功率譜；功率譜丟失了頻譜的相位信息；頻譜不同的信號其功率譜是可能相同的；功率譜是幅度取模后平方，結果是個實數matlab中自功率譜密度直接用psd函數就可以求，按照matlab的說法，psd能實現We

20、lch法估計，即相當于用改進的平均周期圖法來求取隨機信號的功率譜密度估計。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

21、=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(

22、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=102

23、4;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_Px

24、xc);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窗no

25、verlap=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_Px

26、x);pause;figure(2)plot(f,plot_Pxx1);pause;figure(3)plot(f,plot_Pxx2); 第四： 相關性分析1. 首先說說自相關和互相關的概念。    這個是信號分析里的概念，他們分別表示的是兩個時間序列之間和同一個時間序列在任意兩個不同時刻的取值之間的相關程度，即互相關函數是描述隨機信號x(t),y(t)在任意兩個不同時刻t1，t2的取值之間的相關程度，自相關函數是描述隨機信號x(t)在任意兩個不同時刻t1，t2的取值之間的相關程度。    自相關函數是描述隨機信號X(t)在任意兩個不同時刻t

27、1，t2的取值之間的相關程度；互相關函數給出了在頻域內兩個信號是否相關的一個判斷指標，把兩測點之間信號的互譜與各自的自譜聯系了起來。它能用來確定輸出信號有多大程度來自輸入信號，對修正測量中接入噪聲源而產生的誤差非常有效.    事實上，在圖象處理中，自相關和互相關函數的定義如下：設原函數是f(t)，則自相關函數定義為R(u)=f(t)*f(-t)，其中*表示卷積；設兩個函數分別是f(t)和g(t)，則互相關函數定義為R(u)=f(t)*g(-t)，它反映的是兩個函數在不同的相對位置上互相匹配的程度。那么，如何在matlab中實現這兩個相關并用圖像顯示出來呢？dt=.1;t=

28、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)，其中×表示乘法，注：此公式僅表示形式

29、計算，并非實際計算所用的公式。當然也可以直接采用卷積進行計算，但是結果會與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。3. 其他相關問題：（1）相關程度與相關函數的取值有什么聯系？    相關系數只是一個比率，不是等單位量度，無什么單位名稱，也不是相關的百分數，一般取小數點后兩位來表示。相關系數的正負號只表示相關的方向，絕對值表示相關的程度。因為不是等單位的度量，因而不能說相關系數0.7是0.35兩倍，只能說相關系數為0.7的二列變量相關程度比相