arma with next.m

自然激励下建筑结构的模态参数识别

clear all
clc
close all hidden
format long

%%%NExT（自然激励技术）
%采样频率
sf=200;
%输出数据长度
np=512;
load E:\Dome\地脉动（4测点随机减量和ITD）分\wh\wh.txt
%参照点数据
x1=wh(:,3);
x=x1';
%响应点数据
y1=wh(:,1);
y=y1';
t=0:1/sf:(np-1)/sf;
nfft=2^nextpow2(2*np);
p=csd(x,y,nfft);
p(nfft/2+1)=real(p(nfft/2));
p(nfft/2+2:nfft)=conj(p(nfft/2:-1:2));
g=ifft(p);
r=real(g(1:np));
figure(1)
plot(t,r);
xlabel('时间(s)');
ylabel('幅值');
grid on;

%%%ARMA时间序列分析（自回归滑动平均模型）
%模态阶数
mn=3;
%互相关函数
b=r;
%ARMA模型阶数
nm=2*mn;
n=fix(length(b)/2);
h=b(1:2*n,1);
dt=1/sf;
t1=0:dt:(2*n-1)*dt;
%时间序列响应拟合的ARMA参数建模
%A和B分别为ARMA模型传递函数的分子和分母系数向量
[A B]=prony(h,nm,nm);
%多项式求根
V=roots(B);
%计算自振频率
F1=abs(log(V))/(2*pi*dt);
%计算阻尼比
D1=sqrt(1./(((imag(log(V))./real(log(V))).^2)+1));
%计算振型向量（特征向量）
for k=0:(2*n-1)
  Va(k+1,:)=[conj(V').^k];
end
S1=inv(conj(Va')*Va)*conj(Va')*(h);
%计算拟合的脉冲相应函数
h1=real(Va*S1);
%设置时域曲线长度
nn=1:length(t1);
%绘制脉冲响应函数拟合曲线图
figure(2)
plot(t1(nn),h(nn),':',t1(nn),h1(nn));
xlabel('时间(s)');
ylabel('幅值');
legend('实测','拟合');
grid on;
%将自振频率从小到大排序
[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(l);  %振型系数
end
wz=zeros(length(F),3);
wz(:,1)=F';
wz(:,2)=D';
wz(:,3)=S';