FFTW - Difference between FFTW_REDFT00 and FFTW_DHT flags

Question

I am interested in doing 1D FFT with FFTW on real data.

For this, I am using a cosinus signal with a frequency equal to 10 Hz and a sampling frequency of sizex*frequency_signal with sizex the number of sampling points.

With the flag (FFTW_DHT) into fftw_plan_r2r_1d(sizex, Array, Array, FFTW_DHT, FFTW_ESTIMATE);, I get the dirac impulsion at f=10 Hz, here are the output (column 1: k*f_sampling/size x and column 2: X_k) after the forward fft :

0.000000 7.304123e-14
10.000000 5.000000e+01
20.000000 -2.227743e-14
30.000000 -1.300521e-14
40.000000 -3.774757e-15
50.000000 -2.989904e-15
60.000000 -4.879698e-15
70.000000 -2.838093e-15
80.000000 -5.479074e-16
90.000000 1.605429e-15
100.000000 -1.491050e-15
110.000000 -2.587601e-16

...

But with the FFTW_REDFT00, I can't get to have the dirac impulsion at f=10 Hz. In this case, I have the following output :

0.000000 -1.998027e+00
10.000000 2.682414e+00
20.000000 9.843837e+01
30.000000 -1.543229e+00
40.000000 6.493255e-01
50.000000 -3.723752e-01
60.000000 2.449150e-01
70.000000 -1.744771e-01
80.000000 1.310807e-01
90.000000 -1.023168e-01
100.000000 8.221456e-02
110.000000 -6.758738e-02
...

Could I get the dirac at f=10 Hz with FFTW_REDFT00 flag ?

What's exactly the difference between these two flags, i.e how can I find the same results of FFTW_DHT with FFTW_REDFT00 flag.

From fftw DFT doc, I thought that these 2 two flags produced the same results but this is not the case apparently.

I would like just to switch from one to another. if I know how to reverse them, it could help me for a code which uses FFTW_REDFT00 flag.

Answer 1

First, note that in http://www.fftw.org/fftw3_doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html, the REDFT00 misses a factor of 2 before π in the argument of cosine. That's why you saw a peak around 20Hz instead of 10Hz.

Second, REDFT00 is particularly tricky because it requires you to allocate one additional element. That is, Array should contain sizex + 1 elements, and you should create the plan as

fftw_plan_r2r_1d(sizex + 1, Array, Array, FFTW_REDFT00, FFTW_ESTIMATE);

Without the additional 1, the peak is widened as you have seen.

To avoid the peak widening. Follow these rules. If you prepare the signal as

for (i = 0; i <= n; i++) a[i] = cos(2*M_PI*10*i/n);

Then detect it by:

fftw_plan_r2r_1d(n + 1, a, a, FFTW_REDFT00, FFTW_ESTIMATE);

If you prepare the signal as

for (i = 0; i < n; i++) a[i] = cos(2*M_PI*10*(i+0.5)/n);

Then detect it by:

fftw_plan_r2r_1d(n, a, a, FFTW_REDFT10, FFTW_ESTIMATE);

If you prepare the signal as

for (i = 0; i < n; i++) a[i] = cos(2*M_PI*10.25*i/n);

Then detect it by:

fftw_plan_r2r_1d(n, a, a, FFTW_REDFT01, FFTW_ESTIMATE);

If you prepare the signal as

for (i = 0; i < n; i++) a[i] = cos(2*M_PI*10.25*(i + .5)/n);

Then detect it by:

fftw_plan_r2r_1d(n, a, a, FFTW_REDFT11, FFTW_ESTIMATE);