基于随机子空间方法模态参数识别的相关matlab程序（可用哦，我在

Frree · 发表于 2014-12-10 15:59

我也在做这方面谢谢啦

xieniuniu · 发表于 2015-5-20 11:13

确实有些乱，需要整理一下！

wocai · 发表于 2015-12-26 20:34

谢谢楼主！！！！

skyliu · 发表于 2017-3-20 13:58

谢谢楼主！！！！

Agoni · 发表于 2017-3-21 09:00

%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*'); %S1Y用黑色的星号表示
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); %对系统矩阵A进行特征值分解
%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; %阵型矩阵=C*（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))); %abs求绝对值
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]) %y轴的范围
%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); %将Freq2按行升序排列，IX是排序的位置信息
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)); %樊江玲文章中的式（4-48）
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); %find就是找到符号某一条件的数的下标
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 Dpro Dpro2 Dpro3 MAC_shape MODE MODE2...
% M4 IX M2 M3 indicefreq m2 dualfreq fs mm2 n3 m3 m4...
% crtimodeshape MODE3 crtidamping crtifreq X Y
%%%%%%%%%%%%%%%%%%%%%%<函数>模态参数识别绘图函数%%%%%%%%%%%%%
N_MODE=6 %%%%%%%%%绘制一阶振型
MODE_N=MODE4(1:8,N_MODE);
ylabel('振型');
xlabel('节点号');
plot(MODE_N,'b-'); %蓝色实线
ylim([-4 4])
%%%%%%%%%%%%%%%%%%%%%%%%稳定图%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;clear
load b16b1_Waveform.txt
TimeVy=b16b1_Waveform(1:5120,1:9);
%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','2'}; %系统阶数N=2*n
answer=inputdlg(prompt,dlg_title,num_lines,def);
M=str2num(answer{1,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);
freqency_of_Sample=1/(TimeVy(2,1)-TimeVy(1,1));
fs=freqency_of_Sample;
T=1/freqency_of_Sample;
p=1;
%%%%%%%%%%%%%%
N_min=2; %
N_max=40; %
%%%%%%%%%%%%%%
S_freq=-ones((N_max-N_min)/2+1,N_max);
S_damp=-ones((N_max-N_min)/2+1,N_max);
for N_system=N_min:2:N_max
U_N=U(:,1:N_system);
V_N=V(:,1:N_system);
S_N=S(1:N_system,1:N_system);
A_N=S_N^-0.5*(U_N')*Teop2*V_N*S_N^-0.5;
[V_N,D_N]=eig(A_N); % 求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量
[i,j]=size(D_N);
for m=1:i
Eignvalue(1,m)=D_N(m,m); %转化为Eignvalue的行向量
end
Eignvalue_Ac=log(Eignvalue)/T
Freq_N=abs(Eignvalue_Ac)/2/pi;
Damp_N=-real(Eignvalue_Ac)./abs(Eignvalue_Ac);
for q=1:N_system
S_freq(p,q)=Freq_N(1,q);
S_damp(p,q)=Damp_N(1,q);
end
p=p+1
end
%for a=1:N_max
%for b=1:(N_max-N_min)/2
% if abs((S_freq((N_max-N_min)/2,a)-S_freq(b,a))/S_freq((N_max-N_min)/2,a))>=0.05|abs((S_damp((N_max-N_min)/2,a)-S_damp(b,a))/S_damp((N_max-N_min)/2,a))>=0.05 %%%%%%%%%
% S_freq(b,a)=-1
% end
% end
%end
%%%%%%%%%绘制%%%%%%%%%%
Ref=8; %%%%%%%%%%%%%%%%%%%%%%%%%%
y_axis=N_min:2:N_max;
nfft=TimePointNum;
%[Pxy,f]=csd(TimeVy(:,Ref),TimeVy(:,9),nfft,fs,hanning(nfft),nfft/2,'mean');
[Pxx,f]=pwelch(TimeVy(:,Ref),hanning(nfft),nfft/2,nfft,fs); %谱分析函数，pxx和f分别为功率谱密度和频率
[AX,handel_fre,handle_spec]=plotyy(S_freq,y_axis,f,Pxx);
set(get(AX(1),'Ylabel'),'String','Number of Singular Values(N)','fontweight','bold')
set(get(AX(2),'Ylabel'),'String','Power Spectral Density(dB)','fontweight','bold')
set(get(AX(1),'Xlabel'),'String','Frequency,Hz','fontweight','bold')
title('Stability diagram','fontweight','bold')
set(AX(1),'Xlim',[0 fs/2]) %xlim表示x轴的限值
set(AX(2),'Xlim',[0 fs/2])
set(handel_fre,'LineStyle','*');
set(handel_fre,'color','r');
set(AX(2),'Ylim',[0 0.0001])

复制代码

账号		自动登录	找回密码
密码			我要加入

[随机振动] 基于随机子空间方法模态参数识别的相关matlab程序（可用哦，我在