[Octave-bug-tracker] [bug #37297] Meaningless test in fftfilt.m, patch to address this issue as well enhance final touchup of routine and fix rounding bug

[Dan Sebald]

URL:
  <http://savannah.gnu.org/bugs/?37297>

                 Summary: Meaningless test in fftfilt.m, patch to address this
issue as well enhance final touchup of routine and fix rounding bug
                 Project: GNU Octave
            Submitted by: sebald
            Submitted on: Fri 07 Sep 2012 09:08:33 PM GMT
                Category: None
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: dev
        Operating System: Any

    _______________________________________________________

Details:

One of the tests in fftfilt.m failed on my system:

  ***** 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);
!!!!! test failed
assert (fftfilt (b, r * x),r * r * [1, 1, 0, 0, 0, 0, 0, 0, 0, 0],eps)
expected
 Columns 1 through 3:
   0.00000 + 1.00000i   0.00000 + 1.00000i   0.00000 + 0.00000i
 Columns 4 through 6:
   0.00000 + 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i
 Columns 7 through 9:
   0.00000 + 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i
 Column 10:
   0.00000 + 0.00000i
but got
 Columns 1 through 3:
   0.00000 + 1.00000i   0.00000 + 1.00000i   0.00000 - 0.00000i
 Columns 4 through 6:
  -0.00000 + 0.00000i  -0.00000 + 0.00000i   0.00000 + 0.00000i
 Columns 7 through 9:
  -0.00000 + 0.00000i   0.00000 - 0.00000i  -0.00000 + 0.00000i
 Column 10:
  -0.00000 + 0.00000i
maximum absolute error 2.22478e-16 exceeds tolerance 2.22045e-16
shared variables   scalar structure containing the fields:
    b =
       1   1
    x =
       1   0   0   0   0   0   0   0   0   0
    r = [](0x0)

This brought to light the fact that another test is somewhat meaningless:

%!test
%! b  = rand (10, 1);
%! x  = rand (10, 1);
%! y0 = filter (b, 1, x);
%! y  = filter (b, 1, x);
%! assert (y, y0);

And then came forth a bug in the code itself.


I've attached a patch that

1) Relaxes the tolerance of the test that failed from eps to 2*eps. 
(Explained later.)
2) Changes the meaningless test so that convolution and fftfilt are compared
for both real and imaginary inputs, with appropriate tolerance.
3) In the final touchup of the results, handles the cases of b real and x
imaginary, b imaginary and x real, and b and x imaginary.  (Note that the
existing code for integerizing the output works appropriately for both x and b
integer requirement, which is nice.)
4) Adds another test in which both x and b are complex.
5) Moves the transpose to the very end of the routine because rounding was
failing.

Here is greater details on the tests and choice of tolerance:

%!test
%! b = [1 1];
%! x = [1, zeros(1,9)];
%! 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].');

The tolerance of epsilon was removed.  Because the inputs are integers, the
output is cast to integer and these results should be an exact comparison.  (I
do wonder if there should be an option to remove the final touchups.)  NOTE:
THIS CHANGE BROUGHT OUT A BUG IN rounding.  y was transposed before final
touchup so y could be a ROW vector in which case:

    idx = !any (x - fix (x));

is a value of 1.  Thus for a row vector

    y(:, idx) = round (y(:, idx));

only was rounding the first element of the vector.

%! 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  );

This test passes because although b is complex, x is not complex.  So given
that x is a unit pulse (all imaginary components are 0), the mixing is such
that there is really only one multiply summed throughout the whole transform.

%! assert (fftfilt (b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], 2*eps);

This is the test that failed.  Notice it is similar to the previous, but now
both b and r*x are complex.  No longer is the imaginary component all zeros,
with the implication that butterflies contribute "on both wings".  Hence, with
the unit pulse, but complex operations, I'm thinking that works out to two
multiplications within summations, hence 2*eps can be lost.

%! assert (fftfilt (b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps  );

Similar to a couple tests previous, just different orientation.

%! 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 (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);

I simply added all the combinations of real and imaginary inputs to test the
enhancements to the final touchups.

%!test
%! b  = rand (10, 1);
%! x  = rand (10, 1);
%! y0 = filter (b, 1, x);
%! y  = fftfilt (b, x);
%! assert (y0, y, 15*eps);
%! y0 = filter (b*i, 1, x*i);
%! y  = fftfilt (b*i, x*i);
%! assert (y0, y, 20*eps);

This is a test that really checks the integrity of the algorithm.  (The
meaningless test prior to the patch.)  After some initial guesses and trial
and error, I tested 1e6 trials, relaxing tolerance so that they all passed.  A
tolerance of 15*eps worked for real inputs, but 15*eps failed long into tests
for the imaginary inputs.  I tried that test with 15*eps tolerance again and
all 1e6 computations were within tolerance.  So, I bumped that test to 16*eps
for both cases.

%!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);

OK, this test pushes the limit: both b and x are complex, non-integer so there
is no zeroing out either the real or imaginary component.  The 50*eps
tolerance did not pass 1e6 trials, so I changed that to 55*eps and that
passed.



    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Fri 07 Sep 2012 09:08:33 PM GMT  Name: octave-fftfilt-2012sep07.patch 
Size: 3kB   By: sebald

<http://savannah.gnu.org/bugs/download.php?file_id=26501>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?37297>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/