ex09061.m

ARMA模型时间序列分析法简称为时序分析法

%
% ex09061
clear all; clc;
mn=7;
ig=2;
f=0.1953125;
b=load('x09051.txt');       %读入数据
b=b';
fno='out09061.txt';      
%b=fscanf (fid,'%f', Eig, inf]); 
%status=fclose (fid);
rad=pi/180;

%~ ARMA~IR (~:~1~ 2~)
nm=2 * mn;
if ig==1    %~:~tI~
  sf= f;
  n=fix(length(b) /2);
  h=b(1, 1: 2*n)';
  dt=1/sf;
  t=0 : dt: (2*n-1) *dt;
else     %~tRt~
  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));
  h=real(ifft(H));
  t=linspace(0, 1/df, 2*n);
  dt=t(2) -t(1);
end
%~l'l~Jg~ln~a~a ARMA #~j~
[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);
h1=real(Va*S1);
nn=1: length(t);
figure (1)
plot(t(nn), h(nn),'r','linewidth',2); hold on;
plot(t(nn),h1(nn),'k'); hold off;
xlabel('Time (s)');
ylabel('Amplitude');
legend ('Measured' ,'Fitness');
grid on;
if ig>1
  H1=fft(Va* S1);
  nn=1:  n;
  figure(2)
  subplot(2, 1, 1)
   plot(f, real(H(nn)),'r', 'linewidth',2); hold on; 
   plot(f, real(H1(nn)),'k'); hold off;
  xlabel('Frequency (Hz)');
  ylabel('Real Part');
  legend('Measured','Fitness');
  grid on;
  subplot (2, 1, 2);
  plot(f, b(2,:),'r', 'linewidth',2); hold on;
  plot(f, imag(H1(nn)),'k'); hold off;
  xlabel('Frequency (Hz)');
  ylabel('Imag Part');
  legend('Measured','Fitness');
  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);  %[~g
  D(m) =D1(l);   %l~l]~l~
  S(m) =S1(l);  %M~tR
end
fid=fopen(fno,'w');
fprintf (fid,'   Freq(Hz)   DampRate (MM)   Coef    Am    Angle\n');
for k=1: m
  fprintf (fid,'%10.4f  %10.4f  %10.6f  %10.4f  %10.4f\n', F(k), D(k)*100.0, imag(S(k)), 2*abs(S(k)), angle(S(k))/rad);
end
status=fclose (fid);