clear
clc
format?long
fni=?input（‘ARMA輸入文件名‘）;
fid=?fopen（fni‘r‘）;
mn=?fscanf（fid‘%d‘1）;?%模態階數

%定義輸入實測數據類型
ig=?fscanf（fid‘%d‘1）;%ig=1時域數據如沖擊響應、自耦振動、互相關函數、隨機減量法處理結果
???????????????????????%ig=2頻域數據如頻響函數實部虛部數據
f=?fscanf（fid‘%f‘1）;?%ig=1時f為采樣頻率；ig=2時f為采樣間df

fno=?fscanf（fid‘%s‘1）;?%輸出數據文件名
b=?fscanf（fid‘%f‘[iginf]）;%實測時域或頻域數據
status=?fclose（fid）;

%建立ARMA模型的階數（為模態階數的2倍）
nm=?2*mn;
%組織識別計算所用的時域數據及參數
if?ig==1?%實測時域數據
????%取采樣頻率
????sf=?f;
????%取時域數據長度的1/2
????n=?fix（length（b）/2）;
????%將輸入時域數據賦值給列向量h
????h=?b（1?1:2*n）‘;
????%計算時間間隔
????df=?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（01/df2*n）;
????%計算時間間隔
????dt=?t（2）?-t（1）;
end

%時間序列響應擬合的ARMA參數建模
%A和B分別為ARMA模型傳遞函數的分子和分母系數向量
[A?B]?=prony（hnmnm）;
%多項式求根（零點）
V=?roots（B）;
%計算自振頻率
F1=?abs（log（V））?/（2*pi*dt）;
%計算阻尼比
D1=?sqrt（1./?（?（?（imag?（log?（V））./real?（