%1、噪声添加程序 y向加速度
%读取数据文件
clc;clear
[FileName,PathName]=uigetfile('*.txt','Select the TimeVy File'); %显示选择时程文件的对话框
File=strcat(PathName,FileName); %把数横向连接成单个字符串
fip=fopen(File,'r');
TimeVy=load('-ascii',File); %载入时程信息
[TimePointNum,MesNodeTotalNum]=size(TimeVy); %这个矩阵与hankel块方向不同,下面为转置
MesNodeTotalNum=MesNodeTotalNum-1; %第一列为时间序列
TimePointNum=TimePointNum;
Baizaosheng=wgn(TimePointNum,MesNodeTotalNum,1);
%产生TimePointNum行1列的随机数
Zero1=zeros(TimePointNum,1);
Kuozhan=[Zero1,Baizaosheng];
B=0.05*Kuozhan; %加入5%的噪声
C=TimeVy+B;
save TimeVy2.mat C -ascii; %保存白噪声
%2、奇异值分解
%读取数据文件
clc;clear
[FileName,PathName]=uigetfile('*.txt','Select the TimeVy File');
File=strcat(PathName,FileName);
fip=fopen(File,'r');
TimeVy=load('-ascii',File);
[TimePointNum,MesNodeTotalNum]=size(TimeVy);
MesNodeTotalNum=MesNodeTotalNum-1;
TimePointNum=TimePointNum;
%输入计算参数
%创建输入对话框
prompt={'TOEPLITZ矩阵块数M(>=采样频率/(2*结构基频)):'}; %j>20M
dlg_title='Input Parameter of Calculate';
num_lines=1;
def={'100'};
answer=inputdlg(prompt,dlg_title,num_lines,def);
%获取输入的值
M=str2num(answer{1,1});
j=TimePointNum-2*M; %Hankel矩阵最后一个矩阵块对应坐标为(2*i,j),即2*i+j需<=TimePointNum
Hankel=zeros((2*M+1)*MesNodeTotalNum,j);
for ti=1:j
m=0
for tii=1:(2*M+1)
km=ti+tii-1;
m=m+1;
%计算每次赋值在所在行的起始位置
beginy=(m-1)*MesNodeTotalNum+1;
%计算每次赋值在所在行的结束位置
endy=m*MesNodeTotalNum;
Hankel(beginy:endy,ti)=TimeVy(km,2:(MesNodeTotalNum+1))';
end
end
Yp=Hankel(1:M*MesNodeTotalNum,:);
Yf=Hankel(M*MesNodeTotalNum+1:2*M*MesNodeTotalNum,:);
Yf2=Hankel((M+1)*MesNodeTotalNum+1:(2*M+1)*MesNodeTotalNum,:);
Teop1=Yf*Yp'/j;
%组成第一个Toeplitz矩阵
[U,S,V]=svd(Teop1);
S1=diag(S);
%在窗口中绘制矩阵奇异值分解结果图
cla
newplot %为后续的图形命令创建图形窗口和坐标轴,并返回当前坐标轴的句柄
subplot(1,1,1) %作为一个二维曲线绘制函数,它的功能是:将一个窗口分为若干块,在选中的某一块内可以绘制图形
plot(S1,'k*');
ylabel('Singular value');
xlabel('Order of system');
title('Singular value decomposition results ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Baizaosheng=wgn(TimePointNum,MesNodeTotalNum,1);
%产生TimePointNum %行1列的随机数
% Zero1=zeros(TimePointNum,1);
% Kuozhan=[Zero1,Baizaosheng];
% B=0.15*Kuozhan; %加入5%的噪声
% TimeVy=TimeVy+B;
%3、协方差驱动SSI模态参数识别程序; %TimeVy需为偶数
clc;clear
%[FileName,PathName]=uigetfile('*.txt','Select the TimeVy File');
%File=strcat(PathName,FileName);
%fip=fopen(File,'r');
%TimeVy=load('-ascii',File);
load b16z2_Waveform.txt %%%%%%%%%%%%%%%%%%%%%%%%%%
TimeVy=b16z2_Waveform(1:5120,1:16); %%%%%%%%%%%%%%%%%%%%%%%%%%
%yacc=[];
[TimePointNum,MesNodeTotalNum]=size(TimeVy);
MesNodeTotalNum=MesNodeTotalNum-1;
TimePointNum=TimePointNum;
prompt={'TOEPLITZ矩阵块数M(>=采样频率/(2*结构基频)):','系统阶数N(<i*实际测点数):'};
dlg_title='Input Parameter of Calculate';
num_lines=1;
def={'50','30'}; %系统阶数N=2*n
answer=inputdlg(prompt,dlg_title,num_lines,def);
M=str2num(answer{1,1});
N=str2num(answer{2,1});
j=TimePointNum-2*M;
Hankel=zeros((2*M+1)*MesNodeTotalNum,j);
for ti=1:j
m=0
for tii=1:(2*M+1)
km=ti+tii-1;
m=m+1;
%计算每次赋值在所在行的起始位置
beginy=(m-1)*MesNodeTotalNum+1;
%计算每次赋值在所在行的结束位置
endy=m*MesNodeTotalNum;
Hankel(beginy:endy,ti)=(TimeVy(km,2:(MesNodeTotalNum+1)))';
end
end
Yp=Hankel(1:M*MesNodeTotalNum,:);
Yf=Hankel(M*MesNodeTotalNum+1:2*M*MesNodeTotalNum,:);
Yf2=Hankel((M+1)*MesNodeTotalNum+1:(2*M+1)*MesNodeTotalNum,:);
%组成第一个Toeplitz矩阵
Teop1=Yf*Yp'/j;
%组成第二个Toeplitz矩阵
Teop2=Yf2*Yp'/j;
%Hankel=[];
[U,S,V]=svd(Teop1);
%提取前N阶奇异值分解结果
U_de=U(:,1:N);
S_de=S(1:N,1:N);
V_de=V(:,1:N);
%计算C,A矩阵的估计
O_GJ=U_de*S_de^0.5;
C_GJ=O_GJ(1:MesNodeTotalNum,:);
A_GJ=S_de^-0.5*(U_de')*Teop2*V_de*S_de^-0.5;
%U_de=[];
%S_de=[];
%V_de=[];
%计算状态矩阵A的特征值和特征向量
[V_of_A,D_of_A]=eig(A_GJ);
%D_of_A=zeros(N,N);
%V_of_A=zeros(N,N);
%for i=1:N/2
%D_of_A(i,i)=D_of_A1(i*2-1,i*2-1);
%V_of_A(:,i)=V_of_A1(:,i*2-1);
%end
%for i=N/2+1:N
%D_of_A(i,i)=D_of_A1(i*2-N,i*2-N);
%V_of_A(:,i)=V_of_A1(:,i*2-N);
%end
%计算系统的特征值和特征向量
freqency_of_Sample=1/(TimeVy(2,1)-TimeVy(1,1));
fs=freqency_of_Sample;
T=1/freqency_of_Sample;
Total_mode_jieShu=N;
%计算振型
MODE=C_GJ*V_of_A;
%计算频率
for i=1:N
LN_lamuda(i,1)=log(D_of_A(i,i))/T;
Re_Abs(i,1)=abs(real(LN_lamuda(i,1)));
Im_Abs(i,1)=abs(imag(LN_lamuda(i,1)));
Freq(i,1)=abs(LN_lamuda(i,1))/2/pi;
%计算阻尼比
Dpro(i,1)=Re_Abs(i,1)/abs(LN_lamuda(i,1));
end
%将实部和虚步分别归一化
for opp_index=1:Total_mode_jieShu
min_of_real=min(real(MODE(:,opp_index)));
max_of_real=max(real(MODE(:,opp_index)));
min_of_imag=min(imag(MODE(:,opp_index)));
max_of_imag=max(imag(MODE(:,opp_index)));
if max_of_real~=min_of_real
Mode_value_Total_of_real_GYH=2*(real(MODE(:,opp_index))-min_of_real)/(max_of_real-min_of_real)-1;
else
Mode_value_Total_of_real_GYH=zeros(MesNodeTotalNum,1);
end
if max_of_imag~=min_of_imag
Mode_value_Total_of_imag_GYH=2*(imag(MODE(:,opp_index))-min_of_imag)/(max_of_imag-min_of_imag)-1;
else
Mode_value_Total_of_imag_GYH=zeros(MesNodeTotalNum,1);
end
MODE(:,opp_index)=Mode_value_Total_of_real_GYH+Mode_value_Total_of_imag_GYH*sqrt(-1);
if imag(MODE(:,opp_index))~=0;
MODE(:,opp_index)=imag(MODE(:,opp_index));
else
MODE(:,opp_index)=real(MODE(:,opp_index));
end
end
%N_MODE=2 %%%%%%%%%绘制一阶振型
%MODE_N=MODE(:,N_MODE);
%ylabel('振型');
%xlabel('节点号');
%plot(MODE_N,'k*');
%ylim([-4 4])
%clear prompt answer def dlg_title beginy acc6 freqency_of_Sample...
% endy V_de V U_de Yp Yf2 Yf num_lines min_of_real...
% min_of_imag tii ti opp_index km j i max_of_real...
% max_of_imag m U MesNodeTotalNum M Mode_value_Total_of_imag_GYH...
% N Mode_value_Total_of_real_GYH Im_Abs C_GJ A_GJ...
% Hankel TimeVy T Teop1 TimePointNum Teop2 Re_Abs...
% O_GJ S Total_mode_jieShu S_de D_of_A LN_lamuda V_of_A
%%%%%%%%%%%%%%%%%%%%%%真实模态判断%%%%%%%%%%%%%%%%%%%%%%%%
M2=length(Freq);
mm2=1;
for m2=1:M2
if (Freq(m2)<=(0.5*fs))&(Freq(m2)>0)&(Dpro(m2)>0)&(Dpro(m2)<0.1)
Freq2(mm2)=Freq(m2);
Dpro2(mm2)=Dpro(m2);
MODE2(:,mm2)=MODE(:,m2);
mm2=mm2+1;
end
end
[Freq3,IX]=sort(Freq2);
for m3=1:length(Freq3)
Dpro3(m3)=Dpro2(IX(m3));
MODE3(:,m3)=MODE2(:,IX(m3));
end
%%%%%%%%%重频合并%%%%%
%%重频判断系数
crtifreq=1;
crtidamping=5;
crtimodeshape=98; %相同模态判断
M3=length(Freq3);
dualfreq=ones(1,M3);
for m3=1:M3
for n3=m3+1:M3
X=MODE3(:,m3);
Y=MODE3(:,n3);
MAC_shape=(X'*Y)^2/((X'*X)*(Y'*Y));
if (((Freq3(m3)-Freq3(n3))*100/Freq3(m3))<crtifreq)...
&(((Dpro3(m3)-Dpro3(n3))*100/Dpro3(m3))<crtidamping)...
&(MAC_shape*100>crtimodeshape)
dualfreq(n3)=0;
end
end
end
indicefreq=find(dualfreq);
M4=length(indicefreq);
for m4=1:M4
Freq4(m4)=Freq3(indicefreq(m4));
Dpro4(m4)=Dpro3(indicefreq(m4));
MODE4(:,m4)=MODE3(:,indicefreq(m4));
end
%clear Freq2 Freq3 Freq Dp
没有合适的资源?快使用搜索试试~ 我知道了~
matlab (SSI)随机子空间法程序
共1个文件
txt:1个
5星 · 超过95%的资源 需积分: 44 165 下载量 196 浏览量
2017-12-31
09:44:32
上传
评论 17
收藏 4KB ZIP 举报
温馨提示
模态分析程序,提取结构模态自振频率,阻尼比等模态信息
资源推荐
资源详情
资源评论
收起资源包目录
SSI.zip (1个子文件)
SSI.txt 11KB
共 1 条
- 1
资源评论
- ll5704422023-03-07需要导入什么样的文件啊
qq_37452676
- 粉丝: 4
- 资源: 1
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功