clear all
clc
close all hidden
format long

%%%NExT（自然激勵技術）
%采樣頻率
sf=200;
%輸出數(shù)據(jù)長度
np=512;
load E:\Dome\地脈動（4測點隨機減量和ITD）分\wh\wh.txt
%參照點數(shù)據(jù)
x1=wh(:,3);
x=x1';
%響應點數(shù)據(jù)
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時間序列分析（自回歸滑動平均模型）
%模態(tài)階數(shù)
mn=3;
%互相關函數(shù)
b=r;
%ARMA模型階數(shù)
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參數(shù)建模
%A和B分別為ARMA模型傳遞函數(shù)的分子和分母系數(shù)向量
[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);
%計算擬合的脈沖相應函數(shù)
h1=real(Va*S1);
%設置時域曲線長度
nn=1:length(t1);
%繪制脈沖響應函數(shù)擬合曲線圖
figure(2)
plot(t1(nn),h(nn),':',t1(nn),h1(nn));
xlabel('時間(s)');
ylabel('幅值');
legend('實測','擬合');
grid on;
%將自振頻率從小到大排序
[F2,I]=sort(F1);
%剔除方程解中的非模態(tài)項（非共軛根）和共軛項
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);  %振型系數(shù)
end
wz=zeros(length(F),3);
wz(:,1)=F';
wz(:,2)=D';
wz(:,3)=S';