信号检测与处理实验报告

Harbin Institute of Technology 信号检测与处理 实验报告 2016年01月问题：最小二乘估计一次完成算法1问题描述考虑仿真对象 其中，是服从正态分布的白噪声N。输入信号采用4阶M序列（伪随机序列模拟白噪声），幅度为1。试对模型参数进行估计。2问题分析设输入信号的取值是从k =1到k =16的M序列，由最小二乘估计原理可知，待估计参数为=。其中，被估计参数、观测矩阵z L、H L的表达式为 , ， 通过matlab对系统进行仿真，仿真算法程序流程图如图1所示。赋输入信号初值u定义输出观测值的长度并计算系统的输出值画出输入和输出观测值的图形给样本矩阵HL和zL赋值计算参数从中分离出并显示出被辨识参数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); %定义输出观测值的长度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; %横坐标范围是1到16，步长为1plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线subplot(3,1,3) %画三行一列图形窗口中的第三个图形stem(z),grid on %画出输出观测值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) % 给样本矩阵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实验运行结果如下： 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 最小二乘一次完成算法仿真实例中输入信号和输出观测值可以看到，由于仿真模型中未引入噪声，估计结果与实际系统一致。3.对该问题局限性的讨论与改进在对上面问题的讨论中，我们附加了两个条件，一个是虽然不知道其具体参数，但我们知道系统的阶次。另一个假设是系统的输出中没有附加噪声。而实际情况往往不满足这两个已知假设条件。下面讨论在有附加噪声和系统阶次未知的情况下，如何对系统的参数进行估计。最小二乘法估计本身就可以很好的滤除噪声对系统参数的影响。这里我们对系统附加了一个随机误差干扰。问题主要集中在如何对系统的阶次进行估计。常用的系统阶次结构辨识方法有以下几种:根据汉格尔矩阵估计系统阶次设一个可观可控的SISO过程的脉冲响应序列为个g(1),g(2),g(L),可以通过汉格尔（Hankel）矩阵的秩来确定系统的阶次。令Hankel阵为：，其中l决定H(l,k)阵的维数，k可在1至(L-2l+2)间任意选择。则有rankH(l,k)=n0，l=n0。如果l=n0（过程的真实阶次），那么Hankel阵的秩等于n0。因此可以利用Hankel阵的奇异性来确定系统的阶次n0。根据残差平方和估计模型的阶次SISO过程的差分方程模型的输出残差为z(k)，数据长度L，Hn为n阶时的数据矩阵，为n阶时的参数的估计量，n为模型阶次估计值，n0为真实阶次，则残差平方和函数J：残差平方和有这样的性质：当L足够大时，随着n增加J(n)先是显著地下降，当nn0时，J(n)值显著下降的现象就终止。这就是损失函数法来定阶的原理。 根据上述原理，编写程序如下：L=300; % M序列的周期x1=1;x2=1;x3=1;x4=0;x5=1;x6=0; %四个移位积存器的输出初始值for k=1:L; %开始循环，长度为L u(k)=xor(x3,x4); %第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或” x6=x5;x5=x4;x4=x3;x3=x2;x2=x1;x1=u(k);end %大循环结束，产生输入信号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);%用循环产生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阶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; %假设为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; %假设为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; %假设为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; %假设为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)