[SciPy-dev] accuracy problem in filterdesign.py

[Pearu Peterson]

Hi Heiko,

On Sun, 24 Mar 2002, Heiko Henkelmann wrote:

> while writing a test driver for a minimum phase calculation routine I came
> across the following problem. It is causing asymmetriesin the output of
> freqz. In line 107 and 109 in filter_design.py the follwoing is happening:
> 
> 
> >>> N=512
> >>> lastpoint=2*pi
> >>> w1=arange(0,lastpoint,lastpoint/N)
> >>> w2=arange(0,N)*(lastpoint/N)
> >>> lastpoint-w1[511]-w1[1]
> -6.3546390371982397e-014
> >>> lastpoint-w2[511]-w2[1]
> 4.0245584642661925e-016
> >>> w1[511]
> 6.2709134608765646
> >>> w2[511]
> 6.2709134608765007
> >>> w2[511]-w1[511]
> -6.3948846218409017e-014
> >>>

What Numeric are you using? Platform? Compiler? Because I don't see these
rounding errors on Linux, Numeric 20.2.1, gcc-2.95.4:

Python 2.1.2 (#1, Mar 16 2002, 00:56:55) 
[GCC 2.95.4 20011002 (Debian prerelease)] on linux2
Type "copyright", "credits" or "license" for more information.
>>> from Numeric import *
>>> N=512
>>> lastpoint=2*pi
>>> w1=arange(0,lastpoint,lastpoint/N)
>>> w2=arange(0,N)*(lastpoint/N)
>>> lastpoint-w1[511]-w1[1]
4.0245584642661925e-16
>>> lastpoint-w2[511]-w2[1]
4.0245584642661925e-16
>>> w1[511]
6.2709134608765007
>>> w2[511]
6.2709134608765007
>>> w2[511]-w1[511]
0.0

If I recall correctly, then gcc uses longer floats when doing internal
float operations and may be that's why I don't have these errors.

> It appears that arange for floating point values is accumulating an error
> during summation. Do you think it would be worthwhile to tweak
> filter_design.py to work around this problem?

Could you bug report this to Numeric people? The corresponding code
fragment of arange is

    value = start;
    for (i=0; i < length; i++) {
        dbl_descr->cast[type](&value, 0, rptr, 0, 1);
        value += step;
        rptr += elsize;
    }

and suggest the following fix:

    for (i=0; i < length; ++i) {
        value = start + i * step;
        dbl_descr->cast[type](&value, 0, rptr, 0, 1);
        rptr += elsize;
    }

Regards,
	Pearu