CODEWITHSUNDEEP
csignal-processingfftfftw
Krist
Krist

Reputation: 57

FFTW library: verifying fourier transform derivative property

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

enter image description here

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:

enter image description here

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?

Upvotes: 2

Views: 901

Answers (1)

francis
francis

Reputation: 9817

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!

Upvotes: 3

Related Questions