2011-06-23 11 views
1

私はフォワードFFTとIFFT(結果を正規化)を使っていくつかの実際の関数をテストしましたが、これはうまくいきます。FFTW3で二次微分をとる

しかし、私は実際の関数の2次導関数をとることが好きです。簡単にするために、私はsin(2 * pi * t)をテストケースとして扱います。ここで私は(ライブラリ内のFFT機能)を使用し、関連するコードは次のとおりです。

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

は、今ここで問題です:このアルゴリズムはどこ、フォームの罪(2 *パイ*トンの*のK)のいずれかの機能のために完璧に動作しますKは整数ですが、テスト関数sin(3 * pi * t)として取ると、アルゴリズムは失敗します。私はコーディングの間違いを見ることができません。

関数が実際であるため、私はk値の半分を取る必要があることに注意してください。これは問題ではありません。

ありがとうございます。

答えて

1

私の推測では、sin(3*pi*t)はサンプル間隔に整数のサイクル数を与えないため、不連続性が導入されていると思います。ほとんどのFFT関連アプリケーションでは、このような不連続性に対処するためにウィンドウ関数を適用しますが、明らかにこれは微分に誤差項を導入し、これを訂正できるかどうかはわかりません。

0

あなたがこの問題を解決したかどうかわかりません...しかし、大きな問題は、sin(3 Pi)がドメイン[0,1](sin(0)!= sin(3 * Pi))。

FFTが正しく機能しない...