查看: 911|回复: 5

[随机振动] 随机子空间法函数(代码求修正)

[复制链接]
发表于 2014-10-29 16:52 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
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

以上是我根据文献参考,自己编写的随机子空间法函数代码。感觉两种驱动方法识别出的固有频率区别挺大。
想请各位老师和高手指点,代码有哪些要修正的地方?谢谢

回复
分享到:

使用道具 举报

发表于 2021-9-9 20:24 | 显示全部楼层
blackdot 发表于 2014-11-4 14:49
不知正确与否,特别是协方差驱动,感觉有些怪怪的。。。紧接着要尝试稳定图的编写了

是否进行修改,数据驱动识别的结果也存在较大差距
 楼主| 发表于 2014-11-4 14:49 | 显示全部楼层
不知正确与否,特别是协方差驱动,感觉有些怪怪的。。。紧接着要尝试稳定图的编写了
发表于 2014-10-31 16:33 | 显示全部楼层
发表于 2014-10-31 16:32 | 显示全部楼层
好简单的代码啊

点评

反对: 5.0
反对: 5
小建议,多参与讨论,少说543  发表于 2014-11-5 13:21
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表