非参数模型辨识：脉冲响应曲线

1.1 问题描述

1.2 方法思路
首先根据系统脉冲响应函数G(t)计算出51(0-50)个真值数据G,再利用以下公式产生550组带噪声的仿真数据。

最后利用上述仿真矩阵构造相应的数据矩阵，再用最小二乘法求解系统的脉冲响应估计值，其构造公式如下所示。

1.3 实验结果

1.4 实验代码（Matlab）

main.m

%%%采用最小二乘法辨识系统的脉冲响应向量，并与表达式计算出来的真值做比较
%%%第-步：根据系统脉冲响应函数G(t)计算出51(0-50)个真值数据G
for k=1:1:51 
    %%%矩阵下标不能为0，令k=t+1,即0时刻时的下标设为1
    G(k)=0.25*exp(-0.5*(k-1))-exp(-0.25*(k-1))+0.75*exp(-(k-1)/6);%产生第0时刻到50时刻的真实数据
end  

%%%第二步:根据课本上式4-1-8计算带噪声的仿真数据
%%采用M序列作为系统的输入，由于要产生550个仿真结果，为了充分激励系统的模态，选择M序列的周期为2^10-1
u=F2(1023,1);%%利用作业2中代码封装成的函数产生周期为1023的M序列
%根据题目要求产生550组仿真数据，而真实数据为51个，故其构造的U矩阵为550*51
for m=1:1:550  
    %每一行的数据由u(p+m)到u(m)组成
    for p=1:1:51
            U(m,p)=u(51+m-p);%%%u(p+m)到u(m)
    end
end  
%%产生均值为0，方差为0.05的正态分布白噪声V
V=F1(0.05,0,550);%%利用作业1中代码封装成的函数产生白噪声
%%由课本上式4-1-7可得含噪声的仿真输出
y=U*G'+V';
%%%第三步：采用最小二乘一次完成法求解系统的脉冲响应估计值并与真实值比较
G1=inv(U'*U)*U'*y;%最小二乘估计值
%画图
figure 
t=[0:1:50];
plot(t,G,'r','linewidth',2); %%画出真实数据
xlim([-5 55])
hold on
plot(t,G1','b-.','linewidth',2);  %%画出最小二乘估计数据
legend('真值','估计值')
title('脉冲响应曲线比较图')

F1.m

function Y3=F1(v,m,n)
%利用乘同余法和变换抽样法产生均值为m，方差为v的正态分布的白噪声
%公式：Xi=A*Xi*mod(M)
M1=2^11;%为2的方幂
M2=2^17;%为2的方幂
A1=119;%%不能太小
A2=279;%%不能太小
x1=1;%伪随机序列1的初值
x2=11;%伪随机序列2的初值
X1=x1;
X2=x2;
Y1=[];%伪随机序列1
Y2=[];%伪随机序列1
% n=700; %伪随机序列长度
%%%%%采用乘同余法产生伪随机序列%%%%%
for i=1:n
    %产生伪随机序列1
    X1=mod(A1*X1,M1);
   Y1=[Y1 X1/M1];
    %产生伪随机序列2
    X2=mod(A2*X2,M2);
    Y2=[Y2 X2/M2];
end
%正态分布的白噪声
Y3=m+v*sqrt(-2*log(Y1)).*cos(2*pi*Y2);%正态分布的白噪声
figure
plot(Y3,'b');
title('正态分布的白噪声');
end

F2.m

function X=F2(T,a)
%%%采用伪随机信号发生器生成周期为1023的M序列
y=[1 1 0 1 1 0 1 0 1 0];%%%初始化10级移位寄存器，保证其初值不全为0
for i = 1:1:T
    Y(i)=y(10);%%%第n级的输出即是M序列的两种状态
    M(i)=mod(y(10)+y(7),2);%第n级与第N1级模二加结果，此时分别为10和7
    %%%%移位寄存器工作原理%%%
    y(10)=y(9);
    y(9)=y(8);
    y(8)=y(7);
    y(7)=y(6);
    y(6)=y(5);
    y(5)=y(4);
    y(4)=y(3);
    y(3)=y(2);
    y(2)=y(1);
    y(1)= M(i);
end
for i=1:1:T
    if Y(i)==1
        X(i)=a;
    else
        X(i)=-a;
    end
end
stairs(X,'b'); %画出伪随机信号
title('周期为1023、幅值为1的伪随机信号');
ylim([-1.2 1.2])
end