e1081.m

matlab算法集 matlab算法集

%----------------------------------------------------------------
% Example 10.8.1: ARMA model 
%----------------------------------------------------------------

% Initialize

   clear
   clc
   p = 256;
   n1 = 3;
   m1 = 3;
   n2 = 3;
   m2 = 3;
   r1 = n1 + m1 + 1;
   r2 = n2 + m2 + 1;
   fs = 1000;
   theta1 = [-.5 .6 -.3 1. -2. 8. -4.]';
   theta2 = zeros(r2,1);
   x1 = zeros (r1,1);
   x2 = zeros (r2,1);
   A = zeros (p,2);
   phi = zeros (p,2); 
   Y = zeros (p,2);
   t = zeros (p,1);

% Identify system  

   fprintf ('Example 10.8.1: ARMA Model\n'); 
   u = randu (p,1,-1,1);
   t = [0 : p-1]';
   for i = 1 : p
      [x1,Y(i,1)] = arma (u(i),theta1,x1,n1,m1);
   end
   theta2 = getarma (u,Y(:,1),n2,m2);
   show ('theta',theta2)
      
% Compare responses 
   
   disp('Computing frequency responses ... ')
   x1 = zeros (r1,1);
   x2 = zeros (r2,1);
   for i = 1 : p
      if i == 1
         u(i) = 1;
      else
         u(i) = 0;
      end
      [x1,Y(i,1)] = arma (u(i),theta1,x1,n1,m1);
      [x2,Y(i,2)] = arma (u(i),theta2,x2,n2,m2);
   end

% Display results

   for i = 1 : 2
      [A(:,i),phi(:,i),f] = freqrsp (u,Y(:,i),fs);
   end
   graphxy (f(1:p/2),A(1:p/2,:),'Magnitude Responses','f (Hz)','A')
   graphxy (f(1:p/2),phi(1:p/2,:),'Phase Responses',...
            'f (Hz)','phi (deg)')
   j = min([40,p]);
   graphxy (t(1:j),Y(1:j,:),'Time Responses','k','y(k)')
%----------------------------------------------------------------