arma.m

来自「内燃机转子仿真」· M 代码 · 共 83 行

clear all;
clc;
close all;
format long;
%function [wf1 wd1 wv1]=arma(v)
mn=12;
ig=1;
load 'outsjjl';%ig=1，时域振动数据
b=y;
f=2000;%ig=1，f为采样频率，ig=2时为频率间隔
nm=2*mn;
if ig==1
    sf=f;
    n=fix(length(b)/2);%n=500
    h=b(1,1:2*n)';
    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,:)';
    H(n+1)=real(H(n));
    H(n+2:2*n)=conj(H(n:-1:2));
    h=real(ifft(H));
    t=linspace(0,1/df,2*n);
    dt=t(2)-t(1);
end

[A B]=prony(h,mn,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(t);
figure(1);
plot(t(nn),h(nn),':',t(nn),h1(nn));
xlabel('时间(s)');
ylabel('幅值');
legend('实测','拟合');
grid on;
if ig>1
    H1=fft(Va*Sl);
    nn=1:n;
    figure(2);
    subplot(2,1,1);
    plot(f,real(H(nn)),':',f,real(H1(nn)));
    xlabel('频率(Hz)');
    ylabel('实部');
    legend('实测','拟合');
    grid on；
    subplot(2,1,2);
    plot(f(nn),b(2:nn),':',f(nn),imag(H1(nn)));
    xlabel('频率(Hz)');
    ylabel('虚部');
    lelend('实测','拟合');
    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;
    l=I(k);
    F(m)=F1(l);
    D(m)=D1(l);
    S(m)=S1(l);
end