數(shù)字信號(hào)處理實(shí)驗(yàn)

1、實(shí)驗(yàn)一 自適應(yīng)濾波器一、實(shí)驗(yàn)?zāi)康?、掌握功率譜估計(jì)方法2、會(huì)用matlab對(duì)功率譜進(jìn)行仿真二、實(shí)驗(yàn)原理 功率譜估計(jì)方法有很多種，一般分成兩大類，一類是經(jīng)典譜估計(jì)；另一類是現(xiàn)代譜估計(jì)。經(jīng)典譜估計(jì)可以分成兩種，一種是BT法，另一種是周期法；BT法是先估計(jì)自相關(guān)函數(shù)，然后將相關(guān)函數(shù)進(jìn)行傅里葉變換得到功率譜函數(shù)。相應(yīng)公式如下所示：周期圖法是采用功率譜的另一種定義，但與BT法是等價(jià)的，相應(yīng)的功率譜估計(jì)如下所示： 其計(jì)算框圖如下所示：圖1.1周期圖法計(jì)算用功率譜框圖由于觀測(cè)數(shù)據(jù)有限，所以周期圖法估計(jì)分辨率低，估計(jì)誤差大。針對(duì)經(jīng)典譜估計(jì)的缺點(diǎn)，一般有三種改進(jìn)方法：平均周期圖法、窗函數(shù)法和修正的周期圖平均法。

2、三、實(shí)驗(yàn)要求信號(hào)是正弦波加正態(tài)零均值白噪聲，信噪比為10dB，信號(hào)頻率為2kHZ，取樣頻率為100kHZ。四、實(shí)驗(yàn)程序與實(shí)驗(yàn)結(jié)果（1）用周期圖法進(jìn)行譜估計(jì)A、實(shí)驗(yàn)程序：%用周期法進(jìn)行譜估計(jì)clear all;N1=128;%數(shù)據(jù)長(zhǎng)度N2=256;N3=512;N4=1024;f=2;%正弦波頻率，單位為kHZfs=100;%抽樣頻率,單位為kHZn1=0:N1-1;n2=0:N2-1;n3=0:N3-1;n4=0:N4-1;a=sqrt(20);%由信噪比為10dB計(jì)算正弦信號(hào)的幅度wn1=randn(1,N1);xn1=a*sin(2*pi*f*n1./fs)+wn1;Pxx1=10*log

3、10(abs(fft(xn1).2)/N1);%周期法求功率譜f1=(0:length(Pxx1)-1)/length(Pxx1);wn2=randn(1,N2);xn2=a*sin(2*pi*f*n2./fs)+wn2;Pxx2=10*log10(abs(fft(xn2).2)/N2);f2=(0:length(Pxx2)-1)/length(Pxx2);wn3=randn(1,N3);xn3=a*sin(2*pi*f*n3./fs)+wn3;Pxx3=10*log10(abs(fft(xn3).2)/N3);f3=(0:length(Pxx3)-1)/length(Pxx3);wn4=ra

4、ndn(1,N4);xn4=a*sin(2*pi*f*n4./fs)+wn4;Pxx4=10*log10(abs(fft(xn4).2)/N4);f4=(0:length(Pxx4)-1)/length(Pxx4);subplot(2,2,1);plot(f1,Pxx1);xlabel('頻率');ylabel('功率(dB)');title('功率譜Pxx，N=128');subplot(2,2,2);plot(f2,Pxx2);xlabel('頻率');ylabel('功率(dB)');title('功

5、率譜Pxx，N=256');subplot(2,2,3);plot(f3,Pxx3);xlabel('頻率');ylabel('功率(dB)');title('功率譜Pxx，N=512');subplot(2,2,4);plot(f4,Pxx4);xlabel('頻率');ylabel('功率(dB)');title('功率譜Pxx，N=1024');B、實(shí)驗(yàn)仿真結(jié)果：（2）采用漢明窗，分段長(zhǎng)度L=32，用修正的周期圖求平均法進(jìn)行譜估計(jì)A：實(shí)驗(yàn)程序:clear all;N=512;%數(shù)據(jù)長(zhǎng)度

6、Ns=32;%分段長(zhǎng)度f(wàn)1=2;%正弦波頻率，單位為kHZfs=100;%抽樣頻率,單位為kHZn=0:N-1;a=sqrt(20);%由信噪比為10dB計(jì)算正弦信號(hào)的幅度wn=randn(1,N);xn=a*sin(2*pi*f1*n./fs)+wn;w=hamming(32)'%漢明窗Pxx1=abs(fft(w.*xn(1:32),Ns).2)/norm(w)2;Pxx2=abs(fft(w.*xn(33:64),Ns).2)/norm(w)2;Pxx3=abs(fft(w.*xn(65:96),Ns).2)/norm(w)2;Pxx4=abs(fft(w.*xn(97:128)

7、,Ns).2)/norm(w)2;Pxx5=abs(fft(w.*xn(129:160),Ns).2)/norm(w)2;Pxx6=abs(fft(w.*xn(161:192),Ns).2)/norm(w)2;Pxx7=abs(fft(w.*xn(193:224),Ns).2)/norm(w)2;Pxx8=abs(fft(w.*xn(225:256),Ns).2)/norm(w)2;Pxx9=abs(fft(w.*xn(257:288),Ns).2)/norm(w)2;Pxx10=abs(fft(w.*xn(289:320),Ns).2)/norm(w)2;Pxx11=abs(fft(w.*x

8、n(321:352),Ns).2)/norm(w)2;Pxx12=abs(fft(w.*xn(353:384),Ns).2)/norm(w)2;Pxx13=abs(fft(w.*xn(385:416),Ns).2)/norm(w)2;Pxx14=abs(fft(w.*xn(417:448),Ns).2)/norm(w)2;Pxx15=abs(fft(w.*xn(449:480),Ns).2)/norm(w)2;Pxx16=abs(fft(w.*xn(481:512),Ns).2)/norm(w)2;Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+

9、Pxx8+Pxx9+Pxx10+Pxx11+Pxx12+Pxx13+Pxx14+Pxx15+Pxx16)/16);f=(0:length(Pxx)-1)/length(Pxx);plot(f,Pxx);xlabel('頻率');ylabel('功率(dB)');title('加窗平均周期圖法功率譜Pxx，N=512');grid on;B：實(shí)驗(yàn)仿真結(jié)果：五參考文獻(xiàn)：1 丁玉美，闊永紅，高新波.數(shù)字信號(hào)處理-時(shí)域離散隨機(jī)信號(hào)處理M.西安：西安電子科技大學(xué)出版社，2002.2 萬(wàn)建偉，王玲.信號(hào)處理仿真技術(shù)M.長(zhǎng)沙：國(guó)防科技大學(xué)出版社，2008.實(shí)驗(yàn)

10、二 卡爾曼濾波器的設(shè)計(jì)一實(shí)驗(yàn)?zāi)康?.熟悉并掌握卡爾曼濾波、自適應(yīng)濾波和譜估計(jì)的原理。2.可以仿真符合要求的卡爾曼濾波器、自適應(yīng)濾波器和各種譜估計(jì)方法。3.掌握卡爾曼濾波器的遞推公式和仿真方法。4.熟悉matlab的用法。二實(shí)驗(yàn)原理卡爾曼濾波是用狀態(tài)空間法描述系統(tǒng)的，由狀態(tài)方程和測(cè)量方程所組成?？柭鼮V波用前一個(gè)狀態(tài)的估計(jì)值和最近一個(gè)觀測(cè)數(shù)據(jù)來(lái)估計(jì)狀態(tài)變量的當(dāng)前值，并以狀態(tài)變量的估計(jì)值的形式給出。其狀態(tài)方程和量測(cè)方程如下所示：其中，表示時(shí)間，輸入信號(hào)是一白噪聲，輸出信號(hào)的觀測(cè)噪聲也是一個(gè)白噪聲，輸入信號(hào)到狀態(tài)變量的支路增益等于1，即B=1；表示狀態(tài)變量之間的增益矩陣，可隨時(shí)間變化，表示第次迭代的

11、取值，表示狀態(tài)變量與輸出信號(hào)之間的增益矩陣，可隨時(shí)間變化，其信號(hào)模型如圖1.1所示(用代替)。圖1.1 卡爾曼濾波器的信號(hào)模型卡爾曼濾波是采用遞推的算法實(shí)現(xiàn)的，其基本思想是先不考慮輸入信號(hào)和觀測(cè)噪聲的影響，得到狀態(tài)變量和輸出信號(hào)的估計(jì)值，再用輸出信號(hào)的估計(jì)誤差加權(quán)矯正狀態(tài)變量的估計(jì)值，使?fàn)顟B(tài)變量估計(jì)誤差的均方值最小。其遞推公式如下所示： 假設(shè)初始條件已知，其中，那么遞推流程見圖1.2所示。 圖1.2 卡爾曼濾波遞推流程圖三實(shí)驗(yàn)要求一連續(xù)平穩(wěn)的隨機(jī)信號(hào)x(t),自相關(guān)，信號(hào)x(t)為加性噪聲所干擾，噪聲是白噪聲，測(cè)量值的離散值y(k)為已知。Matlab仿真程序如下：%編卡爾曼濾波遞推程序，估計(jì)

12、信號(hào)x(t)的波形clear all;clc;Ak=exp(-0.02); %各系數(shù)由前面確定；Ck=1; Rk=0.1; p(1)=20; %各初值；Qk=1-exp(-0.04); p1(1)=Ak*p(1)*Ak'+Qk; %由p1代表p'；x(1)=0; %設(shè)信號(hào)初值為0；H(1)=p1(1)*Ck'*inv(Ck*p1(1)*Ck'+Rk); zk=-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

13、.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 %zk為測(cè)量出來(lái)的離散值；N=length(zk); %要測(cè)量的點(diǎn)數(shù)；for k=2:N p1(k)=Ak*p(k-1)*Ak'+Qk; %未考慮噪聲時(shí)的均方誤差陣；H(k)=p1(k)*Ck'*inv(Ck*p1(k)*Ck'+Rk); %增益方程；I=eye(size(H(k); %產(chǎn)生和H（k）維數(shù)相同的單位矩陣；p(k)=(I-H(k)*Ck)*p1(k); %濾波的均方誤差陣；x(k)=Ak*x(k-1)+

14、H(k)*(zk(k)-Ck*Ak*x(k-1); %遞推公式；end,x %顯示信號(hào)x(k)的數(shù)據(jù)；m=1:N; n=m*0.02; plot(n,zk,'-r*',n,x,'-bo'); %便于比較zk和x(k)在同一窗口輸出；xlabel('t/s','Fontsize',16); ylabel('z(t),x(t)','fontsize',16); title('卡爾曼濾波遞推x(t)的估計(jì)波形與z(t)波形','fontsize',16) legend('觀測(cè)數(shù)據(jù)z(t)','信號(hào)估計(jì)值x(t)',2); grid;四實(shí)驗(yàn)結(jié)果五實(shí)驗(yàn)小結(jié)通過(guò)卡爾曼濾波估