# High precision Discrete Fourier Transform in c

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. 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.

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. 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

• 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. – Kirill Nov 14 '17 at 16:36
• 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. – Kirill Nov 14 '17 at 20:04
• @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. – Kirill Nov 14 '17 at 23:37
• @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. – Federico Poloni Nov 15 '17 at 8:49
• @RoyLiao If you need 300 digits of accuracy, your best bet is using a multiprecision library such as gmplib.org. – Federico Poloni Nov 15 '17 at 13:23

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.
• @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. – Roy Liao Nov 15 '17 at 8:34