課程名稱：現(xiàn)代信號(hào)分析考試形式：□專題研究報(bào)告□論文大作業(yè)□綜合考試學(xué)生姓名：學(xué)號(hào)：序號(hào)分項(xiàng)類別得分1題目12題目23題目34題目45題目56題目67題目78題目89題目910題目10總分評(píng)閱人：時(shí)間：年月目錄一、FFT算法的特點(diǎn)。 -3-二、平穩(wěn)隨機(jī)過(guò)程 -8-三、功率譜估計(jì) -10-四、自適應(yīng)濾波器（LMS和RLS） -16-五、維納濾波器 -23-六、FIR維納濾波器 -25-七、AR/MA濾波器 -28-八、卡爾曼濾波器 -31-九、小波分析方法及應(yīng)用 -37-十、其他現(xiàn)代信號(hào)分析方法 -43-現(xiàn)代信號(hào)分析48-47-現(xiàn)代信號(hào)分析27-一、FFT算法的特點(diǎn)。題目：請(qǐng)舉例說(shuō)明FFT算法的特點(diǎn)計(jì)算離散傅里葉變換的一種快速算法，簡(jiǎn)稱FFT?？焖俑道锶~變換是1965年由J.W.庫(kù)利和T.W.圖基提出的。采用這種算法能使計(jì)算機(jī)計(jì)算離散傅里葉變換所需要的乘法次數(shù)大為減少，特別是被變換的抽樣點(diǎn)數(shù)N越多，F(xiàn)FT算法計(jì)算量的節(jié)省就越顯著。根據(jù)對(duì)序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法，基本算法是基2DIT和基2DIF。FFT在離散傅里葉反變換、線性卷積和線性相關(guān)等方面也有重要應(yīng)用?？焖俑凳献儞Q（FFT），是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實(shí)等特性，對(duì)離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。有限列長(zhǎng)為N的序列x（n）的DFT變換公式在MATLAB中的表達(dá)式為其逆變換為FFT主要就是利用了以下兩個(gè)特性使長(zhǎng)序列的DFT分解為更小點(diǎn)數(shù)的DFT來(lái)實(shí)現(xiàn)的。利用的對(duì)稱性使DFT運(yùn)算中有些項(xiàng)合并利用的周期性和對(duì)稱性使長(zhǎng)序列的DFT分解為更小點(diǎn)數(shù)的DFTFFT算法在應(yīng)用時(shí)能夠大大減少計(jì)算量，但是也存在著和DFT同樣的一些問題，主要有：頻譜混疊，泄露問題以及柵欄效應(yīng)的問題，下面通過(guò)運(yùn)用MATLAB的FFT相關(guān)函數(shù)對(duì)信號(hào)進(jìn)行變換來(lái)進(jìn)行實(shí)際分析。頻譜混疊奈奎斯特定理已被眾所周知了，所以幾乎所有人的都知道為了不讓頻譜混疊，理論上采樣頻譜大于等于信號(hào)的最高頻率。采樣周期的倒數(shù)是頻譜分辨率，最高頻率的倒數(shù)是采樣周期。設(shè)定采樣點(diǎn)數(shù)為N，采樣頻率fs，最高頻率fh，故頻譜分辨率f=fs/N,而fs>=2fh，所以可以看出最高頻率與頻譜分辨率是相互矛盾的，提高頻譜分辨率f的同時(shí)，在N確定的情況下必定會(huì)導(dǎo)致最高頻率fh的減??；同樣的，提高最高頻率fh的同時(shí)必會(huì)引起f的增大，即分辨率變大,下面舉例說(shuō)明。設(shè)輸入信號(hào)為，即它的主要頻率為120Hz和150Hz。分別取采樣頻率為400Hz，300Hz，200Hz，100Hz，程序如下：L=1024;%采樣點(diǎn)數(shù)NFFT=2^nextpow2(L);%Nextpowerof2fromlengthofyFs1=400;%采樣頻率t1=(0:L-1)/Fs1;s1=3*cos(2*pi*120*t1)+5*sin(2*pi*150*t1);Y1=fft(s1,NFFT)/L;f1=Fs1/2*linspace(0,1,NFFT/2+1);Fs2=300;t2=(0:L-1)/Fs2;s2=3*cos(2*pi*120*t2)+5*sin(2*pi*150*t2);Y2=fft(s2,NFFT)/L;f2=Fs2/2*linspace(0,1,NFFT/2+1);Fs3=200;t3=(0:L-1)/Fs3;s3=3*cos(2*pi*120*t3)+5*sin(2*pi*150*t3);Y3=fft(s3,NFFT)/L;f3=Fs3/2*linspace(0,1,NFFT/2+1);Fs4=100;t4=(0:L-1)/Fs4;s4=3*cos(2*pi*120*t4)+5*sin(2*pi*150*t4);Y3=fft(s4,NFFT)/L;f4=Fs4/2*linspace(0,1,NFFT/2+1);subplot(2,2,1);plot(f1,2*abs(Y1(1:NFFT/2+1)));title('采樣頻率為400Hz');gridon;xlabel('頻率(Hz)');subplot(2,2,2);plot(f2,2*abs(Y2(1:NFFT/2+1)));title('采樣頻率為300Hz');gridon;xlabel('頻率(Hz)');subplot(2,2,3);plot(f3,2*abs(Y3(1:NFFT/2+1)));title('采樣頻率為200Hz');gridon;xlabel('頻率(Hz)');subplot(2,2,4);plot(f4,2*abs(Y4(1:NFFT/2+1)));title('采樣頻率為100Hz');gridon;xlabel('頻率(Hz)');gridon;結(jié)果如圖1-1所示，可見當(dāng)采樣頻率為300Hz時(shí)，沒有發(fā)生頻率混疊，但當(dāng)頻率小于兩倍的信號(hào)頻率時(shí)，經(jīng)過(guò)采樣就會(huì)發(fā)生混疊現(xiàn)象，輸入信號(hào)中的某些頻率就無(wú)法辨識(shí)出。圖1-1FFT的混淆現(xiàn)象泄漏現(xiàn)象實(shí)際中的信號(hào)序列往往很長(zhǎng)，甚至是無(wú)限長(zhǎng)序列。為了方便，我們往往用截短的序列來(lái)近似它們。這樣可以使用較短的DFT來(lái)對(duì)信號(hào)進(jìn)行頻譜分析。這種截短等價(jià)于給原信號(hào)序列乘以一個(gè)窗函數(shù)。而窗函數(shù)的頻譜不是有限帶寬的，從而它和原信號(hào)的頻譜進(jìn)行卷積以后會(huì)擴(kuò)展原信號(hào)的頻譜。值得一提的是，泄漏是不能和混疊完全分離開的，因?yàn)樾孤秾?dǎo)致頻譜的擴(kuò)展，從而造成混疊。為了減小泄漏的影響，可以選擇適當(dāng)?shù)拇昂瘮?shù)使頻譜的擴(kuò)散減到最小。下面舉例說(shuō)明。設(shè)輸入信號(hào)為，即它的主要頻率為120Hz和150Hz。，用MATLAB自帶的信號(hào)處理工具箱（通過(guò)命令sptool調(diào)出）分別對(duì)信號(hào)加矩形窗（boxcar），海明窗（Hamming），漢寧窗（Hanning）和布萊克曼窗（Blackman），分析不同的窗函數(shù)對(duì)減小泄露現(xiàn)象的影響。具體操作如下：先生成信號(hào)s，然后在sptool工具箱的file菜單中import信號(hào)s，F(xiàn)s為400Hz，如圖1-2所示。圖1-2sptool工具箱界面然后點(diǎn)create按鈕，在patameters欄選擇welch方法，Nfft點(diǎn)數(shù)1024，Nwind為256，分別加布萊克曼窗,矩形窗，海明窗，和三角窗，得如下PSD圖。(a)blackman(b)boxcar(c)hamming(d)triang圖1-3FFT的泄露現(xiàn)象由以上四圖可見，布萊克曼窗主瓣寬，旁瓣小，頻率識(shí)別精度最低，但幅值識(shí)別精度最高。矩形窗主瓣窄，旁瓣大，頻率識(shí)別精度最高，幅值識(shí)別精度最低；海明窗是余弦窗的一種，又稱改進(jìn)的升余弦窗，海明窗與漢寧窗都是余弦窗，只是加權(quán)系數(shù)不同。海明窗加權(quán)的系數(shù)能使旁瓣達(dá)到更小，海明窗的頻譜是由3個(gè)矩形時(shí)窗的頻譜合成，但其旁瓣衰減速度為20dB／（10oct）。三角窗與矩形窗比較,主瓣寬約等于矩形窗的兩倍,但旁瓣小,而且無(wú)負(fù)旁瓣。所以對(duì)于窗函數(shù)的選擇，應(yīng)考慮被分析信號(hào)的性質(zhì)與處理要求。如果僅要求精確讀出主瓣頻率，而不考慮幅值精度，則可選用主瓣寬度比較窄而便于分辨的矩形窗，如果分析窄帶信號(hào)，且有較強(qiáng)的干擾噪聲，則應(yīng)選用旁瓣幅度小的窗函數(shù)，如漢寧窗。柵欄效應(yīng)柵欄效應(yīng)是因?yàn)橛糜?jì)算頻譜只限制為基頻的整數(shù)倍而不可能將頻譜視為一連續(xù)函數(shù)而產(chǎn)生的。設(shè)信號(hào)的最高頻率為，要求的頻率分辨率為F（即在頻率軸上所能得到的最小頻率間隔），則采樣點(diǎn)數(shù)N需要滿足以下條件：若不滿足，就會(huì)產(chǎn)生柵欄效應(yīng)。減小柵欄效應(yīng)通常有兩個(gè)方法：一個(gè)方法是在原紀(jì)錄末端添加一些零值點(diǎn)來(lái)變動(dòng)時(shí)間周期內(nèi)的點(diǎn)數(shù)，并保持記錄不變。從而在保持原有頻譜連續(xù)形式不變的情況下，變更了譜線的位置。這樣，原來(lái)看不到的頻譜分量就能移到可見的位置上了。但是，補(bǔ)零并沒有增加序列的有效長(zhǎng)度，所以并不能提高分辨率。但補(bǔ)零同時(shí)可以使數(shù)據(jù)長(zhǎng)度N為2的整數(shù)冪，以便于應(yīng)用FFT。另一個(gè)方法是增加采樣點(diǎn)數(shù)，從而減小頻率間隔，使原來(lái)被擋住的一些頻譜的峰點(diǎn)或谷點(diǎn)顯露出來(lái)。注意，這時(shí)候每根譜線多對(duì)應(yīng)的頻率和原來(lái)的已經(jīng)不相同了。下面通過(guò)實(shí)際的函數(shù)來(lái)分析兩種方法的效果。程序如下：L=90;%采樣點(diǎn)數(shù)NFFT=2^nextpow2(L);%Nextpowerof2fromlengthofyFs=400;%2采樣頻率t=(0:L-1)/Fs;s=cos(200*pi*t)+2*cos(204*pi*t);Y1=fft(s,NFFT)/L;f1=Fs/2*linspace(0,1,NFFT/2+1);plot(f1,2*abs(Y1(1:NFFT/2+1)));xlabel('頻率');%原序列末端補(bǔ)零L3=1024;NFFT=2^nextpow2(L3);s3=[szeros(1,L3-L)];Y3=fft(s3,NFFT)/L3;f3=Fs/2*linspace(0,1,NFFT/2+1);figure;plot(f3,2*abs(Y3(1:NFFT/2+1)));xlabel('頻率');('頻率');%增加采樣點(diǎn)數(shù)L2=1024;NFFT=2^nextpow2(L2);t=(0:L2-1)/Fs;s2=cos(200*pi*t)+2*cos(204*pi*t);Y2=fft(s2,NFFT)/L2;f2=Fs/2*linspace(0,1,NFFT/2+1);figure;plot(f2,2*abs(Y2(1:NFFT/2+1)));xlabel('頻率');設(shè)輸入信號(hào)為，采樣頻率為400Hz。可見，信號(hào)含有2個(gè)頻率分量，分別為100Hz和102Hz，則要求的頻率分辨率為2Hz，可算得N為202，即至少采樣202個(gè)點(diǎn)，才不會(huì)發(fā)生柵欄現(xiàn)象。首先取采樣點(diǎn)數(shù)為90，對(duì)其進(jìn)行FFT，得到的圖形如圖1-4所示。圖1-4FFT的柵欄效應(yīng)可見，雖然采樣頻率滿足香農(nóng)采樣定理，但由于采樣點(diǎn)數(shù)較少，頻率為20Hz和20.5Hz的信號(hào)無(wú)法分辨。下面采用補(bǔ)零的方法，將N補(bǔ)滿1024，得到的FFT圖形如圖1-5所示。圖1-5信號(hào)末端補(bǔ)零改善柵欄效應(yīng)由上圖可以看出，補(bǔ)零對(duì)分辨率沒有影響，但是對(duì)頻譜起到了平滑的作用。同時(shí)，由于增加了大量零值，使得信號(hào)的平均幅頻響應(yīng)大大減小。然后采用增加采樣點(diǎn)的方法，將采樣點(diǎn)數(shù)增加到1024，F(xiàn)FT后的圖形如圖1-6所示。圖1-6增加采樣點(diǎn)數(shù)改善柵欄效應(yīng)由上圖可以看出，增加采樣點(diǎn)數(shù)后就能分辨出100Hz和102Hz的信號(hào)，有效地改善了柵欄效應(yīng)。從時(shí)域看，這種方法相當(dāng)于對(duì)信號(hào)進(jìn)行整周期采樣，實(shí)際中常用此方法來(lái)提高周期信號(hào)的頻譜分析精度。二、平穩(wěn)隨機(jī)過(guò)程題目：什么是寬平穩(wěn)隨機(jī)過(guò)程？什么是嚴(yán)平穩(wěn)隨機(jī)過(guò)程？它們之間有什么聯(lián)系？（5分）平穩(wěn)隨機(jī)過(guò)程是指它的統(tǒng)計(jì)特性不隨時(shí)間的推移而變化，在譜估計(jì)中可以認(rèn)為平穩(wěn)隨機(jī)信號(hào)是由一個(gè)白噪聲激勵(lì)一個(gè)線性時(shí)不變系統(tǒng)產(chǎn)生的。2.1寬平穩(wěn)隨機(jī)過(guò)程設(shè)有一個(gè)二階隨機(jī)過(guò)程，它的均值為常數(shù)，自相關(guān)函數(shù)僅是時(shí)間差的函數(shù)，則稱它為寬平穩(wěn)隨機(jī)過(guò)程，用符號(hào)化語(yǔ)言表示，給定一個(gè)二階矩過(guò)程，如果對(duì)任意，有：（1）（常數(shù)）（2）自相關(guān)函數(shù)則稱為寬平穩(wěn)隨機(jī)過(guò)程或廣義平穩(wěn)隨機(jī)過(guò)程。2.2嚴(yán)平穩(wěn)隨機(jī)過(guò)程在數(shù)學(xué)中，平穩(wěn)隨機(jī)過(guò)程（Stationaryrandomprocess）或者嚴(yán)平穩(wěn)隨機(jī)過(guò)程（Strictly-sensestationaryrandomprocess），又稱狹義平穩(wěn)隨機(jī)過(guò)程，是指隨機(jī)過(guò)程在某一固定時(shí)間和位置的概率分布與所有時(shí)間和位置的概率分布相同，即隨機(jī)過(guò)程的統(tǒng)計(jì)特性不隨時(shí)間的推移而變化。這樣，數(shù)學(xué)期望和方差這些參數(shù)也不隨時(shí)間和位置變化。用符號(hào)化語(yǔ)言表示出來(lái)，設(shè)隨機(jī)過(guò)程，對(duì)于任意的n（n=1,2,···）和任意選定的以及為任意值，當(dāng)時(shí)，n維隨機(jī)變量和具有相同的概率分布，即：則稱隨機(jī)過(guò)程具有平穩(wěn)性，稱此過(guò)程為嚴(yán)平穩(wěn)隨機(jī)過(guò)程，簡(jiǎn)稱隨機(jī)過(guò)程。2.3嚴(yán)平穩(wěn)隨機(jī)過(guò)程與寬平穩(wěn)隨機(jī)過(guò)程區(qū)別聯(lián)系（1）一個(gè)寬平穩(wěn)過(guò)程不一定是嚴(yán)平穩(wěn)過(guò)程，一個(gè)嚴(yán)平穩(wěn)過(guò)程也不一定是寬平穩(wěn)過(guò)程。例1：，其中w服從U(0,2π)，隨機(jī)過(guò)程是寬平穩(wěn)過(guò)程，但不是嚴(yán)平穩(wěn)過(guò)程；例2：服從柯西分布的隨機(jī)變量序列是嚴(yán)平穩(wěn)隨機(jī)過(guò)程，但不是寬平穩(wěn)隨機(jī)過(guò)程。（2）寬平穩(wěn)過(guò)程定只涉及與一維、二維分布有關(guān)的數(shù)字特征，所以一個(gè)嚴(yán)平穩(wěn)過(guò)程只要二階矩存在（或有界），則必定是寬平穩(wěn)隨機(jī)過(guò)程。但反過(guò)來(lái)一般不成立。注：二階矩過(guò)程定義，如果隨機(jī)過(guò)程對(duì)第一個(gè)，二階矩都存在，那么稱其為二階矩過(guò)程。（3）正態(tài)過(guò)程是一個(gè)重要特例，一個(gè)寬平穩(wěn)的正態(tài)過(guò)程必定是嚴(yán)平穩(wěn)的。這是因?yàn)椋赫龖B(tài)過(guò)程的概率密度是由均值函數(shù)和自相關(guān)函數(shù)完全確定的，因而如果均值函數(shù)和自相關(guān)函數(shù)不隨時(shí)間的推移而變化，則概率密度函數(shù)也不隨時(shí)間的推移發(fā)生變化。三、功率譜估計(jì)題目：簡(jiǎn)述經(jīng)典功率譜估計(jì)與現(xiàn)代功率譜估計(jì)的差別，并計(jì)算下題：（14分）設(shè)隨機(jī)過(guò)程是單位方差白噪聲激勵(lì)如下的系統(tǒng)而產(chǎn)生的：取序列的長(zhǎng)度分別為128和1024數(shù)據(jù)段，并用經(jīng)典和現(xiàn)代譜估計(jì)方法進(jìn)行譜估計(jì)，分析不同參數(shù)對(duì)最終結(jié)果的影響。解答：3.1經(jīng)典功率譜估計(jì)和現(xiàn)代功率譜估計(jì)的差別信號(hào)的頻譜分析是研究信號(hào)特性的重要手段之一，通常是求其功率譜來(lái)進(jìn)行頻譜分析。功率譜反映了隨機(jī)信號(hào)各頻率成份功率能量的分布情況，可以提示信號(hào)中隱含的周期性及靠得很近的譜峰等有用信息，在許多領(lǐng)域都發(fā)揮了重要作用。然而，實(shí)際應(yīng)用中的平穩(wěn)隨機(jī)信號(hào)通常是有限長(zhǎng)的，只能根據(jù)有限長(zhǎng)信號(hào)估計(jì)原信號(hào)的真實(shí)功率譜。因此，功率譜估計(jì)就是基于有限的數(shù)據(jù)尋找信號(hào)、隨機(jī)過(guò)程或系統(tǒng)的頻率成分。一般功率譜估計(jì)分為經(jīng)典譜估計(jì)和現(xiàn)代譜估計(jì)。其各自的特點(diǎn)及相互差別如下。根據(jù)所用理論的不同，通常將基于相關(guān)函數(shù)的傅里葉變換的估計(jì)方法稱為經(jīng)典功率譜估計(jì)，而將參數(shù)模型估計(jì)方法和基于相關(guān)矩陣特征值分解的信號(hào)頻率估計(jì)方法，稱為現(xiàn)代功率譜估計(jì)方法。經(jīng)典譜估計(jì)為線性估計(jì)方法，是建立在傳統(tǒng)的傅里葉變換的基礎(chǔ)之上，它首先由給定的數(shù)據(jù)估計(jì)自相關(guān)序列，然后對(duì)估計(jì)出的自相關(guān)序列進(jìn)行傅里葉變換得到功率譜估計(jì)，沒有將信號(hào)的可用信息結(jié)合到估計(jì)過(guò)程中，有計(jì)算效率高、估計(jì)值正比于正弦波信號(hào)的功率等優(yōu)勢(shì)。但性能依賴于數(shù)據(jù)序列的長(zhǎng)度，頻率分辨率低，且需在方差和分辨率之間做出權(quán)衡，不適用于短時(shí)數(shù)據(jù)的情況，典型代表有Blackman和Tukey提出的自相關(guān)譜估計(jì)法（簡(jiǎn)稱BT法），周期圖法以及周期圖法的改進(jìn)方法——Barlett法和Welth法；而現(xiàn)代功率譜估計(jì)是非線性估計(jì)方法，估計(jì)性能依賴于參數(shù)模型，最終可獲得方差小、分辨率高的譜估計(jì)，但所用模型必須適合于所分析信號(hào)，否則譜估計(jì)將是錯(cuò)誤或者不準(zhǔn)確的。有AR模型，MA模型，ARMA模型等；非參數(shù)模型譜估計(jì)有最小方差法和MUSIC法等。3.2功率譜計(jì)算由上式的輸入輸出關(guān)系可以看出這是一個(gè)AR模型，由已知數(shù)據(jù)得到：本題中的隨機(jī)過(guò)程是由單位方差白噪聲激勵(lì)一個(gè)差分系統(tǒng)方程而產(chǎn)生。首先利用Matlab產(chǎn)生該信號(hào)源：n=1:1:1024;wn=randn(1,length(n));%產(chǎn)生隨機(jī)單位白噪聲序列sn=filter([1],[1,-2.2377,3.7476,-2.6293,0.9224],wn);subplot;plot(sn,'r','LineWidth',1.4);title('s(n)','fontsize',12);產(chǎn)生的信號(hào)源如下圖3-1所示。對(duì)信號(hào)源信號(hào)求FFT變換，得到信號(hào)源的頻譜圖，作為估計(jì)算法的參考標(biāo)準(zhǔn)，如圖3-2所示。圖3-1Matlab產(chǎn)生的信號(hào)源圖3-2信號(hào)源的頻譜圖1）經(jīng)典譜估計(jì)下面首先采用直接法（周期圖法）進(jìn)行功率譜估計(jì)，程序如下：%nfft=128window=boxcar(length(sn));nfft=128;[Pxx,f]=periodogram(sn,window,nfft);subplot(211);plot(f,10*log10(Pxx));xlabel('頻率');ylabel('相對(duì)功率譜密度dB');title('周期圖法點(diǎn)數(shù)128');%nfft=1024window=boxcar(length(sn));%矩形窗nfft=1024;[Pxx,f]=periodogram(sn,window,nfft);subplot(212);;plot(f,10*log10(Pxx));xlabel('頻率');ylabel('相對(duì)功率譜密度dB');title('周期圖法點(diǎn)數(shù)1024');改變采樣點(diǎn)數(shù)，獲得的周期圖法譜估計(jì)如下圖所示：圖3-3不同序列長(zhǎng)度周期法譜估計(jì)運(yùn)用周期圖法通過(guò)對(duì)采樣數(shù)據(jù)長(zhǎng)度的調(diào)整，即在采樣數(shù)據(jù)為128、1024時(shí)比較周期圖法的效果，可知采集數(shù)據(jù)序列長(zhǎng)度越大，譜分辨率越高，但譜曲線起伏也越大。因此需要改進(jìn)，這里的改進(jìn)指的是數(shù)據(jù)方差的改進(jìn)。當(dāng)前主要有兩種改進(jìn)方法：一種是周期圖的平滑，即采用間接法估計(jì)功率譜，即BT法；另一種是平均法，它將長(zhǎng)度為N的數(shù)據(jù)X(n)分成L段，分別求出每段的功率譜，然后加以平均，如Barlett法和Welch法。這里我們就采用以上3種方法對(duì)直接法進(jìn)行改進(jìn)。Matlab程序如下：%BT法nfft=1024;%以點(diǎn)數(shù)1024為例cxn=xcorr(sn,'unbiased');%求序列自相關(guān)函數(shù)P1=fft(cxn,nfft);subplot(311);plot(f,10*log10(P1(1:257)),'r','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('BT法','fontsize',12)%Bartlett法window=boxcar(length(sn));%矩形窗noverlap=0;%數(shù)據(jù)無(wú)重疊p=0.9;%置信概率p2=psd(sn,nfft,1,window,noverlap,p);%計(jì)算PSDsubplot(312);plot(f,10*log10(p2(1:257)),'r','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('Bartlett法','fontsize',12);%Welch法window=boxcar(length(sn));%矩形窗noverlap=0;%數(shù)據(jù)無(wú)重疊range='onesided';%頻率間隔為[0,1/2]，只計(jì)算一半頻率[p3,w]=pwelch(sn,window,noverlap,nfft,range);subplot(313);plot(f,10*log10(p3),'b','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('Welch法','fontsize',12);各個(gè)經(jīng)典功率譜估計(jì)方法的結(jié)果如圖3-4所示。圖3-4經(jīng)典譜估計(jì)方法的估計(jì)結(jié)果為了了解使用不同種窗函數(shù)對(duì)功率譜的影響，以Welth法為例，進(jìn)行了應(yīng)用各種常見窗函數(shù)的仿真實(shí)驗(yàn)。程序如下:%矩形窗window=boxcar(length(sn));%矩形窗noverlap=0;%數(shù)據(jù)無(wú)重疊range='onesided';%頻率間隔為[0,1/2]，只計(jì)算一半頻率[p3,w]=pwelch(sn,window,noverlap,nfft,range);subplot(221);plot(f,10*log10(p3),'b','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('矩形窗','fontsize',12);%漢明窗window=hamming(length(sn));%矩形窗noverlap=0;%數(shù)據(jù)無(wú)重疊range='onesided';%頻率間隔為[0,1/2]，只計(jì)算一半頻率[p3,w]=pwelch(sn,window,noverlap,nfft,range);subplot(222);plot(f,10*log10(p3),'b','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('漢明窗','fontsize',12);%布萊克曼窗window=blackman(length(sn));%矩形窗noverlap=0;%數(shù)據(jù)無(wú)重疊range='onesided';%頻率間隔為[0,1/2]，只計(jì)算一半頻率[p3,w]=pwelch(sn,window,noverlap,nfft,range);subplot(223);plot(f,10*log10(p3),'b','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('布萊克曼窗','fontsize',12);%三角窗window=triang(length(sn));%矩形窗noverlap=0;%數(shù)據(jù)無(wú)重疊range='onesided';%頻率間隔為[0,1/2]，只計(jì)算一半頻率[p3,w]=pwelch(sn,window,noverlap,nfft,range);subplot(224);plot(f,10*log10(p3),'b','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('三角窗','fontsize',12);估計(jì)結(jié)果如圖3-5所示：圖3-5不同窗函數(shù)對(duì)功率譜估計(jì)結(jié)果的影響綜合以上的波形，可以對(duì)經(jīng)典功率譜估計(jì)作如下總結(jié)：經(jīng)典功率譜估計(jì)，不論是周期圖法還是BT法都可以用FFT快速計(jì)算，且計(jì)算量較小，物理概念明確，因而仍是目前較常用的譜估計(jì)方法。功率譜的分辨率正比于數(shù)據(jù)長(zhǎng)度，當(dāng)數(shù)據(jù)長(zhǎng)度較大時(shí)，譜的曲線起伏加劇。BT法是對(duì)周期圖法的改進(jìn)，但其譜的平滑是以犧牲分辨率為代價(jià)的。Welth方法是對(duì)BT法德改進(jìn)，它能使用更多的窗函數(shù)，并能更好的控制功率譜密度的估計(jì)的方差特性。在實(shí)際應(yīng)用中，在方差、偏差和分辨率之間存在著矛盾，只能根據(jù)需要作折中處理。2）現(xiàn)代譜估計(jì)下面，根據(jù)AR模型法和burg算法對(duì)信號(hào)源進(jìn)行現(xiàn)代功率譜估計(jì)。Matlab程序如下：%現(xiàn)代功率譜分析%n=1:1:1024;wn=randn(1,length(n));%產(chǎn)生隨機(jī)單位白噪聲序列sn=filter([1],[1,-2.7377,3.7476,-2.6293,0.9224],wn);[Pxx,f]=periodogram(sn,window,nfft);nfft=1024;%AR模型法[p4,w]=pyulear(sn,100,nfft);figure(3);subplot(2,1,1);plot(f,10*log10(abs(p4)),'b','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('AR模型法','fontsize',12);%burg算法[p5,w]=pburg(sn,100,nfft);subplot(2,1,2);plot(f,10*log10(abs(p5)),'b','LineWidth',1.4);xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('burg算法','fontsize',12);估計(jì)結(jié)果如圖3-6所示：圖3-6現(xiàn)代功率譜估計(jì)方法的估計(jì)結(jié)果為了分析現(xiàn)代功率譜中模型階數(shù)對(duì)譜估計(jì)的影響，現(xiàn)利用AR模型法，改變模型階數(shù)。分別用10、20、100作為模型階數(shù)在采樣點(diǎn)數(shù)分別為128和1024下進(jìn)行分析，結(jié)果如下：%=現(xiàn)代功率譜分析=%%AR模型法n=1:1:1024;wn=randn(1,length(n));%產(chǎn)生隨機(jī)單位白噪聲序列sn=filter([1],[1,-2.7377,3.7476,-2.6293,0.9224],wn);[P1,w4]=pmcov(sn,10,128);[P2,w5]=pmcov(sn,50,128);[P3,w6]=pmcov(sn,100,128);[P4,w1]=pmcov(sn,10,1024);[P5,w2]=pmcov(sn,50,1024);[P6,w3]=pmcov(sn,100,1024);subplot(3,2,1);plot(w4/pi,10*log10(P1));xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('12810階AR模型功率譜圖');gridon;subplot(3,2,3);plot(w5/pi,10*log10(P2));xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('12850階AR模型功率譜圖');gridon;subplot(3,2,5);plot(w6/pi,10*log10(P3));xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('128100階AR模型功率譜圖');gridon;subplot(3,2,2);plot(w1/pi,10*log10(P4));xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('102410階AR模型功率譜圖');gridon;subplot(3,2,4);plot(w2/pi,10*log10(P5));xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('102450階AR模型功率譜圖');gridon;subplot(3,2,6);plot(w3/pi,10*log10(P6));xlabel('頻率','fontsize',12);ylabel('功率譜/dB','fontsize',12);title('1024100階AR模型功率譜圖');gridon;圖3-7不同階數(shù)下AR模型法進(jìn)行功率譜估計(jì)可以看到，在相同采樣數(shù)據(jù)長(zhǎng)度下階數(shù)過(guò)小會(huì)導(dǎo)致頻率成分的辨識(shí)出現(xiàn)偏差，階數(shù)過(guò)大會(huì)使隨機(jī)起伏變化密集和加劇，方差增大。在50階的情況下效果較為理想，但階數(shù)越高的情況下計(jì)算量同時(shí)也在增加，應(yīng)該適當(dāng)?shù)倪M(jìn)行選擇。在相同AR模型階數(shù)的條件下，采樣數(shù)據(jù)越多，頻率的分辨率也越高，但會(huì)出現(xiàn)波動(dòng)。由上述分析可以看出，AR模型法、burg算法估計(jì)出的功率譜相對(duì)經(jīng)典的功率譜估計(jì)方法得出的結(jié)果方差較小，頻率的準(zhǔn)確性也相對(duì)較高?，F(xiàn)代譜估計(jì)的分辨率可以不受樣本長(zhǎng)度的限制。這是因?yàn)閷?duì)于給定的N點(diǎn)有限長(zhǎng)序列x(n),雖然其估計(jì)出的自相關(guān)函數(shù)也是有限長(zhǎng)的,但是現(xiàn)代譜估計(jì)的一些隱含著數(shù)據(jù)和自相關(guān)函數(shù)的外推,使其可能的長(zhǎng)度超過(guò)給定的長(zhǎng)度,不像經(jīng)典譜估計(jì)那樣受窗函數(shù)的影響。與經(jīng)典譜估計(jì)結(jié)果比較，現(xiàn)代譜估計(jì)在分辨率和方差上都優(yōu)于經(jīng)典譜，而且現(xiàn)代譜線要平滑得多，且在采樣數(shù)據(jù)短的情況下也非常適用。四、自適應(yīng)濾波器（LMS和RLS）題目：請(qǐng)簡(jiǎn)述一下自適應(yīng)濾波器（LMS和RLS），并任選下面的一題，（2）、（3）數(shù)據(jù)可從網(wǎng)上下載（MIT-BIH數(shù)據(jù)庫(kù)：/physiobank/database/mitdb），分析步長(zhǎng)及指數(shù)加權(quán)因子對(duì)濾波結(jié)果的影響。同時(shí)給出迭代次數(shù)與濾波器系數(shù)，迭代次數(shù)與均方誤差之間的關(guān)系曲線。（13分）（1）利用自適應(yīng)濾波器從寬帶信號(hào)中提取單頻信號(hào)。設(shè)，是寬頻信號(hào)，A,B,,任選，要求提取兩個(gè)單頻信號(hào)；同時(shí)對(duì)不同的信噪比的情況比較一下。（2）消除心電圖的電源干擾（3）胎兒心電監(jiān)護(hù)濾波技術(shù)是指從復(fù)雜的含有噪聲的信號(hào)中提取出有用的所需信號(hào)，而濾波器是一個(gè)選頻系統(tǒng)，讓有用的所需信號(hào)以盡量小的衰減通過(guò)，其實(shí)我們可以把大部分的信號(hào)處理系統(tǒng)看作濾波器系統(tǒng)。濾波器可以分為經(jīng)典濾波器和現(xiàn)代濾波器，也可以分為模擬濾波器和數(shù)字濾波器，而自適應(yīng)濾波器是屬于現(xiàn)代濾波器，它是相對(duì)固定濾波器而言。本題首先介紹自適應(yīng)濾波器的概念和工作原理，然后對(duì)LMS和RLS兩種算法進(jìn)行介紹，隨后以胎兒心電信號(hào)（FECG）的提?。}3）為例來(lái)介紹自適應(yīng)濾波器的應(yīng)用，并分析步長(zhǎng)及指數(shù)加權(quán)因子對(duì)濾波結(jié)果的影響。1.自適應(yīng)濾波器簡(jiǎn)介自適應(yīng)濾波器是指能夠根據(jù)輸入信號(hào)自動(dòng)調(diào)整性能進(jìn)行數(shù)字信號(hào)處理的數(shù)字濾波器。自適應(yīng)濾波器在許多領(lǐng)域都有應(yīng)用，如自適應(yīng)天線系統(tǒng)、數(shù)字通信接收機(jī)、自適應(yīng)噪聲對(duì)消技術(shù)、系統(tǒng)建模等，本題的例子便是應(yīng)用自適應(yīng)濾波器進(jìn)行自適應(yīng)噪聲對(duì)消技術(shù)。自適應(yīng)濾波器通常由兩個(gè)部分組成，一部分是數(shù)字濾波器，用來(lái)完成希望的信號(hào)處理；另一部分是自適應(yīng)算法，用來(lái)調(diào)節(jié)濾波器系數(shù)（或加權(quán)），如圖3-1所示。-參數(shù)可調(diào)數(shù)字濾波器-參數(shù)可調(diào)數(shù)字濾波器自適應(yīng)濾波器算法+∑d(n)x(n)e(n)y(n)圖4-1自適應(yīng)濾波器的結(jié)構(gòu)框圖如圖，參數(shù)可調(diào)數(shù)字濾波器可以是FIR數(shù)字濾波器或IIR數(shù)字濾波器，也可以是格型數(shù)字濾波器。輸入信號(hào)x(n)通過(guò)參數(shù)可調(diào)數(shù)字濾波器后產(chǎn)生輸出信號(hào)(或響應(yīng))y(n)，將其與參考信號(hào)(或稱期望響應(yīng))d(n)進(jìn)行比較，形成誤差信號(hào)e(n)。e(n)(有時(shí)還要利用x(n))通過(guò)某種自適應(yīng)算法對(duì)濾波器參數(shù)進(jìn)行調(diào)整，最終使e(n)的均方值最小。根據(jù)算法的不同，可以將自適應(yīng)濾波器分為：遞推最小二乘（RLS）濾波器、最小均方差（LMS）濾波器、格型濾波器、無(wú)限沖激響應(yīng)（IIR）濾波器等。（1）LMS濾波器最小均方誤差(LMS)算法是由Widrow和Hoff提出的，具有計(jì)算量小、易于實(shí)現(xiàn)等優(yōu)點(diǎn)在實(shí)踐中廣泛應(yīng)用。LMS算法的基本思想：調(diào)整濾波器自身的參數(shù)，使濾波器的輸出信號(hào)與期望輸出信號(hào)之間的均方誤差最小，系統(tǒng)輸出為有用信號(hào)的最佳估計(jì)。實(shí)質(zhì)上LMS可以看成是一種隨機(jī)梯度或者隨機(jī)逼近算法，可以寫成如下的基本迭代方程：其中μ為步長(zhǎng)因子，是控制穩(wěn)定性和收斂速度的參量，X(k)表示第k時(shí)刻參考信號(hào)矢量，，k為迭代次數(shù)，M為濾波器的階數(shù)。d(k)表示第k時(shí)刻的輸入信號(hào)矢量，y(k)、e(k)分別表示第k時(shí)刻的輸出信號(hào)與輸出誤差，W(k)表示k時(shí)刻權(quán)系數(shù)矢量，W(k)=[W(k，0)，W(k，1)…W(k，M-1)]。從上式可以看出其結(jié)構(gòu)簡(jiǎn)單、計(jì)算量小且穩(wěn)定性好等優(yōu)點(diǎn)，但固定步長(zhǎng)的LMS算法在收斂速度、跟蹤速率及權(quán)失調(diào)噪聲之間的要求是相互矛盾的，為了克服其缺點(diǎn)，人們提出了各種變步長(zhǎng)的LMS改進(jìn)算法，主要是采用減小均方誤差或者以某種規(guī)則基于時(shí)變步長(zhǎng)因子跟蹤信號(hào)的時(shí)變，其中有正規(guī)LMS算法(NLMS)、梯度自適應(yīng)步長(zhǎng)算法、自動(dòng)增益控制自適應(yīng)算法、符號(hào)-誤差LMS算法、符號(hào)-數(shù)據(jù)LMS算法、數(shù)據(jù)復(fù)用LMS算法等。自適應(yīng)濾波器收斂的條件是。其中是輸入信號(hào)的自相關(guān)矩陣R的最大特征值。μ的選取必須在收斂速度和失調(diào)之間取得較好的折中，既要具有較快的收斂速度，又要使穩(wěn)態(tài)誤差最小。它控制了算法穩(wěn)定性和自適應(yīng)速度，如果μ很小，算法的自適應(yīng)速度會(huì)很慢；如果μ很大，算法會(huì)變得不穩(wěn)定。（2）RLS濾波器RLS（RecursiveLeast-Square）算法是通過(guò)矩陣求逆引理而導(dǎo)出的，它對(duì)輸入信號(hào)的自相關(guān)矩陣的逆矩陣進(jìn)行遞推估計(jì)更新，其迭代公式如下：上式中，n是遞推次數(shù)，輸入為x(n)，P(n)為信號(hào)自相關(guān)矩陣的逆，g(n)為增益矢量，誤差為e(n)，期望輸出為d(n)，權(quán)重因子為λ。RLS算法能實(shí)現(xiàn)快速收斂，且其收斂速率不隨輸入向量x(n)集平均相關(guān)矩陣R特征值的擴(kuò)散度而改變。當(dāng)工作于時(shí)變環(huán)境中時(shí)，這類算法具有極好的性能，但其實(shí)現(xiàn)都以增加計(jì)算復(fù)雜度和穩(wěn)定性為代價(jià)。2.胎兒心電監(jiān)護(hù)胎兒心電圖(FetalElectrocardiogram，F(xiàn)ECG)是研究胎兒心臟電生理活動(dòng)的一項(xiàng)客觀指標(biāo)，反映了胎兒在孕期中的成長(zhǎng)和健康狀況。胎兒心電一般是從母親腹部采集獲得，不可避免地會(huì)含有許多噪聲。首先，加在母體上的工頻干擾；其次，母體心電（MaternalElectrocardiogram，MECG），而且與胎兒心電相比，它會(huì)強(qiáng)得多；最后，還有夾雜的除心電以外的其他生物電信號(hào)。雖然采集獲得的信號(hào)含有許多種類的噪聲，但經(jīng)前期處理均可濾掉，或者得到很好的抑制，最后剩下的主要是母體心電。因此，我們采用自適應(yīng)濾波消噪算法把混合信號(hào)中的母體心電信號(hào)去除。自適應(yīng)噪聲抵消濾波器的結(jié)構(gòu)框圖如圖3-2所示：圖4-2自適應(yīng)噪聲抵消濾波器上圖中抵消器的“信號(hào)源”為s+n0，其中s為有用信號(hào)，n0為一個(gè)與信號(hào)s不相關(guān)的噪聲，抵消器的“噪聲源”為噪聲n1，n1與信號(hào)s不相關(guān)，卻以某種未知的方式與噪聲n0相關(guān)，由圖可以看出噪聲n1經(jīng)自適應(yīng)濾波器輸出y，再?gòu)脑驾斎雜+n0中減去該輸出，產(chǎn)生了系統(tǒng)的輸出e=s+n0-y。本題中，信號(hào)源為母體腹部心電信號(hào)，噪聲源為母體胸部心電信號(hào)MECG，系統(tǒng)輸出e為胎兒心電信號(hào)FECG。數(shù)據(jù)是從MIT-BIH數(shù)據(jù)庫(kù)下載，進(jìn)入網(wǎng)址/cgi-bin/atm/ATM，在左側(cè)Input選項(xiàng)里選擇Database為胎兒心電圖選項(xiàng)，Record項(xiàng)任選一個(gè)，Signals項(xiàng)選擇為all，右側(cè)Toolbox選擇為Exportsignalsas.mat。如圖4.3所示。圖4-3數(shù)據(jù)下載設(shè)置完上面步驟后，會(huì)找到網(wǎng)頁(yè)里有一個(gè)下載項(xiàng)，ecgca323_edfm.mat，點(diǎn)擊下載。這就是對(duì)我們有用的數(shù)據(jù)，MATLAB打開后為一個(gè)6*10000的數(shù)據(jù)，我們只要對(duì)第一行和第三行的數(shù)據(jù)進(jìn)行處理，分別為母親胸部和腹部的數(shù)據(jù)信號(hào)。我們分別用LMS算法和RLS算法得到胎兒心電。胎兒心電監(jiān)護(hù)原理：監(jiān)測(cè)胎兒心電，需要將母親的心音及背景干擾去除。母親的心音強(qiáng)度通常是胎兒心音的2倍到10倍，肌肉活動(dòng)及胎兒運(yùn)動(dòng)產(chǎn)生的背景干擾也常常大于胎兒心音的強(qiáng)度。采用自適應(yīng)對(duì)消的方法來(lái)監(jiān)護(hù)胎兒心電圖，將母親胸部的信號(hào)作為參考輸入，主要包括母親心音和背景噪聲，也有胎兒的心音信號(hào)，如圖4.4所示（數(shù)據(jù)為MIT-BIH數(shù)據(jù)庫(kù)下載）；從母親腹部取出的信號(hào)作為原始輸入，包括母親心音、胎兒心音和背景干擾，如圖4.5所示（數(shù)據(jù)為下載數(shù)據(jù)），可以看出，這個(gè)系統(tǒng)是一個(gè)有信號(hào)分量泄露到參考輸入端的自適應(yīng)對(duì)消系統(tǒng)。我們希望的結(jié)果是對(duì)消后母親的心電信號(hào)的強(qiáng)度小于胎兒的心電信號(hào)強(qiáng)度，且胎兒的心電信號(hào)要比較清楚。程序如下：load('ecgca102_edfm.mat');E=val(1,:)./10000;E1=val(3,:)./10000;thorax=E(1:3000);%胸部信號(hào)，參考信號(hào)dn=thorax.';abdomen=E1(1:3000);%腹部信號(hào)，輸入信號(hào)xn=abdomen.';figure(1);plot((1:3000),dn);xlabel('采樣點(diǎn)數(shù)');ylabel('幅度');title('母親心電圖');figure(2);plot((1:3000),xn);xlabel('采樣點(diǎn)數(shù)');ylabel('幅度');title('胎兒心電圖');母親胎兒心電圖：圖4-4母親胸部的心電圖，作為期望信號(hào)圖4-5胎兒的心電圖，作為輸入信號(hào)3.1用LMS算法得到胎兒心電圖比較圖4.4和圖4.5可以發(fā)現(xiàn)，胎兒的心電信號(hào)緊隨著母親的心電信號(hào)，這里采樣點(diǎn)數(shù)選取了3000個(gè)，階數(shù)采用20階（即權(quán)系數(shù)的個(gè)數(shù)），步長(zhǎng)因子為改變的參數(shù)，分別設(shè)置為0.005，0.05，0.5。另外由于設(shè)置的階數(shù)為20階，所以權(quán)矢量就有20個(gè)，如果查看每一個(gè)與迭代次數(shù)的關(guān)系曲線，就要占很大篇幅，由于權(quán)系數(shù)（濾波器系數(shù)）更新的方程是一樣的，所以它們的規(guī)律是相似的，這里只顯示前三個(gè)權(quán)系數(shù)與迭代次數(shù)的關(guān)系曲線圖。其用LMS算法得到胎兒心電圖的MATLAB源代碼如下（此處以步長(zhǎng)為0.005為例）：load('ecgca102_edfm.mat');E=val(1,:)./10000;E1=val(3,:)./10000;thorax=E(1:3000);%胸部信號(hào)，參考信號(hào)dn=thorax.';abdomen=E1(1:3000);%腹部信號(hào)，輸入信號(hào)xn=abdomen.';M=20;%階數(shù)itr=3000;%迭代次數(shù)en=zeros(itr,1);%誤差序列W=zeros(M,itr);%權(quán)系數(shù)mu=0.5;%步長(zhǎng)因子s=zeros(itr,1);%均方誤差%迭代計(jì)算fork=M:itrx=xn(k:-1:k-M+1);%濾波器M個(gè)抽頭的輸入y=W(:,k-1).'*x;%濾波器輸出en(k)=dn(k)-y;%第k次迭代的誤差W(:,k)=W(:,k-1)+2*mu*en(k)*x;%權(quán)系數(shù)更新s(k)=s(k-1)+abs(en(k)).^2;c(k)=sqrt(s(k)/k);%均方誤差end%求最優(yōu)時(shí)濾波器的輸出序列yn=inf*ones(size(xn));fork=M:length(xn)x=xn(k:-1:k-M+1);yn(k)=W(:,end).'*x;endfigure;subplot(311),plot((1:3000),xn-yn,'r',(1:3000),dn,'b');legend('胎兒心電圖','母親心電圖');xlabel('采樣點(diǎn)數(shù)');ylabel('幅度');subplot(312);plot(c);xlabel('迭代次數(shù)');ylabel('均方誤差');subplot(313);plot(W(1,:),'r');holdon;plot(W(2,:));plot(W(3,:),'k');xlabel('迭代次數(shù)');ylabel('濾波器系數(shù)');legend('權(quán)系數(shù)1','權(quán)系數(shù)2','權(quán)系數(shù)3');修改步長(zhǎng)參數(shù)，即mu，分別取0.005，0.05，0.5。運(yùn)行程序，得到的仿真圖形分別如圖4-6、圖4-7、圖4-8所示。圖4-6步長(zhǎng)為0.005時(shí)的LMS算法圖4-7步長(zhǎng)為0.05時(shí)的LMS算法圖4-8步長(zhǎng)為0.5時(shí)的LMS算法比較上面3幅圖，可以很明顯的發(fā)現(xiàn)，隨著步長(zhǎng)因子的增加，LMS算法的收斂速度，失調(diào)量也隨著變化，當(dāng)步長(zhǎng)因子增大，系統(tǒng)的收斂速度變快，但是相應(yīng)的失調(diào)量也越來(lái)越大（如圖4.7所示），當(dāng)步長(zhǎng)因子大到某個(gè)數(shù)值時(shí)，系統(tǒng)不再收斂，也就是系統(tǒng)不再穩(wěn)定（如圖4.8所示）；但是這并不意味著越小越好，因?yàn)椴介L(zhǎng)因子越小，則系統(tǒng)收斂速度變慢，這在實(shí)際應(yīng)用中也是不容許的。因此參數(shù)的選擇應(yīng)從整個(gè)系統(tǒng)要求出發(fā),在滿足精度要求的前提下,盡量減少自適應(yīng)時(shí)間。3.2用RLS算法得到胎兒心電圖遞歸最小二乘（RLS）算法，其實(shí)現(xiàn)步驟有點(diǎn)類似卡爾曼濾波器，也是通過(guò)幾個(gè)遞歸公式來(lái)循環(huán)。采樣點(diǎn)數(shù)3000，濾波器階數(shù)設(shè)置為80，初始化條件設(shè)置為：是小的正數(shù)，取0.001。。選擇改變的參數(shù)為遺忘因子，分別取0.97、0.7、0.5，觀察濾波效果。同樣，由于濾波器系數(shù)過(guò)多，這里只觀測(cè)前三個(gè)濾波器系數(shù)與迭代次數(shù)的曲線。RLS算法得到胎兒心電圖的MATLAB源代碼如下（此處以遺忘因子=0.97為例）：%迭代計(jì)算fork=M:itrx=xn(k:-1:k-M+1);%濾波器M個(gè)抽頭的輸入K(:,k)=(Rxx*x)./(la+x.'*Rxx*x);%更新濾波器增益e(k)=dn(k)-W(:,k-1).'*x;%誤差序列W(:,k)=W(:,k-1)+K(:,k)*e(k);%更新權(quán)系數(shù)Rxx=(Rxx-K(:,k)*x.'*Rxx)/la;%更新自相關(guān)矩陣逆更新y(k)=W(:,k).'*x;%濾波器輸出e1(k)=dn(k)-y(k);s(k)=s(k-1)+abs(e1(k)).^2;c(k)=sqrt(s(k)/k);%均方誤差endfigure;subplot(311),plot((1:3000),xn-y,'r',(1:3000),dn,'b');legend('胎兒心電圖','母親心電圖');xlabel('采樣點(diǎn)數(shù)');ylabel('幅度');subplot(312);plot(c);xlabel('迭代次數(shù)');ylabel('均方誤差');subplot(313);plot(W(1,:),'r');holdon;plot(W(2,:));plot(W(3,:),'k');xlabel('迭代次數(shù)');ylabel('濾波器系數(shù)');legend('權(quán)系數(shù)1','權(quán)系數(shù)2','權(quán)系數(shù)3');仿真波形分別如圖4-9、圖4-10、圖4-11所示。圖4-9遺忘因子為0.97時(shí)的RLS算法圖4-10遺忘因子為0.7時(shí)的RLS算法圖4-11遺忘因子為0.5時(shí)的RLS算法比較上面3幅圖，根據(jù)均方誤差與迭代次數(shù)的關(guān)系曲線圖可以發(fā)現(xiàn)當(dāng)=0.9時(shí)，系統(tǒng)在1800點(diǎn)處收斂；而當(dāng)=0.5時(shí)，系統(tǒng)在500點(diǎn)就收斂了，所以得出結(jié)論：遺忘因子越小，系統(tǒng)的跟蹤能力越強(qiáng)，但同時(shí)可以發(fā)現(xiàn)其對(duì)噪聲也越敏感；相反，遺忘因子值越大，系統(tǒng)跟蹤能力減弱，但對(duì)噪聲不敏感，收斂時(shí)估計(jì)誤差也越小。所以，在用普通遞推RLS算法時(shí)，一定要對(duì)有個(gè)準(zhǔn)確的取值，才能保證系統(tǒng)性能的最佳狀態(tài)，所以一般文獻(xiàn)要求的取值范圍為。五、維納濾波器題目：觀測(cè)信號(hào)為,其中有信號(hào)是為恒量平穩(wěn)序列，其統(tǒng)計(jì)特征已求得為噪聲是零均值白噪聲，且與有用信號(hào)不相關(guān)，即求維納濾波器？（10分）維納濾波器實(shí)際上就是在最小均方誤差條件下探確定濾波器的沖擊響應(yīng)或系統(tǒng)函數(shù)，也就是求解維納-霍夫方程的問題。因噪聲是零均值白噪聲，且與有用信號(hào)不相關(guān),根據(jù)已知條件及可求得（為恒量平穩(wěn)信號(hào)）則維納霍夫方程為。以一階維納濾波器為例，則維納霍夫方程可寫為求解該方程得維納濾波器為：均方誤差六、FIR維納濾波器題目：假定下圖所需的信號(hào)的d(n)是一個(gè)正弦序列d(n)=，噪聲序列和都是AR(1)過(guò)程，分別由如下的一階差分方程產(chǎn)生：式中,是零均值，單位方差的白噪聲，與d(n)不相關(guān)，取不同的值，如0.1，0.9，試采,2個(gè)不同階數(shù)的FIR維納濾波器對(duì)信號(hào)進(jìn)行濾波。并分析其結(jié)果。(12分)維納濾波器（Wienerfilter）是由數(shù)學(xué)家維納（RorbertWiener）提出的一種以最小平方為最優(yōu)準(zhǔn)則的線性濾波器。其數(shù)字系統(tǒng)描述如下圖所示，它是要從噪聲觀測(cè)x(n)=d(n)+v(n)中恢復(fù)信號(hào)d(n)。圖6-1維納濾波器的一般框圖題目中所示系統(tǒng)的Wiener-Hopf方程可以推到如下：將v2(n)作為維納濾波器的輸入，以估計(jì)噪聲v1(n)，則Wiener-Hopf方程為：，其中是v2(n)的自相關(guān)矩陣，是期望信號(hào)v1(n)和濾波器輸入v2(n)之間的互相關(guān)矢量。對(duì)互相關(guān)量，可推導(dǎo)rv1v2=E{v1(n)v2*(n-k)}=E{[x(n)-d(n)]v2*(n-k)}=E{x(n)v2*(n-k)}-E{d(n)v2*(n-k)}題目中v2(n)與d(n)不相關(guān)，則上式第二項(xiàng)為零，上式簡(jiǎn)化為rv1v2=E{x(n)v2*(n-k)}因此，Wiener-Hopf方程變?yōu)椋篟v2ω=rxv2式中Rv2是v2(n)的自相關(guān)矩陣，ω是維納濾波器系數(shù)矢量，rxv2是x(n)和v2(n)的互相關(guān)矢量。因此，要求解出維納濾波器，需要先求出維納濾波器系數(shù)矢量ω。根據(jù)上述分析，在matlab中，首先利用白噪聲產(chǎn)生v1(n)和v2（n）序列，之后求解出v2(n)的自相關(guān)矩陣和x(n)和v2(n)的互相關(guān)矢量，然后求出維納濾波器系數(shù)矢量ω。再利用已求解的維納濾波器進(jìn)行題目所示的噪聲濾除工作。取ω0=0.05，具體求解過(guò)程如下程序所示：clear;clc;N=512;n=1:1:N;en=randn(1,length(n));%白噪聲d=sin(0.05*n*pi+pi/3);k=2;%維納濾波器的階數(shù)v1=filter([1],[1-0.1],en);v2=filter([1],[10.6],en);x=d+v1;Rv2=xcorr(v2,'unbiased');%計(jì)算序列的自相關(guān)函數(shù)Rv1v2=xcorr(v1,v2,'unbiased');%計(jì)算序列的互相關(guān)函數(shù)fori=1:k;%構(gòu)造自相關(guān)矩陣和互相關(guān)矩陣,Rv21(i,:)=Rv2(N-i+1:N-i+k);Rv1v21(i,:)=Rv1v2(N+i-1);end;h=inv(Rv21)*Rv1v21;v11=filter(h',1,v2);v1_1=filter(h',1,v2);e=x-v11;subplot(4,1,1),plot(n,d,'b'),title('輸入信號(hào)波形')subplot(4,1,2),plot(n,x,'b'),title('信號(hào)加白噪聲')subplot(4,1,3),plot(n,v1_1,'r'),title('白噪聲估計(jì)')subplot(4,1,4),plot(n,e,'g'),title('恢復(fù)后信號(hào)波形')e1=e-d;a=sum(e1.^2)當(dāng)=0.1，維納濾波器階數(shù)為2階時(shí)，濾波結(jié)果如下圖6-2所示。 圖6-2=0.1，2階維納濾波器濾波效果 當(dāng)=0.1，維納濾波器階數(shù)為5階時(shí)，濾波結(jié)果如下圖6-3所示。圖6-3=0.1，5階維納濾波器濾波效果當(dāng)=0.9，維納濾波器階數(shù)為20階時(shí)，濾波結(jié)果如下圖6-4所示。圖6-4=0.9，20階維納濾波器濾波效果當(dāng)=0.9，維納濾波器階數(shù)為50階時(shí)，濾波結(jié)果如下圖6-5所示。圖6-5=0.9，50階維納濾波器濾波效果根據(jù)上述分析和結(jié)果圖進(jìn)行對(duì)比可以得出，當(dāng)=0.1時(shí)，維納濾波器階數(shù)為5時(shí)，維納濾波器濾波得到的信號(hào)幅值與相位與原信號(hào)就有較好的一致性。但當(dāng)當(dāng)=0.9時(shí)，維納濾波器階數(shù)為20時(shí)，維納濾波器濾波得到的信號(hào)幅值與相位與原信號(hào)都沒有太好的一致性，只有階數(shù)為50時(shí)，濾波得到的信號(hào)幅值與相位與原信號(hào)才有較好的一致性?？梢钥闯鼍S納濾波器階數(shù)越高，濾波效果更好一些，在選擇50階濾波時(shí)，曲線平滑且接近原波形曲線，但實(shí)際情況下，并不是階數(shù)越高越好，當(dāng)濾波階數(shù)過(guò)大時(shí)，維納濾波器在仿真計(jì)算過(guò)程中會(huì)引入舍入誤差，當(dāng)舍入誤差逐漸積累，會(huì)導(dǎo)致高計(jì)算量的高階情況下均方誤差增大，引起濾波不理想。但當(dāng)超過(guò)一定階數(shù)時(shí)，階數(shù)的增加對(duì)濾波效果的改善已經(jīng)基本可以忽略，此時(shí)單純靠提高階數(shù)的辦法不能改進(jìn)濾波效果了。說(shuō)明此時(shí)維納濾波器的效能已經(jīng)達(dá)到最佳，若要再改善效果，必須選用結(jié)構(gòu)更先進(jìn)的濾波器。七、AR/MA濾波器題目：用雷達(dá)測(cè)量地球和月球之間的距離，測(cè)量過(guò)程用下列方程描述其中是均值為零，方差為的白噪聲序列，它表示測(cè)量誤差。為了提高測(cè)量精度，現(xiàn)采用以下兩種濾波器分別對(duì)進(jìn)行處理，試比較其方差的大小。濾波器１：濾波器２：式中。（10分）由于測(cè)量過(guò)程為測(cè)量真值與零均值、方差為的白噪聲序列，因此，在濾波前，的均值為，方差為白噪聲序列方差，即：本題是對(duì)測(cè)量過(guò)程經(jīng)過(guò)兩種濾波器后信號(hào)序列的方差的比較。題中所給出的濾波器1和濾波器2，分別為AR和MA模型，兩個(gè)濾波器中值未確定，但給出了相應(yīng)范圍，即，為了對(duì)這個(gè)范圍內(nèi)濾波器1和濾波器2的方差進(jìn)行比較，在MATLAB程序中，需要設(shè)定一個(gè)參數(shù)，其范圍為0至1，采樣頻率為1000，并將其存入一個(gè)數(shù)組中，以比較不同值下的方差大小，即兩種濾波器對(duì)原信號(hào)的濾波效果。另外為了得出白噪聲序列方差對(duì)濾波后方差的影響，可以通過(guò)輸入白噪聲不同的輸入值獲得，本題取方差為0.04觀測(cè)其效果，Matlab程序見如下：clear;clc;N=1000;var0=input('輸入白噪聲方差：');init=2055615866;randn('seed',init);w=sqrt(var0)*randn(1,N);d=3800;%距離x=d+w;%輸入信號(hào)y1(1)=d;y2(1)=d;fori=1:Na=i/N;forj=2:Ny1(j)=a*y1(j-1)+(1-a)*x(j);y2(j)=a*x(j)+(1-a)*x(j-1);endv1(i)=var(x);%輸入信號(hào)方差v2(i)=var(y1);%濾波器1后的方差v3(i)=var(y2);%濾波器2后的方差endt=(1:N)/N;plot(t,v1,'r',t,v2,'b',t,v3,'g');xlabel('a值');ylabel('方差');legend('未濾波','濾波器1','濾波器2');運(yùn)行結(jié)果如下圖所示：圖7-1=0.04時(shí)不同值的方差比較如圖7-1可知，運(yùn)用濾波器濾波后的信號(hào)方差小于濾波前，即提高測(cè)量精度；濾波器1（AR）的信號(hào)方差在值在小于0.2時(shí)信號(hào)方差和濾波器2（MA）得到的大致相同；但當(dāng)值慢慢增大時(shí)，大概到點(diǎn)0.53時(shí)，濾波器2的信號(hào)方差也突然隨之增大，而濾波器1（AR）的方差值卻始終在減小，因此可知濾波器1的性能要比濾波器2好。因此在選擇用濾波器2（MA）時(shí)應(yīng)該選擇值盡量小于其方差開始增大的點(diǎn)，即0.5左右。下面給出值為0.3和0.9的情況下，對(duì)白噪聲序列方差不同時(shí)對(duì)濾波效果進(jìn)行了比較，具體程序Matlab程序如下：clear;clc;N=1000;fori=1:Nvar0=i/N;init=2055615866;randn('seed',init);w=sqrt(var0)*randn(1,N);d=3800;%距離x=d+w;%輸入信號(hào)y1(1)=d;y2(1)=d;K=900;a=K/N;%確定a值forj=2:Ny1(j)=a*y1(j-1)+(1-a)*x(j);y2(j)=a*x(j)+(1-a)*x(j-1);endv1(i)=var(x);%輸入信號(hào)方差v2(i)=var(y1);%濾波器1后的方差v3(i)=var(y2);%濾波器2后的方差endt=(1:N)/N;plot(t,v1,'r',t,v2,'b',t,v3,'g');xlabel('白噪聲方差');ylabel('方差');legend('未濾波','濾波器1','濾波器2運(yùn)行結(jié)果見圖7-2。（a）=0.3（b）=0.9圖7-2不同值白噪聲方差對(duì)濾波結(jié)果的影響由圖可以看出，兩種濾波器下信號(hào)方差的變化趨勢(shì)都為線性增大的狀態(tài)，當(dāng)=0.9時(shí)，濾波器1（AR）性能優(yōu)于=0.3，且明顯好于=0.9時(shí)濾波器2（MA），在這種情況下選擇濾波器時(shí)優(yōu)先考慮濾波器1（AR）。八、卡爾曼濾波器題目：簡(jiǎn)述卡爾曼濾波器。并計(jì)算下面的例子：（15分）一連續(xù)平穩(wěn)的隨機(jī)信號(hào)，自相關(guān)函數(shù)，信號(hào)為加性噪聲所干擾，噪聲是白噪聲，測(cè)量值的離散值為已知，=0.01s,-3.2,-0.8,-14,-16,-17,-18,-3.3,-2.4,-18,-0.3,-0.4,-0.8,-19,-2.0,-1.2,-11,-14,-0.9,0.8,10,0.2,0.5,-0.5,2.4,-0.5,0.5,-13,0.5,10,-12,0.5,-0.6,-15,-0.7,15,0.5,-0.7,-2.0,-19,-17,-11,-14.利用卡爾曼濾波，估計(jì)信號(hào)的波形。并給出不同的量測(cè)噪聲對(duì)卡爾曼濾波的影響，同時(shí)給出迭代次數(shù)與均方誤差的關(guān)系曲線。1．卡爾曼濾波器簡(jiǎn)介卡爾曼濾波是一種高效率的遞歸濾波器(自回歸濾波器),它能夠從一系列的不完全包含噪聲的測(cè)量中，估計(jì)動(dòng)態(tài)系統(tǒng)的狀態(tài)。簡(jiǎn)單來(lái)說(shuō)，卡爾曼濾波器是一個(gè)“optimalrecursivedataprocessingalgorithm（最優(yōu)化自回歸數(shù)據(jù)處理算法）”。對(duì)于解決很大部分的問題，他是最優(yōu)，效率最高甚至是最有用的。他的廣泛應(yīng)用已經(jīng)超過(guò)30年，包括機(jī)器人導(dǎo)航，控制，傳感器數(shù)據(jù)融合甚至在軍事方面的雷達(dá)系統(tǒng)以及導(dǎo)彈追蹤等等。近年來(lái)更被應(yīng)用于計(jì)算機(jī)圖像處理，例如頭臉識(shí)別，圖像分割，圖像邊緣檢測(cè)等等。首先，我們先要引入一個(gè)離散控制過(guò)程的系統(tǒng)。該系統(tǒng)可用一個(gè)線性隨機(jī)微分方程（LinearStochasticDifferenceequation）來(lái)描述：（8-1）再加上系統(tǒng)的測(cè)量值：（8-2）上面兩式子中，是時(shí)刻的系統(tǒng)狀態(tài)，是時(shí)刻對(duì)系統(tǒng)的控制量。A和B是系統(tǒng)參數(shù)，對(duì)于多模型系統(tǒng)，他們?yōu)榫仃?。是時(shí)刻的測(cè)量值，H是測(cè)量系統(tǒng)的參數(shù)，對(duì)于多測(cè)量系統(tǒng)，為矩陣。和分別表示過(guò)程和測(cè)量的噪聲。他們被假設(shè)成高斯白噪聲(WhiteGaussianNoise)，他們的協(xié)方差矩陣分別是Q，R（這里我們假設(shè)他們不隨系統(tǒng)狀態(tài)變化而變化）。對(duì)于滿足上面的條件(線性隨機(jī)微分系統(tǒng)，過(guò)程和測(cè)量都是高斯白噪聲)，卡爾曼濾波器是最優(yōu)的信息處理器。下面我們來(lái)用他們結(jié)合他們的協(xié)方差矩陣來(lái)估算系統(tǒng)的最優(yōu)化輸出。首先我們要利用系統(tǒng)的過(guò)程模型，來(lái)預(yù)測(cè)下一狀態(tài)的系統(tǒng)。假設(shè)現(xiàn)在的系統(tǒng)狀態(tài)是，根據(jù)系統(tǒng)的模型，可以基于系統(tǒng)的上一狀態(tài)而預(yù)測(cè)出現(xiàn)在狀態(tài)：(8-3)式(8-3)中，是利用上一狀態(tài)預(yù)測(cè)的結(jié)果，是上一狀態(tài)最優(yōu)的結(jié)果，為現(xiàn)在狀態(tài)的控制量，如果沒有控制量，它可以為0。到現(xiàn)在為止，我們的系統(tǒng)結(jié)果已經(jīng)更新了，可是，對(duì)應(yīng)于的協(xié)方差矩陣還沒更新。我們用表示協(xié)方差矩陣：(8-4)式(5-4)中，是對(duì)應(yīng)的協(xié)方差矩陣，是對(duì)應(yīng)的協(xié)方差矩陣，表示的轉(zhuǎn)置矩陣，是系統(tǒng)過(guò)程的協(xié)方差矩陣。式子8-1，8-2就是卡爾曼濾波器5個(gè)公式當(dāng)中的前兩個(gè)，也就是對(duì)系統(tǒng)的預(yù)測(cè)?，F(xiàn)在我們有了現(xiàn)在狀態(tài)的預(yù)測(cè)結(jié)果，然后我們?cè)偈占F(xiàn)在狀態(tài)的測(cè)量值。結(jié)合預(yù)測(cè)值和測(cè)量值，我們可以得到現(xiàn)在狀態(tài)的最優(yōu)化估算值：(8-5)其中為卡爾曼增益(KalmanGain)：(8-6)到現(xiàn)在為止，我們已經(jīng)得到了狀態(tài)下最優(yōu)的估算值。但是為了要另卡爾曼濾波器不斷的運(yùn)行下去直到系統(tǒng)過(guò)程結(jié)束，我們還要更新狀態(tài)下的協(xié)方差矩陣：(8-7)其中為1的矩陣，對(duì)于單模型單測(cè)量，。2.利用卡爾曼濾波，估計(jì)信號(hào)的波形（1）求過(guò)程的功率譜密度將積分按＋和－分成兩部分進(jìn)行，信號(hào)的自相關(guān)函數(shù)是，則其功率譜為（2）用復(fù)頻率s來(lái)表示功率譜密度（）（3）譜分解與Z變換根據(jù)沖激響應(yīng)不變?cè)?，將連續(xù)系統(tǒng)變成離散系統(tǒng)，即其中，根據(jù)等式得變換到時(shí)域得：,因此，所以得出狀態(tài)方程，即B=1，又因?yàn)榱繙y(cè)方程，所以，D(k)=0，T=0.01s，。（4）下面給出的matlab程序，包含自編卡爾曼濾波算法和matlab自帶卡爾曼函數(shù)。在用matlab自帶卡爾曼函數(shù)前，先調(diào)用sys=ss（A,B,C,D,Ts）函數(shù)，這里采樣周期Ts=0.02s，A,B,C,D均已求出，得到離散卡爾曼狀態(tài)模型，然后調(diào)用函數(shù)[kalmf,L,P,H,E]=kalman(sys,Q,R),設(shè)計(jì)離散卡爾曼濾波器。H和E表示系統(tǒng)穩(wěn)態(tài)的修正矩陣和均方誤差，把它們代入估計(jì)值的更新公式（如下）就可以根據(jù)觀測(cè)信號(hào)得到卡爾曼濾波的估計(jì)值了。其中，S(k)是卡爾曼濾波估計(jì)值，Y(k)是已知的測(cè)量值。clearall;Y=[-3.2,-0.8,-14,-16,-17,-18,-3.3,-2.4,-18,-0.3,-0.4,-0.8,-19,-2.0,-1.2,-11,-14,-0.9,0.8,10,0.2,0.5,-0.5,2.4,-0.5,0.5,-13,0.5,10,-12,0.5,-0.6,-15,-0.7,15,0.5,-0.7,-2.0,-19,-17,-11,-14];N=length(Y);%卡爾曼函數(shù)A=0.99;B=1;C=1;D=0;T=0.01;sys=ss(A,B,C,D,T);%求出離散卡爾曼狀態(tài)模型v=std(randn(1,N));R=v^2;Q=0.02;[kalmf,L,P,H,E]=kalman(sys,Q,R);%H：修正矩陣E:最小均方誤差S=zeros(1,N);fork=2:NS(k)=A*S(k-1)+H*(Y(k)-C*A*S(k-1));%卡爾曼濾波后的估計(jì)值endk=1:1:42;plot(k*0.02,S,'b-',k*0.02,Y,'r');legend('估計(jì)值','測(cè)量值');ylabel('信號(hào)幅值');title('卡爾曼函數(shù)濾波結(jié)果');gridon;運(yùn)行結(jié)果如圖8-1所示：圖8-1卡爾曼濾波估計(jì)圖上圖是直接用卡爾曼濾波函數(shù)進(jìn)行濾波的效果。在Matlab當(dāng)中，可以根據(jù)狀態(tài)方程和觀測(cè)方程，調(diào)用dlqe函數(shù)，直接求解其最佳遞推估計(jì)方程，可以利用該方程對(duì)輸入信號(hào)進(jìn)行濾波，比較兩種方式的結(jié)果，進(jìn)行結(jié)果驗(yàn)證。結(jié)果驗(yàn)證的Matlab程序如下：%采用dlqe函數(shù)求解kalman濾波器的穩(wěn)態(tài)參數(shù)a=0.99;G=1;Qk=0.02;Rk=1;c=1;[L,P1,P]=dlqe(a,G,c,Qk,Rk);%%%%%%%%%根據(jù)kalman濾波器的穩(wěn)態(tài)參數(shù)設(shè)計(jì)的濾波器fori=1:size_x(1)xw(i)=X(i);endfori=2:size_x(1)xw(i)=0.99*xw(i-1)+L*(X(i)-0.99*xw(i-1));end%結(jié)果顯示subplot(312);plot(k*0.02,xk_s(k),'b-',k*0.02,xw(k),'r');legend('Kalman','Wiener');xlabel('時(shí)間Ts');ylabel('信號(hào)幅值');title('Kalman和穩(wěn)態(tài)濾波比較');gridon;本題中，信號(hào)為加性噪聲所干擾，噪聲是白噪聲，但是白噪聲的方差沒有給出，以為例，Matlab程序得到的結(jié)果如圖8-2所示。圖8-2卡爾曼濾波結(jié)果和穩(wěn)態(tài)濾波結(jié)果比較穩(wěn)態(tài)濾波結(jié)果是采用dlqe函數(shù)編寫的濾波函數(shù)和kalman，從圖可以看出，卡爾曼濾波在0.3s以前和穩(wěn)態(tài)濾