|
楼主 |
发表于 2010-4-30 09:54
|
显示全部楼层
随机子空间法程序共享
以下是我写的随机子空间法模态识别法的程序,信号是15通道长度为1024的动应变响应,识别出拱的应变模态和频率及阻尼比,得到的结果和实际的差异很大,得到的应变振型是一个150*150的复数矩阵,而不是实数矩阵,我现在不知道怎么来提取应变模态振型,请高手指点,谢谢!!
% 2n1是hankel矩阵的行数数, n2是hankel矩阵的列数, d是采样频率
load('strainsig4.mat'); %strainsig是1024长度的15通道动应变响应
h=strainsig4; [r,c]=size(h);
% r=1024; c=15 ;
n1=10; n2=1000; d=2000;
for j=1:n2, for i=1:2*n1
hankel1(c*(i-1)+1:c*i,j)=h(i+j-1,:)'; % generate Hankel matrix
end; end
hankel1= hankel1/n2^0.5;
% generate Hankel-matrix from R-matrix
yp=hankel1(1:c*n1,:); yf=hankel1(c*n1+1:c*2*n1,:);
yp1=hankel1(1:c*(n1-1),:); yf1=hankel1(c*(n1-1)+1:c*2*n1,:);
pref=yf*yp'*pinv(yp*yp')*yp; pref1=yf1*yp1'*pinv(yp1*yp1')*yp1;
[u s v]=svd(pref); [ud sd vd]=svd(pref1); [r1 c1]=size(s);
for i=1:r1
ss(i)=s(i,i); ssd(i)=sd(i,i); n=find(ss); nd=find(ssd); n=max(n); nd=max(nd);
end
u1=u(:,1:n); s1=s(1:n,1:n); u1d=ud(:,1:nd); s1d=sd(1:nd,1:nd);
v1=v(:,1:n); oi=u1*s1^0.5; oii=u1d*s1d^0.5;
xi=pinv(oi)*pref; xii=pinv(oii)*pref1;
%compute A C
A=xii*pinv(xi); Yii=hankel1(n1*c+1:(n1+1)*c,:); C=Yii*pinv(xi);
%模态识别
[v,z]=eig(A); dis=C*v;
Ph=180*(abs(angle(dis))-pi/2)/pi; %计算相位角
Mdis=abs(dis).*sign(Ph); %应变振型矩阵
z1=diag(z); z2=log(z1);
fr=abs(z2)*d; %系统圆频率
damp=real(z2)./fr*d; %系统阻尼比
fr=fr/(2*pi); %频率
figure(1); plot(ss,'*'); xlabel('阶数');ylabel('奇异值')
[ 本帖最后由 ChaChing 于 2010-5-28 23:41 编辑 ] |
评分
-
2
查看全部评分
-
|