From bug-octave-request at bevo dot che dot wisc dot edu Wed Jan 16 14:30:29 2002
Subject: Adapted fftfilt that handles matrices
From: Roberto Hernandez <robher at adinet dot com dot uy>
To: bug-octave at bevo dot che dot wisc dot edu
Date: Wed, 16 Jan 2002 17:32:31 -0300

This is a multi-part message in MIME format.
--------------000604000500060007000402
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit

Hello,

I've been using the "fftfilt" function to do some Filterbank 
time-frequency analysis. I found that Paul Kienzle has a modified 
version of "fftfilt" in Octave-Forge which handles the filtering of many 
signals simultaneously. Each signal is one column in a matrix, which is 
passed to fftfilt instead of a single-signal vector.

There was one small bug when using complex matrices. I've made Paul 
aware of it and he's fixed it.

Anyway, I thought it would be a positive change in the Octave sources. 
I'm sending in the patch against version 1.17 of "fftfilt.m" at:

http://www.octave.org/cgi-bin/cvsweb.cgi/~checkout~/octave/scripts/signal/fftfilt.m

Hope it's as useful to everyone as it's been to me.

Roberto.

--------------000604000500060007000402
Content-Type: text/plain;
 name="fftfilt.patch"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="fftfilt.patch"

*** fftfilt.m	Wed Jan 16 17:03:50 2002
--- fftfilt.m.old	Wed Jan 16 17:15:55 2002
***************
*** 25,43 ****
  ##
  ## Given the optional third argument,  at var{n}, @code{fftfilt} uses the
  ## overlap-add method to filter  at var{x} with @var{b} using an N-point FFT.
- ##
- ## If  at var{x} is a matrix, filter each column of the matrix.
  ##  at end deftypefn
  
  ## Author: KH <Kurt dot Hornik at ci dot tuwien dot ac dot at>
  ## Created: 3 September 1994
  ## Adapted-By: jwe
  
- ## 2000-04-04 Paul Kienzle <pkienzle at users dot sf dot net>
- ## * handle matrices
- ## 2001-01-15 RH <robher at adinet dot com dot uy>
- ## * test for real/integer using matrix rather than column-wise ops
- 
  function y = fftfilt (b, x, N)
  
    ## If N is not specified explicitly, we do not use the overlap-add
--- 25,36 ----
***************
*** 50,71 ****
      usage (" fftfilt (b, x, N)");
    endif
  
-   transpose = ( rows (x) == 1 );
-   if transpose, x = x.'; endif
    [r_x, c_x] = size (x);
    [r_b, c_b] = size (b);
!   if min ([r_b, c_b]) != 1
!     error ("fftfilt: b should be a vector");
    endif
    l_b  = r_b * c_b;
!   b = reshape (b, l_b, 1);
  
    if (nargin == 2)
      ## Use FFT with the smallest power of 2 which is >= length (x) +
      ## length (b) - 1 as number of points ...
!     N    = 2^(ceil (log (r_x + l_b - 1) / log(2)));
!     B = fft (b, N);
!     y = ifft (fft (x, N) .* B (:,ones(1,c_x)));
    else
      ## Use overlap-add method ...
      if (! (is_scalar (N)))
--- 43,69 ----
      usage (" fftfilt (b, x, N)");
    endif
  
    [r_x, c_x] = size (x);
    [r_b, c_b] = size (b);
!   if (! (min ([r_x, c_x]) == 1 || min ([r_b, c_b]) == 1))
!     error ("fftfilt: both x and b should be vectors");
    endif
+   l_x  = r_x * c_x;
    l_b  = r_b * c_b;
! 
!   if ((l_x == 1) && (l_b == 1))
!     y = b * x;
!     return;
!   endif
! 
!   x = reshape (x, 1, l_x);
!   b = reshape (b, 1, l_b);
  
    if (nargin == 2)
      ## Use FFT with the smallest power of 2 which is >= length (x) +
      ## length (b) - 1 as number of points ...
!     N    = 2^(ceil (log (l_x + l_b - 1) / log(2)));
!     y    = ifft (fft (x, N) .* fft(b, N));
    else
      ## Use overlap-add method ...
      if (! (is_scalar (N)))
***************
*** 74,106 ****
      N = 2^(ceil (log (max ([N, l_b])) / log(2)));
      L = N - l_b + 1;
      B = fft (b, N);
!     B = B (:, ones (c_x,1));
!     R = ceil (r_x / L);
!     y = zeros (r_x, c_x);
      for r = 1:R;
        lo  = (r - 1) * L + 1;
!       hi  = min (r * L, r_x);
!       tmp = zeros (N, c_x);
!       tmp (1:(hi-lo+1), :) = x (lo:hi,:);
!       tmp = ifft (fft (tmp) .* B);
!       hi  = min (lo+N-1, r_x);
!       y (lo:hi, :) = y (lo:hi, :) + tmp (1:(hi-lo+1), :);
      endfor
    endif
  
!   y = y(1:r_x,:);
!   if transpose, y=y.'; endif
  
    ## Final cleanups: if both x and b are real respectively integer, y
!   ## should also be; note that this doesn't handle the case where x is
!   ## mixed real/complex
  
!   if ( isreal (b) && isreal (x) )
      y = real (y);
    endif
!   if ( !any (b - round (b)) )
!     idx = !any (x - round (x));
!     y (:, idx) = round (y (:, idx));
    endif
  
  endfunction
--- 72,98 ----
      N = 2^(ceil (log (max ([N, l_b])) / log(2)));
      L = N - l_b + 1;
      B = fft (b, N);
!     R = ceil (l_x / L);
!     y = zeros (1, l_x);
      for r = 1:R;
        lo  = (r - 1) * L + 1;
!       hi  = min (r * L, l_x);
!       tmp = ifft (fft (x(lo:hi), N) .* B);
!       hi  = min (lo+N-1, l_x);
!       y(lo:hi) = y(lo:hi) + tmp(1:(hi-lo+1));
      endfor
    endif
  
!   y = reshape (y(1:l_x), r_x, c_x);
  
    ## Final cleanups: if both x and b are real respectively integer, y
!   ## should also be
  
!   if (! (any (imag (x)) || any (imag (b))))
      y = real (y);
    endif
!   if (! (any (x - round (x)) || any (b - round (b))))
!     y = round (y);
    endif
  
  endfunction

--------------000604000500060007000402--




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------