MATLAB處理信號得到頻譜、相譜、功率譜

1、第一：頻譜一.調(diào)用方法X=FFT(x)；X=FFT(x，N)；x=IFFT(X);x=IFFT(X,N)用MATLAB進(jìn)行譜分析時注意：（1）函數(shù)FFT返回值的數(shù)據(jù)結(jié)構(gòu)具有對稱性。例：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.2929iXk與xn的維數(shù)相同，共有8個元素。Xk的第一個數(shù)對應(yīng)于直流分量，即頻率值為0。（2）做FFT分析時，幅值

2、大小與FFT選擇的點(diǎn)數(shù)有關(guān)，但不影響分析結(jié)果。在IFFT時已經(jīng)做了處理。要得到真實(shí)的振幅值的大小，只要將得到的變換后結(jié)果乘以2除以N即可。二.FFT應(yīng)用舉例例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采樣頻率fs=100Hz，分別繪制N=128、1024點(diǎn)幅頻圖。clf;fs=100;N=128; %采樣頻率和數(shù)據(jù)點(diǎn)數(shù)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); %對信號進(jìn)行快速Fourier變換mag=abs(y); %求得Fourier變換后的振幅f

3、=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;%對信號采樣數(shù)據(jù)為1024點(diǎn)的處理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=ff

4、t(x,N); %對信號進(jì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;運(yùn)行結(jié)果： fs=100Hz，Nyquist頻率為fs/2=50Hz。整個頻譜圖是以Ny

5、quist頻率為對稱軸的。并且可以明顯識別出信號中含有兩種頻率成分：15Hz和40Hz。由此可以知道FFT變換數(shù)據(jù)的對稱性。因此用FFT對信號做譜分析，只需考察0Nyquist頻率范圍內(nèi)的福頻特性。若沒有給出采樣頻率和采樣間隔，則分析通常對歸一化頻率01進(jìn)行。另外，振幅的大小與所用采樣點(diǎn)數(shù)有關(guān)，采用128點(diǎn)和1024點(diǎn)的相同頻率的振幅是有不同的表現(xiàn)值，但在同一幅圖中，40Hz與15Hz振動幅值之比均為4：1，與真實(shí)振幅0.5：2是一致的。為了與真實(shí)振幅對應(yīng)，需要將變換后結(jié)果乘以2除以N。例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,繪制：（

6、1）數(shù)據(jù)個數(shù)N=32，F(xiàn)FT所用的采樣點(diǎn)數(shù)NFFT=32；（2）N=32，NFFT=128；（3）N=136，NFFT=128；（4）N=136，NFFT=512。clf;fs=100; %采樣頻率Ndata=32; %數(shù)據(jù)長度N=32; %FFT的數(shù)據(jù)長度n=0:Ndata-1;t=n/fs; %數(shù)據(jù)對應(yīng)的時間序列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; %真實(shí)頻率subplot(2,2,1),plot(f(1:N/2),mag(

7、1:N/2)*2/N); %繪出Nyquist頻率之前的振幅xlabel(頻率/Hz);ylabel(振幅);title(Ndata=32 Nfft=32);grid on;Ndata=32; %數(shù)據(jù)個數(shù)N=128; %FFT采用的數(shù)據(jù)長度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; %真實(shí)頻率subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率之前的振幅xlabel(頻率/H

8、z);ylabel(振幅);title(Ndata=32 Nfft=128);grid on;Ndata=136; %數(shù)據(jù)個數(shù)N=128; %FFT采用的數(shù)據(jù)個數(shù)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; %真實(shí)頻率subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率之前的振幅xlabel(頻率/Hz);ylabel(振幅);title(Ndata=136 Nfft=128)

9、;grid on;Ndata=136; %數(shù)據(jù)個數(shù)N=512; %FFT所用的數(shù)據(jù)個數(shù)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; %真實(shí)頻率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;結(jié)論：（1）當(dāng)數(shù)據(jù)個數(shù)和FFT采用的數(shù)據(jù)個數(shù)均為32時，頻

10、率分辨率較低，但沒有由于添零而導(dǎo)致的其他頻率成分。（2）由于在時間域內(nèi)信號加零，致使振幅譜中出現(xiàn)很多其他成分，這是加零造成的。其振幅由于加了多個零而明顯減小。（3）FFT程序?qū)?shù)據(jù)截斷，這時分辨率較高。（4）也是在數(shù)據(jù)的末尾補(bǔ)零，但由于含有信號的數(shù)據(jù)個數(shù)足夠多，F(xiàn)FT振幅譜也基本不受影響。 對信號進(jìn)行頻譜分析時，數(shù)據(jù)樣本應(yīng)有足夠的長度，一般FFT程序中所用數(shù)據(jù)點(diǎn)數(shù)與原含有信號數(shù)據(jù)點(diǎn)數(shù)相同，這樣的頻譜圖具有較高的質(zhì)量，可減小因補(bǔ)零或截斷而產(chǎn)生的影響。例3：x=cos(2*pi*0.24*n)+cos(2*pi*0.26*n)（1）數(shù)據(jù)點(diǎn)過少，幾乎無法看出有關(guān)信號頻譜的詳細(xì)信息；（2）中間的圖是將

11、x(n)補(bǔ)90個零，幅度頻譜的數(shù)據(jù)相當(dāng)密，稱為高密度頻譜圖。但從圖中很難看出信號的頻譜成分。（3）信號的有效數(shù)據(jù)很長，可以清楚地看出信號的頻率成分，一個是0.24Hz，一個是0.26Hz，稱為高分辨率頻譜。 可見，采樣數(shù)據(jù)過少，運(yùn)用FFT變換不能分辨出其中的頻率成分。添加零后可增加頻譜中的數(shù)據(jù)個數(shù)，譜的密度增高了，但仍不能分辨其中的頻率成分，即譜的分辨率沒有提高。只有數(shù)據(jù)點(diǎn)數(shù)足夠多時才能分辨其中的頻率成分。第二： 相譜（相位譜和頻率普是回事兒，想著把頻譜中的幅值部分換成相角就可以了）由于沒有找到具體的理論，就舉幾個例子說明一下。 比如要求y=sin(2*pi*60*t)的相位譜，程序如下：fs

12、=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函數(shù)求出每個點(diǎn)的角度，并由弧度轉(zhuǎn)化成角度！angle函數(shù)解釋：Phase angleSyntaxP = angle(Z)DescriptionP = angle(Z

13、) 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 the statementZ = R.*exp(i*theta)converts back to the original complex Z.ExamplesZ = 1 - 1i 2 + 1i 3 - 1i

14、 4 + 1i 1 + 2i 2 - 2i 3 + 2i 4 - 2i 1 - 3i 2 + 3i 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 -1.1071 0.9273 -0.7854AlgorithmsThe angle function can be expressed as angle(z) = imag(log(z) = atan2(imag(z),real(z).第三：功率譜

15、matlab實(shí)現(xiàn)經(jīng)典功率譜估計fft做出來是頻譜，psd做出來是功率譜；功率譜丟失了頻譜的相位信息；頻譜不同的信號其功率譜是可能相同的；功率譜是幅度取模后平方，結(jié)果是個實(shí)數(shù)matlab中自功率譜密度直接用psd函數(shù)就可以求，按照matlab的說法，psd能實(shí)現(xiàn)Welch法估計，即相當(dāng)于用改進(jìn)的平均周期圖法來求取隨機(jī)信號的功率譜密度估計。psd求出的結(jié)果應(yīng)該更光滑吧。1、直接法：直接法又稱周期圖法，它是把隨機(jī)序列x(n)的N個觀測數(shù)據(jù)視為一能量有限的序列，直接計算x(n)的離散傅立葉變換，得X(k)，然后再取其幅值的平方，并除以N，作為序列x(n)真實(shí)功率譜的估計。Matlab代碼示例：clea

16、r;Fs=1000; %采樣頻率n=0:1/Fs:1;%產(chǎn)生含有噪聲的序列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)估計出自相關(guān)函數(shù)R(n)，然后對R(n)進(jìn)行傅立葉變換，便得到x(n)的功率譜估計。Matlab代碼示例：clear;Fs=1000; %采樣頻率n=0:1/Fs:1;%產(chǎn)生含有噪聲的序

17、列xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n);nfft=1024;cxn=xcorr(xn,unbiased); %計算序列的自相關(guān)函數(shù)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、改進(jìn)的直接法：對于直接法的功率譜估計，當(dāng)數(shù)據(jù)長度N太大時，譜曲線起伏加劇，若N太小，譜的分辨率又不好，因此需要改進(jìn)。3.1、Bartlett法Bartlett平均周期

18、圖的方法是將N點(diǎ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; %數(shù)據(jù)無重疊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);plo

19、t_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法進(jìn)行了兩方面的修正，一是選擇適當(dāng)?shù)拇昂瘮?shù)w(n)，并再周期圖計算前直接加進(jìn)去，加窗的優(yōu)點(diǎn)是無論什么樣的窗函數(shù)均可使譜估計非負(fù)。二是在分段時，可使各段之間有重疊，這樣會使方差減小。Matlab代碼示例：clear;Fs=1000;n=0:1/Fs:1;xn=cos(2*pi*40*n)+3*cos

20、(2*pi*100*n)+randn(size(n);nfft=1024;window=boxcar(100); %矩形窗window1=hamming(100); %海明窗window2=blackman(100); %blackman窗noverlap=20; %數(shù)據(jù)無重疊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,noverla

21、p,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);第四： 相關(guān)性分析1. 首先說說自相關(guān)和互相關(guān)的概念。 這個是信號分析里的概念，他們分別表示的是兩個時間序列之間和同一個時間序列在任意兩個不同時刻的取值之間的相關(guān)程度，即互相關(guān)函數(shù)是描述隨機(jī)信號x(t),y(t)在任意兩個不同時刻t1，

22、t2的取值之間的相關(guān)程度，自相關(guān)函數(shù)是描述隨機(jī)信號x(t)在任意兩個不同時刻t1，t2的取值之間的相關(guān)程度。 自相關(guān)函數(shù)是描述隨機(jī)信號X(t)在任意兩個不同時刻t1，t2的取值之間的相關(guān)程度；互相關(guān)函數(shù)給出了在頻域內(nèi)兩個信號是否相關(guān)的一個判斷指標(biāo)，把兩測點(diǎn)之間信號的互譜與各自的自譜聯(lián)系了起來。它能用來確定輸出信號有多大程度來自輸入信號，對修正測量中接入噪聲源而產(chǎn)生的誤差非常有效. 事實(shí)上，在圖象處理中，自相關(guān)和互相關(guān)函數(shù)的定義如下：設(shè)原函數(shù)是f(t)，則自相關(guān)函數(shù)定義為R(u)=f(t)*f(-t)，其中*表示卷積；設(shè)兩個函數(shù)分別是f(t)和g(t)，則互相關(guān)函數(shù)定義為R(u)=f(t)*g(

23、-t)，它反映的是兩個函數(shù)在不同的相對位置上互相匹配的程度。那么，如何在matlab中實(shí)現(xiàn)這兩個相關(guān)并用圖像顯示出來呢？dt=.1;t=0:dt:100;x=cos(t);a,b=xcorr(x,unbiased);plot(b*dt,a)上面代碼是求自相關(guān)函數(shù)并作圖，對于互相關(guān)函數(shù)，稍微修改一下就可以了，即把a(bǔ),b=xcorr(x,unbiased);改為a,b=xcorr(x,y,unbiased);便可。2. 實(shí)現(xiàn)過程： 在Matalb中，求解xcorr的過程事實(shí)上是利用Fourier變換中的卷積定理進(jìn)行的，即R(u)=ifft(fft(f)fft(g)，其中表示乘法，注：此公式僅表示形

24、式計算，并非實(shí)際計算所用的公式。當(dāng)然也可以直接采用卷積進(jìn)行計算，但是結(jié)果會與xcorr的不同。事實(shí)上，兩者既然有定理保證，那么結(jié)果一定是相同的，只是沒有用對公式而已。下面是檢驗兩者結(jié)果相同的代碼：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. 其他相關(guān)問題：（1）相關(guān)程度與相關(guān)函數(shù)的取值有什么聯(lián)系？ 相關(guān)系數(shù)只是一個比率，不是等單位量度，無什么單位名稱，也不是相關(guān)的百分?jǐn)?shù)，一般取小數(shù)點(diǎn)后兩位來表示。相關(guān)系數(shù)的正負(fù)號只表示相關(guān)的方向，絕對值表示相關(guān)的程度。因為不是等單位的度量，因而不能說相關(guān)系數(shù)0.7是0.35兩倍，只能說相關(guān)系數(shù)為0.7的二列變量相關(guān)程度比相關(guān)系數(shù)為0.35的二列變量相關(guān)程度更為密切和更高。也不能說相關(guān)系數(shù)從0.70到0.80與相關(guān)系數(shù)從0.30到0.4