Фильтр частотной области в FFTW - PullRequest
Новая
       48

Фильтр частотной области в FFTW

0 голосов
/ 12 мая 2019

Я хочу применить полупроизводный фильтр sqrt(d/dt) к моему сигналу. Мы знаем, что это трудно реализовать во временной области, но это просто умножить на sqrt(iw) (w - угловая частота) в частотной области. И сделать IFFT, чтобы вернуть его во временную область. Но я не могу получить правильный ответ. 

    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];
        }            
    }

Я использую простую функцию синуса, чтобы проверить это. На выходе все нули. Кто-нибудь может понять это?

...