Fourier transform rounding error

Question

I'm messing around with Fourier transformations. Now I've created a class that does an implementation of the DFT (not doing anything like FFT atm). This is the implementation I've used:

public static Complex[] Dft(double[] data)
    {
        int length = data.Length;
        Complex[] result = new Complex[length];

        for (int k = 1; k <= length; k++)
        {
            Complex c = Complex.Zero;
            for (int n = 1; n <= length; n++)
            {
                c += Complex.FromPolarCoordinates(data[n-1], (-2 * Math.PI * n * k) / length);
            }
            result[k-1] =  1 / Math.Sqrt(length) * c;
        }
        return result;
    }

And these are the results I get from Dft({2,3,4})

Well it seems pretty okay, since those are the values I expect. There is only one thing I find confusing. And it all has to do with the rounding of doubles.

First of all, why are the first two numbers not exactly the same (0,8660..443 8 ) vs (0,8660..443). And why can't it calculate a zero, where you'd expect it. I know 2.8E-15 is pretty close to zero, but well it's not.

Anyone know how these, marginal, errors occur and if I can and want to do something about it.

It might seem that there's not a real problem, because it's just small errors. However, how do you deal with these rounding errors if you're for example comparing 2 values.

5,2 + 0i != 5,1961524 + i2.828107*10^-15

Cheers

Answer 1

I think you've already explained it to yourself - limited precision means limited precision. End of story.

If you want to clean up the results, you can do some rounding of your own to a more reasonable number of siginificant digits - then your zeros will show up where you want them.

To answer the question raised by your comment, don't try to compare floating point numbers directly - use a range:

if (Math.Abs(float1 - float2) < 0.001) {
  // they're the same!
}

The comp.lang.c FAQ has a lot of questions & answers about floating point, which you might be interested in reading.

Answer 2

From http://support.microsoft.com/kb/125056

Emphasis mine.

There are many situations in which precision, rounding, and accuracy in floating-point calculations can work to generate results that are surprising to the programmer. There are four general rules that should be followed:

1. In a calculation involving both single and double precision, the result will not usually be any more accurate than single precision. If double precision is required, be certain all terms in the calculation, including constants, are specified in double precision.

2. Never assume that a simple numeric value is accurately represented in the computer. Most floating-point values can't be precisely represented as a finite binary value. For example .1 is .0001100110011... in binary (it repeats forever), so it can't be represented with complete accuracy on a computer using binary arithmetic, which includes all PCs.

3. Never assume that the result is accurate to the last decimal place. There are always small differences between the "true" answer and what can be calculated with the finite precision of any floating point processing unit.

4. Never compare two floating-point values to see if they are equal or not- equal. This is a corollary to rule 3. There are almost always going to be small differences between numbers that "should" be equal. Instead, always check to see if the numbers are nearly equal. In other words, check to see if the difference between them is very small or insignificant.

---

Note that although I referenced a microsoft document, this is not a windows problem. It's a problem with using binary and is in the CPU itself.

And, as a second side note, I tend to use the Decimal datatype instead of double: See this related SO question: decimal vs double! - Which one should I use and when?

Answer 3

In C# you'll want to use the 'decimal' type, not double for accuracy with decimal points.

As to the 'why'... repsensenting fractions in different base systems gives different answers. For example 1/3 in a base 10 system is 0.33333 recurring, but in a base 3 system is 0.1.

The double is a binary value, at base 2. When converting to base 10 decimal you can expect to have these rounding errors.