利用相关分析法辨识脉冲响应

利用相关分析法辨识脉冲响应

自1205 刘彬 41251141

1 实验方案设计

1.1 生成输入数据和噪声

用M 序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。 生成白噪声时，首先利用乘同余法生成U[0,1]均匀分布的随机数，再利用U[0,1]均匀分布的随机数生成标准正态分布的白噪声。 1.2 过程仿真

模拟过程传递函数)(s G ，获得输出数据y(k)。)(s G 采取串联传递函数仿真，

2

12111

11)(T s T s T T K s G ++=

，用M 序列作为辨识的输入信号。

1.3 计算互相关函数

∑++=-=

p

p N r N i p

Mz i z k i u rN k R )1(1

)()(1

)(

其中r 为周期数，1+=p N i 表示计算互相关函数所用的数据是从第二个周期开始的，目的是等过程仿真数据进入平稳状态。

1.4 计算脉冲响应估计值、脉冲响应理论值、脉冲响应估计误差

脉冲响应估计值[]

)1()()1()(?2

--?+=p Mz Mz

p

p

N R k R

t

a N N k g

脉冲响应理论值[]

21//2

10)(T t k T t k e e T T K

k g ?-?---=

脉冲响应估计误差

()()

∑∑==-=

p

p

N k N k g k g

k g

k g

1

2

1

2

)()(?)(δ 1.5 计算噪信比

信噪比()()2

2

)()(v k v y k y --=η

2 编程说明

M 序列中，M 序列循环周期取

63

126=-=p N ，时钟节拍t ?=1Sec ，幅度1=a ，

特征多项式为1)(56⊕⊕=s s s F 。白噪声循环周期为32768215=。

)(s G 采样时间0T 设为1Sec ，Sec 2.6 Sec,3.8 ,12021===T T K

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); %计算噪声方差

oy=fangcha(y); %计算信号方差

yita=sqrt(oy/ov); %计算信噪比End

3.5 画图函数1

%mainPlot.m

figure(1)

for n=4:40

[og yita]=main(n);

y1(n)=og;

end

y1=y1(4:40);

plot([4:40],y1);

xlabel('周期数');

ylabel('脉冲响应估计误差');

figure(2)

for n=4:40

[og yita]=main(n);

y2(n)=yita;

end

y2=y2(4:40);

plot([4:40],y2);

xlabel('周期数');

ylabel('噪信比');

3.5 画图函数2

%mainPlot2.m

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序列驱动，生成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

figure(1)

plot([1:252],y,[1:252],z);

Legend('不含噪声的输出序列','含噪声的输出序列');

figure(2)

plot([1:63],g_k,[1:63],Eg);

Legend('脉冲响应估计值','脉冲响应理论值');

4 数据记录

表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 脉冲响应理论值 6.07 5.94 5.74 5.49 5.21 4.91 4.60 t 15 16 17 18 19 20 21

脉冲响应估计值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 t 36 37 38 39 40 41 42

脉冲响应估计值-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.03

5 曲线打印

图1 信噪比随着周期数增大的变化

图2 脉冲响应计算误差随着周期数增大的变化

图3 加入噪声前后的输出序列比较

图4 脉冲响应理论值与估计值的比较

6 结果分析

6.1 信噪比脉冲响应计算误差随周期的变化

随着周期数的增加，信噪比减小，说明噪声随着周期数的增加变得更强烈，而计算误差的减小表示周期数的增加使得不确定因素的影响减小，使得计算结果与实际更接近。

6.2 加入噪声前后的输出序列比较

加入噪声前后的变化并不大，说明噪声对输出序列影响不大，在第二个周期之后输出序列变得稳定，具有周期变化。

6.3 脉冲响应理论值与估计值的比较

随着时间的增加，脉冲响应理论值与估计值全部变小，且差值变小，与实验前的理论推导一致。

7 实验体会

通过本次试验，我不仅更深层次的了解了系统辨识的内容，而且也复习和运用了matlab，为以后的工作实践打了一些基础。