# HG changeset patch
# User Daniel J Sebald
# Date 1365201250 18000
# Node ID d7ee444b4bf0edd139c03b2599bfabb5e0c79746
# Parent a971d8bdaadc7799a4b1c772f9449f4f4e447734
Modify fftfilt to be more user friendly (patch #7999).
* fftfilt.m (fftfilt): Remove line that modifies N to something other than
user input. Error if N is too short. Warn if N is so short that there is too
much loop iteration. Move real and imaginary tests inside filter conditional.
diff r a971d8bdaadc r d7ee444b4bf0 scripts/signal/fftfilt.m
 a/scripts/signal/fftfilt.m Fri Apr 05 03:04:56 2013 0400
+++ b/scripts/signal/fftfilt.m Fri Apr 05 17:34:10 2013 0500
@@ 20,10 +20,15 @@
## @deftypefn {Function File} {} fftfilt (@var{b}, @var{x}, @var{n})
##
## With two arguments, @code{fftfilt} filters @var{x} with the FIR filter
## @var{b} using the FFT.
+## @var{b} using an FFT of both sequences, multiplication and inverse FFT
+## of the product.
##
## Given the optional third argument, @var{n}, @code{fftfilt} uses the
## overlapadd method to filter @var{x} with @var{b} using an Npoint FFT.
+## It is good practice to choose @var{N} sufficiently large (e.g., twice the
+## length of @var{b}) in order to achieve efficiency improvement. It is also
+## good practice to choose N a power of 2 to improve the efficiency of the
+## FFT and inverse FFT algorithms.
##
## If @var{x} is a matrix, filter each column of the matrix.
## @seealso{filter, filter2}
@@ 36,10 +41,10 @@
function y = fftfilt (b, x, n)
## If N is not specified explicitly, we do not use the overlapadd
 ## method at all because loops are really slow. Otherwise, we only
 ## ensure that the number of points in the FFT is the smallest power
 ## of two larger than N and length(b). This could result in length
 ## one blocks, but if the user knows better ...
+ ## method at all because loops are really slow. If N is specified,
+ ## an error occurs if the wraparound effect covers the whole FFT
+ ## and a warning occurs if the wraparound is so great that it makes
+ ## overlapadd extremely inefficient.
if (nargin < 2  nargin > 3)
print_usage ();
@@ 76,8 +81,13 @@
if (! (isscalar (n)))
error ("fftfilt: N has to be a scalar");
endif
 n = 2 ^ nextpow2 (max ([n, l_b]));
+ if (n < l_b)
+ error ("fftfilt: N not long enough to create block free of wraparound");
+ endif
L = n  l_b + 1;
+ if (L < l_b)
+ warning ("fftfilt: N is small relative to length of filter, may run slowly");
+ end
B = fft (b, n);
B = B(:, ones (c_x,1));
R = ceil (r_x / L);
@@ 101,13 +111,14 @@
##  If b is real and x is imaginary, y should be imaginary.
##  If b is imaginary and x is real, y should be imaginary.
##  If both b and x are imaginary, y should be real.
 xisreal = all (imag (x) == 0);
 xisimag = all (real (x) == 0);

if (all (imag (b) == 0))
+ xisreal = all (imag (x) == 0);
+ xisimag = all (real (x) == 0);
y (:,xisreal) = real (y (:,xisreal));
y (:,xisimag) = complex (real (y (:,xisimag)) * 0, imag (y (:,xisimag)));
elseif (all (real (b) == 0))
+ xisreal = all (imag (x) == 0);
+ xisimag = all (real (x) == 0);
y (:,xisreal) = complex (real (y (:,xisreal)) * 0, imag (y (:,xisreal)));
y (:,xisimag) = real (y (:,xisimag));
endif