MATLAB與現(xiàn)代信號(hào)處理

第8章參數(shù)化建模 38.1時(shí)間域建模 38.1.1線性預(yù)報(bào)法(AR模型) 38.1.2Prony法（ARMA模型） 48.1.3Steiglity-McBride法(ARMA模型) 48.2頻率域建模 58.2.1模擬濾波器的s域建模 68.2.2數(shù)字濾波器的z域內(nèi)建模 6第9章隨機(jī)信號(hào)分析 89.1隨機(jī)信號(hào)的數(shù)字特征 89.1.1均值、均方值、方差 89.1.2離散隨機(jī)信號(hào) 99.1.3估計(jì) 99.2相關(guān)函數(shù)和協(xié)方差 109.2.1相關(guān)函數(shù)和自協(xié)方差 109.2.2互相關(guān)函數(shù)和互協(xié)方差 119.2.3計(jì)算相關(guān)函數(shù)和協(xié)方差的MATLAB函數(shù) 119.3功率譜估計(jì) 169.3.1功率譜密度 169.3.2周期圖法 169.3.3多窗口法 269.3.4最大熵法（Maxmumentropymethod，MEM法） 289.3.5多信號(hào)分類法 309.4傳遞函數(shù)估計(jì) 319.5相干函數(shù) 33 運(yùn)用功率譜提取地球自由振蕩信息 359.6.1資料 35方法 36結(jié)果 37第10章 數(shù)字信號(hào)處理的幾個(gè)前沿課題 4210.1時(shí)譜（倒譜）分析 4210.2地震觀測(cè)系統(tǒng)的仿真和地面運(yùn)動(dòng)的恢復(fù) 44 小波分析舉例 5710.3.1小波變換與突變 5810.3.2資料 6210.3.3本尼奧夫應(yīng)變資料的小波分析 6311900~2001年中國(guó)大陸地震資料的小波分析 6321500~2001年華北地區(qū)地震資料的小波分析 64第8章參數(shù)化建模由前幾章的討論可見，無論是模擬濾波器還是數(shù)字濾波器，均是根據(jù)輸入信號(hào)和系統(tǒng)的傳遞函數(shù)求得系統(tǒng)的輸出信號(hào)。反過來，能否通過系統(tǒng)的輸入信號(hào)和輸出信號(hào)，或?yàn)V波器的頻率響應(yīng)或脈沖響應(yīng)求得系統(tǒng)的傳遞函數(shù)呢？這在一定的假設(shè)下是可以得到的。這種技術(shù)涉及本章要講的參數(shù)化建模。所謂參數(shù)化建模就是根據(jù)未知系統(tǒng)的某些信息（如脈沖響應(yīng)、頻率響應(yīng)或輸入輸出序列）建立該系統(tǒng)的有理傳遞函數(shù)模型。這項(xiàng)技術(shù)用于求一個(gè)信號(hào)、系統(tǒng)或過程的模型參數(shù)，并廣泛應(yīng)用于語(yǔ)言分析、數(shù)據(jù)壓縮、高分辨譜估計(jì)、信號(hào)處理等領(lǐng)域。參數(shù)化建模分為時(shí)間域建模和頻率域建模兩類，本章將分別介紹。8.1時(shí)間域建模時(shí)間域建模就是給定系統(tǒng)時(shí)間域內(nèi)的信息，如脈沖響應(yīng)或系統(tǒng)輸入和輸出時(shí)間序列，求數(shù)字濾波器有理傳遞函數(shù)分子和分母多項(xiàng)式系數(shù)。MATLAB中提供了三種時(shí)間域建模函數(shù)，下面分別敘述。8.1.1線性預(yù)報(bào)法(AR模型)如果一個(gè)信號(hào)的各個(gè)采樣值不是獨(dú)立的，其任一時(shí)刻k的采樣值x(k)是它過去的n個(gè)采樣值及一個(gè)白噪聲序列k時(shí)刻的值u(k)線性組合而成（線性預(yù)報(bào)）,即x(k)=-a(2)x(k-1)-a(3)x(k-2)…-a(n)x(k-n-1)-a(n+1)x(k-n)+u(k)（8-1）如果知道參數(shù)個(gè)數(shù)n，根據(jù)這個(gè)信號(hào)序列可以得到這n個(gè)系數(shù)的值。數(shù)字信號(hào)處理中稱這個(gè)信號(hào)x符合n階自回歸(automaticregression,AR)線性模型，它可用一個(gè)全極點(diǎn)IIR濾波器脈沖響應(yīng)來逼近，全極點(diǎn)濾波器的傳遞函數(shù)具有如下形式：（8-2）MATLAB信號(hào)處理工具箱函數(shù)lpc利用自回歸(AR)模型的相關(guān)法來求出AR模型系數(shù)a，也是全極點(diǎn)濾波器模型。其調(diào)用格式為：[a,g]=lpc(x,n)其中，x為信號(hào)，是一個(gè)實(shí)時(shí)間序列；n為AR模型階次；a為AR模型系數(shù)，a=[1,a(2),…,a(n+1);g為AR模型的增益。這就是說，要建立起AR模型，必須知道其階次。函數(shù)lpc的算法可簡(jiǎn)要地分為幾個(gè)步驟：（1）調(diào)用函數(shù)xcorr求信號(hào)x的相關(guān)估計(jì)。（2）調(diào)用函數(shù)levinson實(shí)現(xiàn)Levinson-Durbin遞推算法，求出系數(shù)a。函數(shù)levinson的調(diào)用格式為：a=levinson(r,n)其中，r為x的自相關(guān)序列；n為AR模型階次。【例8-1】設(shè)信號(hào)x是一個(gè)帶白噪聲的4階IIR濾波器的脈沖響應(yīng)，IIR濾波器的傳遞函數(shù)無零點(diǎn)，其分母多項(xiàng)式系數(shù)為a=[10.10.10.10.1]，用線性預(yù)報(bào)法建立其AR模型。%Samp8_1randn('state',0);%設(shè)置隨機(jī)函數(shù)的狀態(tài)a=[10.10.10.10.1];%濾波器分母多項(xiàng)式系數(shù)x=impz(1,a,10)+randn(10,1)/20;%求得帶噪聲的濾波器脈沖響應(yīng)%采用第一種方法r=xcorr(x);%求信號(hào)x的自相關(guān)r，參考第9章的該函數(shù)用法。r(1:length(x)-1)=[];%該序列的前面部分(相當(dāng)于x的長(zhǎng)度-1)受到邊界的影響，扣除aa=levinson(r,4)%用Levison-Durbin遞推算法采用自相關(guān)序列求出系數(shù)aa%采用第二種方法a1=lpc(x,4)%直接調(diào)用lpc（LinearPredictorCoefficients）函數(shù)求得結(jié)果該程序的運(yùn)行結(jié)果為：程序中兩種方法求得的結(jié)果完全相同。若信號(hào)不包含隨機(jī)噪聲，用上面的程序求得信號(hào)AR模型和所采用全極點(diǎn)濾波器模型完全相同。即使加入隨機(jī)噪聲，得到結(jié)果也與原模型參數(shù)a接近。這說明，如果一個(gè)信號(hào)x是IIR濾波器的脈沖響應(yīng)，可以利用函數(shù)lpc建立的向量a作為濾波器參數(shù)化模型。8.1.2Prony法（ARMA模型）如果已知某濾波器（或系統(tǒng)）的脈沖響應(yīng)序列，可以采用Prony方法求得該濾波器（IIR）傳遞函數(shù)的系數(shù)。MATLAB信號(hào)處理工具箱提供了Prony法建模的函數(shù)prony,其調(diào)用格式為：[b,a]=prony(h,nb,na)式中，h為時(shí)間域脈沖響應(yīng)序列；nb,na分別為濾波器分子和分母多項(xiàng)式階數(shù)；b,a分別為濾波器分子和分母多項(xiàng)式系數(shù)向量。這就是說，要想建立起ARMA模型，必須知道傳遞函數(shù)分子和分母多項(xiàng)式的最高階次。濾波器傳遞函數(shù)表示為：（8-3）Prony法是用一個(gè)IIR濾波器脈沖響應(yīng)逼近信號(hào)h。【例8-2】建立一個(gè)4階Butterworth數(shù)字濾波器，歸一化截止頻率為0.2，求其脈沖響應(yīng)，再用脈沖響應(yīng)系數(shù)和prony函數(shù)求濾波器模型,并與原始模型進(jìn)行比較。%Samp8_2[b,a]=butter(4,0.2)%設(shè)計(jì)4階截止頻率為0.2(歸一化)的Butterworth濾波器h=impz(b,a,26);%求26個(gè)點(diǎn)組成的濾波器的脈沖響應(yīng)[bb,aa]=prony(h,4,4)%運(yùn)用脈沖響應(yīng)h建立4階系統(tǒng)的分子和分母多項(xiàng)式系數(shù)運(yùn)行結(jié)果為：aa=程序運(yùn)行結(jié)果表明，用函數(shù)prony建模和原始濾波器模型完全相同，顯然函數(shù)prony所建立的新濾波器的脈沖響應(yīng)和信號(hào)h是完全吻合的。因此，若已知一個(gè)系統(tǒng)的脈沖響應(yīng)和該系統(tǒng)分子分母多項(xiàng)式階數(shù)，可以運(yùn)用該方法求得該系統(tǒng)的傳遞函數(shù)。8.1.3Steiglity-McBride法(ARMA模型)已知某濾波器（或系統(tǒng)）的脈沖響應(yīng)序列建立濾波器傳遞函數(shù)的另一種方法是Steiglity-McBride法。這種方法利用Steiglity-McBride迭代建模，使系統(tǒng)b(z)/a(z)的脈沖響應(yīng)x’和輸出信號(hào)x的均方差最小，即（8-4）MATLAB信號(hào)處理工具箱函數(shù)stmcb實(shí)現(xiàn)這種方法建模，可用于濾波器設(shè)計(jì)和系統(tǒng)辨識(shí)。其調(diào)用格式為：[b,a]=stmcb(x,nb,na[,niter[,ai]])式中,x為系統(tǒng)的脈沖響應(yīng)；nb,na分別為系統(tǒng)傳遞函數(shù)b(z)/a(z)分子和分母多項(xiàng)式的階次；niter為迭代計(jì)算次數(shù)；缺省值為5；ai為分母系數(shù)的初步估計(jì)，可缺省；b,a分別為系統(tǒng)（IIR型濾波器）分子和分母多項(xiàng)式系數(shù)向量。系統(tǒng)(IIR濾波器)傳遞函數(shù)具有與（8-3）式相同的形式。函數(shù)的另一種調(diào)用形式適用于系統(tǒng)的輸入和輸出信號(hào)作為函數(shù)輸入：[b,a]=stmcb(x,u,nb,na[,niter[,ai]])式中，x和u分別為系統(tǒng)（濾波器）輸出和輸入信號(hào)，其余參數(shù)同上。函數(shù)的這種用法可以根據(jù)濾波器的輸入和輸出（而不是前面的脈沖響應(yīng)）求得濾波器傳遞函數(shù)?！纠?-3】建立一個(gè)6階Butterworth濾波器，截止頻率為0.2，求其脈沖響應(yīng)，并用脈沖響應(yīng)數(shù)據(jù)和函數(shù)stmcb建立新濾波器模型，比較兩個(gè)濾波器的頻率特性。%Samp8_3h=filter(b,a,[1zeros(1,100)]);%采用輸入為脈沖函數(shù)的方法求得系統(tǒng)的脈沖響應(yīng)[bb,aa]=stmcb(h,6,6)%采用stmcb方法求得系統(tǒng)傳遞函數(shù)多項(xiàng)式系數(shù)運(yùn)行結(jié)果為：結(jié)果表明，函數(shù)stmcb自濾波器脈沖響應(yīng)建立的模型和原濾波器模型完全相同。為了比較lpc，prony和stmcb這三個(gè)函數(shù),下例給出了三種方法建模的誤差【例8-4】采用例8-1中的濾波器，用階數(shù)3模擬比較函數(shù)lpc、prony和stmcb三種建模方法的誤差。%Samp8_4randn('seed',0);%設(shè)置隨機(jī)數(shù)種子x=impz(1,[10.10.10.1],10)+randn(10,1)/20;%a1=lpc(x,3);%運(yùn)用lpc函數(shù)建立傳遞函數(shù)系數(shù)[b2,a2]=prony(x,3,3);%運(yùn)用prony函數(shù)建立傳遞函數(shù)系數(shù)[b3,a3]=stmcb(x,3,3);%運(yùn)用stmcb函數(shù)建立傳遞函數(shù)系數(shù)%計(jì)算誤差disp('函數(shù)LPC輸出:')err1=x-impz(1,a1,10);%x與lpc函數(shù)確定模型脈沖響應(yīng)的每個(gè)數(shù)的誤差sumerr1=sum(err1.^2)%x與lpc函數(shù)確定模型脈沖響應(yīng)的總誤差disp('函數(shù)PRONY輸出:')err2=x-impz(b2,a2,10);%x與prony函數(shù)確定模型脈沖響應(yīng)的每個(gè)數(shù)誤差sumerr2=sum(err2.^2)%x與prony函數(shù)確定模型脈沖響應(yīng)的總誤差disp('函數(shù)STMCB輸出:')err3=x-impz(b3,a3,10);%x與stmcb函數(shù)確定模型脈沖響應(yīng)的每個(gè)數(shù)誤差sumerr3=sum(err3.^2)%x與stmcbc函數(shù)確定模型脈沖響應(yīng)的總誤差程序輸出結(jié)果為：函數(shù)LPC輸出:sumerr1=0.0219;函數(shù)PRONY輸出:sumerr2=0.0147;該例表明：在相同的階數(shù)條件下，函數(shù)stmcb精度最好，prony次之，lpc較差。8.2頻率域建模前面一節(jié)我們根據(jù)系統(tǒng)（濾波器）的時(shí)間域信號(hào)（脈沖響應(yīng)或輸入、輸出信號(hào)）建立系統(tǒng)傳遞函數(shù)。但如果已知某系統(tǒng)（濾波器）的頻率響應(yīng)，能否建立起該系統(tǒng)（濾波器）的模型呢？這就涉及到頻率域建模。頻率域建模就是由濾波器（系統(tǒng)）的頻率響應(yīng)建立其模型。根據(jù)頻率域類型分為面向模擬濾波器的s域內(nèi)建模和面向數(shù)字濾波器的z域內(nèi)建模。8.2.1模擬濾波器的s域建模模擬濾波器傳遞函數(shù)的一般形式為：（8-5）MATLAB信號(hào)處理工具箱函數(shù)invfreqs用于由給定的頻率響應(yīng)數(shù)據(jù)求形如上式的模擬濾波器傳遞函數(shù)。調(diào)用格式為：[b,a]=invfreqs(h,w,nb,na[,wt[niter]])式中，h為復(fù)頻率響應(yīng)向量；w為對(duì)應(yīng)的頻率向量；nb,na分別為模擬濾波器分子和分母多項(xiàng)式階次；wt為與w同維數(shù)的權(quán)系數(shù)向量，調(diào)整對(duì)頻率的擬合誤差；niter為迭代次數(shù)。若已知復(fù)頻率響應(yīng)的幅值向量mag和相位向量phase，在MATLAB中要采用h=mag.*exp(j*phase)將其復(fù)合為復(fù)數(shù)形式?！纠?-5】已知濾波器原始傳遞函數(shù)中b=[1232],a=[12321],試求其頻率響應(yīng)，并利用頻率響應(yīng)數(shù)據(jù)求原始濾波器的傳遞函數(shù)。%Samp8_5disp('濾波器初始模型:');b=[1232]a=[12321][h,w]=freqs(b,a,64);%計(jì)算64個(gè)點(diǎn)的模擬濾波器頻率響應(yīng)disp('INVFREQS函數(shù)得到的模型:')[bb,aa]=invfreqs(h,w,3,4)%運(yùn)用invfreqs函數(shù)求濾波器系數(shù)程序的運(yùn)行結(jié)果為：濾波器初始模型:b=1232a=12321INVFREQS函數(shù)得到的模型aa=1.00002.00003.00002.00001.0000運(yùn)行結(jié)果表明，由函數(shù)invfreqs從頻率響應(yīng)數(shù)據(jù)辨識(shí)的濾波器模型和原始模型完全相同。8.2.2數(shù)字濾波器的z域內(nèi)建模數(shù)字濾波器傳遞函數(shù)的一般形式為:（8-6）MATLAB信號(hào)處理工具箱函數(shù)invfreqz用于由給定的頻率響應(yīng)數(shù)據(jù)求形如上式的數(shù)字濾波器的傳遞函數(shù)。調(diào)用格式為：[b,a]=invfreqz(h,w,nb,na[,wt[,niter]])式中，h為復(fù)頻率響應(yīng)向量；w為對(duì)應(yīng)的頻率向量；nb,na分別為數(shù)字濾波器分子和分母多項(xiàng)式階次；wt為與w同維數(shù)的權(quán)系數(shù)向量，調(diào)整對(duì)頻率的擬合誤差；niter為迭代次數(shù)。函數(shù)invfreqs和invfreqz利用等誤差方法由給定的頻率響應(yīng)數(shù)據(jù)來辨識(shí)最佳模型，當(dāng)niter缺省時(shí)，函數(shù)利用下面的算法：（8-7）式中，和分別表示多項(xiàng)式a,b在頻率w(k)處的Fourier變換，為權(quán)系數(shù)向量，n為頻率點(diǎn)數(shù)（h和w的長(zhǎng)度）。這種算法不能保證辨識(shí)模型的穩(wěn)定性。當(dāng)?shù)螖?shù)niter確定后，函數(shù)用輸出誤差迭代算法進(jìn)行運(yùn)算：（8-8）這樣可保證辨識(shí)模型的穩(wěn)定性，且提高辨識(shí)模型的精度?！纠?-6】已知原始濾波器為4階Butterworth低通數(shù)字濾波器，歸一化截止頻率為0.4。用其頻率響應(yīng)數(shù)據(jù)產(chǎn)生一個(gè)三階濾波器。試比較采用等誤差和輸出誤差迭代算法的結(jié)果。%Samp8_6[b,a]=butter(4,0.4);%產(chǎn)生4階Butterworth濾波器[h,w]=freqz(b,a,64);%計(jì)算頻率響應(yīng)[bb,aa]=invfreqz(h,w,3,3);%采用invfreqz求解濾波器的系數(shù)wt=ones(size(w));%產(chǎn)生權(quán)向量niter=30;%迭代次數(shù)[bbb,aaa]=invfreqz(h,w,3,3,wt,niter);%采用權(quán)向量和迭代次數(shù)求濾波器系數(shù)[h1,w1]=freqz(bb,aa,w);%計(jì)算模型1的頻率響應(yīng)[h2,w2]=freqz(bbb,aaa,w);%計(jì)算模型2的頻率響應(yīng)plot(w/pi,abs(h),'k',w1/pi,abs(h1),'k:',w2/pi,abs(h2),'b.');%繪制頻率響應(yīng)xlabel('歸一化頻率');ylabel('振幅')legend('理想模型','未采用權(quán)向量和迭代次數(shù)','采用權(quán)向量和迭代次數(shù)');gridonsumerr1=sum(abs(h-h1).^2)%輸出模型1誤差sumerr2=sum(abs(h-h2).^2)%輸出模型2誤差圖8-1例8-6的兩種算法估計(jì)模型的幅頻特性的比較程序命令窗口的輸出為：sumerr1=0.0200；sumerr2=0.0096。程序運(yùn)行的圖形顯示為圖8-1。程序運(yùn)行結(jié)果表明，輸出誤差迭代法（預(yù)估算法）具有較好的擬合精度。第八章習(xí)題用MATLAB函數(shù)產(chǎn)生一個(gè)5階低通ChebyshevI型數(shù)字濾波器（歸一化截止頻率為0.4）的脈沖響應(yīng)序列，再用Prony法建模并與原濾波器模型進(jìn)行比較。用MATLAB函數(shù)產(chǎn)生一個(gè)5階低通Butterworth數(shù)字濾波器（歸一化截止頻率為0.4）并計(jì)算其脈沖響應(yīng)序列，再用Steiglitz-McBride法建模并與原濾波器模型進(jìn)行比較。用MATLAB函數(shù)產(chǎn)生一個(gè)低通ChebyshevI型高通數(shù)字濾波器，設(shè)計(jì)指標(biāo)為：通帶截止頻率Fc=2.5kHz,通帶波紋為1dB，阻帶衰減大于20dB。試求（1）高通數(shù)字濾波器的系統(tǒng)函數(shù)；（2）濾波器的頻率響應(yīng)；（3）由頻率響應(yīng)數(shù)據(jù)辨識(shí)濾波器模型并與原濾波器模型進(jìn)行比較。第9章隨機(jī)信號(hào)分析隨機(jī)信號(hào)和確定信號(hào)是兩類性質(zhì)完全不同的信號(hào)，對(duì)它們的描述、分析和處理方法也不相同。隨機(jī)信號(hào)是一種不能用確定數(shù)學(xué)關(guān)系式來描述的，無法預(yù)測(cè)未來某時(shí)刻精確值的信號(hào)，也無法用實(shí)驗(yàn)的方法重復(fù)再現(xiàn)。隨機(jī)信號(hào)分為平穩(wěn)和非平穩(wěn)兩類。平穩(wěn)隨機(jī)信號(hào)又分為各態(tài)歷經(jīng)和非各態(tài)歷經(jīng)。本章所討論的隨機(jī)信號(hào)是平穩(wěn)的且是各態(tài)歷經(jīng)的。在研究無限長(zhǎng)信號(hào)時(shí)，總是取某段有限長(zhǎng)信號(hào)作分析。這一有限長(zhǎng)信號(hào)稱為一個(gè)樣本（或稱子集），而無限長(zhǎng)信號(hào)x(t)稱為隨機(jī)信號(hào)總體（或稱集）。各態(tài)歷經(jīng)的平穩(wěn)隨機(jī)過程中的一個(gè)樣本的時(shí)間均值和集的平均值相等。因此一個(gè)樣本的統(tǒng)計(jì)特征代表了隨機(jī)信號(hào)總體，這使得研究大大簡(jiǎn)化。工程上的隨機(jī)信號(hào)一般均按各態(tài)歷經(jīng)平穩(wěn)隨機(jī)過程來處理。僅在離散時(shí)間點(diǎn)上給出定義的隨機(jī)信號(hào)稱為離散時(shí)間隨機(jī)信號(hào)，即隨機(jī)信號(hào)序列。隨機(jī)信號(hào)序列可以是連續(xù)隨機(jī)信號(hào)的采樣結(jié)果，也可以是自然界里實(shí)際存在的物理現(xiàn)象，即它們本身就是離散的。平穩(wěn)隨機(jī)過程在時(shí)間上是無始無終的，即其能量是無限的，本身的Fourier變換也是不存在的；但功率是有限的。通常用功率譜密度來描述隨機(jī)信號(hào)的頻域特征，這是一個(gè)統(tǒng)計(jì)平均的頻譜特性。平穩(wěn)隨機(jī)過程統(tǒng)計(jì)特征的計(jì)算要求信號(hào)x(n)無限長(zhǎng)，而實(shí)際上這是不可能的，只能用一個(gè)樣本，即有限長(zhǎng)序列來計(jì)算。因此得到的計(jì)算值不是隨機(jī)信號(hào)真正的統(tǒng)計(jì)值，而僅僅是一種估計(jì)。本章首先介紹隨機(jī)信號(hào)的數(shù)字特征，旨在使大家熟悉描述隨機(jī)信號(hào)的常用特征量。然后介紹描述信號(hào)之間關(guān)系的相關(guān)函數(shù)和協(xié)方差。這些是數(shù)字信號(hào)時(shí)間域內(nèi)的描述。在頻率域內(nèi)，本章介紹功率譜及其估計(jì)方法，并給出了功率譜在傳遞函數(shù)估計(jì)方面的應(yīng)用。最后介紹描述頻率域信號(hào)之間關(guān)系的函數(shù)相干函數(shù)。9.1隨機(jī)信號(hào)的數(shù)字特征9.1.1均值、均方值、方差若連續(xù)隨機(jī)信號(hào)x(t)是各態(tài)歷經(jīng)的，則隨機(jī)信號(hào)x(t)均值可表示為：（9-1）均值描述了隨機(jī)信號(hào)的靜態(tài)（直流）分量，它不隨時(shí)間而變化。隨機(jī)信號(hào)x(t)的均方值表達(dá)式為：（9-2）均方值表示了信號(hào)的強(qiáng)度或功率。隨機(jī)信號(hào)x(t)的均方根值表示為：（9-3）也是信號(hào)能量或強(qiáng)度的一種描述。隨機(jī)信號(hào)x(t)的方差表達(dá)式為：（9-4）是信號(hào)幅值相對(duì)于均值分散程度的一種表示，也是信號(hào)純波動(dòng)（交流）分量大小的反映。隨機(jī)信號(hào)x(t)的均方差（標(biāo)準(zhǔn)差）可表示為：（9-5）其意義與相同。9.1.2離散隨機(jī)信號(hào)若x(n)是離散的各態(tài)歷經(jīng)的平穩(wěn)隨機(jī)信號(hào)序列。類似連續(xù)隨機(jī)信號(hào)，其數(shù)字特征可由下面式子來表示：均值：（9-6）均方值：（9-7）方差：（9-8）式中，為均方根值，為均方差值。9.1.3估計(jì)以上計(jì)算都是對(duì)無限長(zhǎng)信號(hào)而言，而工程實(shí)際中所取得信號(hào)是有限長(zhǎng)的，計(jì)算中均無法取或。對(duì)于有限長(zhǎng)模擬隨機(jī)信號(hào)，計(jì)算均值時(shí)將（9-1）式改寫為：（9-9）式中，均值僅僅是對(duì)的估計(jì)。當(dāng)T足夠長(zhǎng)時(shí)，均值估計(jì)能精確逼近真實(shí)均值。對(duì)于周期信號(hào)，常取T為一個(gè)周期，估計(jì)均值就能完全代表真實(shí)均值。對(duì)于有限長(zhǎng)隨機(jī)信號(hào)序列，計(jì)算其均值估計(jì)可將（9-6）式改寫為：（9-10）當(dāng)序列長(zhǎng)度足夠長(zhǎng)時(shí)，均值估計(jì)能精確逼近真實(shí)均值。類似地，可以寫出均方值和方差估計(jì)表達(dá)式。在MATLAB工具箱中，沒有專門函數(shù)來計(jì)算均值、均方值和方差。但隨機(jī)信號(hào)的統(tǒng)計(jì)數(shù)字特征值計(jì)算都可以通過MATLAB編程實(shí)現(xiàn)。在數(shù)值計(jì)算中，常將連續(xù)隨機(jī)信號(hào)離散化，當(dāng)作隨機(jī)序列來處理?！纠?-1】計(jì)算以長(zhǎng)度N=100000的正態(tài)分布高斯隨機(jī)信號(hào)的均值、均方值、均方根值、方差和均方差。%Samp9_1N=100000;%數(shù)據(jù)個(gè)數(shù)randn('state',0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的狀態(tài)y=randn(1,N);%產(chǎn)生一個(gè)隨機(jī)序列disp('平均值:');yMean=mean(y)%計(jì)算隨機(jī)序列的均值disp('均方值:');y2p=y*y'/N%計(jì)算其均方值,這里利用了矩陣相乘的算法disp('均方根:');ysq=sqrt(y2p)%計(jì)算其均方根值disp('標(biāo)準(zhǔn)差:');ystd=std(y,1)%相當(dāng)于ystd=sqrt(sum((y-yMean).^2)/(N-1))disp('方差:');yd=ystd.*ystd該程序的運(yùn)行結(jié)果為：均方根標(biāo)準(zhǔn)差方差注意，函數(shù)s=std(x,flag)計(jì)算標(biāo)準(zhǔn)差。x為向量或矩陣；s為標(biāo)準(zhǔn)差；flag為控制符，用來控制標(biāo)準(zhǔn)算法。當(dāng)flag=0(或缺省)時(shí)，按下式計(jì)算無偏標(biāo)準(zhǔn)差：（9-11）當(dāng)flag=1時(shí)，按下式計(jì)算有偏標(biāo)準(zhǔn)差：（9-12）9.2相關(guān)函數(shù)和協(xié)方差9.2.1相關(guān)函數(shù)和自協(xié)方差對(duì)于隨機(jī)信號(hào)x(t)，自相關(guān)函數(shù)為：（9-13）式中，為時(shí)移。若去掉x(t)的均值部分，則相應(yīng)的自相關(guān)函數(shù)稱為自協(xié)方差，即：（9-14）（9-15）對(duì)于離散隨機(jī)信號(hào)序列，x(n)的自相關(guān)函數(shù)和自協(xié)方差為：（9-16）（9-17）式中，m為延遲。9.2.2互相關(guān)函數(shù)和互協(xié)方差對(duì)于兩個(gè)不同隨機(jī)信號(hào)x(t),y(t)，互相關(guān)函數(shù)為：（9-18）式中，為時(shí)移?；f(xié)方差為：（9-19）對(duì)于互協(xié)方差有：（9-20）對(duì)于離散隨機(jī)序列x(n)和y(n)，互相關(guān)函數(shù)和互協(xié)方差為：（9-21）（9-22）注意到上面的公式中要求或，而在工程中信號(hào)是有限長(zhǎng)的，因此只能得到相關(guān)函數(shù)和協(xié)方差的估計(jì)值，當(dāng)t或N足夠長(zhǎng)時(shí)，估計(jì)值能精確地逼近真實(shí)值。9.2.3計(jì)算相關(guān)函數(shù)和協(xié)方差的MATLAB函數(shù)MATLAB信號(hào)處理工具箱提供了計(jì)算隨機(jī)信號(hào)相關(guān)函數(shù)xcorr。函數(shù)xcorr用于計(jì)算隨機(jī)序列自相關(guān)和互相關(guān)函數(shù)。調(diào)用格式為：[c,lags]=xcorr(x,y[,maxlags,’option’])式中，x,y為兩個(gè)獨(dú)立的隨機(jī)信號(hào)序列，長(zhǎng)度均為N；c為x,y的互相關(guān)估計(jì);lags為相關(guān)估計(jì)c的序號(hào)向量，其范圍為[-maxlags:maxlags]。option缺省或’none’時(shí)，函數(shù)xcorr按下式執(zhí)行非歸一化的相關(guān)：（9-23）上角加星號(hào)的變量表示其共軛。option為’biased’,計(jì)算有偏互相關(guān)函數(shù)估計(jì)（9-24）option為’unbiased’,計(jì)算無偏互相關(guān)函數(shù)估計(jì)（9-25）option為’coeff’,序列歸一化，使零延遲的自相關(guān)函數(shù)為1。maxlags為x和y的最大延遲,該項(xiàng)缺省時(shí)，函數(shù)返回值c的長(zhǎng)度為2N-1;該項(xiàng)不缺省時(shí)，函數(shù)返回值c的長(zhǎng)度為2*maxlags+1。該函數(shù)也可用于求一個(gè)隨機(jī)信號(hào)序列x(n)的自相關(guān)函數(shù)，調(diào)用格式為：c=xcorr(x,maxlags)MATLAB信號(hào)處理工具箱還提供了計(jì)算協(xié)方差函數(shù)xcov。若x為一個(gè)向量，調(diào)用c=xcov(x)函數(shù)，則輸出c為x序列的方差；若x為一個(gè)矩陣，則c=xcov(x)輸出矩陣x的協(xié)方差矩陣，其對(duì)角線元素為每一列數(shù)據(jù)的方差，非對(duì)角線元素表示對(duì)應(yīng)行和對(duì)應(yīng)列的協(xié)方差。調(diào)用c=xcov(x,y)，其中x,y具有相同的長(zhǎng)度，則c為式（9-22）所對(duì)應(yīng)的互協(xié)方差?！纠?-2】求帶有白噪聲干擾的頻率為10Hz的正弦信號(hào)和白噪聲信號(hào)的自相關(guān)函數(shù)并進(jìn)行比較。%Samp9_2clf;N=1000;Fs=500;%數(shù)據(jù)長(zhǎng)度和采樣頻率n=0:N-1;t=n/Fs;%時(shí)間序列Lag=100;%延遲樣點(diǎn)數(shù)randn('state',0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)x=sin(2*pi*10*t)+0.6*randn(1,length(t));%原始信號(hào)[c,lags]=xcorr(x,Lag,'unbiased');%對(duì)原始信號(hào)進(jìn)行無偏自相關(guān)估計(jì)subplot(2,2,1),plot(t,x);%繪原始信號(hào)xxlabel('時(shí)間/s');ylabel('x(t)');title('帶噪聲周期信號(hào)');gridon;subplot(2,2,2);plot(lags/Fs,c);%繪x信號(hào)自相關(guān),lags/Fs為時(shí)間序列xlabel('時(shí)間/s');ylabel('Rx(t)');title('帶噪聲周期信號(hào)的自相關(guān)');gridon;%信號(hào)x1x1=randn(1,length(x));%產(chǎn)生一與x長(zhǎng)度一致的隨機(jī)信號(hào)[c,lags]=xcorr(x1,Lag,'unbiased');%求隨機(jī)信號(hào)x1的無偏自相關(guān)subplot(2,2,3),plot(t,x1);%繪制隨機(jī)信號(hào)x1xlabel('時(shí)間/s');ylabel('x1(t)');title('噪聲信號(hào)');gridon;subplot(2,2,4);plot(lags/Fs,c);%繪制隨機(jī)信號(hào)x1的無偏自相關(guān)xlabel('時(shí)間/s');ylabel('Rx1(t)');title('噪聲信號(hào)的自相關(guān)');gridon圖9-1例9-2求得帶有白噪聲的周期信號(hào)和噪聲信號(hào)的自相關(guān)程序運(yùn)行結(jié)果見圖9-1。可以看出，含有周期成分和噪聲成分的自相關(guān)函數(shù)在時(shí)具有最大值。且在較大時(shí)仍具有明顯的周期性，其頻率和周期信號(hào)的頻率相同；而不含周期成分的純?cè)肼曅盘?hào)在時(shí)具有最大值。且在稍大時(shí)很快衰減至零。自相關(guān)的這一性質(zhì)被用來識(shí)別隨機(jī)信號(hào)中是否含有周期成分并用于確定所含周期成分的頻率?！纠?-3】已知兩個(gè)周期信號(hào)：，其中f=10Hz，求互相關(guān)函數(shù)。%Samp9_3clf;N=1000;Fs=500;%數(shù)據(jù)長(zhǎng)度和采樣頻率n=0:N-1;t=n/Fs;%時(shí)間序列Lag=200;%最大延遲單位x=sin(2*pi*10*t);%周期信號(hào)xy=0.5*sin(2*pi*10*t+90*pi/180);%與x有90o相移的信號(hào)y[c,lags]=xcorr(x,y,Lag,'unbiased');%求無偏互相關(guān)subplot(2,1,1);plot(t,x,'r');%繪制x信號(hào)holdon;plot(t,y,':');%在同一幅圖中繪y信號(hào)legend('x信號(hào)','y信號(hào)')xlabel('時(shí)間/s');ylabel('x(t)y(t)');title('原始信號(hào)');gridon;holdoffsubplot(2,1,2),plot(lags/Fs,c,'r');%繪制x,y的互相關(guān)xlabel('時(shí)間/s');ylabel('Rxy(t)');title('信號(hào)x和y的相關(guān)');gridon圖9-2例9-3所給的原始信號(hào)及其互相關(guān)程序運(yùn)行結(jié)果為圖9-2。由結(jié)果可見，也是周期為10Hz的余弦信號(hào)（注意兩個(gè)圖形其中的時(shí)間尺度有改變），其幅值為0.25，初始相位為90o。由此可得到互相關(guān)函數(shù)的一個(gè)重要特性：兩個(gè)均值為零且具有相同頻率的周期信號(hào)，其互相關(guān)函數(shù)保留原信號(hào)頻率、相位差和幅值信息。下面例子給出互相關(guān)函數(shù)的另一個(gè)工程應(yīng)用實(shí)例。若信號(hào)x(t)和信號(hào)y(t)是同類信號(hào)且有時(shí)移，用互相關(guān)函數(shù)可以準(zhǔn)確地計(jì)算出兩個(gè)信號(hào)時(shí)移大小。這種特性使得互相關(guān)函數(shù)廣泛地應(yīng)用于工程測(cè)量技術(shù)中。【例9-4】?jī)蓚€(gè)sinc信號(hào)有0.2s的時(shí)移，用互相關(guān)函數(shù)計(jì)算時(shí)移的大小。%Samp9_4clfN=1000;n=0:N-1;Fs=500;t=n/Fs;%數(shù)據(jù)個(gè)數(shù)采樣頻率和時(shí)間序列Lag=200;%最大延遲單位數(shù)[c,lags]=xcorr(x1,y1,Lag,'unbiased');%計(jì)算兩個(gè)函數(shù)的互相關(guān)subplot(2,1,1),plot(t,x1,'r');%繪第一個(gè)信號(hào)holdon;plot(t,y1,'b:');%在同一幅圖中繪第二個(gè)信號(hào)legend('信號(hào)x','信號(hào)y');%繪制圖例xlabel('時(shí)間/s');ylabel('x(t)y(t)');title('信號(hào)x和y');holdoffsubplot(2,1,2),plot(lags/Fs,c,'r');%繪制互相關(guān)信號(hào)xlabel('時(shí)間/s');ylabel('Rxy(t)');title('信號(hào)x和y的相關(guān)');程序運(yùn)行結(jié)果見圖9-3。可以清楚地看到第二個(gè)信號(hào)相對(duì)于第一個(gè)信號(hào)延遲了0.2s，即在-0.2s處出現(xiàn)相關(guān)極大值。因此可以采用該項(xiàng)技術(shù)檢測(cè)延遲信號(hào)。圖9-3例9-4兩個(gè)sinc信號(hào)的互相關(guān)9.3功率譜估計(jì)功率譜估計(jì)的目的是根據(jù)有限數(shù)據(jù)給出信號(hào)、隨機(jī)過程的頻率成分分布的描述。上節(jié)所述的相關(guān)分析是時(shí)域內(nèi)在噪聲背景下提取有用信息的途徑，而功率譜是頻域內(nèi)提取淹沒在噪聲中有用信息的分析方法。9.3.1功率譜密度假如隨機(jī)信號(hào)x(t)的自相關(guān)函數(shù)為，的Fourier變換為：（9-26）則定義Sx(f)為x(t)的自功率譜密度或稱為自功率譜。因?yàn)镾x(f)可解釋為x(t)的平均功率相對(duì)于頻率的分布函數(shù)。自功率譜Sx(f)包含的全部信息。如果噪聲信號(hào)中含有某種頻率成分，可以從自功率譜中看出。若隨機(jī)信號(hào)x(t)的自功率譜為Sx(f)，則根據(jù)Fourier逆變換可得：（9-27）和自功率譜類似，兩個(gè)隨機(jī)信號(hào)x(t)和y(t)互相關(guān)的頻率特性可用互功率譜密度來描述。互功率譜密度和互相關(guān)函數(shù)也是一Fourier變換對(duì)。（9-28）（9-29）對(duì)于離散隨機(jī)序列x(n)，自功率譜密度Sx(f)和自相關(guān)函數(shù)Rx(m)的關(guān)系為：（9-30）這里，Ts為數(shù)據(jù)采樣間隔。對(duì)于離散隨機(jī)序列x(n)和y(n)，互功率譜密度Sxy(f)和互相關(guān)函數(shù)Rxy(m)關(guān)系為：（9-31）且有（9-32）實(shí)際工程中隨機(jī)序列長(zhǎng)度均為有限長(zhǎng)，因此利用有限長(zhǎng)隨機(jī)序列計(jì)算的自功率譜密度和互功率譜密度只是真實(shí)值的一種估計(jì)。功率譜估計(jì)方法一般可分為參數(shù)估計(jì)和非參數(shù)估計(jì)兩類。MATLAB信號(hào)處理工具箱介紹的功率譜密度非參數(shù)估計(jì)方法有WELCH法、MTM法和MUSIC法；參數(shù)估計(jì)方法有MEM法，下面先介紹求取功率譜估計(jì)的基本方法，然后逐一介紹信號(hào)處理工具箱中給出的上述方法。9.3.2周期圖法1.基本方法周期圖法是直接將信號(hào)的采樣數(shù)據(jù)x(n)進(jìn)行Fourier變換求取功率譜密度估計(jì)的方法。假定有限長(zhǎng)隨機(jī)信號(hào)序列為x(n)。它的Fourier變換和功率譜密度估計(jì)存在下面的關(guān)系：（9-33）式中，N為隨機(jī)信號(hào)序列x(n)的長(zhǎng)度。在離散的頻率點(diǎn),有：（9-34）其中，F(xiàn)FT[x(n)]為對(duì)序列x(n)的Fourier變換，由于FFT[x(n)]的周期為N，求得的功率譜估計(jì)以N為周期，因此這種方法稱為周期圖法。下面用例子說明如何采用這種方法進(jìn)行功率譜估計(jì)。【例9-5】用Fourier變換算法求信號(hào)的功率譜。其中f1=50Hz,f2=120Hz,w(t)為白噪聲，采樣頻率Fs=1000Hz。（1）信號(hào)長(zhǎng)度N=256；（2）信號(hào)長(zhǎng)度N=1024。%Samp9_5clf;Fs=1000;%第一種情況:N=256N=256;Nfft=256;%數(shù)據(jù)長(zhǎng)度和FFT所用的數(shù)據(jù)長(zhǎng)度n=0:N-1;t=n/Fs;%采用的時(shí)間序列xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);%帶有噪聲的信號(hào)Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅譜平方的平均值，并轉(zhuǎn)換為dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%給出頻率序列subplot(2,1,1),plot(f,Pxx);%繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('周期圖N=256');gridon;%第二種情況:N=1024Fs=1000;N=1024;Nfft=1024;%采樣頻率、數(shù)據(jù)長(zhǎng)度和FFT所用數(shù)據(jù)長(zhǎng)度n=0:N-1;t=n/Fs;%時(shí)間序列xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);%帶有噪聲原始信號(hào)Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅譜平方的平均值，并轉(zhuǎn)換為dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%頻率序列subplot(2,1,2),plot(f,Pxx);%繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('周期圖N=1024');gridon程序的運(yùn)行結(jié)果為圖9-4?？梢钥闯觯陬l率50Hz和120Hz處，功率譜有兩個(gè)峰值，說明信號(hào)中含有50Hz和120Hz兩種頻率成分。但功率譜密度在很大范圍內(nèi)波動(dòng)，而且并沒有因信號(hào)取樣點(diǎn)數(shù)N的增加而有明顯改進(jìn)。注意本例中500Hz為Nyquist頻率，由第3章介紹的Fourier變換的知識(shí)可知，只觀看500Hz之前的功率譜可以推知以后的功率譜。在本章后面的例題中我們將只研究Nyquist頻率之前的功率譜。圖9-4例9-5中含有噪聲信號(hào)的功率譜用有限長(zhǎng)樣本序列的Fourier變換來表示隨機(jī)序列的功率譜，只是一種估計(jì)或近似，不可避免存在誤差。為了減少誤差，使功率譜估計(jì)更加平滑，可采用分段平均周期圖法（Bartlett法）、加窗平均周期圖法（Welch法）等方法加以改進(jìn)。2.分段平均周期圖法（Bartlett法）將信號(hào)序列x(n)，n=0,1,…,N-1，分成互不重疊的P個(gè)小段，每小段由m個(gè)采樣值，則P*m=N。對(duì)每個(gè)小段信號(hào)序列進(jìn)行功率譜估計(jì)，然后再取平均作為整個(gè)序列x(n)的功率譜估計(jì)。平均周期圖法還可以對(duì)信號(hào)x(n)進(jìn)行重疊分段，如按2:1重疊分段，即前一段信號(hào)和后一段信號(hào)有一半是重疊的。對(duì)每一小段信號(hào)序列進(jìn)行功率譜估計(jì)，然后再取平均值作為整個(gè)序列x(n)的功率譜估計(jì)。這兩種方法都稱為平均周期圖法，一般后者比前者好?！纠?-6】對(duì)例9-5中的信號(hào)序列，用兩種平均周期圖法求信號(hào)的功率譜密度估計(jì)。%Samp9_6clf;Fs=1000;%數(shù)據(jù)采樣點(diǎn)數(shù)%運(yùn)用信號(hào)不重疊分段估計(jì)功率譜N=1024;Nsec=256;n=0:N-1;t=n/Fs;%數(shù)據(jù)點(diǎn)數(shù)，分段間隔，時(shí)間序列randn('state',0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);%帶噪聲的原始信號(hào)pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;%第一段功率譜pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第二段功率譜pxx3=abs(fft(xn(515:768),Nsec).^2)/Nsec;%第三段功率譜pxx4=abs(fft(xn(769:1024),Nsec).^2)/Nsec;%第四段功率譜Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);%平均得到整個(gè)序列功率譜f=(0:length(Pxx)-1)*Fs/length(Pxx);%給出功率譜對(duì)應(yīng)的頻率subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('平均周期圖(無重疊)N=4*256');gridon%運(yùn)用信號(hào)重疊分段估計(jì)功率譜pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;%第一段功率譜pxx2=abs(fft(xn(129:384),Nsec).^2)/Nsec;%第二段功率譜pxx3=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第三段功率譜pxx4=abs(fft(xn(385:640),Nsec).^2)/Nsec;%第四段功率譜pxx5=abs(fft(xn(513:768),Nsec).^2)/Nsec;%第五段功率譜pxx6=abs(fft(xn(641:896),Nsec).^2)/Nsec;%第六段功率譜pxx7=abs(fft(xn(769:1024),Nsec).^2)/Nsec;%第七段功率譜Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7);%功率譜平均并轉(zhuǎn)化為dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%頻率序列subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('平均周期圖(重疊一半)N=1024');gridon圖9-5例9-6中重疊和不重疊分段計(jì)算的功率譜比較程序運(yùn)行結(jié)果為圖9-5，上圖采用不重疊分段法的功率譜估計(jì)，下圖為2:1重疊分段的功率譜估計(jì)，可見后者估計(jì)曲線較為平滑。與上例比較，平均周期圖法功率譜估計(jì)具有明顯效果（漲落曲線靠近0dB）。3.加窗平均周期圖法加窗平均周期圖法是對(duì)分段平均周期圖法的改進(jìn)。在信號(hào)序列x(n)分段后，用非矩形窗口對(duì)每一小段信號(hào)序列進(jìn)行預(yù)處理，再采用前述分段平均周期圖法進(jìn)行整個(gè)信號(hào)序列x(n)的功率譜估計(jì)。由窗函數(shù)的基本知識(shí)（第7章）可知，采用合適的非矩形窗口對(duì)信號(hào)進(jìn)行處理可減小“頻譜泄露”，同時(shí)可增加頻峰的寬度，從而提高頻譜分辨率?！纠?-7】對(duì)例9-5中的信號(hào)序列，用加窗平均周期圖法進(jìn)行功率譜密度估計(jì)。%Samp9_7clf;Fs=1000;%采樣頻率N=1024;Nsec=256;n=0:N-1;t=n/Fs;%數(shù)據(jù)長(zhǎng)度，分段數(shù)據(jù)長(zhǎng)度、時(shí)間序列w=hanning(256)';%采用的窗口數(shù)據(jù)randn('state',0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);%帶噪聲的信號(hào)%采用不重疊加窗方法的功率譜估計(jì)pxx1=abs(fft(w.*xn(1:256),Nsec).^2)/norm(w)^2;%第一段加窗振幅譜平方pxx2=abs(fft(w.*xn(257:512),Nsec).^2)/norm(w)^2;%第二段加窗振幅譜平方pxx3=abs(fft(w.*xn(513:768),Nsec).^2)/norm(w)^2;%第三段加窗振幅譜平方pxx4=abs(fft(w.*xn(769:1024),Nsec).^2)/norm(w)^2;%第四段加窗振幅譜平方Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);%求得平均功率譜，轉(zhuǎn)換為dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%求得頻率序列subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('加窗平均周期圖(無重疊)N=4*256');gridon%采用重疊加窗方法的功率譜估計(jì)pxx1=abs(fft(w.*xn(1:256),Nsec).^2)/norm(w)^2;%第一段加窗振幅譜平方pxx2=abs(fft(w.*xn(129:384),Nsec).^2)/norm(w)^2;%第二段加窗振幅譜平方pxx3=abs(fft(w.*xn(257:512),Nsec).^2)/norm(w)^2;%第三段加窗振幅譜平方pxx4=abs(fft(w.*xn(385:640),Nsec).^2)/norm(w)^2;%第四段加窗振幅譜平方pxx5=abs(fft(w.*xn(513:768),Nsec).^2)/norm(w)^2;%第五段加窗振幅譜平方pxx6=abs(fft(w.*xn(641:896),Nsec).^2)/norm(w)^2;%第六段加窗振幅譜平方pxx7=abs(fft(w.*xn(769:1024),Nsec).^2)/norm(w)^2;%第七段加窗振幅譜平方Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7);%平均功率譜轉(zhuǎn)換為dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%頻率序列subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%繪制功率譜曲線xlabel('頻率/Hz');ylabel('功率譜/dB');title('加窗平均周期圖(重疊一半)N=1024');gridon程序運(yùn)行結(jié)果為圖9-6。其中上圖采用無重疊數(shù)據(jù)分段的加窗平均周期圖法進(jìn)行功率譜估計(jì)，而下圖采用重疊數(shù)據(jù)分段的加窗平均周期圖法進(jìn)行功率譜估計(jì)，顯然后者是更佳的，信號(hào)譜峰加寬，而噪聲譜均在0dB附近，更為平坦（注意采用無重疊數(shù)據(jù)分段噪聲的最大的下降分貝數(shù)大于5dB，而重疊數(shù)據(jù)分段周期圖法噪聲的最大下降分貝數(shù)小于5dB）。圖9-6例9-7給出的不重疊加窗和重疊加窗周期圖法得到的功率譜曲線4.Welch法估計(jì)及其MATLAB函數(shù)Welch功率譜密度就是用改進(jìn)的平均周期圖法來求取隨機(jī)信號(hào)的功率譜密度估計(jì)的。Welch法采用信號(hào)重疊分段、加窗函數(shù)和FFT算法等計(jì)算一個(gè)信號(hào)序列的自功率譜估計(jì)(PSD如上例中的下半部分的求法)和兩個(gè)信號(hào)序列的互功率譜估計(jì)（CSD）。MATLAB信號(hào)處理工具箱函數(shù)提供了專門的函數(shù)PSD和CSD自動(dòng)實(shí)現(xiàn)Welch法估計(jì)，而不需要像以上例子一樣自己編程。函數(shù)psd利用Welch法估計(jì)一個(gè)信號(hào)自功率譜密度，函數(shù)調(diào)用格式為：[Pxx[,f]]=psd(x[,Nfft,Fs,window,Noverlap,’dflag’])式中，x為信號(hào)序列；Nfft為采用的FFT長(zhǎng)度。這一值決定了功率譜估計(jì)速度，當(dāng)Nfft采用2的冪時(shí)，程序采用快速算法；Fs為采樣頻率；Window定義窗函數(shù)和x分段序列的長(zhǎng)度。窗函數(shù)長(zhǎng)度必須小于或等于Nfft，否則會(huì)給出錯(cuò)誤信息；Noverlap為分段序列重疊的采樣點(diǎn)數(shù)（長(zhǎng)度），它應(yīng)小于Nfft;dflag為去除信號(hào)趨勢(shì)分量的選擇項(xiàng)：’linear’,去除線性趨勢(shì)分量，’mean’去除均值分量，’none’不做去除趨勢(shì)處理。Pxx為信號(hào)x的自功率譜密度估計(jì)。f為返回的頻率向量，它和Pxx對(duì)應(yīng)，并且有相同長(zhǎng)度。在psd函數(shù)調(diào)用格式中，缺省值為：Nfft=min(256,length(x)),Fs=2Hz,window=hanning(Nfft),noverlap=0.若x是實(shí)序列，函數(shù)psd僅計(jì)算頻率為正的功率譜估計(jì)?！纠?-8】對(duì)例9-5中的信號(hào)序列，用函數(shù)psd繪制自功率譜密度估計(jì)曲線。%Samp9_8clf;Fs=1000;%采樣頻率N=1024;Nfft=256;n=0:N-1;t=n/Fs;%數(shù)據(jù)長(zhǎng)度、時(shí)間序列window=hanning(256);%選用的窗口noverlap=128;%分段序列重疊的采樣點(diǎn)數(shù)（長(zhǎng)度）dflag='none';%不做趨勢(shì)處理randn('state',0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);%帶噪聲原始信號(hào)Pxx=psd(xn,Nfft,Fs,window,noverlap,dflag);%功率譜估計(jì)f=(0:Nfft/2)*Fs/Nfft;%求得對(duì)應(yīng)的頻率向量plot(f,10*log10(Pxx));%繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');title('PSD—Welch方法');gridon程序運(yùn)行結(jié)果見圖9-7?？梢钥吹剑捎肳elch法求得的功率譜與前面所述方法相比效果最好。使信號(hào)得以突出，而抑制了噪聲。圖9-7Welch法得到的功率譜估計(jì)注意程序前半部分中頻率向量f的創(chuàng)建方法。它與函數(shù)psd的輸出Pxx長(zhǎng)度的關(guān)系如下：若x為實(shí)序列，當(dāng)Nfft為奇數(shù)時(shí)，f=(0:(Nfft+1)/2-1)/Nfft;當(dāng)Nfft為偶數(shù)時(shí)，f=(0:Nfft/2)/Nfft。函數(shù)還有一種缺省返回值的調(diào)用格式，用于直接繪制信號(hào)序列x的功率譜估計(jì)曲線。函數(shù)還可以計(jì)算帶有置信區(qū)間的功率譜估計(jì)，調(diào)用格式為：[Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,Noverlap,p)式中，p為置信區(qū)間，?！纠?-9】設(shè)計(jì)一個(gè)歸一化頻率為0.2的FIR濾波器，對(duì)一個(gè)白噪聲信號(hào)序列進(jìn)行濾波，對(duì)濾波后的信號(hào)繪制置信區(qū)間為0.95的功率譜估計(jì)曲線。%Samp9_9r=randn(16384,1);%產(chǎn)生白噪聲隨機(jī)向量x=filter(h,1,r);%對(duì)隨機(jī)向量濾波psd(x,1024,10000,kaiser(512,5),0,0.95);%采用Kaiser窗進(jìn)行功率譜估計(jì)xlabel('頻率/Hz');ylabel('功率譜/dB')程序運(yùn)行結(jié)果見圖9-8。程序中首先產(chǎn)生一個(gè)FIR濾波器，其截止頻率為0.2，相對(duì)于本例中10000Hz的采樣頻率，5000Hz對(duì)應(yīng)的歸一化頻率為1，0.2的歸一化頻率對(duì)應(yīng)于1000Hz。噪聲為白噪聲信號(hào)，含有多種頻率成分。因此采用上述濾波器對(duì)噪聲進(jìn)行濾波后，截止頻率1000Hz以后的信號(hào)大大衰減，具有豐富頻率成分的白噪聲經(jīng)過濾波器濾波后輸出信號(hào)頻譜相當(dāng)于濾波器的幅頻特性曲線。其中上面曲線和下面曲線分別對(duì)應(yīng)于95%的置信區(qū)間的上下限。圖9-8例9-9設(shè)計(jì)濾波器對(duì)白噪聲濾波后的功率譜由此可知，濾波器輸入白噪聲序列的輸出信號(hào)的功率譜或自相關(guān)可以確定濾波器的頻率特性。（2）函數(shù)csd利用welch法估計(jì)兩個(gè)信號(hào)的互功率譜密度，函數(shù)調(diào)用格式為：[Pxy[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,’dflag’)[Pxy,Pxyc[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,p)這里，x,y為兩個(gè)信號(hào)序列；Pxy為x,y的互功率譜估計(jì)；其他參數(shù)的意義同自功率譜函數(shù)psd?！纠?-10】產(chǎn)生兩個(gè)長(zhǎng)度為16384的白噪聲信號(hào)和一個(gè)帶有白噪聲的1000Hz的周期信號(hào)，求兩個(gè)白噪聲信號(hào)及白噪聲與帶有噪聲的周期信號(hào)的互相關(guān)譜，并繪制互相關(guān)譜曲線。FFT所采用的長(zhǎng)度為1024，采用500個(gè)點(diǎn)的三角形窗，并且沒有重疊，采樣頻率Fs=10000Hz。%Samp9_10Fs=10000;%信號(hào)采樣頻率randn('state',0);%設(shè)置隨機(jī)狀態(tài)初始值x=randn(16384,1);%產(chǎn)生第一個(gè)白噪聲信號(hào)y=randn(16384,1);%產(chǎn)生第二個(gè)白噪聲信號(hào)z=2*sin(2*pi*1000*[1:16384]'/Fs)+randn(16384,1);%產(chǎn)生含周期成分的噪聲信號(hào)[Pxy,f]=csd(x,y,1024,Fs,triang(500),0);%白噪聲信號(hào)互相關(guān)譜,無重疊三角窗[Pxz,f]=csd(x,z,1024,Fs,triang(500),0);%噪聲與周期信號(hào)互相關(guān)譜無重疊三角窗subplot(2,1,1),plot(f,10*log10(Pxy)),gridon;%繪制x,y信號(hào)互相關(guān)譜ylabel('互相關(guān)譜/dB'),title('噪聲信號(hào)的互相關(guān)譜')subplot(2,1,2),plot(f,10*log10(Pxz)),gridon;%繪制x,z信號(hào)的互相關(guān)譜ylabel('互相關(guān)譜/dB'),xlabel('頻率/Hz')title('帶噪聲周期信號(hào)的互相關(guān)譜')程序運(yùn)行結(jié)果為圖9-9?？梢钥吹?，兩個(gè)白噪聲信號(hào)的互功率譜（上圖）雜亂無章，看不出周期成分，大部分功率譜在-5dB以下。然而白噪聲與帶有噪聲的周期信號(hào)的功率譜在其周期（頻率為1000Hz）處有一峰值，清楚地表明了周期信號(hào)的周期或頻率。因此，利用未知信號(hào)與白噪聲信號(hào)的互功率譜也可以檢測(cè)未知信號(hào)中所含有的頻率成分。圖9-9例9-10的兩個(gè)噪聲信號(hào)及噪聲信號(hào)和含噪聲周期信號(hào)的互功率譜9.3.3多窗口法多窗口法(Multitapermethod,簡(jiǎn)稱MTM法)利用多個(gè)正交窗口（Tapers）獲得各自獨(dú)立的近似功率譜估計(jì)，然后綜合這些估計(jì)得到一個(gè)序列的功率譜估計(jì)。相對(duì)于普通的周期圖法，這種功率譜估計(jì)具有更大的自由度，并在估計(jì)精度和估計(jì)波動(dòng)方面均有較好的效果。普通的功率譜估計(jì)只利用單一窗口，因此在序列始端和末端均會(huì)丟失相關(guān)信息，而且無法找回。而MTM法估計(jì)增加窗口使得丟失的信息盡量減少。MTM法簡(jiǎn)單地采用一個(gè)參數(shù)：時(shí)間帶寬積(Time-bandwidthproduct)NW，這個(gè)參數(shù)用以定義計(jì)算功率譜所用窗的數(shù)目，為2*NW-1。NW越大，功率譜計(jì)算次數(shù)越多，時(shí)間域分辨率越高，而頻率域分辨率降低，使得功率譜估計(jì)的波動(dòng)減小。隨著NW增大，每次估計(jì)中譜泄漏增多，總功率譜估計(jì)的偏差增大。對(duì)于每一個(gè)數(shù)據(jù)組，通常有一個(gè)最優(yōu)的NW使得在估計(jì)偏差和估計(jì)波動(dòng)兩方面求得折中，這需要在程序中反復(fù)調(diào)試來獲得。MATLAB信號(hào)處理工具箱中函數(shù)PMTM就是采用MTM法估計(jì)功率譜密度。函數(shù)調(diào)用格式為：[Pxx[,f]]=pmtm(x[,nw,Nfft,Fs])式中，x為信號(hào)序列；nw為時(shí)間帶寬積，缺省值為4。通?？扇?，5/2，3，7/2；Nfft為FFT長(zhǎng)度；Fs為采樣頻率。上面的函數(shù)還可以通過無返回值而繪出置信區(qū)間，如pmtm(x,nw,Nfft,Fs,’option’,p)繪制帶置信區(qū)間的功率譜密度估計(jì)曲線，?！纠?-11】用多窗口法(MTM)，分別采用NW=4和NW=2，估計(jì)例9-5的含有噪聲和50Hz和120Hz周期信號(hào)的功率譜密度。%Samp9_11clf;Fs=1000;%采樣頻率N=1024;Nfft=256;n=0:N-1;t=n/Fs;%數(shù)據(jù)長(zhǎng)度、分段數(shù)據(jù)長(zhǎng)度，時(shí)間序列randn('state',0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);%原始信號(hào)[Pxx1,f]=pmtm(xn,4,Nfft,Fs);%用多窗口法(NW=4)估計(jì)功率譜subplot(2,1,1),plot(f,10*log10(Pxx1));%繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');title('多窗口法(MTM)nw=4');gridon[Pxx,f]=pmtm(xn,2,Nfft,Fs);%用多窗口法(NW=2)估計(jì)功率譜subplot(2,1,2),plot(f,10*log10(Pxx));%繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');title('多窗口法(MTM)nw=2');gridon程序運(yùn)行結(jié)果見圖9-10。可以看到，NW=4和2均得到了比較好的功率譜估計(jì)，圖9-10上圖中采用NW=4估計(jì)的功率譜采用較多的窗口，使得每個(gè)窗口的數(shù)據(jù)點(diǎn)數(shù)減少，因此使得頻率域分辨率降低，因此上圖功率譜的波動(dòng)較少。圖9-10下圖為NW=2的情況，采用較少的窗口，每個(gè)窗口的數(shù)據(jù)點(diǎn)增多，使得頻率域的分辨率提高，但噪聲譜的分辨率也隨之增高，使得功率譜具有相對(duì)較大的波動(dòng)。圖9-10例9-11采用多窗口法估計(jì)含有周期信號(hào)的功率譜9.3.4最大熵法（Maxmumentropymethod，MEM法）如上所述，周期圖法功率譜估計(jì)需要對(duì)信號(hào)序列“截?cái)唷被蚣哟疤幚?，其結(jié)果是使估計(jì)的功率譜密度為信號(hào)序列真實(shí)譜和窗譜的卷積，導(dǎo)致誤差的產(chǎn)生。最大熵功率譜估計(jì)的目的是最大限度地保留截?cái)嗪髞G失的“窗口”以外信號(hào)的信息，使估計(jì)譜的熵最大。主要方法是以已知的自相關(guān)序列rxx(0),rxx(1),…,rxx(p)為基礎(chǔ)，外推自相關(guān)序列rxx(p+1),rxx(p+2),…，保證信息熵最大。最大熵功率譜估計(jì)法假定隨機(jī)過程是平穩(wěn)高斯過程，可以證明，隨機(jī)信號(hào)的最大熵譜與AR自回歸（全極點(diǎn)濾波器）模型譜（參看第8章）是等價(jià)的。MATLAB信號(hào)處理工具箱提供最大熵功率譜估計(jì)函數(shù)pmem，其調(diào)用格式為：[Pxx,f,a]=pmem(x,p,Nfft,Fs,’xcorr’)式中，x為輸入信號(hào)序列或輸入相關(guān)矩陣；p為全極點(diǎn)濾波器階次；a為全極點(diǎn)濾波器模型系數(shù)向量；’xcorr’是把x認(rèn)為是相關(guān)矩陣。【例9-12】用最大熵法估計(jì)例9-5的含有噪聲和50Hz和120Hz周期信號(hào)的功率譜密度，并與Welch法得到的結(jié)果進(jìn)行比較。%Samp9_12clf;Fs=1000;%采樣頻率N=1024;Nfft=256;n=0:N-1;t=n/Fs;%數(shù)據(jù)長(zhǎng)度、分段長(zhǎng)度和時(shí)間序列window=hanning(256);%采用窗口randn('state',0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);%原始信號(hào)[Pxx1,f]=pmem(xn,14,Nfft,Fs);%采用最大熵法，采用濾波器階數(shù)14，估計(jì)功率譜subplot(2,1,1),plot(f,10*log10(Pxx1));%繪制功率譜xlabel('頻率/Hz');ylabel('功率譜/dB');title('最大熵法(MTM)Order=14');gridon%采用Welch方法估計(jì)功率譜noverlap=128;%重疊數(shù)據(jù)dflag='none';subplot(2,1,2)psd(xn,Nfft,Fs,window,noverlap,dflag);%采用Welch方法估計(jì)功率譜xlabel('頻率/Hz');ylabel('功率譜/dB')title('Welch方法估計(jì)的功率譜');gridon程序運(yùn)行結(jié)果為圖9-11。比較最大熵功率譜估計(jì)(MEM)和改進(jìn)的平均周期圖功率譜估計(jì)，可見，MEM法估計(jì)的功率譜曲線較光滑。在這一方法中，MEM法選定全極點(diǎn)濾波器的階數(shù)取得越大，能夠獲得的窗口外的信息越多，但計(jì)算量也越大，需要根據(jù)情況折中考慮。圖9-11最大熵法與Welch法估計(jì)的功率譜的比較9.3.5多信號(hào)分類法MATLAB信號(hào)處理工具箱還提供另一種功率譜估計(jì)函數(shù)pmusic。該函數(shù)執(zhí)行多信號(hào)分類法（multiplesignalclassification，Music法)。將數(shù)據(jù)自相關(guān)矩陣看成由信號(hào)自相關(guān)矩陣和噪聲自相關(guān)矩陣兩部分組成，即數(shù)據(jù)自相關(guān)矩陣R包含有兩個(gè)子空間信息：信號(hào)子空間和噪聲子空間。這樣，矩陣特征值向量（Eigenvector）也可分為兩個(gè)子空間：信號(hào)子空間和噪聲子空間。為了求得功率譜估計(jì)，函數(shù)pmusic計(jì)算信號(hào)子空間和噪聲子空間的特征值向量函數(shù)，使得在周期信號(hào)頻率處函數(shù)值最大，功率譜估計(jì)出現(xiàn)峰值，而在其他頻率處函數(shù)值最小。其調(diào)用格式為：[Pxx,f,a]=pmusic(x,p[[,thresh],Nfft,Fs,window,Noverlap])式中,x為輸入信號(hào)的向量或矩陣；p為信號(hào)子空間維數(shù)；thresh為閾值，其他參數(shù)的意義與函數(shù)psd相同。【例9-13】用多信號(hào)分類法，采用7個(gè)窗口，估計(jì)例9-5的含有噪聲和50Hz和120Hz周期信號(hào)的功率譜密度。%Samp9_13clf;Fs=1000;%采樣頻率N=1024;Nfft=256;n=0:N-1;t=n/Fs;%數(shù)據(jù)長(zhǎng)度、分段數(shù)據(jù)和時(shí)間序列randn('state',0)%設(shè)置產(chǎn)生隨機(jī)數(shù)的初始狀態(tài)xn=sin(2*pi*100*t)+2*sin(2*pi*200*t)+randn(1,N);%原始信號(hào)pmusic(xn,[7,1.1],Nfft,Fs,32,16);%采用多信號(hào)分類法估計(jì)功率譜xlabel('頻率/Hz');ylabel('功率譜/dB')title('通過MUSIC法估計(jì)的偽譜')程序運(yùn)行結(jié)果為圖9-12?？梢妶D形較為清楚地識(shí)別出了信號(hào)中所含的頻率成分，并且具有較高的分辨率。但要注意，函數(shù)pmusic參數(shù)的選擇對(duì)估計(jì)的功率譜影響較大。圖9-12多信號(hào)分類法估計(jì)的功率譜9.4傳遞函數(shù)估計(jì)在第8章我們學(xué)過，如果已知一個(gè)系統(tǒng)（濾波器）的輸入和輸出，我們可以估計(jì)該系統(tǒng)的傳遞函數(shù)。然而濾波過程是存在誤差的，即信號(hào)中含有噪聲。利用含有噪聲的輸入和輸出信號(hào)對(duì)系統(tǒng)（濾波器）傳遞函數(shù)的估計(jì)用功率譜密度比較合適。一個(gè)線性時(shí)不變系統(tǒng)，頻率特性為H()，x(n)和y(n)分別為系統(tǒng)的輸入和輸出?？梢宰C明：（9-35）式中，Pxx為x(n)的自功率譜密度；Pxy為x(n)和y(n)的互功率譜密度。由此還可以證明，輸入x(n)和輸出y(n)之間的傳遞函數(shù)估計(jì)為：（9-36）用MATLAB編程很容易實(shí)現(xiàn)上面的估計(jì)。當(dāng)x(n)，y(n)已知時(shí)，由函數(shù)psd求得，由函數(shù)csd求得，再由上式求得系統(tǒng)的傳遞函數(shù)估計(jì)。MATLAB信號(hào)處理工具箱還提供函數(shù)tfe可直接完成上面的運(yùn)算。函數(shù)tfe用于從系統(tǒng)輸入和輸出的功率譜估計(jì)中求取傳遞函數(shù)的估計(jì)。函數(shù)的調(diào)用格式為：[Txy[,f]=]tfe(x,y[[,Nfft],window,noverlap])式中，x和y分別為系統(tǒng)輸入和輸出信號(hào)向量；Txy為系統(tǒng)傳遞函數(shù)的復(fù)數(shù)頻率響應(yīng)。Nfft=min(256,length(x)),window=hanning(Nfft),noverlap=0,其意義與csd函數(shù)的意義相同?！纠?-14】設(shè)FIR濾波器的脈沖響應(yīng)為元素均為1/10、長(zhǎng)度為10的信號(hào)。產(chǎn)生一個(gè)長(zhǎng)度為1024的帶有白噪聲的周期信號(hào)，并采用上述濾波器進(jìn)行濾波，采用輸入信號(hào)和濾波后的輸出數(shù)據(jù)估計(jì)該濾波器的幅頻響應(yīng)并與實(shí)際幅頻響應(yīng)作比較。%Samp9_14clf;Fs=1000;%信號(hào)采樣頻率N=1024;Nfft=256;n=0:N-1;t=n/Fs;%數(shù)據(jù)長(zhǎng)度、分段長(zhǎng)度和時(shí)間序列window=hanning(256);%采用的窗口noverlap=128;%重疊數(shù)據(jù)個(gè)數(shù)dflag='none';randn('state',0);%設(shè)置產(chǎn)生隨機(jī)數(shù)的狀態(tài)xn=sin(2*pi*50*t)+randn(1,N);%輸入信號(hào)h=ones(1,10)/10;%設(shè)計(jì)一個(gè)濾波器的脈沖響應(yīng)yn=filter(h,1,xn);%濾波后信號(hào)%采用自功率譜和互功率譜估計(jì)的方法[Pxx,f1]=psd(xn,Nfft,Fs,window,noverlap,dflag);%估計(jì)自功率譜[Pxy,f1]=csd(xn,yn,Nfft,Fs,window,noverlap,dflag);%估計(jì)互功率譜He=Pxy./Pxx;%求得系統(tǒng)傳遞函數(shù)[HEST,f]=tfe(xn,yn,Nfft,Fs,window,noverlap,dflag);%采用tfe函數(shù)估計(jì)傳遞函數(shù)%比較結(jié)果H=freqz(h,1,f,Fs);%求得系統(tǒng)的頻率響應(yīng)subplot(3,1,1),plot(f,abs(H));%繪制原濾波器幅頻特性ylabel('振幅');title('實(shí)際模型的幅頻響應(yīng)');axis([050001]);gridon;subplot(3,1,2),plot(f,abs(HEST));%繪制用tfe函數(shù)估計(jì)的傳遞函數(shù)幅頻特性ylabel('振幅');title('TFE函數(shù)估計(jì)的幅頻響應(yīng)');axis([050001]);gridonsubplot(3,1,3),plot(f1,abs(He));%繪制用功率譜估計(jì)的傳遞函數(shù)幅頻特性xlabel('頻率/Hz');ylabel('振幅');title('功率譜估計(jì)的幅頻響應(yīng)');axis([050001]);gridon程序的運(yùn)行結(jié)果為圖9-13。圖中上圖為真實(shí)濾波器模型的幅頻響應(yīng)，中圖為用tfe函數(shù)計(jì)算得到的濾波器幅頻響應(yīng)，下圖為采用自功率譜和互功率譜得到的濾波器幅頻響應(yīng)。可以看到估計(jì)的濾波器幅頻響應(yīng)與真實(shí)濾波器幅頻響應(yīng)具有較好的一致性。另外，采用自功率譜和互功率譜計(jì)算得到的濾波器幅頻響應(yīng)和采用tfe函數(shù)得到的濾波器幅頻響應(yīng)一致表明了處理方法的一致性。圖9-13例9-14真實(shí)FIR濾波器幅頻特性和估計(jì)幅頻特性的比較9.5相干函數(shù)兩個(gè)信號(hào)x(t)和y(t)之間的相干函數(shù)（Coherencefunction）為：（9-37）式中，為x(t)和y(t)的互功率譜密度，，分別為x(t),y(t)的自功率譜密度。相干函數(shù)為0~1的實(shí)數(shù)，它用來檢測(cè)信號(hào)x(t)和y(t)在頻域內(nèi)的相關(guān)程度。若y(t)為x(t)的線性響應(yīng)，則=1，若x(t)和y(t)完全互不相