北京工商大学
《系统辨识》课程
实验报告
（2014-2015 1学期）
课程名称：系统辨识
题目：利用相关分析法辨识脉冲响应
专业班级：控制工程
学生姓名：
指导教师：刘刘
成绩：
2015年1月18日
一、实验目的
通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。
二、实验内容
图1为本实验的原理框图。
过程传递函数为)
(s
G,其中
Sec
2
6
T
Sec,
3
8
120
2
1
.
.
,=
=
=T
K；)
(
)
(k
z
k
u和分别为过程的输入和输出变量；)
(k
v
为过程测量白噪声，服从正态分布，均值为零，方差为2
v
σ，记作)
,
(
~
)
(2
v
N
k
vσ；
)
(k
g
为过程的脉冲响应理论值，)
(
ˆ
k
g为过程脉冲响应估计值，)
(~k
g为过程脉冲响应估计误差。
过程的输入驱动采用M序列，输出受到白噪声)
(k
v的污染。
根据过程的输入和输出数据{})(
),
(k
z
k
u，利用相关分析算法根据输出过程的脉冲响应值)
(
ˆ
k
g，并与过程脉冲响应理论值)
(k
g
比较，得到过程脉冲响应估计误差值)
(~k
g，当∞
→
k
时，应该有0
→
)
(~k
g。
图1 相关分析法辨识脉冲响应原理框图
三、实验要求
进行方案设计，模拟过程传递函数，获得输出数据，用M序列作为辨识的输入信号，噪声采用标准正态分布的白噪声，计算互相关函数，不同λ值的脉冲响应估计值、脉冲响应理论值和脉冲响应估计误差，计算信噪比，画出实验流程图，用MATLAB编程实现。
四、实验原理
相关分析法
v(k)
u(k) z(k)
)1
)(
1
(
)
(
2
1
+
+
=
s
T
s
T
K
s
G
y(k)
1、采用串联传递函数)
(s
G仿真
2
1
2
1
1
1
1
1
T
s
T
s
T
T
K
s
G
/
/
)
(
+
+
=
令
2
1
1T
T
K
K＝，则)
(s
G的表达框图为：
2、一个单输入单输出线性定常系统的动态特性可用它的脉冲响应函数g(σ)来描述。
这样，只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数g(τ)。
而在系统有正常输入的情形下，辨识脉冲响应的原理图如下图所示。
()()()
y t g x t d
σσσ
∞
=-
⎰
则
000
()
11
lim()()(){lim()()}
T T
T T
x t
y t x t dt g x t x t dt d
T T
τ
τσστσ
∞
→∞→∞
-
-=--
⎰⎰⎰
上式两端同乘，进而取时间均值，有
()()()
xy x
R g R d
τστσσ
∞
=-
-
⎰
则
这就是著名的维纳霍夫积分方程。
()
()(),()()
()()()()
()
()
x x
xy x
xy
x t
R k R k
R g R d kg
R
g
k
τδττσδτσ
τστσστ
τ
τ
∞
=-=-
-
=-=
=
⎰
如果输入是，这时的自相关函数为
则根据维纳霍夫积分方程可得
或者
白噪声
1
1
/1T
s
K
+
u(k) x(k)
2
1
1
T
s/
+
y(k)
五、实验框图
六、实验代码
function ex2
clc;
clear all;
close all;
%创建M序列
Np=63;%循环周期
delta_T = 1;%时钟节拍
a=1;%幅度
M(1)=1;
M(2)=0;
M(3)=0;
M(4)=1;
M(5)=1;
M(6)=0;%初始化M序列
M_XuLie(Np) = 0;
for n = 1 : Np
temp = xor(M(6), M(5));
if(temp == 0)
M_XuLie(n) = a;
else
M_XuLie(n) = -a;
end
M(6) = M(5);
M(5) = M(4);
M(4) = M(3);
M(3) = M(2);
M(2) = M(1);
M(1) = temp;
end
%生成M序列完毕
r=3;%周期数
u=repmat(M_XuLie,1,r+1);%将M序列赋给输入，作为输入信号%第一步，从u(k)得到x(k),y(k)
K = 120;
T0 = 1; % 采样时间
T1 = 8.3;
T2 = 6.2;
K1=K/(T1*T2);
%初始化X(k),Y(k)为0
K2=1
x(63)=0;
y(63)=0
for k = 2 : 63*4
%取得x(k)序列
x(k)=exp(-T0/T1)*x(k-1)+T1*K1*(1-exp(-T0/T1))*u(k-1)+T1*K1... *(T1*(exp(-T0/T1)-1)+T0)*(u(k)-u(k-1))/T0
%取得y(k)序列
y(k)=exp(-T0/T2)*y(k-1)+T2*K2*(1-exp(-T0/T2))*x(k-1)+T2*K2... *(T2*(exp(-T0/T2)-1)+T0)*(x(k)-x(k-1))/T0
end
%获取没有白噪声时候输出完毕
%作图
figure(1);
plot(u,'r');
hold on;
plot(x,'k');
plot(y,'b');
legend('u(k)','x(k)','y(k)');
%第二步，将白噪声添加入输出信号
%产生白噪声信号v
fangcha = 0.5;%随意指定的方差
v = fangcha * randn(1,63*4);
%信号叠加,输出实际信号z(k)
z = y + v;
figure(2);
%打印无白噪声污染信号
plot(y,'b');
hold on;
%打印白噪声信号
plot(v,'m');
%打印白噪声污染后的信号
plot(z,'k');
legend('y(k)','v(k)','z(k)');
%计算Rmz(k)
for k = 1 : Np
Rmz(k)=0;%初始化为0
for i = (Np + 1) : ((r+1)*Np)
Rmz(k)=Rmz(k) + u(i-k)*z(i);
end
Rmz(k)=Rmz(k)/(r*Np);
end
%计算c
c=-Rmz(Np - 1);
%计算脉冲响应估计值g1
g1=Np*(Rmz+c)/((Np+1)*a^2*delta_T);
%计算理论脉冲g0
for k = 1: Np
g0(k)=K/(T1-T2)*(exp(-k*delta_T/T1)-exp(-k*delta_T/T2));
end
%计算脉冲响应估计误差delta_g
delta_g=sqrt(sum((g0-g1).^2)/sum(g0.^2));
figure(3);
plot(g0,'k');
hold on;
plot(g1,'r');
%axis([0,100,0,10]);
legend('脉冲响应理论值g0(k)','脉冲响应估计值g1');
七、实验结果
1、输入u(k),中间输入x(k)，无干扰输入(k）
2、白噪声标准差为1.5时，理想输出y(k)，带干扰的输出z(k)，干扰v(k)
3、输入白噪声标准差为1.5,周期数r为3时，脉冲响应理论值与估计值:
脉冲响应估计误差: 0.0467
八、实验结论
1、根据维纳-霍夫积分方程，只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数；
2、通过仿真，看到白噪声方差越大，实际输出结果的偏差也就越大；
3、周期数越大脉冲响应的估计值与理论值越接近，同时会增大数据量。
可以证明当k很大时，误差趋于0。