CODEWITHSUNDEEP
cfftwderivative
Guest
Guest

Reputation: 11

Taking a second derivative in FFTW3

I have tested my code for some real functions using a forward FFT and IFFT (normalized the result), this works fine.

I would, however, like to take a second derivative of a real function. For simplicity sake, I take sin(2*pi*t) as a test case. Here is the relevant code I use (FFT functions in a library):

int main(void)
{
    int i;
    int nyh = (N/2) + 1;

    double result_array[nyh][2];

    double x_k[nyh][2];
    double x_r[N];

    FILE* psit;
    psit=fopen("psitest.txt","w");  

    init();

    fft(x, result_array);  //function in a library, this has been tested
    psi(result_array, x_k);  
    ifft(x_k, x_r);        //function in a library, this has been tested

    for(i=0;i<N;i++)
    {
        fprintf(psit, "%f\n", x_r[i]);
    }


    fclose(psit);

    return 0;
}

void psi(double array[nyh][2], double out[nyh][2])
{
    int i;

    for ( i = 0; i < N/2; i++ )
    {
        out[i][0] = -4.0*pi*pi*i*i*array[i][0];
        out[i][1] = -4.0*pi*pi*i*i*array[i][1];
    }   
    out[N/2][0]=0.0;
    out[N/2][1]=0.0;
}

void init()
{
    int i;

    for(i=0;i<N;i++) 
    {
        x[i] = sin(2.0*pi*i/N);
    }
}

Now here is the problem: This algorithm works perfectly for any function of the form sin( 2*pi*t*K), where K is an integer, but if I take as a test function sin(3*pi*t), the algorithm fails. I am not able to see the mistake in my coding.

Please note that because the function is real, I only need to take half of the k values. This is not the problem.

Upvotes: 1

Views: 751

Answers (2)

masaishi
masaishi

Reputation: 21

I don't know if you have fixed this problem... But I guess the major problem is that sin(3 Pi t) is not periodic in the domain [0,1](sin(0) != sin (3 * Pi)).

FFT could not work properly...

Upvotes: 0

Paul R
Paul R

Reputation: 213080

My guess is that sin(3*pi*t) introduces a discontinuity, since it does not give an integer number of cycles in your sample interval. For most FFT-related applications you would apply a window function to deal with such discontinuities, but obviously this will introduce an error term into your derivative and I'm not sure whether you will be able to correct for this.

Upvotes: 1

Related Questions