home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / DECIMATE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.1 KB  |  112 lines
  1. function odata = decimate(idata,r,nfilt,option)
  2. %DECIMATE Resample data at a lower rate after lowpass filtering.
  3. %    Y = DECIMATE(X,R) resamples the sequence in vector X at 1/R
  4. %    times the original sample rate.  The resulting resampled
  5. %    vector Y is R times shorter, LENGTH(Y) = LENGTH(X)/R.
  6. %    DECIMATE filters the data with an eighth order Chebyshev
  7. %    type I lowpass filter with cutoff frequency .8*(Fs/2)/R,
  8. %    before resampling.
  9. %    Y = DECIMATE(X,R,N) uses an N'th order Chebyshev filter.
  10. %    Y = DECIMATE(X,R,'FIR') uses the 30 point FIR filter
  11. %    generated by FIR1(30,1/R) to filter the data.
  12. %    Y = DECIMATE(X,R,N,'FIR') uses the N-point FIR filter.
  13. %    See also INTERP.
  14.  
  15. %    L. Shure 6-4-87
  16. %    Revised 12-11-88 LS
  17. %    Copyright (c) 1987-88 by the MathWorks, Inc.
  18. %    All rights reserved.
  19.  
  20. %    References:
  21. %      "Programs for Digital Signal Processing", IEEE Press
  22. %      John Wiley & Sons, 1979, Chap. 8.3.
  23.  
  24. if nargin < 2
  25.    error('Not enough input arguments.')
  26. end
  27. if abs(r-fix(r)) > eps
  28.    error('Resampling rate R must be a positive integer.')
  29. end
  30. if fix(r) == 1
  31.     odata = idata;
  32.     return
  33. end
  34. if r <= 0
  35.     error('Resampling rate R must be a positive integer.')
  36. end
  37. fopt = 1;
  38. if nargin == 2
  39.     nfilt = 8;
  40. else
  41.     if isstr(nfilt)
  42.         if nfilt(1) == 'f' | nfilt(1) == 'F'
  43.             fopt = 0;
  44.         elseif nfilt(1) == 'i' | nfilt(1) == 'I'
  45.             fopt = 1;
  46.         else
  47.             error('Filter options are FIR and IIR.')
  48.         end
  49.         if nargin == 4
  50.             nfilt = option;
  51.         else
  52.             nfilt = 8*fopt + 30*(1-fopt);
  53.         end
  54.     else
  55.         if nargin == 4
  56.             if option(1) == 'f' | option(1) == 'F'
  57.                 fopt = 0;
  58.             elseif option(1) == 'i' | option(1) == 'I'
  59.                 fopt = 1;
  60.             else
  61.                 error('Filter options are FIR and IIR.')
  62.             end
  63.         end
  64.     end
  65. end
  66. if fopt == 1 & nfilt > 13
  67.     disp('Warning: IIR filters above order 13 may be unreliable.')
  68. end
  69.  
  70. nd = max(size(idata));
  71. nout = ceil(nd/r);
  72. [m,n]=size(idata);
  73.  
  74. if fopt == 0    % FIR filter
  75.     b = fir1(nfilt,1/r);
  76. % prepend sequence with mirror image of first points so that transients
  77. % can be eliminated. then do same thing at right end of data sequence.
  78.     nfilt = nfilt+1;
  79.     itemp = 2*idata(1) - idata((nfilt+1):-1:2);
  80.     [itemp,zi]=filter(b,1,itemp);
  81.     [odata,zf] = filter(b,1,idata,zi);
  82.     if m == 1    % row data
  83.         itemp = zeros(1,2*nfilt);
  84.     else    % column data
  85.         itemp = zeros(2*nfilt,1);
  86.     end
  87.     itemp(:) = [2*idata(nd)-idata((nd-1):-1:(nd-2*nfilt))];
  88.     itemp = filter(b,1,itemp,zf);
  89. % finally, select only every r'th point from the interior of the lowpass
  90. % filtered sequence
  91.     gd = grpdelay(b,1,8);
  92.     list = round(gd(1)+1.25):r:nd;
  93.     odata = odata(list);
  94.     lod = length(odata);
  95.     nlen = nout - lod;
  96.     nbeg = r - (nd - list(length(list)));
  97.     if m == 1    % row data
  98.         odata = [odata itemp(nbeg:r:nbeg+nlen*r-1)];
  99.     else
  100.         odata = [odata; itemp(nbeg:r:nbeg+nlen*r-1)];
  101.     end
  102. else    % IIR filter
  103.     rip = .05;    % passband ripple in dB
  104.     [b,a] = cheby1(nfilt, rip, .8/r);
  105. % be sure to filter in both directions to make sure the filtered data has zero phase
  106. % make a data vector properly pre- and ap- pended to filter forwards and back
  107. % so end effects can be obliterated.
  108.     odata = filtfilt(b,a,idata);
  109.     nbeg = r - (r*nout - nd);
  110.     odata = odata(nbeg:r:nd);
  111. end
  112.