數(shù)字信號處理實驗

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

2、三、實驗要求信號是正弦波加正態(tài)零均值白噪聲，信噪比為10dB，信號頻率為2kHZ，取樣頻率為100kHZ。四、實驗程序與實驗結(jié)果（1）用周期圖法進(jìn)行譜估計A、實驗程序：%用周期法進(jìn)行譜估計clear all;N1=128;%數(shù)據(jù)長度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計算正弦信號的幅度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、實驗仿真結(jié)果：（2）采用漢明窗，分段長度L=32，用修正的周期圖求平均法進(jìn)行譜估計A：實驗程序:clear all;N=512;%數(shù)據(jù)長度

6、Ns=32;%分段長度f1=2;%正弦波頻率，單位為kHZfs=100;%抽樣頻率,單位為kHZn=0:N-1;a=sqrt(20);%由信噪比為10dB計算正弦信號的幅度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：實驗仿真結(jié)果：五參考文獻(xiàn)：1 丁玉美，闊永紅，高新波.數(shù)字信號處理-時域離散隨機(jī)信號處理M.西安：西安電子科技大學(xué)出版社，2002.2 萬建偉，王玲.信號處理仿真技術(shù)M.長沙：國防科技大學(xué)出版社，2008.實驗

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

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

12、信號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è)信號初值為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為測量出來的離散值；N=length(zk); %要測量的點(diǎn)數(shù)；for k=2:N p1(k)=Ak*p(k-1)*Ak'+Qk; %未考慮噪聲時的均方誤差陣；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 %顯示信號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)的估計波形與z(t)波形','fontsize',16) legend('觀測數(shù)據(jù)z(t)','信號估計值x(t)',2); grid;四實驗結(jié)果五實驗小結(jié)通過卡爾曼濾波估