octave_packages/m/signal/fftfilt.m

   1 ## Copyright (C) 1994-2012 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
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License 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, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} fftfilt (@var{b}, @var{x}, @var{n})
  21 ##
  22 ## With two arguments, @code{fftfilt} filters @var{x} with the FIR filter
  23 ## @var{b} using the FFT.
  24 ##
  25 ## Given the optional third argument, @var{n}, @code{fftfilt} uses the
  26 ## overlap-add method to filter @var{x} with @var{b} using an N-point FFT.
  27 ##
  28 ## If @var{x} is a matrix, filter each column of the matrix.
  29 ## @seealso{filter, filter2}
  30 ## @end deftypefn
  31
  32 ## Author: Kurt Hornik <Kurt.Hornik@wu-wien.ac.at>
  33 ## Created: 3 September 1994
  34 ## Adapted-By: jwe
  35
  36 function y = fftfilt (b, x, n)
  37
  38   ## If N is not specified explicitly, we do not use the overlap-add
  39   ## method at all because loops are really slow.  Otherwise, we only
  40   ## ensure that the number of points in the FFT is the smallest power
  41   ## of two larger than N and length(b).  This could result in length
  42   ## one blocks, but if the user knows better ...
  43
  44   if (nargin < 2 || nargin > 3)
  45     print_usage ();
  46   endif
  47
  48   transpose = (rows (x) == 1);
  49
  50   if (transpose)
  51     x = x.';
  52   endif
  53
  54   [r_x, c_x] = size (x);
  55   [r_b, c_b] = size (b);
  56
  57   if (! isvector (b))
  58     error ("fftfilt: B must be a vector");
  59   endif
  60
  61   if (ndims (x) != 2)
  62     error ("fftfilt: X must be a 1-D or 2-D array");
  63   endif
  64
  65   l_b = r_b * c_b;
  66   b = reshape (b, l_b, 1);
  67
  68   if (nargin == 2)
  69     ## Use FFT with the smallest power of 2 which is >= length (x) +
  70     ## length (b) - 1 as number of points ...
  71     n = 2 ^ nextpow2 (r_x + l_b - 1);
  72     B = fft (b, n);
  73     y = ifft (fft (x, n) .* B(:, ones (1, c_x)));
  74   else
  75     ## Use overlap-add method ...
  76     if (! (isscalar (n)))
  77       error ("fftfilt: N has to be a scalar");
  78     endif
  79     n = 2 ^ nextpow2 (max ([n, l_b]));
  80     L = n - l_b + 1;
  81     B = fft (b, n);
  82     B = B(:, ones (c_x,1));
  83     R = ceil (r_x / L);
  84     y = zeros (r_x, c_x);
  85     for r = 1:R;
  86       lo = (r - 1) * L + 1;
  87       hi = min (r * L, r_x);
  88       tmp = zeros (n, c_x);
  89       tmp(1:(hi-lo+1),:) = x(lo:hi,:);
  90       tmp = ifft (fft (tmp) .* B);
  91       hi  = min (lo+n-1, r_x);
  92       y(lo:hi,:) = y(lo:hi,:) + tmp(1:(hi-lo+1),:);
  93     endfor
  94   endif
  95
  96   y = y(1:r_x, :);
  97   if (transpose)
  98     y = y.';
  99   endif
 100
 101   ## Final cleanups: If both x and b are real, y should be real.
 102   ## If both x and b are integer, y should be integer.
 103
 104   if (isreal (b) && isreal (x))
 105     y = real (y);
 106   endif
 107   if (! any (b - fix (b)))
 108     idx = !any (x - fix (x));
 109     y(:, idx) = round (y(:, idx));
 110   endif
 111
 112 endfunction
 113
 114
 115 %!shared b, x, r
 116 %!test
 117 %!  b = [1 1];
 118 %!  x = [1, zeros(1,9)];
 119 %!  assert(fftfilt(b,  x  ), [1 1 0 0 0 0 0 0 0 0]  , eps);
 120 %!  assert(fftfilt(b,  x.'), [1 1 0 0 0 0 0 0 0 0].', eps);
 121 %!  assert(fftfilt(b.',x  ), [1 1 0 0 0 0 0 0 0 0]  , eps);
 122 %!  assert(fftfilt(b.',x.'), [1 1 0 0 0 0 0 0 0 0].', eps);
 123
 124 %!test
 125 %!  r = sqrt(1/2) * (1+i);
 126 %!  b = b*r;
 127 %!  assert(fftfilt(b, x  ), r*[1 1 0 0 0 0 0 0 0 0]  , eps);
 128 %!  assert(fftfilt(b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], eps);
 129 %!  assert(fftfilt(b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps);
 130
 131 %!test
 132 %!  b = [1 1];
 133 %!  x = zeros (10,3); x(1,1)=-1; x(1,2)=1;
 134 %!  y0 = zeros (10,3); y0(1:2,1)=-1; y0(1:2,2)=1;
 135 %!  y = fftfilt (b, x);
 136 %!  assert (y,y0);
 137
 138 %!test
 139 %!  b  = rand (10, 1);
 140 %!  x  = rand (10, 1);
 141 %!  y0 = filter (b, 1, x);
 142 %!  y  = filter (b, 1, x);
 143 %!  assert (y, y0);
 144
 145 %% Test input validation
 146 %!error fftfilt (1)
 147 %!error fftfilt (1, 2, 3, 4)
 148 %!error fftfilt (ones (2), 1)
 149 %!error fftfilt (2, ones (3,3,3))
 150 %!error fftfilt (2, 1, ones (2))
 151