利用相关分析法辨识脉冲响应1 生成输入数据和噪声用M 序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。
生成白噪声时，首先利用乘同余法生成U[0,1]均匀分布的随机数，再利用U[0,1]均匀分布的随机数生成标准正态分布的白噪声。
白噪声循环周期为15232768=。
2 过程仿真模拟过程传递函数)(s G ，获得输出数据y(k)。
)(s G 采取串联传递函数仿真，21211111)(T s T s T T K s G ++=，用M 序列作为辨识的输入信号。
()G s 采样时间0T 设为1Sec ，12120, 8.3Sec, 6.2Sec K T T ===(1) 惯性环节其中，T 为惯性环节的时间常数，K 为惯性环节的静态放大倍数。
若采样时间记作0T ，则惯性环节的输出可写成：[]011111000T k u k u T e T TK k u e TK k y e k y T T T T T T )()())()()()()(///--+-+--+-=---(2) 传递函数()G s 仿真（串联）21211111T s T s T T K s G //)(++=令112KK T T ＝，则()G s 的表达框图为：[]011111000T k u k u T e T TK k u e TK k y e k y TT T T T T )()())()()()()(///--+-+--+-=---编程语句可写成：[][][][]};);()();()();()(;/)()(*)(**)(*)(*)(*)(;/)()(*)(***)(*)(**)(*)({);;(;)(;)();/();/();*/(k y k y k x k x k u k u T k x k x T E T T k x E T k y E k y T k u k u T E T K T k u E K T k x E k x k k k y x T T E T T E T T K K =-=-=---+-+--+-=--+-+--+-=++<===-=-==11111111111112521for 0000EXP EXP 0022222200111111112021012113、白噪声生成● 利用乘同余法生成U[0,1]均匀分布的随机数)(,)(mod )(,,],[~)(mod ,奇数循环周期其中118317923276821002151=====⎪⎩⎪⎨⎧==-+x A M U Mx M Ax x k i i i i ξ● 利用U[0,1]均匀分布的随机数生成正态分布的白噪声),(~)(212106v i i v N k v σξσ⎪⎭⎫⎝⎛-=∑= 其中，标准差v σ分别取0,0.1,0.5。
● 编程语句};);.(*};);();;*{);;(;{),;(06FLOAT ,MOD(121for 02521for -=+===++<==++<=ksai Sigm a v(k)xi/M ksai ksai M xi xi xi A xi i i i ksai k k k 4、M 序列生成● 用M 序列作为辨识的输入信号，N 序列的循环周期取62163P N =-=,时钟节拍1Sec t ∆=，幅度1a =,逻辑“0”为a,逻辑“1”为-a ，特征多项式自选，如65()1F s s s =⊕⊕。
● 生成M 序列的结构图● 编程语句};;)()(;)()(};);()({);;(;)()();()()({),;(a k u M a k u M i M i M i i P i for M M M M M k k k -======-=-->====+=++<=then 10if then 00if 1000then 20if 2102521for 5、互相关函数的计算∑++=-=PP N r N i PMz i z k i u rN k R )()()()(111其中，r 为周期数，1P i N =+表示计算互相关函数所用的数据是从第二个周期开始的，目的是等过程仿真数据进入平稳状态。
6、c 的补偿补偿量c 应取(1)Mz P R N -－，不能取()Mz P R N －。
因为()Mz R k 是周期函数，则有()(0)Mz P Mz R N R ＝，故不能取()Mz P R N －。
7、计算脉冲响应估计值● 脉冲响应估计值 []2ˆ()()(1)PMz P N gk R k c N a t=++∆ ● 脉冲响应理论值[]21//210)(T t k T t k e e T T Kk g ∆-∆---=● 脉冲响应估计误差()()∑∑==-=ppN k N k g k g k gk g12012)()(ˆ)(δ8 源程序清单8.1 均匀分布随机数生成函数function sita=U(N)%生成N个[0 1]均匀分布随机数A=179; x0=11; M=2^15;for k=1:Nx2=A*x0;x1=mod(x2,M);v1=x1/(M+1);v(:,k)=v1;x0=x1;endsita=v;end8.2 正态分布白噪声生成函数function v=noise(aipi)%生成正态分布N(0,sigma)sigma=1; %标准差for k=1:length(aipi)ksai=0;for i=1:12temp=mod(i+k,length(aipi))+1;ksai=ksai+aipi(temp);endv(k)=sigma*(ksai-6);endend8.3 M序列生成函数function [Np r M]=createM(n,a)%生成长度为n的M序列,周期为Np，周期数为r x=[1 1 1 1 1 1]; %初始化初态for i=1:ny=x;x(2:6)=y(1:5);x(1)=xor(y(5),y(6));U(i)=y(6);endM=U*a;lenx=length(x);Np=2^lenx-1;r=n/Np;end8.4 过程仿真函数function y=createy(u,K,T1,T2,T0)n=length(u);K1=K/(T1*T2);E1=exp(-T0/T1);E2=exp(-T0/T2);x(1)=0;y(1)=0;for k=2:nx(k)=E1*x(k-1)+T1*K1*(1-E1)*u(k-1)...+T1*K1*(T1*(E1-1)+T0)*(u(k)-u(k-1))/T0;y(k)=E2*y(k-1)+T2*(1-E2)*x(k-1)...+T2*(T2*(E1-1)+T0)*(x(k)-x(k-1))/T0;u(k-1)=u(k);x(k-1)=x(k);y(k-1)=y(k);endend8.5 相关函数计算函数function R_Mz=RMz(Np,r,u,z)r=r-1;y=zeros(1,Np);for k=1:Npy(k)=0;for i=Np+1:(r+1)*Npy(k)=y(k)+u(i-k)*z(i);endy(k)=y(k)/(r*Np);endR_Mz=y;end8.6 主函数N=252;K=120; T1=8.3; T2=6.2; T0=1; a=1;sita=U(N); %生成[0 1]均匀分布随机数v=noise(sita); %利用aipi生成正态分布白噪声[Np r u]=createM(N,a); %生成长度为N的M序列y=createy(u,K,T1,T2,T0); %利用M序列驱动，生成yz=y+v;R_Mz=RMz(Np,r,u,z); %计算相关函数% 计算脉冲响应估计值g_k=zeros(1,Np);for k=1:Npg_k(1,k)=(R_Mz(1,k)-R_Mz(Np-1))*Np/((Np+1)*a*a*T0);end% 计算脉冲响应理论值Eg=zeros(1,Np);for k=1:NpEg(1,k)=K/(T1-T2)*(exp(-k*T0/T1)-exp(-k*T0/T2));endfigure(1)plot([1:252],y,[1:252],z);legend('不含噪声的输出序列','含噪声的输出序列');figure(2)plot([1:63],g_k,[1:63],Eg);legend('脉冲响应估计值','脉冲响应理论值');9 数据记录表1脉冲响应估计值与脉冲响应理论值的比较t 1 2 3 4 5 6 7脉冲响应估计值0.79 0.92 1.02 1.04 1.05 1.01 0.92 脉冲响应理论值 2.03 3.52 4.59 5.32 5.77 6.02 6.11 t 8 9 10 11 12 13 14脉冲响应估计值0.87 0.80 0.74 0.65 0.57 0.50 0.42脉冲响应估计值0.33 0.23 0.17 0.10 0.05 -0.01 -0.06 脉冲响应理论值 4.29 3.99 3.69 3.40 3.12 2.86 2.62 t 22 23 24 25 26 27 28脉冲响应估计值-0.10 -0.16 -0.19 -0.22 -0.25 -0.29 -0.28 脉冲响应理论值 2.39 2.18 1.98 1.80 1.63 1.48 1.33 t 29 30 31 32 33 34 35脉冲响应估计值-0.30 -0.31 -0.32 -0.36 -0.37 -0.39 -0.41 脉冲响应理论值 1.20 1.09 0.98 0.88 0.79 0.71 0.64脉冲响应估计值-0.44 -0.46 -0.47 -0.46 -0.49 -0.51 -0.52 脉冲响应理论值0.58 0.52 0.46 0.41 0.37 0.33 0.30 t 43 44 45 46 47 48 49脉冲响应估计值-0.53 -0.54 -0.55 -0.55 -0.56 -0.54 -0.56 脉冲响应理论值0.27 0.24 0.21 0.19 0.17 0.15 0.13 t 50 51 52 53 54 55 56脉冲响应估计值-0.57 -0.57 -0.56 -0.57 -0.57 -0.56 -0.55 脉冲响应理论值0.12 0.11 0.10 0.09 0.08 0.07 0.06 t 57 58 59 60 61 62 63脉冲响应估计值-0.53 -0.52 -0.53 -0.52 -0.53 0.00 0.61 脉冲响应理论值0.05 0.05 0.04 0.04 0.03 0.03 0.0310 曲线打印图1 加入噪声前后的输出序列比较图2 脉冲响应理论值与估计值的比较作业：如图所示一阶系统，采用5级移位寄存器产生M序列作为输入信号，辨识该系统的脉冲响应。