GSL (GNU Scientific Library) を用いた FFT

前準備

GSL のインストール

Windows での GSL のインストール（ソースコードを使用）（MSYS2，configure，make を利用）: 別ページ »で説明

実数列，半複素数列の FFT

1. MSYS2 MINGW64 を実行する

   スタートメニューで，「MSYS2」の下の「MSYS2 MINGW64」を選ぶ．

2. プログラムの準備

   #include<stdio.h>
   #include<math.h>
   #include<gsl/gsl_fft_real.h>
   #include<gsl/gsl_fft_halfcomplex.h>
        
   #define N 18
   int main (int argc, char **argv)
   {
     int i; 
     double data[N];
        
     gsl_fft_real_wavetable * real;
     gsl_fft_halfcomplex_wavetable * hc;
     gsl_fft_real_workspace * w;
        
     data[0]  = 0.0; 
     data[1]  = 0.0; 
     data[2]  = 0.0; 
     data[3]  = 0.0; 
     data[4]  = 0.0; 
     data[5]  = 0.0; 
     data[6]  = 1.0; 
     data[7]  = 1.0; 
     data[8]  = 1.0; 
     data[9]  = 1.0; 
     data[10] = 1.0; 
     data[11] = 1.0; 
     data[12] = 0.0; 
     data[13] = 0.0; 
     data[14] = 0.0; 
     data[15] = 0.0; 
     data[16] = 0.0; 
     data[17] = 0.0; 
   
     /* 元データの表示 */
     for (i = 0; i < N; i++) {
         printf ("%d: %e\n", i, data[i]);
     }
        
     w = gsl_fft_real_workspace_alloc(N);
     real = gsl_fft_real_wavetable_alloc(N);
     gsl_fft_real_transform (data, /* stride */ 1, N, real, w);
     gsl_fft_real_wavetable_free (real);
        
     /* 変換結果の表示 */
     for (i = 0; i < N; i++) {
         printf ("%d: %e\n", i, data[i]);
     }
   
     hc = gsl_fft_halfcomplex_wavetable_alloc(N);
     gsl_fft_halfcomplex_inverse (data, 1, N, hc, w);
     gsl_fft_halfcomplex_wavetable_free(hc);
        
     /* 逆変換結果の表示 */
     for (i = 0; i < N; i++) { 
       printf ("%d: %e\n", i, data[i]);
     }
        
     gsl_fft_real_workspace_free (w);
   
     return 0;
   }
   

3. ビルドして実行

   gcc -I"C:\gsl\include" -c -o a.o a.c
   gcc -L"C:\gsl\lib" -o a.exe a.o -lgsl -lgslcblas
   ./a.exe
   

   

機能まとめ

• 実数，半複素数データ FFT用の作業領域の確保（長さが2の累乗 (radix 2)でない場合）

  長さ N の実数データの FFT のための作業領域の確保

  gsl_fft_real_workspace *w = gsl_fft_real_workspace_alloc(N);
  

• 実数データ FFT用の波形表の確保（長さが2の累乗 (radix 2)でない場合）

  gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(N);
  

• 半複素数 (half-complex) データ FFT用の波形表の確保（長さが2の累乗 (radix 2)でない場合）

  gsl_fft_halfcomplex_wavetable *hc = = gsl_fft_halfcomplex_wavetable_alloc(N);
  

• 実数データの FFT（長さが2の累乗 (radix 2)でない場合）

  作業領域 w, 波形表 real，波形データ data，長さ N，とび数(stride) が１．

  gsl_fft_real_transform(data, 1, N, real, w);
  

• 半複素数 (half-complex) データの FFT（長さが2の累乗 (radix 2)でない場合）

  作業領域 w, 波形表 hc，波形データ data，長さ N，とび数(stride) が１．

  gsl_fft_halfcomplex_inverse(data, 1, N, hc, w);
  

• 実数データ FFT用の波形表の解放（長さが2の累乗 (radix 2)でない場合）

  gsl_fft_real_wavetable_free(r);
  

• 半複素数 (half-complex) データ FFT用の波形表の解放（長さが2の累乗 (radix 2)でない場合）

  gsl_fft_halfcomplex_wavetable_free(hc);
  

• 実数，半複素数データ FFT用の作業領域の解放（長さが2の累乗 (radix 2)でない場合）

  gsl_fft_real_workspace_free(w);