main.m

随机减量方法

clear;
clc;
load data.txt;%%%采用1号传感器的数据
mean_data=mean(data);
data=data-mean_data;
for i=1:6
    temp(:,i)=data((i-1)*10240+1:i*10240);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num=10240;
noverlap=2048;
fss=20;nfft=10240;window=hanning(num);
[px_qian,f]=pwelch(temp(:,1),window,noverlap,nfft,fss);
wp=0.7/(20/2);
[b,a]=butter(4,wp,'low');%巴特沃斯低通滤波器
x=filter(b,a,temp(:,1));
[px_hou,f]=pwelch(x,window,noverlap,nfft,fss);
%%%%%%%%%%计算自由振动信号
[t,y]=rdt(x);%%分别采用样本点数为900、1024和1200三种不同的计算条件，最终选择待确定

mn=30;%模态阶数
nm=2*mn;%建立特征方程矩阵的阶数（为模态阶次的2倍）
n=fix(length(y)/2);
h=y;
dt=1/20;
t=0:dt:(2*n-1)*dt;
%计算自由振动响应矩阵
L=length(h);
M=L/2;
for k=1:nm
    x1(k,:)=h(k:L-(nm-k+1));
    x2(k,:)=h(k+1:L-(nm-k));
end

B=x1\x2;%B=x2*x1'*inv(x1*x1');%用最小二乘法求解特征方程矩阵
[A,V]=eig(B);%计算特征值和特征向量
%变换特征值对角阵为以向量
for k=1:nm
    U(k)=V(k,k);
end
F1=abs(log(U'))./(2*pi*dt);%计算模态频率向量
D1=sqrt(1./(((imag(log(U'))./real(log(U'))).^2)+1));%计算阻尼比向量
%计算振型系数向量
l=1;
for k=0:(2*n-1)
    Va(k+1,:)=[conj(U).^k];
end
%S1=(inv(conj(Va')*Va)*conj(Va')*h);
[F2,I]=sort(F1);%将自振频率从小到大排序
%%提出方程解中的非模态项（非共轭根）和共轭根（重复项）
m=0;
for k=1:nm-1
    if F2(k)~=F2(k+1)
        continue;
    end
    m=m+1;
    l=I(k);
    F(m)=F1(l);
    D(m)=D1(l);
    %S(m)=S1(1);
end