信号检测与处理实验报告

HarbinInstituteofTechnology信号检测与处理实验报告2016年01月问题：最小二乘估计一次完成算法1．问题描述考虑仿真对象)()2(5.0)1()2(7.0)1(5.1)(kvkukukzkzkz其中，)(kv是服从正态分布的白噪声N)1,0(。输入信号采用4阶M序列（伪随机序列模拟白噪声），幅度为1。试对模型参数进行估计。2．问题分析设输入信号的取值是从k=1到k=16的M序列，由最小二乘估计原理可知，待估计参数LSθˆ为LSθˆ=LτL1LτLzH)HH(。其中，被估计参数LSθˆ、观测矩阵zL、HL的表达式为2121ˆbbaaLSθ,)16()4()3(zzzLz，)14()2()1()15()3()2()14()2()1()15()3()2(uuuuuuzzzzzzLH通过matlab对系统进行仿真，仿真算法程序流程图如图1所示。赋输入信号初值u定义输出观测值的长度并计算系统的输出值画出输入和输出观测值的图形给样本矩阵HL和zL赋值计算参数LSθˆ从LSθˆ中分离出并显示出被辨识参数a1,a2,b1,b2停机图1最小二乘一次完成算法程序框图程序代码如下：%二阶系统的最小二乘一次完成算法估计程序u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1];%系统估计的输入信号为一个周期的M序列z=zeros(1,16);%定义输出观测值的长度fork=3:16z(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;%横坐标范围是1到16，步长为1plot(i,z)%图形的横坐标是采样时刻i,纵坐标是输出观测值z,图形格式为连续曲线subplot(3,1,3)%画三行一列图形窗口中的第三个图形stem(z),gridon%画出输出观测值z的径线图形，并显示坐标网格u,z%显示输入信号和输出观测信号%L=14%数据长度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)]%给样本矩阵zL赋值%CalculatingParametersc1=HL'*HL;c2=inv(c1);c3=HL'*ZL;c=c2*c3%计算并显示%DisplayParametersa1=c(1),a2=c(2),b1=c(3),b2=c(4)%从中分离出并显示a1、a2、b1、b2%End实验运行结果如下：u=[-1，1，-1，1，1，1，1，-1，-1，-1，1，-1，-1，1，1]z=[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.5185]HL=001.0000-1.0000-0.50000-1.00001.0000-0.2500-0.50001.0000-1.0000-0.5250-0.25001.00001.0000-2.1125-0.52501.00001.0000-4.3012-2.11251.00001.0000-6.4731-4.3012-1.00001.0000-6.1988-6.4731-1.0000-1.0000-3.2670-6.1988-1.0000-1.00000.9386-3.26701.0000-1.00003.19490.9386-1.00001.00004.63523.1949-1.0000-1.00006.21654.63521.0000-1.00005.58006.21651.00001.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.5185]Tc=[-1.5000，0.7000，1.0000，0.5000]Ta1=-1.5000a2=0.7000b1=1.0000b2=0.5000输入信号与输出观测值波形如图2所示。可以看到，由于仿真模型中未引入噪声，估计结果与实际系统一致。051015-1010246810121416-100100246810121416-10010STEMPLOTOFINPUTMSERIESWAVEOFMEASUREMENTSOFOUTPUTSTEMPLOTOFMEASUREMENTSOFOUTPUT图2最小二乘一次完成算法仿真实例中输入信号和输出观测值3.对该问题局限性的讨论与改进在对上面问题的讨论中，我们附加了两个条件，一个是虽然不知道其具体参数，但我们知道系统的阶次。另一个假设是系统的输出中没有附加噪声。而实际情况往往不满足这两个已知假设条件。下面讨论在有附加噪声和系统阶次未知的情况下，如何对系统的参数进行估计。最小二乘法估计本身就可以很好的滤除噪声对系统参数的影响。这里我们对系统附加了一个随机误差干扰。问题主要集中在如何对系统的阶次进行估计。常用的系统阶次结构辨识方法有以下几种:根据汉格尔矩阵估计系统阶次设一个可观可控的SISO过程的脉冲响应序列为{个g(1),g(2),……g(L)},可以通过汉格尔（Hankel）矩阵的秩来确定系统的阶次。令Hankel阵为：)22()1()1()()2()1()1()1()(),(lkgkglkglkgkgkglkgkgkgklH，其中l决定H(l,k)阵的维数，k可在1至(L-2l+2)间任意选择。则有rank[H(l,k)]=n0，l=n0。如果l=n0（过程的真实阶次），那么Hankel阵的秩等于n0。因此可以利用Hankel阵的奇异性来确定系统的阶次n0。根据残差平方和估计模型的阶次SISO过程的差分方程模型的输出残差为z(k)，数据长度L，Hn为n阶时的数据矩阵，nˆˆ为n阶时的参数的估计量，n为模型阶次估计值，n0为真实阶次，则残差平方和函数J：)(~1)ˆ()ˆ(1~~1)ˆ(12ˆˆˆˆˆˆ00kzLHzHzLzzLnJLnnknnnTnnnnTn残差平方和有这样的性质：当L足够大时，随着n增加J(n))ˆ(nJ先是显著地下降，当nn0时，J(n)值显著下降的现象就终止。这就是损失函数法来定阶的原理。根据上述原理，编写程序如下：L=300;%M序列的周期x1=1;x2=1;x3=1;x4=0;x5=1;x6=0;%四个移位积存器的输出初始值fork=1:L;%开始循环，长度为Lu(k)=xor(x3,x4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或”x6=x5;x5=x4;x4=x3;x3=x2;x2=x1;x1=u(k);end%大循环结束，产生输入信号uplot(u)%绘图M序列v=randn(300,1);%随机误差干扰z=zeros(1,300);fork=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”矩阵fori=4:300H(i,:)=[-z(i-1)-z(i-2)-z(i-3)u(i-1)u(i-2)u(i-3)];%用循环产生H矩阵z1(i,:)=[z(i)];%用循环产生z矩阵end%计算参数%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)%辨识出参数%系统阶次辨识AIC算法%bb=zeros(5,1);n=1;%假设为1阶fori=2:300H1(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;%假设为2阶fori=3:300H2(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;%假设为3阶fori=4:300H3(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;%假设为4阶fori=5:300H4(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;%假设为5阶fori=6:300H5(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)fori=1:5if(AIC(i)==x)N=i%所辨识出的阶次Nendendplot(1:5,[AIC(1)AIC(2)AIC(3)AIC(4)AIC(5)])程序运行结果如下a1=-1.4961a2=0.6941a3=0.0027b1=0.9995b2=0.5041b3=0.0017其中a3，b3接近于0，辨识出的系统参数与实际系统参数之间的误差不超过2%，可见该方法能够较好的完成对系统阶次和参数的辨识。