ssi程序有点问题，请高手检查！！！

luyunjun023

%hankel矩阵
clear
clc
close all hidden
d=100     %采样频率
n1=20     %行(2*n1)
n2=40     %列
load Y.txt %响应数据数列
y=Y';
yy=zeros(2*n1,n2);    %空hankel矩阵
for i=1:n2
    yy(:,i)=y(i:i+2*n1-1) %hankel矩阵第i列
end
hankel_a=yy/((n2)^0.5)     %最后的hankel矩阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求奇异值
yp=hankel_a(1:n1,:);
yf=hankel_a(n1+1:2*n1,:);
pref=yf*yp'*pinv(yp*yp')*yp;
[r,c]=size(pref);
pi=pref(1:r-1,:);   
[p0 d0 q0]=svd(pi); %奇异值分解得d0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%定阶
%假如阶数为i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%参数识别
pr=p0(:,1:i);    %U1
dd=d0(1:i,1:i);  %S1
dr=dd^0.5;       %S1^0.5
qr=q0(:,1:i);    %V1’

oi=pr*dr;        %U1*S1^0.5
xi=pinv(oi)*pi;

%compuet xi1
clear yf p0 d0 q0 pr dd dr qr

pi=pref(2:r,:);   
[p0 d0 q0]=svd(pi);
xi1=pinv(oi)*pi;

a=xi1*pinv(xi);  %A

[v,z]=eig(a);
z1=diag(z);  
z2=log(z1);  
fr=abs(z2)*d;
damp=real(z2)./fr*d;

luyunjun023

算出来结果和实际值差别很大

luyunjun023

怎么没有人来了？？？？都回家过年了吗

luyunjun023

高手，你在哪里？？？

huazi071783

楼主不知道你这个信号Y是几个通道的数据，看你的语句Y就是单通道的信号了，单个通道数据怎么能识别出阵型呢？