M序列的产生

M 序列的产生
1. 对象或参数
数学模型如下：
)()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z +-+-=-+--
其中，)(k v 是白噪声N （0,1）（服从正态分布），输入信号采用4阶M 序列，幅度为1. 辨识模型如下：
)()2()1()2()1()(2121k v k u b k b k z a k z a k z +-+-=-+--
设输入信号的取值为k=1到k=16的M 序列，则待辨识参数L L L L LS z H H H ττθ1-）（=
∧。

⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=∧4321a a a a LS θ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)16(...)4()3(z z z z L ⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡---=)14()15()14()15(......)2()3(z(2)-)3(-)1()2()1()2(-u u z z u u z u u z z H L 编制仿真程序，，获取输入输出数据，运用最小二乘法对系统的参数进行辨识，并将辨识结果与实际参数进行对比。

2. 程序框图
3. 程序实现
function [] = JM-232()
%--------------------------------------------------------实验题目及初始化定义
disp('最小二乘法的实现')
disp(' ')
disp('数学模型为：z(k)-1.5z(k-1)+0.7z(k-2)=u(k-1)+0.5u(k-2)+v(k).')
disp(' 所选择的辨识模型为：z(k)+a1*z(k-1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k).') disp(' 故实际参数a1=-1.5, a2=0.7, b1=1, b2=0.5. ')
disp(' ')
disp('系统输入采用4 阶M 序列且其幅度为1，输出受到白噪声序列v(k)干扰.') ds = input(' 假设白噪声序列v(k)服从正态分布，均值为0，请输入方差ds = '); disp(' ')
disp(' ')
a0=65539;M=2147483647;x=123456;b=10000;
U=[];V=[];Z=[];T=[];
P=zeros(20,4);H=zeros(14,4);LS=zeros(4,1);ZL=zeros(14,1);
%------------------------------------------------产生4 阶M 序列作为输入信号u(k)
fori=1:4
P(1,i)=1;
U(1)=-2*(P(1,4)-0.5);
end
for j=2:16
fori=2:4
P(j,i)=P(j-1,i-1);
end
P(j,1)=mod(P(j-1,3)+P(j-1,4),2);
U(j)=-2*(P(j,4)-0.5);
end
disp('** 通过仿真获得：')
disp(' ')
disp(' 系统输入u（k）取k=1 到k=16 图形如figure 1 所示；')
disp(' ')
figure(1)
plot(U);hold on
axis([0 17 -1.5 1.5])
title('系统输入u（k）图形')
%--------------------------------------------------产生白噪声作为干扰信号v(k)
fori=1:192
x=mod(a0*x+b,M);
T(i)=x/M;
end
aT=mean(T);vT=var(T);
fori=1:16
tempsum=0;
for j=1:12
tempsum=tempsum+T(12*i-j+1);
end
V(i)=sqrt(ds)*(tempsum-12*aT)/sqrt(12*vT);
end
aV=mean(V);vV=var(V);
disp(' 系统白噪声干扰v（k）取k=1 到k=16 图形如figure 2 所示；') disp([' 该系统白噪声v（k）均值为' num2str(aV)])
disp([' 方差为' num2str(vV)])
disp(' ')
figure(2)
plot(V);hold on
axis([0 17 -3 3])
title('系统白噪声v（k）图形')
%--------------------------------计算输出z(k)，同时由系统输出构成观测矩阵ZL 与H
Z(1)=0;Z(2)=0;
for k=3:16
Z(k)=1.5*Z(k-1)-0.7*Z(k-2)+U(k-1)+0.5*U(k-2)+V(k);
ZL(k-2,1)=Z(k);
end
fori=1:14
H(i,1)=-Z(i+1);H(i,2)=-Z(i);H(i,3)=U(i+1);H(i,4)=U(i);
end
disp(' 系统输出z（k）图形如figure 3 所示；')
disp(' ')
figure(3)
plot(Z)
title('系统输出z（k）图形')
%---------------------------------------------------------最小二乘法参数辨识
LS=inv(H'*H)*H'*ZL;
disp(' 运用最小二乘法对这一系统参数进行辨识，计算得到：')
disp([' a1(实际值) = -1.5，a1(辨识值)= ' num2str(LS(1,1))])
disp([' a2(实际值) = 0.7，a2(辨识值)= ' num2str(LS(2,1))])
disp([' b1(实际值) = 1，b1(辨识值)= ' num2str(LS(3,1))])
disp([' b2(实际值) = 0.5，b2(辨识值)= ' num2str(LS(4,1))])
disp(' ')
end
4. 仿真结果
本系统输入采用4 阶M 序列且其幅度为1，选取k=1 到k=16 数据，如图1所示系统输入图形。

图1 系统输入图形
输出受到白噪声序列v(k)污染，已知白噪声序列v(k)服从正态分布，均值为0，ds=2的白噪声序列图形如图2所示。

通过仿真，得到系统输出z(k)如图3所示。

图2系统白噪声图形（ds=2）图3 系统输出图形（ds=2）
运用最小二乘法对系统参数进行辨识，得到如图4所示的辨识结果。

图4 辨识结果截图（ds=2）
ds=1的白噪声序列图形如图5所示。

通过仿真，得到系统输出z(k)如图6所示。

图5系统白噪声图形（ds=1）图6 系统输出图形（ds=1）
图7 辨识结果截图（ds=1）
ds=1的白噪声序列图形如图8所示。

通过仿真，得到系统输出z(k)如图9所示。

图8 系统白噪声图形（ds=0）图9 系统输出图形（ds=0）
图10 辨识结果截图（ds=0）
5.结论
在研究系统辨识问题时，将把待辨识的系统看做“黑箱”，只考虑系统的输入输出特性，而不强调系统的内部机理。

本文即采用最小二乘法根据系统输入与输出特性对系统参数进行辨识，得到的辨识结果与实际值接近。

从ds=2和ds=1的辨识结果来看，a1、a2吻合的较好，b1、b2的误差较大。

选取白噪声序列方差ds=0，再次仿真得到如图10所示的辨识结果。

不难发现，当方差ds=0时，即没有系统白噪声干扰时，系统辨识参数与实际参数完全一致，进而可得系统辨识的效果与系统噪声有关。

噪声对参数估计有影响，噪声方差增大，估计值的偏差也增大。