ARMA模型预测

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 （