home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / FILTFILT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.9 KB  |  64 lines
  1. function y = filtfilt(b,a,x)
  2. %FILTFILT Zero-phase forward and reverse digital filtering.
  3. %    Y = FILTFILT(B, A, X) filters the data in vector X with the
  4. %    filter described by vectors A and B to create the filtered
  5. %    data Y.  The filter is described by the difference equation:
  6. %
  7. %      y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
  8. %                       - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
  9. %
  10. %    After filtering in the forward direction, the filtered
  11. %    sequence is then reversed and run back through the filter.
  12. %    The resulting sequence has precisely zero-phase distortion
  13. %    and double the filter order.  Care is taken to minimize
  14. %    startup and ending transients by matching initial conditions.
  15. %    See also FILTER.
  16.  
  17. %    L. Shure 5-17-88
  18. %    (c) Copyright 1988, by The MathWorks, Inc.
  19.  
  20. [m,n] = size(x);
  21. b = b(:).';
  22. a = a(:).';
  23. nb = length(b);
  24. na = length(a);
  25. nfilt = max(nb,na);
  26. % set up initial conditions to remove dc offset problems at the beginning of
  27. % the sequence
  28. if nb == na
  29.     zi = cumsum(b(nfilt:-1:1)-a(nfilt:-1:1)); zi = zi(nfilt-1:-1:1);
  30. else
  31.     if na < nfilt
  32.         bb = b;
  33.         aa = a;
  34.         aa(nfilt) = 0;
  35.     else    % nb < nfilt
  36.         aa = a;
  37.         bb = b;
  38.         bb(nfilt) = 0;
  39.     end
  40.     zi = cumsum(bb(nfilt:-1:1)-aa(nfilt:-1:1)); zi = zi(nfilt-1:-1:1);
  41.     clear bb aa 
  42. end
  43. % reflect beginning sequence of data so slopes match and filter this
  44. % precursory signal so end effects are reduced
  45. nfact = 3*(nfilt-1);
  46. if m == 1    % row data
  47.     y = [2*x(1)-x((nfact+1):-1:2) x];
  48. else     % column data
  49.     y = [2*x(1)-x((nfact+1):-1:2);x];
  50. end
  51. y = filter(b,a,y,[zi*y(1)]);
  52. % reverse data, prepend and filter again
  53. y(:) = y(length(y):-1:1);
  54. if m == 1    % row data
  55.     y = [2*y(1)-y((nfact+1):-1:2) y];
  56. else     % column temp
  57.     y = [2*y(1)-y((nfact+1):-1:2);y];
  58. end
  59. y = filter(b,a,y,[zi*y(1)]);
  60. % remove pieces of y that were tacked on to eliminate end effects
  61. y([1:nfact max(m,n)+nfact+(1:nfact)]) = [];
  62. % reverse the final sequence to be in original direction
  63. y = y(length(y):-1:1);
  64.