pulse.m

根据脉冲序列

function   Gs=Pulse(G,T)  %根据脉冲响应数据，辨识脉冲传递函数

clear all
G=[0 0.3025 0.4183 0.3008 0.1226 -0.086];
T=1;
syms x s
k=max(size(G));
n=k/2;
n1=n;
n2=n;
n3=n;
n4=n;
n5=n;

i1=1;
i3=1;
i4=1;
i5=1;

j1=1;
j2=1;
j3=0;

%构造矩阵H
for i1=1:n1
    for j1=j2:n1
       H(i1,j1-j3)=G(j1+1);
    end 
    j2=j2+1;
    j3=j3+1;
    n1=n1+1;
end


%构造矩阵g1
for i2=1:n2
     g1(i2)=G(i2);
end
%确定待定系数a(i)
a=-inv(H)*g1';

%求特征方程的单根x(i)
i2=1:n3;
for i3=1:n3
    x1(i3)=x^(i2(i3)*T);
end  
%求传递函数的极点s
x2=double(solve(x1*a+1));
for i4=1:n4
    s1(i4)=log10(x2(i4));
end
%求系数c
s2=vander(s1);
s3=rot90(s2);
b=inv(s3)*g1';
c=b;


%确定传递函数
for i5=1:n5
    fs6(i5)=c(i5)/(s-s1(i5))
    
end


Gs=sum(fs6)