FFTW library: verifying fourier transform derivative property

Question

I am trying to verify this relation with the fftw library:

Therefore, I chose f to be a Gaussian, computed the Fourier Transform of its derivative and compared it with the one of the Fourier Transform of the Gaussian multiplied by ik. Here is what I get:

It is very strange, especially because the plot of the Fourier Transform of the derivative of the Gaussian (i.e the red one) is not 0 at the origin, while it should be (I checked with the analytic one).

The code seems ok to me, anyhow here it is (I am using C):

int main() {

    int i, N = 100;
    double v[N], x[N], k[N/2+1], vd[N];
    double dx = 2*pi/N, dk=2*pi/(N*dx), tmp;
    fftw_complex *out;
    fftw_plan forward,inverse;

    out = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex )*( N/2 + 1 )); 


    forward = fftw_plan_dft_r2c_1d(N, v, out, FFTW_ESTIMATE);
    inverse = fftw_plan_dft_c2r_1d(N, out, vd, FFTW_ESTIMATE);

    //Initialise arrays 

    for( i = 0; i < N; i++ ) {   
            x[i] = dx*i;
            v[i] = -2*x[i]*exp( -pow( x[i], 2) );

            printf( " %le %le \n ", x[i], v[i] );
    }  


    for( i = 0; i < N; i++ ) {
       k[i]=i*dk;
    }

    k[N/2]=0.;

    //Compute fft
    fftw_execute( forward ); 

    //Print the results

    for( i = 0; i < N/2 + 1 ; i++ ) {
        printf( "%le %le %le \n", i*dk, out[i][0], out[i][1] );
    }

    //Multiply by ik 

    for( i = 0; i < ( N/2 + 1 ); i++ ) {

            tmp=out[i][0];
            out[i][0]=-k[i]*out[i][1];
            out[i][1]=k[i]*tmp;

            printf( "%le %le %le \n", i*dk, out[i][0], out[i][1] );

    }

    fftw_destroy_plan(forward);
    fftw_destroy_plan(inverse);
    fftw_free(out);

    return 0;  
}

Could anyone please tell me what I am doing wrong?

Please post the (relevant) code directly in your question, not links to images of code. — Oliver Charlesworth, Feb 15 '18 at 11:04
Note that whilst that relationship holds for the Fourier Transform, it isn't quite that straightforward for the Discrete Fourier Transform. — Oliver Charlesworth, Feb 15 '18 at 11:07
Yes, in principle you are right, but I checked (for example here: http://fourier.eng.hmc.edu/e101/lectures/handout3/node2.html ) and it should still be valid. — Krist, Feb 15 '18 at 11:13
Sorry, wrong link, I meant here: http://fourier.eng.hmc.edu/e101/lectures/handout3/node7.html — Krist, Feb 15 '18 at 11:18
One obvious bug is that `k` only has size `N/2+1`, yet you initialise it as if it had size `N` (`N` should probably `N/2` in your `for` loop.) — Paul R, Feb 15 '18 at 11:41

score 3 · Accepted Answer · answered Feb 16 '18 at 07:58

3

The discretized signal of the derived gaussian is supposed to be:

        x[i] = dx*i;
        v[i] = -2*x[i]*exp( -pow( x[i], 2) );

Nevertheless, the Discrete Fourier Transform corresponds to the Fourier Transform of periodic signals. Indeed, the underlining discretized function is written as an infinite weighted sum of sine waves.

Therefore, as DFT is applied, the discretized signal above corresponds to a periodized half of the derivative of a gaussian. Indeed, its average (zero frequency) is not zero, as all the values are negative.

To mimic the derivative of a gaussian (or any other signal of "finite" range ), the whole range of the signal must be covered. Hence, dx must be chosen such that dx*N>>sigma, where sigma is the standard deviation of the gaussian. And all the support of the function must be covered, including the positive side of the derivative.

Could you try something like:

        double sigma=dx*N*0.1;
        x[i] = dx*i-dx*(N/2);
        v[i] = -2*x[i]*exp( -pow( x[i]/sigma, 2) );

There must be a scaling left due to the value of the standard deviation.

The DFT can still be useful for non-periodic function, but the non-periodic functions are to be mapped to a periodic function by using a window. What is observed here is that applying a rectangular window, periodizing and then deriving is not the same as deriving, applying a rectangular window and finally periodizing. While all linear, these operators do not commute!

answered Feb 16 '18 at 07:58

francis

9,525
2
25
41

Setting `sigma=dx=1` should be fine, makes the computations easier. I get a nice DFT out of that. – Cris Luengo Feb 16 '18 at 08:06
I agree with you about the periodicity argument, indeed I should transform a function which has the same values at the extremes of the "box". The problem is that, as far as I understood, FFTW works in such a way that the first element of the output array is associated with the origin of frequencies, so whenever I transform a signal which is not starting from the origin, I get problems ( by the way if you knew ho to solve that it would be great). I tried with a gaussian which starts from the origin and is centered in the half of the box, but still those plots are very different. – Krist Feb 16 '18 at 10:10
@Mike : Indeed, this is the idea behind more refined window than the rectangular one. Most windows features values going to zero near the edges so as to recover a continuous periodized signal and avoid introducing an artificial discontinuity which pollutes the estimate of the derivative. In fact, some windows feature a null derivative at the edge so that even the first derivative of the periodized signal remain smooth. You can try to apply a Tukey window! – francis Feb 16 '18 at 10:18
@MikeM: You don't need a windowing function here. You need to make sure that the full Gaussian is in your `v[i]`, with the origin at `v[0]`, and with `v[i]` periodic: `v[i+size]` should be filled for `i` from `-1` to `-size/2`. This way you can fit your whole kernel in the signal, and still have the origin where it needs to be for the DFT. If `v` is large enough, the Gaussian will not introduce any discontinuities anywhere. – Cris Luengo Feb 18 '18 at 07:30

FFTW library: verifying fourier transform derivative property

1 Answers1