function [a,gamma]=mywarma(y,n,m,M)%% The modified Yule-Walker ARMA method % given by equation (3.7.9) with the AR coefficients estimated% using the overdetermined set of equation (3.7.4), where W is% the identity matrix.%% [a,gamma]=mywarma(y,n,m,M)% %      y     -> the data vector%      n     -> AR model order%      m     -> MA model order%      M     -> the constant which determine the amount of %               overdetermination%      a     <- the AR coefficients vector%      gamma <- the onesided gamma in equation (3.7.8)% Copyright 1996 by R. Mosesy=y(:);N=length(y);             % data length%-------------------- estimate the AR coefficients --------------------% compute the standard biased ACS estimate [r(0) r(1) r(2) ...r(n)]if (N <= m+M)   disp('Error: m+M must be < the data length.');   returnelseif (M<n)   disp('Error: M must be >= n');   returnendr=zeros(m+M+1,1);for i = 0 : m+M,   r(i+1)=y(1:N-i)'*y(i+1:N)/N;end% form the r vector and R matrix in equation (3.7.1)r1=r(m+2:m+M+1);if ((m-n+1)>=0)   R1=toeplitz(r(m+1:m+M),r(m+1:-1:m-n+2).');else   R1=toeplitz(r(m+1:m+M),[r(m+1:-1:1);conj(r(2:abs(m-n+1)+1))].');end% compute the AR coffecientsa=-R1\r1;a=[1;a];%-------------------- estimate the gamma coefficients --------------------gamma=zeros(m+1,1);for k = 0 : m,   for j = 0 : n,      for p = 0 : n,         ind=k+p-j;         if ind >=0            gamma(k+1)=gamma(k+1)+a(j+1)*a(p+1)'*r(ind+1);         else            gamma(k+1)=gamma(k+1)+a(j+1)*a(p+1)'*r(abs(ind)+1)';         end      end   endend