实验一利用相关函数辨识脉冲响应分析解析

%取得 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);
与过程脉冲响应理论值 g 0 (k ) 比较，得到过程脉冲响应估计误差值
~g ( k) ，当
k
时，应该有 g~(k) 0 。
v(k)
u(k)
K
y(k)
z(k)
G ( s)
( T1s 1)( T 2 s 1)
相关分析法
三、实验要求
图 1 相关分析法辨识脉冲响应原理框图
进行方案设计，模拟过程传递函数，获得输出数据，用 M序列作为辨识的输
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);
%第一步，从 u(k) 得到 x(k),y(k) K = 120;
T0 = 1; T1 = 8.3; T2 = 6.2;
二、实验内容
图 1 为 本 实 验 的 原 理 框 图 。 过 程 传 递 函 数 为 G( s) , 其 中
K 120, T1 8.3S e cT2, 6.2S e ；c u(k )和 z(k) 分别为过程的输入和输出变量； v(k )
为过程测量白噪声， 服从正态分布， 均值为零，方差为
2 v
，记作
% 采样时间
K1=K/(T1*T2);
%初始化 X(k),Y(k) 为 0 K2=1
x(63)=0;
y(63)=0
%将 M序列赋给输入，作为输入信号
5
北京工商大学
计算机与信息工程学院
2014-2015 1
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
1T
1T
lim
y(t) x(t )dt g ( ){lim
x(t ) x(t )dt} d
T T0
0
T T0
则
Rxy ( )
g(
0
) Rx (
)d
这就是著名的维纳 霍夫积分方程。
如果输入是 白噪声 ，这时 x(t)的自相关函数为 Rx ( ) k ( ), Rx( ) k ( )
则根据维纳 霍夫积分方程可得
T1T2 s 1 / T1 s 1 / T2
令 K1＝ K ，则 G ( s) 的表达框图为：
T1T2
u(k)
K1
x(k)
1
y(k)
s 1/ T1
s 1/ T2
2、一个单输入单输出线性定常系统的动态特性可用它的脉冲响应函数 来描述。
g( σ )
则 y(t)
g( ) x(t )d
0
上式两端同乘 x(t )，进而取时间均值，有
%打印无白噪声污染信号 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;
北京工商大学
《系统辨识》课程 实验 报告
（2014-2015 1 学期）
课程名称 ：
系统辨识
题 目 ： 利用相关分析法辨识脉冲响应
专业班级 ：
控制工程
学生姓名 ：
指导教师 ： 刘 刘
成 绩：
2015 年 1 月 18 日
一、实验目的
北京工商大学
计算机与信息工程学院
2014-2015 1
学期
通过仿真实验掌握利用相关分析法辨识脉冲响应的原理和方法。
学期
五、实验框图
4
北京工商大学
计算机与信息工程学院
2014-2015 1
学期
六、ቤተ መጻሕፍቲ ባይዱ验代码
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;
或者
Rxy ( )
g( ) Rx (
0
) d kg( )
这样，只要记录
g( ) Rxy ( ) k
x(t) 、 y(t) 的值，并计算它们的互相关函数，即可求得脉
冲响应函数 g( τ ) 。
而在系统有正常输入的情形下，辨识脉冲响应的原理图如下图所示。
3
北京工商大学
计算机与信息工程学院
2014-2015 1
入信号，噪声采用标准正态分布的白噪声，计算互相关函数，不同
值的脉冲响
应估计值、脉冲响应理论值和脉冲响应估计误差， 计算信噪比， 画出实验流程图，
用 MATLA编B 程实现。
四、实验原理
2
北京工商大学
计算机与信息工程学院
2014-2015 1
学期
1、采用串联传递函数 G( s) 仿真
K
1
1
G ( s)
七、实验结果
1、输入 u(k) , 中间输入 x(k) ，无干扰输入 (k ）
2、白噪声标准差为 1.5 时，理想输出 y(k) ，带干扰的输出 z(k) ，干扰 v(k)
7
北京工商大学
计算机与信息工程学院
2014-2015 1
学期
3、输入白噪声标准差为 1.5, 周期数 r 为 3 时，脉冲响应理论值与估计值 :
v(k )
~
N ( 0,
2 v
)
；
g0 (k ) 为过程的脉冲响应理论值， g?(k ) 为过程脉冲响应估计值， g~( k ) 为过程脉
冲响应估计误差。
过程的输入驱动采用 M序列，输出受到白噪声 v(k ) 的污染。根据过程的输入 和输出数据 u( k), z(k ) ，利用相关分析算法根据输出过程的脉冲响应值 g?( k) ，并
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);
脉冲响应估计误差 : 0.0467
八、实验结论
1、根据维纳 - 霍夫积分方程，只要记录 x(t) 、y(t) 的值，并计算它们的互 相关函数，即可求得脉冲响应函数；
2、通过仿真，看到白噪声方差越大，实际输出结果的偏差也就越大； 3、 周期数越大脉冲响应的估计值与理论值越接近，同时会增大数据量。可 以证明当 k 很大时，误差趋于 0。
8
%初始化为 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
6
学期 ... ...
北京工商大学
计算机与信息工程学院
2014-2015 1
学期
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' );