利用相关分析法辨识脉冲响应实验报告..docVIP

实验一 利用相关分析法辨识脉冲响应 1.1 生成输入数据和噪声 用M序列作为辨识的输入信号，获得输出数据y(k)。采取串联传递函数仿真，，用M序列作为辨识的输入信号 其中r为周期数，表示计算互相关函数所用的数据是从第二个周期开始的，目的是等过程仿真数据进入平稳状态。 1.4 计算脉冲响应估计值、脉冲响应理论值、脉冲响应估计误差 脉冲响应估计值 脉冲响应理论值 脉冲响应估计误差 1.5 计算噪信比 信噪比 2 编程说明 M序列中，M序列循环周期取，时钟节拍=1Sec，幅度，特征多项式为。白噪声循环周期为。采样时间设为1Sec， 3 源程序清单 3.1 均匀分布随机数生成函数 function sita=U(N) %生成N个[0 1]均匀分布随机数 A=179; x0=11; M=2^15; for k=1:N x2=A*x0; x1=mod(x2,M); v1=x1/(M+1); v(:,k)=v1; x0=x1; end sita=v; end 3.2 正态分布白噪声生成函数 function v=noise(aipi) %生成正态分布N(0,sigma) sigma=1; %标准差 for k=1:length(aipi) ksai=0; for i=1:12 temp=mod(i+k,length(aipi))+1; ksai=ksai+aipi(temp); end v(k)=sigma*(ksai-6); end end 3.3 M序列生成函数 function [Np r M]=createM(n,a) %生成长度为n的M序列,周期为Np，周期数为r x=[1 1 1 1 1 1]; %初始化初态 for i=1:n y=x; x(2:6)=y(1:5); x(1)=xor(y(5),y(6)); U(i)=y(6); end M=U*a; lenx=length(x); Np=2^lenx-1; r=n/Np; end 3.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:n x(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); end end 3.5 相关函数计算函数 function R_Mz=RMz(Np,r,u,z) r=r-1; y=zeros(1,Np); for k=1:Np y(k)=0; for i=Np+1:(r+1)*Np y(k)=y(k)+u(i-k)*z(i); end y(k)=y(k)/(r*Np); end R_Mz=y; end 3.5 主函数 function [og yita]=main(time) % 脉冲响应估计误差 og % 噪信比yita N=time*63; 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序列驱动，生成y z=y+v; R_Mz=RMz(Np,r,u,z); %计算相关函数 % 计算脉冲响应估计值 g_k=zeros(1,Np); for k=1:Np g_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:Np Eg(1,k)=K/(T1-T2)*(exp(-k*T0/T1)-exp(-k*T0/T2)); end % 计算脉冲响应估计误差 og=sqrt(norm(Eg-g_k)^2/norm(Eg)^2); ov=fangcha(v); %计算噪声