利用相關分析法辨識脈沖響應

1、精品利用相關分析法辨識脈沖響應自 1205劉彬412511411 實驗方案設計1.1生成輸入數據和噪聲用M 序列作為辨識的輸入信號，噪聲采用標準正態(tài)分布的白噪聲。生成白噪聲時，首先利用乘同余法生成U0,1 均勻分布的隨機數，再利用U0,1 均勻分布的隨機數生成標準正態(tài)分布的白噪聲。1.2過程仿真模擬過程傳遞函數G(s) ，獲得輸出數據 y(k) 。 G( s) 采取串聯(lián)傳遞函數仿真，G (s)K11，用 M 序列作為辨識的輸入信號。T1T2 s1 s1T1T21.3計算互相關函數( r 1) N p1u(i k) z(i )RMz (k)rN p i N p1其中 r為周期數， iN p 1

2、表示計算互相關函數所用的數據是從第二個周期開始的，目的是等過程仿真數據進入平穩(wěn)狀態(tài)。1.4計算脈沖響應估計值、脈沖響應理論值、脈沖響應估計誤差?N pRMz (k)RMz ( N p 1)脈沖響應估計值 g (k)1)a2( N pt-可編輯 -精品脈沖響應理論值 g0 (k)Kek t / T1e k t /T2脈沖響應估計誤差T1T2N pg0 (k )g?(k )2gk1N pg0 ( k)k 121.5計算噪信比信噪比y(k)yv( k)v222 編程說明M 序列中， M 序列循環(huán)周期取N p 26 163 ，時鐘節(jié)拍t =1Sec ，幅度a 1 ，特征多項式為 F( s) s6s51

3、 。白噪聲循環(huán)周期為 21532768 。 G(s) 采樣時間 T0 設為 1Sec ， K120, T18.3Sec, T26.2Sec3 源程序清單3.1均勻分布隨機數生成函數function sita=U(N)%生成 N 個0 1 均勻分布隨機數A=179; x0=11; M=215;for k=1:Nx2=A*x0;x1=mod(x2,M);-可編輯 -精品v1=x1/(M+1);v(:,k)=v1;x0=x1;endsita=v;end3.2正態(tài)分布白噪聲生成函數function v=noise(aipi)%生成正態(tài)分布 N(0,sigma)sigma=1; % 標準差for k=1

4、:length(aipi)ksai=0;for i=1:12temp=mod(i+k,length(aipi)+1;ksai=ksai+aipi(temp);endv(k)=sigma*(ksai-6);endend-可編輯 -精品3.3 M序列生成函數function Np r M=createM(n,a)%生成長度為 n的 M 序列 ,周期為 Np ，周期數為 rx=1 1 1 1 1 1; %初始化初態(tài)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=2lenx-1;r=n

5、/Np;end3.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(

6、k-1)=y(k);endend3.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);end-可編輯 -精品y(k)=y(k)/(r*Np);endR_Mz=y;end3.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=

7、noise(sita); %利用 aipi 生成正態(tài)分布白噪聲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*T

8、0/T1)-exp(-k*T0/T2);end% 計算脈沖響應估計誤差og=sqrt(norm(Eg-g_k)2/norm(Eg)2);ov=fangcha(v); %計算噪聲方差oy=fangcha(y); %計算信號方差yita=sqrt(oy/ov); %計算信噪比 End3.5畫圖函數 1%mainPlot.mfigure(1)for n=4:40og yita=main(n);y1(n)=og;endy1=y1(4:40);plot(4:40,y1);xlabel(' 周期數 ');ylabel(' 脈沖響應估計誤差 ');-可編輯 -精品figure

9、(2)for n=4:40og yita=main(n);y2(n)=yita;endy2=y2(4:40);plot(4:40,y2);xlabel(' 周期數 ');ylabel(' 噪信比 ');3.5畫圖函數 2%mainPlot2.mN=252;K=120; T1=8.3; T2=6.2; T0=1; a=1;sita=U(N);%生成 0 1 均勻分布隨機數v=noise(sita); %利用 aipi 生成正態(tài)分布白噪聲Np r u=createM(N,a); %生成長度為 N 的 M 序列y=createy(u,K,T1,T2,T0); %利用

10、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(' 不含噪聲的輸出序列 ',' 含噪聲的輸出序

11、列 ');figure(2)plot(1:63,g_k,1:63,Eg);Legend(' 脈沖響應估計值 ',' 脈沖響應理論值 ');4 數據記錄表 1 脈沖響應估計值與脈沖響應理論值的比較t1234567脈沖響應估計值0.790.921.021.041.051.010.92-可編輯 -精品脈沖響應理論值2.033.524.595.325.776.026.11t891011121314脈沖響應估計值0.870.800.740.650.570.500.42脈沖響應理論值6.075.945.745.495.214.914.60t1516171819202

12、1脈沖響應估計值0.330.230.170.100.05-0.01-0.06脈沖響應理論值4.293.993.693.403.122.862.62t22232425262728脈沖響應估計值-0.10-0.16-0.19-0.22-0.25-0.29-0.28脈沖響應理論值2.392.181.981.801.631.481.33t29303132333435脈沖響應估計值-0.30-0.31-0.32-0.36-0.37-0.39-0.41脈沖響應理論值1.201.090.980.880.790.710.64t36373839404142脈沖響應估計值-0.44-0.46-0.47-0.46-

13、0.49-0.51-0.52脈沖響應理論值0.580.520.460.410.370.330.30t43444546474849脈沖響應估計值-0.53-0.54-0.55-0.55-0.56-0.54-0.56脈沖響應理論值0.270.240.210.190.170.150.13t50515253545556脈沖響應估計值-0.57-0.57-0.56-0.57-0.57-0.56-0.55-可編輯 -精品脈沖響應理論值0.120.110.100.090.080.070.06t57585960616263脈沖響應估計值-0.53-0.52-0.53-0.52-0.530.000.61脈沖響應理論值0.050.050.040.040.030.030.035 曲線打印圖1 信噪比隨著周期數增大的變化-可編輯 -精品圖2 脈沖響應計算誤差隨著周期數增大的變化圖3 加入噪聲前后的輸出序列比較-可編輯 -精品圖 4 脈沖響應理論值與估計值的比較6 結果分析6.1信噪比脈沖響應計算誤差隨周期的變化隨著周期數的增加，信噪比減小，說明噪聲隨著周期數的增加變得更強烈，而計算誤差的減小表示周期數的增加使得不確定因素的影響減小，使得計算結果