4

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?

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
clduan
  • 41
  • 1
  • 4
    The differences in Imag are in the range of 10 to the power of -15. So not really that different. – Kami Kaze Jun 24 '19 at 14:50
  • 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 Jun 24 '19 at 15:35
  • 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 Jun 24 '19 at 16:16
  • 1
    "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 Jun 24 '19 at 16:25
  • I can also see that you created your Gaussian signal wrong, which is why your DFT doesn't look like a Gaussian. The spatial domain has the origin in the first element as well, you need to circularly shift your signal so that the center of the Gaussian is at the first sample, with the left side of the Gaussian coming in at the right end of the signal. – Cris Luengo Jun 24 '19 at 16:26
  • Yes, if I do inverse transform in both cases, the results are the same. But if I multiply a filter to the Fourier series and inverse transform, there will be slightly different due to not exact conjugate of the imaginary part. – clduan Jun 24 '19 at 16:27
  • 1
    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 Jun 24 '19 at 16:47
  • @schnaader This is a good way. Thx. – clduan Jun 24 '19 at 17:00

2 Answers2

2

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

Kami Kaze
  • 2,069
  • 15
  • 27
  • "Quite small" in such a case can usefully be judged only relative to something else. The absolute magnitude is not sufficient. In this case, however, any of several characteristics of the real signal (maximum, average, absolute range) provides an appropriate basis for comparison for concluding that yes, the imaginary signal is negligible in this case. – John Bollinger Jun 24 '19 at 15:01
  • @JohnBollinger Totally correct, I just infered this in regard to the real part of the spektrum. Made a chance to make this clearer. – Kami Kaze Jun 24 '19 at 15:07
  • 1
    *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.* Not only that, but FFTW can produce different results depending on the exact solution (FFTW's "plan") method chosen. For example, computing the DFT of a 40-element data array will come out slightly different if it's split into 4x10 or 5x8 subelements. Run it on different computers or just happen to get a different memory layout, and you can get different results. To get consistent* results from FFTW requires saving the plans and reusing them. – Andrew Henle Jun 24 '19 at 15:14
  • * - consistent as in "the exact same inputs always produce the exact same outputs down to the last bit". – Andrew Henle Jun 24 '19 at 15:15
  • @AndrewHenle This is good to know I never worked with FFTW, so I can not touch on this matter, but this generally falls under what I said about different algorithm implementation as FFTW seems to have a few to chose from.This somebody might not expect this kind of behaviour, so thanks for the addition. – Kami Kaze Jun 24 '19 at 15:16
0

Rounded to the nearest 0.001% of full scale (real), notice that the imaginary values are all zero.

hotpaw2
  • 70,107
  • 14
  • 90
  • 153