home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / YULEWALK.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.2 KB  |  134 lines
  1. function [B,A] = yulewalk(na, ff, aa, npt, lap)
  2. %YULEWALK Recursive filter design using a least-squares method.
  3. %    [B,A] = YULEWALK(N,F,M) finds the N-th order recursive filter
  4. %    coefficients B and A such that the filter:
  5. %                              -1             -(n-1) 
  6. %           B(z)   b(1) + b(2)z + .... + b(n)z
  7. %           ---- = ---------------------------
  8. %                              -1             -(n-1)
  9. %           A(z)    1   + a(1)z + .... + a(n)z
  10. %
  11. %    matches the magnitude frequency response given by vectors F and M.
  12. %    Vectors F and M specify the frequency and magnitude breakpoints for
  13. %    the filter such that PLOT(F,M) would show a plot of the desired
  14. %    frequency response. The frequencies in F must be between 0.0 and 1.0,
  15. %    with 1.0 corresponding to half the sample rate. They must be in
  16. %    increasing order and start with 0.0 and end with 1.0. 
  17. %
  18. %    See also FIR1, BUTTER, CHEBY, FREQZ and FILTER.
  19.  
  20.  
  21. %   The YULEWALK function performs a least squares fit in the time domain.   
  22. %   The denominator coefficients {a(1),...,a(NA)} are computed by the
  23. %   so called "modified Yule Walker" equations, using NR correlation
  24. %   coefficients computed by inverse Fourier transformation of the specified
  25. %   frequency response H.
  26. %   The numerator is computed by a four step procedure. First, a numerator
  27. %   polynomial corresponding to an additive decomposition of the power 
  28. %   frequency response is computed. Next, the complete frequency response
  29. %   corresponding to the numerator and denominator polynomials is evaluated.
  30. %   Then a spectral factorization technique is used to obtain the 
  31. %   impulse response of the filter. Finally, the numerator polynomial
  32. %   is obtained by a least squares fit to this impulse response.
  33. %   For a more detailed explanation of the algorithm see 
  34. %   B. Friedlander and B. Porat, "The Modified Yule-Walker Method
  35. %   of ARMA Spectral Estimation," IEEE Transactions on Aerospace
  36. %   Electronic Systems, Vol. AES-20, No. 2, pp. 158-173, March 1984.
  37. %
  38. %   Perhaps the best way to get familiar with the proper use of YULEWALK
  39. %   is to examine carefully the demonstration file FILTDEMO.M, which
  40. %   provides a detailed example.
  41.  
  42. %    B. Friedlander 7-16-85
  43. %    Copyright (c) 1985, 1986 by the MathWorks, Inc.
  44. %    Modified 3-5-87 LS
  45. %    Repaired 7-26-90 Charles R. Denham.
  46.  
  47. if (nargin < 3 | nargin > 5)
  48.    error('Wrong number of input parameters.')
  49. end
  50. if (nargin > 3)
  51.    if round(2 .^ round(log(npt)/log(2))) ~= npt
  52.     % NPT is not an even power of two
  53.       npt = round(2.^ceil(log(npt)/log(2)));
  54.    end 
  55. end
  56. if (nargin < 4)
  57.    npt = 512;
  58. end
  59. if (nargin < 5)
  60.    lap = fix(npt/25);
  61. end
  62.    
  63. [mf,nf] = size(ff);
  64. [mm,nn] = size(aa);
  65. if mm ~= mf | nn ~= nf
  66.    error('You must specify the same number of frequencies and amplitudes.')
  67. end
  68. nbrk = max(mf,nf);
  69. if mf < nf
  70.    ff = ff';
  71.    aa = aa';
  72. end
  73.  
  74. if abs(ff(1)) > eps | abs(ff(nbrk) - 1) > eps
  75.    error('The first frequency must be 0 and the last 1.')
  76. end
  77.  
  78. % interpolate breakpoints onto large grid
  79.  
  80. npt = npt + 1;  % For [dc 1 2 ... nyquist].
  81. Ht = zeros(1,npt);
  82.  
  83. nint=nbrk-1;
  84. df = diff(ff');
  85. if (any(df < 0))
  86.    error('Frequencies must be non-decreasing.')
  87. end
  88.  
  89. nb = 1;
  90. Ht(1)=aa(1);
  91. for i=1:nint
  92.     if df(i) == 0
  93.        nb = nb - lap/2;
  94.        ne = nb + lap;
  95.     else
  96.        ne = fix(ff(i+1)*npt);
  97.     end
  98.     if (nb < 0 | ne > npt)
  99.        error('Too abrupt an amplitude change near end of frequency interval.')
  100.     end
  101.     j=nb:ne;
  102.     if ne == nb
  103.         inc = 0;
  104.     else
  105.         inc = (j-nb)/(ne-nb);
  106.     end
  107.     Ht(nb:ne) = inc*aa(i+1) + (1 - inc)*aa(i);
  108.     nb = ne + 1;
  109. end
  110. Ht = [Ht Ht(npt-1:-1:2)];
  111. n = max(size(Ht));
  112. n2 = fix((n+1)/2);
  113. nb = na;
  114. nr = 4*na;
  115. nt = 0:1:nr-1;
  116.  
  117. % compute correlation function of magnitude squared response
  118.  
  119. R = real(ifft(Ht .* Ht));
  120.  
  121. R  = R(1:nr).*(0.54+0.46*cos(pi*nt/(nr-1)));     % pick NR correlations 
  122.  
  123. % Form window to be used in extracting the right "wing" of two-sided
  124. % covariance sequence.
  125. Rwindow = [1/2 ones(1,n2-1) zeros(1,n-n2)]; 
  126.  
  127. A = polystab(denf(R,na));                % compute denominator
  128.  
  129. Qh = numf([R(1)/2,R(2:nr)],A,na);    % compute additive decomposition
  130.  
  131. Ss = 2*real(freqz(Qh,A,n,'whole'))';    % compute impulse response
  132. hh = ifft(exp(fft(Rwindow.*ifft(log(Ss)))));
  133. B  = real(numf(hh(1:nr),A,nb));
  134.