第三部分：相关改进

基本采用（二）中的做法，只是在y的求解摒弃了原理中介绍的方法，直接用matlab自带的lsim函数求解，使本例更有通用性。

源代码：

function mytemp
clear all
close all
clc

%创建M序列
Np=63;%循环周期
delta_T = 1;%时钟节拍
a=1;%幅度
%初始化M序列
M(1)=1;
M(2)=0;
M(3)=0;
M(4)=1;
M(5)=1;
M(6)=0;
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
figure(4)
stairs(M_XuLie)
ylim([-2 2])
%生成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
T = (0:length(u)-1)*delta_T;

num = [K];
den = conv([T1,1],[T2,1]);
y = lsim(num,den,u,T);
y = y';

%获取没有白噪声时候输出完毕
%作图
figure(1);
plot(u,'r');
hold on;
plot(y,'b');
legend('u(k)','y(k)');
%第二步，将白噪声添加入输出信号
%产生白噪声信号v
fangcha = 0.5;%随意指定的方差
v = fangcha * randn(1,63*4);
%信号叠加,输出实际信号z(k)
z = y + v;
%figure(2);

%plot(v,'m');
%hold on;

%plot(z,'y');

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');
figure(3)
plot(g1,'r');

运行结果：

http://s1/middle/84024a4a4cd0b6d8934a0&690

http://s2/middle/84024a4a4cd0b6d97e4b1&690

http://s16/middle/84024a4a4cd0b6da7264f&690