home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-bin / lib / octave / 1.1.1 / m / signal / fftfilt.m < prev    next >
Encoding:
Text File  |  1996-10-12  |  2.8 KB  |  101 lines
  1. # Copyright (C) 1995 John W. Eaton
  2. # 
  3. # This file is part of Octave.
  4. # 
  5. # Octave is free software; you can redistribute it and/or modify it
  6. # under the terms of the GNU General Public License as published by the
  7. # Free Software Foundation; either version 2, or (at your option) any
  8. # later version.
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. # for more details.
  14. # 
  15. # You should have received a copy of the GNU General Public License
  16. # along with Octave; see the file COPYING.  If not, write to the Free
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19. function y = fftfilt (b, x, N)
  20.  
  21. # usage:  fftfilt (b, x [, N])
  22. #
  23. # y = fftfilt (b, x) filters x with the FIR filter b using the FFT.
  24. # y = fftfilt (b, x, N) uses the overlap-add method to filter x with
  25. # b using an N-point FFT.
  26.  
  27. # Written by KH (Kurt.Hornik@ci.tuwien.ac.at) on Sep 3, 1994
  28.  
  29. # Reference:  Oppenheim & Schafer (1989).  Discrete-time Signal
  30. # Processing (Chapter 8).  Prentice-Hall.
  31.  
  32. # If N is not specified explicitly, we do not use the overlap-add 
  33. # method at all because loops are really slow.  Otherwise, we only
  34. # ensure that the number of points in the FFT is the smallest power
  35. # of two larger than N and length(b).  This could result in length
  36. # one blocks, but if the user knows better ...
  37.   
  38.   if (nargin < 2 || nargin > 3)
  39.     usage (" fftfilt (b, x [, N])");
  40.   endif
  41.   
  42.   [r_x, c_x] = size (x);
  43.   [r_b, c_b] = size (b);
  44.   if (! (min ([r_x c_x]) == 1 || min ([r_b c_b]) == 1))
  45.     error ("fftfilt: both x and b should be vectors");
  46.   endif
  47.   l_x  = r_x * c_x;
  48.   l_b  = r_b * c_b;
  49.  
  50.   if ((l_x == 1) && (l_b == 1))
  51.     y = b * x;  
  52.     return;
  53.   endif
  54.   
  55.   x = reshape (x, 1, l_x);
  56.   b = reshape (b, 1, l_b);
  57.   
  58.   if (nargin == 2)
  59. # Use FFT with the smallest power of 2 which is >= length (x) +
  60. # length (b) - 1 as number of points ...
  61.     N    = 2^(ceil (log (l_x + l_b - 1) / log(2)));
  62.     y    = ifft (fft (x, N) .* fft(b, N));
  63.   else
  64. # Use overlap-add method ...
  65.     if !(is_scalar (N))
  66.       error ("fftfilt: N has to be a scalar");
  67.     endif
  68.     N = 2^(ceil (log (max ([N l_b])) / log(2)));
  69.     L = N - l_b + 1;
  70.     B = fft (b, N);
  71.     R = ceil (l_x / L);
  72.     y = zeros (1, l_x);
  73.     for r = 1:R;
  74.       lo  = (r - 1) * L + 1;
  75.       hi  = min (r * L, l_x);
  76.       tmp = ifft (fft (x(lo:hi), N) .* B);
  77.       hi  = min (lo+N-1, l_x);
  78.       y(lo:hi) = y(lo:hi) + tmp(1:(hi-lo+1));
  79.     endfor  
  80.   endif
  81.     
  82.   y = reshape (y(1:l_x), r_x, c_x);
  83.  
  84. # Final cleanups:  if both x and b are real respectively integer, y
  85. # should also be
  86.  
  87.   if (! (any (imag (x)) || any (imag (b))))
  88.     y = real (y);
  89.   endif
  90.   if (! (any (x - round (x)) || any (b - round (b))))
  91.     y = round (y);
  92.   endif
  93.  
  94. endfunction
  95.  
  96.  
  97.  
  98.  
  99.  
  100.  
  101.