[vsipl++] [PATCH] Implement Fir<>, native C++ version

[Jules Bergmann]

Nathan (Jasper) Myers wrote:
> The following patch has been committed.
> 
> It implements vsip::Fir<> using native C++ code, and a comprehensive
> test of all its modes.  
> 
> Note that a few bits of the test are commented out; it uses 
> vsip::Convolution<> to generate the reference output, and that has 
> a little bug I haven't got to tracking down yet.
> 

Nathan,

For your failing cases:

> +  // FIXME: this exposes a bug in vsip::Convolution.
> +  test_fir<float,vsip::sym_even_len_even>(1,32,256);

> +  test_fir<float,vsip::sym_even_len_even>(2,32,256);

> +  test_fir<float,vsip::sym_even_len_even>(3,32,256);

> +  test_fir<float,vsip::sym_even_len_odd>(1,32,256);

> +  test_fir<float,vsip::sym_even_len_odd>(2,32,256);

> +  test_fir<float,vsip::sym_even_len_odd>(3,32,256);

There are two causes for differences between the Fir and Convolution 
results.

The first has to do with exceeding the dynamic range of float:

These kernels have 63 to 64 coefficients (due to symmetry) ranging from 
1 to 63 with an average value of 32.  This gives the filter a gain of 
approx 64 * 32 = 2048 or 2^11.  Floating point numbers have 24 bits of 
precision.  As soon as the input values start to have magnitudes in the 
2^13 range, the filter output magnitude will be in the 2^24 range.  At 
this point, small preturbations in input values (~ 2^1) will be outside 
the dynamic range of the floating point value.  Two computations 
performing the same operations but in a different order will likely 
exceed the dynamic range at different points, resulting in different 
rounding and different answers.

Experimentally, for your first case above (even, D=1, M=32, N=256), when 
running the non-IPP convolution the failures start to occur around n= 
8240, approximate 2^13.  The output value of the filter is 1.681 * 10^7, 
just a bit over 2^24.

By changing the precision to double (which as 53 bits of precision), 
these miscompare go away.  However, if you made your input longer, or if 
you scaled the values appropriately, you could recreate the same type of 
miscompare.


The second cause has to do with algorithm choice.

When using IPP to perform convolution, I see a different type of 
miscompares for smaller values that should be within the dynamic range 
of float.  I suspect this may be do to IPP using a different algorithm 
underneath (such as FFT based convolution) where the size of the FFT 
spreads noise over a wider range than the size of the convolution kernel.


Either way, you shouldn't use '==' to compare floating point values in 
tests.  In most element-wise cases you should use 'equal' (which uses 
almost_equal for floats and doubles, which chekcs relative and absolute 
errors against some semi-arbitrary bounds).  For comparing results of 
non-elementwise operations, its a bit more tricky.  In this case, I 
would recommend using the 'error_db' comparison to compare the results, 
perhaps with a maximum value derived from the entire view to account for 
the noise from the FFT based convolution.

				-- Jules