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所示系统输入图形。