2010-12-12 11 views
5

FFTで何らかのフィルタリングをしようとしています。私はr2r_1dプランを使用していますが、逆変換を行う方法がわかりません。FFTWライブラリで実数FFTを逆実数化する方法

void PerformFiltering(double* data, int n) 
    { 
        /* FFT */ 
     double* spectrum = new double[n]; 

     fftw_plan plan; 

     plan = fftw_plan_r2r_1d(n, data, spectrum, FFTW_REDFT00, FFTW_ESTIMATE); 

     fftw_execute(plan); // signal to spectrum 
     fftw_destroy_plan(plan); 


        /* some filtering here */ 


        /* Inverse FFT */ 
     plan = fftw_plan_r2r_1d(n, spectrum, data, FFTW_REDFT00, FFTW_ESTIMATE); 
     fftw_execute(plan); // spectrum to signal (inverse FFT) 
     fftw_destroy_plan(plan); 

} 

私はすべてのことを正しくやっていますか?私は混乱しているので、FFTWの複雑なDFTの場合、次のようなフラグで変換方向を設定できます。
p = fftw_plan_dft_1d(N、in、out、FFTW_FORWARD、FFTW_ESTIMATE);
または
p = fftw_plan_dft_1d(N、in、out、FFTW_BACKWARD、FFTW_ESTIMATE);

答えて

1

あなたは正しくそれをしました。 FFTW_REDFT00はコサイン変換を意味し、それ自体の逆数です。したがって、「前方」と「後方」を区別する必要はありません。ただし、配列のサイズには注意してください。頻度が10で、データに意味のある点が100個含まれている場合、配列dataはのデータポイントを保持し、n = 101を100の代わりに設定する必要があります。正規化は2*(n-1)である必要があります。以下の例を参考にしてコンパイルしてください。gcc a.c -lfftw3

#include <stdio.h> 
#include <math.h> 
#include <fftw3.h> 
#define n 101 /* note the 1 */ 
int main(void) { 
    double in[n], in2[n], out[n]; 
    fftw_plan p, q; 
    int i; 
    p = fftw_plan_r2r_1d(n, in, out, FFTW_REDFT00, FFTW_ESTIMATE); 
    for (i = 0; i < n; i++) in[i] = cos(2*M_PI*10*i/(n - 1)); /* n - 1 instead of n */ 
    fftw_execute(p); 
    q = fftw_plan_r2r_1d(n, out, in2, FFTW_REDFT00, FFTW_ESTIMATE); 
    fftw_execute(q); 
    for (i = 0; i < n; i++) 
    printf("%3d %9.5f %9.5f\n", i, in[i], in2[i]/(2*(n - 1))); /* n - 1 instead of n */ 
    fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup(); 
    return 0; 
}