GSL (GNU Scientific Library) を用いた FFT

GSL の FFT に関する機能について，一部分を紹介する．

GNU Scientific Library -- Reference Manual の 15 章「Fast Fourier Transforms (FFTs)」

前準備

GSL のインストール

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

実数列，半複素数列の FFT

【要点】

混合基数法は，任意のデータ長に対して使うことができる（長さが2の累乗 (radix 2)でなくても構わないということ）． GSL に実装されている混合基数法は，パウル・シュヴァルツラウバーのFFTPACK ライブラリを実装しなおしたものである．
実数列に FFT を行うと，半複素数列が得られる．この関数が gsl_fft_real();
gsl_fft_real_transform の出力である半複素数列を，一般的な複素数配列に変換するには，次の関数を使う
```
int gsl_fft_halfcomplex_unpack (const double halfcomplex_coefficient [],
gsl_complex_packed_array complex_coefficient, size_t stride, size_t n)
```
データ長 N = 5 半複素数列
半複素数列に逆 FFT を行うと実数列が得られる．この関数が gsl_fft_halfcomplex();

下の2行は必須

#include<gsl/gsl_fft_real.h>
#include<gsl/gsl_fft_halfcomplex.h>

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

プログラムの準備

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

ビルドして実行

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