随机子空间法函数（代码求修正）
function =SSI_fun(h,n1,n2,type)=size(h);
if r<c
h=h';
=size(h);
end
%n1取值为阶数，形成hankel矩阵
%n2取值为取用数据点数
switch type
case 1 %数据驱动子空间法
for i=1:2*n1
hankely(c*(i-1)+1:c*i,:)=h(i:i+n2-1,:)';
end
hankely=hankely/n2^0.5;
Yp=hankely(1:c*n1,:); Yf=hankely(c*n1+1:c*2*n1,:);Pref=Yf*Yp'*pinv(Yp*Yp')*Yp;
Yp1=hankely(1:c*(n1-1),:);Yf1=hankely(c*(n1-1)+1:c*2*n1,:); Pref1=Yf1*Yp1'*pinv(Yp1*Yp1')*Yp1;
=svd(Pref); ss=diag(s); =size(s); n=max(find(ss));u1=u(:,1:n); s1=s(1:n,1:n);Oi=u1*s1^0.5;
=svd(Pref1); ssd=diag(sd); nd=max(find(ssd(1:r1)));ud1=ud(:,1:nd); sd1=sd(1:nd,1:nd);Oii=ud1*sd1^0.5;
xi=pinv(Oi)*Pref; xi1=pinv(Oii)*Pref1;
A=xi1*pinv(xi);Yi=hankely(n1*c+1:(n1+1)*c,:);C=Yi*pinv(xi);
case 2 %协方差驱动子空间法
m=n2; p=n1;q=ceil(n2/c);
for i=1:p
for j=1:q
k1=i+j-1; k2=i+j;
r1=((h(1:m,:))'*h(k1+1:m+k1,:))/m;
hankelr1((i-1)*c+1:i*c,(j-1)*c+1:j*c)=r1;
end
end
=svd(hankelr1); s11=diag(S1); n1=max(find(s11));
U1=U1(:,1:n1); S1=S1(1:n1,1:n1);
Oi=U1*S1^0.5; Oi1=Oi(c+1:end,:); Oi2=Oi(1:end-c,:); A=pinv(Oi2)*Oi1; C=Oi2(1:c,:);
otherwise
return
end
以上是我根据文献参考，自己编写的随机子空间法函数代码。感觉两种驱动方法识别出的固有频率区别挺大。
想请各位老师和高手指点，代码有哪些要修正的地方？谢谢
blackdot 发表于 2014-11-4 14:49
不知正确与否，特别是协方差驱动，感觉有些怪怪的。。。紧接着要尝试稳定图的编写了
是否进行修改，数据驱动识别的结果也存在较大差距 不知正确与否，特别是协方差驱动，感觉有些怪怪的。。。紧接着要尝试稳定图的编写了 {:{41}:}{:{26}:}{:{13}:}{:{13}:} 好简单的代码啊{:{26}:}
页:
[1]