Frequency domain filter in FFTW

Asked by At

I want to apply a half-derivative filter sqrt(d/dt) to my signal. We know it's difficult to implement in time-domain, but it's just multiply by sqrt(iw) (w is the angular frequency) in frequency domain. And do IFFT to get it back to time domain. But I cannot get the correct answer.

    in = (fftw_complex *) fftw_malloc(nt*sizeof(fftw_complex));  
    out = (fftw_complex *) fftw_malloc(nt*sizeof(fftw_complex)); 
    out_ps = (fftw_complex *) fftw_malloc(nt*sizeof(fftw_complex)); 

    for(k=0;k<ng;k++){
        /* FORWARD */
        for (j=0;j<nt;j++){
            in[j][0]=dat[k][j];
            in[j][1]=0.0;
        }

        plan1 = fftw_plan_dft_1d(nt, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

        fftw_execute(plan1);
        fftw_destroy_plan(plan1);

        /* 45 Phase shift in frequency domain */
        for (j=0;j<nt;j++){
               w = j/nt;
               out_ps[j][0] = sqrt(w/2.0)*(out[j][0]-out[j][1]);
               out_ps[j][1] = sqrt(w/2.0)*(out[j][0]+out[j][1]);
        }

        /* INVERSE */
        plan2 = fftw_plan_dft_1d(nt, out_ps, in, FFTW_BACKWARD, FFTW_ESTIMATE);
        fftw_execute(plan2);
        fftw_destroy_plan(plan2);
        for (j=0;j<nt;j++){
                 inv_dat[k][j] = in[j][0];
        }            
    }

I use a simple sine function to test it. The output are all zeros. Can anyvody figure it out?

0 Answers