信號檢測與處理實驗報告

Harbin Institute of Technology 信號檢測與處理 實驗報告 2016年01月問題：最小二乘估計一次完成算法1問題描述考慮仿真對象 其中，是服從正態(tài)分布的白噪聲N。輸入信號采用4階M序列（偽隨機序列模擬白噪聲），幅度為1。試對模型參數(shù)進行估計。2問題分析設(shè)輸入信號的取值是從k =1到k =16的M序列，由最小二乘估計原理可知，待估計參數(shù)為=。其中，被估計參數(shù)、觀測矩陣z L、H L的表達式為 , ， 通過matlab對系統(tǒng)進行仿真，仿真算法程序流程圖如圖1所示。賦輸入信號初值u定義輸出觀測值的長度并計算系統(tǒng)的輸出值畫出輸入和輸出觀測值的圖形給樣本矩陣HL和zL賦值計算參數(shù)從中分離出并顯示出被辨識參數(shù)a1, a2, b1, b2停機圖1 最小二乘一次完成算法程序框圖程序代碼如下：%二階系統(tǒng)的最小二乘一次完成算法估計程序u=-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1; %系統(tǒng)估計的輸入信號為一個周期的M序列z=zeros(1,16); %定義輸出觀測值的長度for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想輸出值作為觀測值endsubplot(3,1,1) %畫三行一列圖形窗口中的第一個圖形stem(u) %畫輸入信號u的徑線圖形subplot(3,1,2) %畫三行一列圖形窗口中的第二個圖形i=1:1:16; %橫坐標(biāo)范圍是1到16，步長為1plot(i,z) %圖形的橫坐標(biāo)是采樣時刻i, 縱坐標(biāo)是輸出觀測值z, 圖形格式為連續(xù)曲線subplot(3,1,3) %畫三行一列圖形窗口中的第三個圖形stem(z),grid on %畫出輸出觀測值z的徑線圖形，并顯示坐標(biāo)網(wǎng)格u,z %顯示輸入信號和輸出觀測信號%L=14 %數(shù)據(jù)長度HL=-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14) %給樣本矩陣HL賦值ZL=z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16) % 給樣本矩陣z L賦值%Calculating Parameters c1=HL*HL; c2=inv(c1); c3=HL*ZL; c=c2*c3 %計算并顯示 %Display Parameters a1=c(1), a2=c(2), b1=c(3),b2=c(4) %從 中分離出并顯示a1 、a2、 b1、 b2%End實驗運行結(jié)果如下： u = -1，1，-1，1，1，1，1，-1，-1，-1，1，-1，-1，1，1z = 0，0，0.5000，0.2500，0.5250，2.1125， 4.3012，6.4731，6.1988，3.2670，-0.9386， -3.1949，-4.6352，6.2165，-5.5800，-2.5185HL = 0 0 1.0000 -1.0000 -0.5000 0 -1.0000 1.0000 -0.2500 -0.5000 1.0000 -1.0000 -0.5250 -0.2500 1.0000 1.0000 -2.1125 -0.5250 1.0000 1.0000 -4.3012 -2.1125 1.0000 1.0000 -6.4731 -4.3012 -1.0000 1.0000 -6.1988 -6.4731 -1.0000 -1.0000 -3.2670 -6.1988 -1.0000 -1.0000 0.9386 -3.2670 1.0000 -1.0000 3.1949 0.9386 -1.0000 1.0000 4.6352 3.1949 -1.0000 -1.0000 6.2165 4.6352 1.0000 -1.0000 5.5800 6.2165 1.0000 1.0000ZL = 0.5000，0.2500，0.5250，2.1125，4.3012，6.4731，6.1988，3.2670，-0.9386，-3.1949， -4.6352，-6.2165，-5.5800，-2.5185Tc = -1.5000，0.7000，1.0000，0.5000Ta1 = -1.5000a2 = 0.7000b1 = 1.0000b2 =0.5000輸入信號與輸出觀測值波形如圖2所示。圖2 最小二乘一次完成算法仿真實例中輸入信號和輸出觀測值可以看到，由于仿真模型中未引入噪聲，估計結(jié)果與實際系統(tǒng)一致。3.對該問題局限性的討論與改進在對上面問題的討論中，我們附加了兩個條件，一個是雖然不知道其具體參數(shù)，但我們知道系統(tǒng)的階次。另一個假設(shè)是系統(tǒng)的輸出中沒有附加噪聲。而實際情況往往不滿足這兩個已知假設(shè)條件。下面討論在有附加噪聲和系統(tǒng)階次未知的情況下，如何對系統(tǒng)的參數(shù)進行估計。最小二乘法估計本身就可以很好的濾除噪聲對系統(tǒng)參數(shù)的影響。這里我們對系統(tǒng)附加了一個隨機誤差干擾。問題主要集中在如何對系統(tǒng)的階次進行估計。常用的系統(tǒng)階次結(jié)構(gòu)辨識方法有以下幾種:根據(jù)漢格爾矩陣估計系統(tǒng)階次設(shè)一個可觀可控的SISO過程的脈沖響應(yīng)序列為個g(1),g(2),g(L),可以通過漢格爾（Hankel）矩陣的秩來確定系統(tǒng)的階次。令Hankel陣為：，其中l(wèi)決定H(l,k)陣的維數(shù)，k可在1至(L-2l+2)間任意選擇。則有rankH(l,k)=n0，l=n0。如果l=n0（過程的真實階次），那么Hankel陣的秩等于n0。因此可以利用Hankel陣的奇異性來確定系統(tǒng)的階次n0。根據(jù)殘差平方和估計模型的階次SISO過程的差分方程模型的輸出殘差為z(k)，數(shù)據(jù)長度L，Hn為n階時的數(shù)據(jù)矩陣，為n階時的參數(shù)的估計量，n為模型階次估計值，n0為真實階次，則殘差平方和函數(shù)J：殘差平方和有這樣的性質(zhì)：當(dāng)L足夠大時，隨著n增加J(n)先是顯著地下降，當(dāng)nn0時，J(n)值顯著下降的現(xiàn)象就終止。這就是損失函數(shù)法來定階的原理。 根據(jù)上述原理，編寫程序如下：L=300; % M序列的周期x1=1;x2=1;x3=1;x4=0;x5=1;x6=0; %四個移位積存器的輸出初始值for k=1:L; %開始循環(huán)，長度為L u(k)=xor(x3,x4); %第一個移位積存器的輸入是第3個與第4個移位積存器的輸出的“或” x6=x5;x5=x4;x4=x3;x3=x2;x2=x1;x1=u(k);end %大循環(huán)結(jié)束，產(chǎn)生輸入信號u plot(u) %繪圖M序列v=randn(300,1); %隨機誤差干擾z=zeros(1,300);for k=4:300z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)/400; endH=zeros(300,6); %定義一個H“0”矩陣for i=4:300 H(i,:)=-z(i-1) -z(i-2) -z(i-3) u(i-1) u(i-2) u(i-3);%用循環(huán)產(chǎn)生H矩陣 z1(i,:)=z(i); %用循環(huán)產(chǎn)生z矩陣end%計算參數(shù)%c=inv(H*H)*H*z1%帶入公式書上3.1.23a1=c(1),a2=c(2),a3=c(3),b1=c(4),b2=c(5),b3=c(6)%辨識出參數(shù)%系統(tǒng)階次辨識AIC算法%bb=zeros(5,1);n=1; %假設(shè)為1階for i=2:300 H1(i,:)=-z(i-1) u(i-1); zz1(i,:)=z(i);endaa1=inv(H1*H1)*H1*zz1bb(1)=(zz1-H1*aa1)*(zz1-H1*aa1)/L;AIC(1)=L*log(bb(1)+4*n;n=2; %假設(shè)為2階for i=3:300 H2(i,:)=-z(i-1) -z(i-2) u(i-1) u(i-2); zz2(i,:)=z(i);endaa2=inv(H2*H2)*H2*zz2bb(2)=(zz2-H2*aa2)*(zz2-H2*aa2)/L;AIC(2)=L*log(bb(2)+4*n;n=3; %假設(shè)為3階for i=4:300 H3(i,:)=-z(i-1) -z(i-2) -z(i-3) u(i-1) u(i-2) u(i-3); zz3(i,:)=z(i);endaa3=inv(H3*H3)*H3*zz3bb(3)=(zz3-H3*aa3)*(zz3-H3*aa3)/L;AIC(3)=L*log(bb(3)+4*n;n=4; %假設(shè)為4階for i=5:300 H4(i,:)=-z(i-1) -z(i-2) -z(i-3) -z(i-4) u(i-1) u(i-2) u(i-3) u(i-4); zz4(i,:)=z(i);endaa4=inv(H4*H4)*H4*zz4bb(4)=(zz4-H4*aa4)*(zz4-H4*aa4)/L;AIC(4)=L*log(bb(4)+4*n;n=5; %假設(shè)為5階for i=6:300 H5(i,:)=-z(i-1) -z(i-2) -z(i-3) -z(i-4) -z(i-5) u(i-1) u(i-2) u(i-3) u(i-4) u(i-5); zz5(i,:)=z(i);endaa5=inv(H5*H5)*H5*zz5bb(5)=(zz5-H5*aa5)*(zz5-H5*aa5)/L;AIC(5)=L*log(bb(5)+4*n;x=min(AIC)