changeset 16429:586972e3ea7a

Correct errors in the fftfilt rounding (bug #37297), add more robust tests. * fftfilt.m (fftfilt): Extend the cleanup of real-only results to that of imaginary-only, and add tests to check proper operation. Replace logical result with find() for rouding index computation. Move transpose to after rounding (bug #37297). Add tests for exact accuracy. Replace y = filter() with y = fftfilt () in second last test. Add random input tests with relaxed tolerance.
author Daniel J Sebald <daniel.sebald@ieee.org>
date Wed, 03 Apr 2013 18:23:36 -0500
parents 610e02bf9579
children 1766d8655006
files scripts/signal/fftfilt.m
diffstat 1 files changed, 69 insertions(+), 23 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/signal/fftfilt.m	Wed Aug 08 20:51:24 2012 -0700
+++ b/scripts/signal/fftfilt.m	Wed Apr 03 18:23:36 2013 -0500
@@ -94,53 +94,99 @@
   endif
 
   y = y(1:r_x, :);
-  if (transpose)
-    y = y.';
+
+  ## 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.
+  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
 
-  ## 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 both x and b are integer in both real and imaginary
+  ##   components, y should be integer.
+  if (! any (b - fix (b)))
+    idx = find (! any (x - fix (x)));
+    y (:, idx) = round (y (:, idx));
+  endif
 
-  if (isreal (b) && isreal (x))
-    y = real (y);
-  endif
-  if (! any (b - fix (b)))
-    idx = !any (x - fix (x));
-    y(:, idx) = round (y(:, idx));
+  ## Transpose after cleanup, otherwise rounding fails.
+  if (transpose)
+    y = y.';
   endif
 
 endfunction
 
 
 %!shared b, x, r
+
 %!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)