0
$\begingroup$

I'm trying to do a high precision discrete fourier transform on a signal. To examine the precision, I use a gaussian function as the signal, because the fourier transform is also a gaussian function. fourier transform of gaussian function using Mathematica

The code is as follows

#include <stdio.h>
#include <tgmath.h>

typedef long double complex cplx;
long double PI;
void _fft(cplx buf[], cplx out[], int n, int step)
{
if (step < n) {
    _fft(out, buf, n, step * 2);
    _fft(out + step, buf + step, n, step * 2);

    for (int i = 0; i < n; i += 2 * step) {
        cplx t = cexpl(- I * PI * i / n) * out[i + step];
        buf[i / 2]     = out[i] + t;
        buf[(i + n)/2] = out[i] - t;
    }
}
}

void fft(cplx buf[], int n)
{
cplx out[n];
for (int i = 0; i < n; i++) out[i] = buf[i];

_fft(buf, out, n, 1);
}


int main()
{
const int nPoints = pow(2, 12);
PI = atan2(1, 1) * 4;
cplx dt = 1e-3;
cplx dOmega = 1 / (dt * nPoints);
long double T[nPoints];
long double DOmega[nPoints];

cplx At[nPoints];
cplx tau = 28.4e-3;
for (int i = 0; i < nPoints; ++i)
{
    T[i] = dt*(i-nPoints/2);
    DOmega[i] = dOmega * (i - nPoints / 2);
    At[i] = cexpl(-T[i]*T[i]/2/(tau*tau));
}

fft(At, nPoints);
FILE* fw;
fw = fopen("fft_01.txt", "w+");

for (int i = 0; i < nPoints; ++i)
{
    fprintf(fw, "%.15Le, %.15Le\n", DOmega[i], fabs(At[i]) );
}

return 0;
}

When I change the defined type from float complex to double complex, the result did show a improvement in precision. However, When I change the data type from double complex to long double complex, there's no improvement while the numbers are far away from the minimal limit.

result with float data result using double data result using long double data

I think that the two wings of the peak should go down like the results by Mathematica, but here it seems stuck around 10^(-35) and can not further get down.

When I further apply a iDFT to the result of DFT, it's more apparent that the long double precision has been lost. The blue one is the original signal, with high precision. The green one is the result after iDFT(DFT(x)), which should be coincide with the blue one, but the minimal numbers are much larger. result by iDFT(DFT(x))

Anyone can tell me where is the problem, thanks in advance.

PS: I am on Mac OS 10.13.1, using gcc as the compiler. I have checked the numerical limits on my platform

storage size for float: 4
minimum float positive value: 1.175494e-38
Maximum float positive value: 3.402823e+38
precision value: 6

storage size for double: 8
minimum double positive value: 2.225074e-308
Maximum double positive value: 1.797693e+308
precision value: 15

storage size for long double: 16
minimum long double positive value: 3.362103e-4932
Maximum long double positive value: 1.189731e+4932
precision value: 18
$\endgroup$
  • 1
    $\begingroup$ The most common reason why long double might give the same answer as double is that there is somewhere a bug where a value is computed as double instead of long double, ruining the subsequent accuracy. I haven't tested this, but PI = atan2(1, 1) * 4; looks like it's always double. Try M_PI. Also, I agree with nicoguaro, the question could be a little easier to read if you made the actual question a bit more explicit. $\endgroup$ – Kirill Nov 14 '17 at 16:36
  • 1
    $\begingroup$ I thought you were asking why double results were the same as long double results, but I guess that's not it. What you're expecting is impossible due to roundoff errors. Also, I suspect M_PI might be double too, try 4*atan2(1.0l,1.0l) to make sure, the long double answer doesn't look right because you'd expect the numbers to be smaller than for double. $\endgroup$ – Kirill Nov 14 '17 at 20:04
  • 2
    $\begingroup$ @RoyLiao Your results are consistent with accuracy being degraded in just such a way, and these bugs are always easy to miss. In fact, there exist compilers that treat long double as a synonym for double, so one really does need to carefully check everything. $\endgroup$ – Kirill Nov 14 '17 at 23:37
  • 3
    $\begingroup$ @RoyLiao The limiting factor is not "the limit of double", but the accuracy with which operations are performed (machine precision), which is 2e-16 for doubles and 2e-34 for long doubles. You should probably stop and learn something more about floating point computations before you go on with your experiments, otherwise there are more suprises awaiting you. $\endgroup$ – Federico Poloni Nov 15 '17 at 8:49
  • 3
    $\begingroup$ @RoyLiao If you need 300 digits of accuracy, your best bet is using a multiprecision library such as gmplib.org. $\endgroup$ – Federico Poloni Nov 15 '17 at 13:23
1
$\begingroup$

In addition to the various comments that point out that most built-in functions such as atan or exp operate on double precision variables, there is also the issue that long double on typical x86-type hardware is not actually much more accurate (if at all) than regular double. In particular, long double is not an arbitrary precision data type such as the ones that Mathematica uses -- Mathematica lets you choose the accuracy of variables however you want it, and it will compute things to this accuracy, but this is not the case for C's long double data type that has a fixed accuracy just like double.

$\endgroup$
  • $\begingroup$ Also, different compilers treat the long double type differently. On some compilers, it's just a synonym for "double" On others you will get the 80 bit extended precision type. Also, the OS might turn on the double precision rounding mode by default. You need to check your particular compiler and OS. See linuxtopia.org/online_books/an_introduction_to_gcc/… $\endgroup$ – Brian Borchers Nov 15 '17 at 5:59
  • $\begingroup$ @Wolfgang Bangerth, I have also tried the exactly double long complex functions for those functions used, like cexpl, atan2l, but it does not make much difference. And I think the problem is, why the result is far away from the double_limit, even using double data type. $\endgroup$ – Roy Liao Nov 15 '17 at 8:34
  • $\begingroup$ @BrianBorchers I have checked the article you referred. Unfortunately, the problem does not lies here. $\endgroup$ – Roy Liao Nov 15 '17 at 8:41
  • 2
    $\begingroup$ It doesn't make much of a difference, or no difference at all? You cannot expect much more accuracy from 80-bit vs the usual 64-bit double precision. $\endgroup$ – Wolfgang Bangerth Nov 15 '17 at 14:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.