%prony法模態參數識別
%%%%%%%%%%%%%%%%%%%%%
%clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('prony模式識別數據文件名:','s');
%fni=out2.signals.values, 'DisplayName', 'out2.signals.values', 'YDataSource', 'out2.signals.values'
fid=fopen(fni,'r');
%mn=9
%fid=out2.signals.values(1:386,1);
mn=fscanf(fid,'%d',1); %模態階數
%[A,mn]=fscanf(fid,'%d',1);
%定義輸出實測數據類型；
%ig=1時域數據如沖擊響應，自由振動，互相關函數，隨即減量法處理結果；
%ig=2頻域數據如頻響函數實部虛部數據；
ig=fscanf(fid,'%d',1);
%ig=1時f為采樣頻率sf，ig=2時，f為頻率間隔df
f=fscanf(fid,'%f',1);
fno=fscanf(fid,'%s',1);%輸出數據文件名
b=fscanf(fid,'%f',[ig,inf]);%實測時域或者頻域數據
status=fclose(fid);
%建立復指數函數的階數
nm=2*mn;
%組織識別計算所用的時域數據及參數
if ig==1%實測時域數據
    %取采樣頻率 
    sf=f;
    %取時域長度1/2的長度 
    n=fix(length(b)/2);           
    %將輸入時域數據賦值給列向量h
    h=b(1,1:2*n)';
   % h=b(1:2*n,1);%%%%%%%
    %計算時間間隔
    dt=1/sf;
    %建立離散時間向量
    t=0:dt:(2*n-1)*dt;
else %實測頻域數據
    %取頻率間隔
    df=f;
    %取實測頻響函數的長度
    n=length(b(1,:));
    f=0:df:(n-1)*df;
    %建立對應正負頻率的實測頻響函數向量
    H=b(1,:)'+b(2,:)'*i;
    H(n+1)=real(H(n));
    H(n+2:2*n)=conj(H(n:-1:2));
    %頻響函數經IFFT并取實部變換成脈沖響應函數
    h=real(ifft(H));
    t=linspace(0,1/df,2*n);
    dt=t(2)-t(1);
end
L=length(h);
M=L/2;
for k=1:nm
   X1(:,k)=h(k:M-1+k);
end;
for k=1:M
  X2(k,:)=-h(nm+k);
end
%用最小二乘法求解待定系數的列向量
%B=X1\X2;
B=inv(X1'*X1)*(X1'*X2);
B(nm+1)=1;
%將B轉換成降冪排列的列向量
B1=B(nm+1:-1:1);
%冪多項式求根
V=roots(B1);%求出特征方程的根
%計算模態頻率向量
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);
%繪制脈沖響應函數擬合曲線圖
figure(1)
plot(t,h,'r:',t,h1);
xlabel('時間’(s)');
ylabel('幅值');
legend('實測','擬合');
grid on ;
%for z=1:1:2*n    
%a=((h(z)-h1(z))/h(z))/100;
%plot(a)
%end
if ig>1
    %計算生成的頻響函數
    H1=fft(Va*S1);
    %繪制頻響函數實部擬合曲線圖
    figure(2)
    nn=1:n;
    subplot(2,1,1);
 %  plot(f,b(1,:),':',f,real(H1(nn)),'r');    
plot(f,real(H(nn)),':',f,real(H1(nn)),'r');      
    xlabel('頻率’(Hz)');
ylabel('實部');
legend('實測','擬合');
grid on;
%繪制頻響函數虛部擬合曲線圖
    subplot(2,1,2);
    plot(f,b(2,:),':',f,imag(H1(nn)),'r');
    xlabel('頻率’(Hz)');
ylabel('虛部');
legend('實測','擬合');
grid on;
end
%將自振頻率從小到大排序
[F2,I]=sort(F1);
%剔除方程解中的非模態項（非共軛根）和共軛相
m=0;
for k=1:nm-1
    if F2(k)~=F2(k+1) %消除其中的共軛相
continue;
    end
m=m+1;
II=I(k);
F(m)=F1(II);%自振頻率
D(m)=D1(II);%阻尼比
S(m)=S1(II);%振型系數
A(m)=abs(S1(II));%振幅
theta(m)=angle(S1(II));%相角
%a=((h-h1)/h)/100
end
%打開文件輸出模態參數數據
fid=fopen(fno,'w');
fprintf(fid,'  頻率(Hz) 阻尼比（%%）振型系數 振幅 相角\n');
for k=1:m
    fprintf(fid,'%10.4f %10.4f %10.6f %10.4f %10.4f,\r\n',F(k),D(k)*100.0,imag(S(k)),A(k),theta(k));
end
status=fclose(fid)