Mercurial > hg > octave-lojdl
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 +++ b/scripts/signal/fftfilt.m @@ -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)