Fri 07 Sep 2012 09:08:33 PM UTC, original submission:
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, rx), rr*[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.22478e16 exceeds tolerance 2.22045e16
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, rx), rr[1 1 0 0 0 0 0 0 0 0], 2eps);
This is the test that failed. Notice it is similar to the previous, but now both b and rx 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 2eps 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 (bi, xi);
%! 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 (bi, 1, xi);
%! y = fftfilt (bi, xi);
%! 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 15eps worked for real inputs, but 15eps failed long into tests for the imaginary inputs. I tried that test with 15eps tolerance again and all 1e6 computations were within tolerance. So, I bumped that test to 16eps 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, noninteger so there is no zeroing out either the real or imaginary component. The 50eps tolerance did not pass 1e6 trials, so I changed that to 55eps and that passed.
