# HG changeset patch
# User Daniel J Sebald <daniel.sebald@ieee.org>
# Date 1350618781 18000
# Node ID 060fcc04c839eafc8bffd2e37790b5d804f21f04
# Parent  9db32cabeacfe74de320c4eddeff9ef1d90150ae
Correct errors in the fftfilt rounding routine and add more robust tests (bug #37297).

* fftfilt.m (fftfilt): Extend the cleanup of real-only results to that of imaginary-only making four cases real/real (real-only), real/imaginary (imaginary-only), imaginary/real (imaginary-only) and imaginary/imaginary (real-only). Add some tests to check proper operation. Correct the integerizing (rounding) index computation using find() rather than the logical result. Move transpose to after the rounding for correct rounding result. Add some tests for proper operation that are exactly accurate rather than within epsilon. Replace meaningless test that had y0 and y equal to the same vector, by definition (bug #37297). This necessitates an increase in tolerance by a couple factors. Add some tests with random inputs to fully exercise the FFT algorithm, and based on some statistical tests noticably increase tolerance values otherwise it isn't uncommon for one of these tests to fail because FFTs have lots of numerical effects due to multiple intermediate stages.

diff -r 9db32cabeacf -r 060fcc04c839 scripts/signal/fftfilt.m
--- a/scripts/signal/fftfilt.m	Wed Oct 17 20:13:19 2012 -0700
+++ b/scripts/signal/fftfilt.m	Thu Oct 18 22:53:01 2012 -0500
@@ -94,21 +94,44 @@
   endif
 
   y = y(1:r_x, :);
+
+  ## Final cleanups:
+  ## - If both b and x are real, y should be real.
+  ## - 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.
+  ## - If both x and b are integer in both real and imaginary
+  ##   components, y should be integer.
+
+  xisreal = zeros (1,size(x,2));
+  xisimag = xisreal;
+  for cx = 1:size(x,2)
+    if (all (imag (x(:,cx)) == 0))
+      xisreal (cx) = 1;
+    elseif (all (real (x (:,cx)) == 0))
+      xisimag (cx) = 1;
+    endif
+  endfor
+  xisreal = find(xisreal);
+  xisimag = find(xisimag);
+  if (all (imag (b) == 0))
+    y (:,xisreal) = real (y (:,xisreal));
+    y (:,xisimag) = complex (real (y (:,xisimag)) * 0, imag (y (:,xisimag)));
+  elseif (all (real (b) == 0))
+    y (:,xisreal) = complex (real (y (:,xisreal)) * 0, imag (y (:,xisreal)));
+    y (:,xisimag) = real (y (:,xisimag));
+  endif
+
+  if (! any (b - fix (b)))
+    idx = find (! any (x - fix (x)));
+    y (:, idx) = round (y (:, idx));
+  endif
+
+  ## Transpose after cleanup, otherwise rounding fails.
   if (transpose)
     y = y.';
   endif
 
-  ## Final cleanups: If both x and b are real, y should be real.
-  ## If both x and b are integer, y should be integer.
-
-  if (isreal (b) && isreal (x))
-    y = real (y);
-  endif
-  if (! any (b - fix (b)))
-    idx = !any (x - fix (x));
-    y(:, idx) = round (y(:, idx));
-  endif
-
 endfunction
 
 
@@ -116,31 +139,53 @@
 %!test
 %! b = [1 1];
 %! x = [1, zeros(1,9)];
-%! assert (fftfilt (b,  x  ), [1 1 0 0 0 0 0 0 0 0]  , eps);
-%! assert (fftfilt (b,  x.'), [1 1 0 0 0 0 0 0 0 0].', eps);
-%! assert (fftfilt (b.',x  ), [1 1 0 0 0 0 0 0 0 0]  , eps);
-%! assert (fftfilt (b.',x.'), [1 1 0 0 0 0 0 0 0 0].', eps);
-
+%! assert (fftfilt (b,  x  ), [1 1 0 0 0 0 0 0 0 0]  );
+%! assert (fftfilt (b,  x.'), [1 1 0 0 0 0 0 0 0 0].');
+%! assert (fftfilt (b.',x  ), [1 1 0 0 0 0 0 0 0 0]  );
+%! assert (fftfilt (b.',x.'), [1 1 0 0 0 0 0 0 0 0].');
+%! assert (fftfilt (b,  [x.' x.']), [1 1 0 0 0 0 0 0 0 0].'*[1 1]);
+%! assert (fftfilt (b,  [x.'+eps x.']) == [1 1 0 0 0 0 0 0 0 0].'*[1 1], [false(10, 1) true(10, 1)]);
+%
 %!test
 %! r = sqrt (1/2) * (1+i);
 %! b = b*r;
-%! assert (fftfilt (b, x  ), r*[1 1 0 0 0 0 0 0 0 0]  , eps);
-%! assert (fftfilt (b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], eps);
-%! assert (fftfilt (b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps);
+%! assert (fftfilt (b, x  ), r*[1 1 0 0 0 0 0 0 0 0]  , eps  );
+%! assert (fftfilt (b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], 2*eps);
+%! assert (fftfilt (b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps  );
 
 %!test
-%! b = [1 1];
-%! x = zeros (10,3); x(1,1)=-1; x(1,2)=1;
+%! b  = [1 1];
+%! x  = zeros (10,3); x(1,1)=-1; x(1,2)=1;
 %! y0 = zeros (10,3); y0(1:2,1)=-1; y0(1:2,2)=1;
-%! y = fftfilt (b, x);
-%! assert (y,y0);
+%! y  = fftfilt (b, x);
+%! assert (y0, y);
+%! y  = fftfilt (b*i, x);
+%! assert (y0*i, y);
+%! y  = fftfilt (b, x*i);
+%! assert (y0*i, y);
+%! y  = fftfilt (b*i, x*i);
+%! assert (-y0, y);
+%! x  = rand (10, 1);
+%! y  = fftfilt (b, [x x*i]);
+%! assert (true, isreal (y(:,1)));
+%! assert (false, any (real (y(:,2))));
 
 %!test
 %! b  = rand (10, 1);
 %! x  = rand (10, 1);
 %! y0 = filter (b, 1, x);
-%! y  = filter (b, 1, x);
-%! assert (y, y0);
+%! y  = fftfilt (b, x);
+%! assert (y0, y, 16*eps);
+%! y0 = filter (b*i, 1, x*i);
+%! y  = fftfilt (b*i, x*i);
+%! assert (y0, y, 16*eps);
+
+%!test
+%! b  = rand (10, 1) + i*rand (10, 1);
+%! x  = rand (10, 1) + i*rand (10, 1);
+%! y0 = filter (b, 1, x);
+%! y  = fftfilt (b, x);
+%! assert (y0, y, 55*eps);
 
 %% Test input validation
 %!error fftfilt (1)
