随机子空间法函数（代码求修正）

blackdot

function [A,C]=SSI_fun(h,n1,n2,type)
[r,c]=size(h);
if r<c
   h=h';
   [r,c]=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;        
        [u s v]=svd(Pref); ss=diag(s); [r1 c1]=size(s); n=max(find(ss));u1=u(:,1:n); s1=s(1:n,1:n);  Oi=u1*s1^0.5;
        [ud sd vd]=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
        [U1 S1 V1]=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
不知正确与否，特别是协方差驱动，感觉有些怪怪的。。。紧接着要尝试稳定图的编写了

是否进行修改，数据驱动识别的结果也存在较大差距

blackdot

不知正确与否，特别是协方差驱动，感觉有些怪怪的。。。紧接着要尝试稳定图的编写了

_gully_frog

_gully_frog

好简单的代码啊