FFTW results differ from FFT in MATLAB

Question

I compare the forward FFT using FFTW and MATLAB fft. The input signal is a Gaussian. Code:

FFTW using C:

float *signal; /*input signal*/
int nt; /*length of the signal*/
fftw_complex *in, *out;
fftw_plan plan1;

in = fftw_malloc(nt*sizeof(fftw_complex)); 
out = fftw_malloc(nt*sizeof(fftw_complex));
for (j=0;j<nt;j++){
        in[j][0]=(double)signal[j];
        in[j][1]=0.0;
}    
plan1 = fftw_plan_dft_1d(nt, in, out, -1, FFTW_ESTIMATE);
fftw_execute(plan1);        
fftw_destroy_plan(plan1);

for (j=0;j<nt;j++){
        real[j]=(float)out[j][0];
        imag[j]=(float)out[j][1];
}

fft function in MATLAB:

fft(signal);

I plot the real and imaginary parts of both results:

fft vs fftw

The real part are almost the same value, while the imaginary part has quite different values. How to fix this problem?

The differences in Imag are in the range of 10 to the power of -15. So not really that different. — Kami Kaze
Another fact is, because my signal is real numbers, but in FFTW I create a complex DFT, which means the signal is complex with zero imaginary values. For real signal, out[i] is the conjugate of out[n-i] as matlab fft result shows. But imag of FFTW does not behave this. — clduan
A good test would be applying the inverse transform in both cases and comparing the results, this will show you if the differences in the imaginary part are relevant (I also think they aren't). — schnaader
"which means the signal is complex with zero imaginary values" The DFT has zero imaginary values because the input is symmetric, not because the input is real-valued. Also, knowing the imaginary values are supposed to be zero, you should have noticed that they are very close to zero. You are just looking at floating-point rounding errors here. Rounding errors do not need to be symmetric. — Cris Luengo
You could solve the "not exact conjugate" part by zeroing all the imaginary values smaller than an epsilon (e.g. 10^-14) as you know it's floating point noise anyway. — schnaader

Kami Kaze Kami Kaze · Accepted Answer · 2019-06-24T14:54:10

You should look at the scale factor of the plot on the left side over the plot of 'Imag'. It says 10^-15. This is quite small in relation to the real signal magnitude (at least the larger parts which is >10^1) so the results are quite similiar.

Floating point algorithms in general tend to not deliver the exact same result as long as they are not implemented exactly in the same way. (And even then they can differ by different options for rounding).

This QA might give some insight: Floating point inaccuracy examples

FFTW results differ from FFT in MATLAB

2 Answers