itd_command.m

随机减量方法

clear;
clc;
clear;
clc;
load data.txt;
mean_data=mean(data);
data=data-mean_data;
for i=1:6
    temp(:,i)=data((i-1)*10*1024+1:i*10*1024);
end
x=temp(:,1);
sf=20;                 %采样频率
np=601;               %输出数据长度
t=0:1/sf:(np-1)/sf;     %建立离散输出时间向量
nt=length(x);
s=1.2*std(x);
m=0;
y=zeros(1,np);
y=y';
for k=2:nt-np
    a=abs(x(k-1)-s);
    b=abs(x(k)-s);
    c=abs(x(k+1)-s);
    if b<a&b<c
        y(1:np)=y(1:np)+x(k:k+np-1);
        m=m+1;
    end
end
y=y./m;
plot(t,y);
mn=30;%模态阶数
nm=2*mn;%建立特征方程矩阵的阶数（为模态阶次的2倍）
n=fix(length(y)/2);
h=y;
dt=1/sf;
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