6. 開発情報

　NumRu/fftライブラリは,次のライブラリを利用する.

入出力ライブラリ

NArray <http://www.ir.isas.ac.jp/~masa/ruby/>

高速フーリエ変換ライブラリ

FFTW <http://www.fftw.org/>
ISPACK(予定) <http://www.gfd-dennou.org/arch/ispack/>

本ライブラリのインストールに際しては,これらのライブラリを先にインストールしておく必要がある.個々のライブラリのインストール方法については,それぞれの付属ドキュメントに譲る.

　FFTWを利用する際には,次の関数を用いている.

fftw_plan fftw_create_plan(int n, fftw_direction dir, int flags);
- fftを実行するための「plan」を作成する初期化関数.
void fftw(fftw_plan plan, int howmany, fftw_complex *in, int istride, int idist, fftw_complex *out, int ostride, int odist);
- 高速フーリエ変換を実行する関数.メモリ的に連続な配列でなくても実行する事ができる.
void fftw_destroy_plan(fftw_plan plan);
- planを解放するための関数.

6.3 ソースの構成

　ソースコードは,以下のような構成になっている.

・fft_fftw.c: メインソース
・na_fft.h: ヘッダファイル
・fft_narray.c: NArrayのサブセット
・extconf.rb: コンパイル用コンフィグスクリプト

6.4 ソースコード

6.4.1 fft_fftw.c

　以下に,fft_fftw.cソースコード全文を示す.

#include "na_fft.h"
#include <fftw.h>
#include <rfftw.h>
#ifdef FFTW_ENABLE_FLOAT
# define NA_FFTW_REAL NA_SFLOAT
# define NA_FFTW_COMP NA_SCOMPLEX
# define COMP_T scomplex
#else
# define NA_FFTW_REAL NA_DFLOAT
# define NA_FFTW_COMP NA_DCOMPLEX
# define COMP_T dcomplex
#endif

typedef enum {
  FFT_NO,
  FFT_YES
} DOING_FFT_DIM_FLAG;

/**** prototype declaration ****/
VALUE fft_fftw(int argc, VALUE *argv, VALUE self);
VALUE ffti_fftw(int argc, VALUE *argv, VALUE self);
static NA *fft_sub(int argc, VALUE *argv, VALUE self, fftw_direction dir);
static NA *fft_DCOMPLEX(int arg_c, int *arg_v, NA *na, fftw_direction dir);
static NA *fft_SCOMPLEX(int arg_c, int *arg_v, NA *na, fftw_direction dir);
static void check_argv(int argc, int *argv, int rank);
/**** end prototype ****/

VALUE 
fft_fftw(int argc, VALUE *argv, VALUE self)
{
  NA *na;
  na=ALLOC(NA);
  na=fft_sub(argc, argv, self, FFTW_FORWARD);
  return Data_Wrap_Struct(cNArray, 0, na_free, na);
}

VALUE 
ffti_fftw(int argc, VALUE *argv, VALUE self)
{
  NA *na;
  na=ALLOC(NA);
  na=fft_sub(argc, argv, self, FFTW_BACKWARD);
  return Data_Wrap_Struct(cNArray, 0, na_free, na);
}

NA *
fft_sub(int argc, VALUE *argv, VALUE self, fftw_direction dir)
{
  NA *na;
  int *arg_v, arg_c, *temp_argv;
  int n;

  /* functions used in this function */
  static NA *fft_DCOMPLEX();
  static NA *fft_SCOMPLEX();
  static void check_argv();

  GetNArray(self, na);

  if(argc){
    temp_argv=ALLOCA_N(int,argc);
    for(n=0 ; n<argc ; n++){
      temp_argv[n]=NUM2INT(argv[n]);
    }
  }

  /* 引数が一個以上のときチェック。以下の条件でエラーをレイズ。
     ・一個で、引数の中身が rank 以上のとき
     ・二個以上で、
       * 負の数を含むとき
       * 引数の中身が rank 以上のとき        */
  check_argv(argc, temp_argv, na->rank);

  /* VALUE argv[] => FFT_OP arg_v[rank] */

  arg_v=ALLOCA_N(int,na->rank);
  for(n=0 ; n<(na->rank) ; n++)
    arg_v[n] = FFT_NO ;    /* initializing int arg_v[rank] */

  if(argc==0){             /* 引数 なし */
    for(n=0 ; n<(na->rank) ; n++){
      arg_v[n] = FFT_YES ;
    }
  }else if(argc==1){
    if(*temp_argv<0){         /* 引数一個 で、負のとき */
      for(n=0 ; n<(na->rank) ; n++){
    arg_v[n] = FFT_YES;
      }
    }else{                  /* 引数一個 で、正のとき */
      arg_v[ *temp_argv ] = FFT_YES ;
    }  
  }
  else{                     /* 引数 二個以上 */
    for(n=0 ; n<argc ;n++)
      arg_v[ temp_argv[n] ] = FFT_YES ;
  }

  switch(na->type){
  case NA_DCOMPLEX:
    return fft_DCOMPLEX(arg_c, arg_v, na, dir);

  case NA_SCOMPLEX:
    return fft_SCOMPLEX(arg_c, arg_v, na, dir);

  default:
    rb_raise(rb_eTypeError, "operand is not valid type");
  }
}

NA *
fft_DCOMPLEX(int arg_c, int *arg_v, NA *na, fftw_direction dir)
{
  fftw_plan p;
  int n, m, rank,total, howmany, dim, length;
  
  dcomplex *ptr;
  char *new_ptr;

  int idist, istride;
  NA *ret_na;

  rank    = na->rank;
  total   = na->total;
  ptr     = (dcomplex *)na->ptr;

  new_ptr = ALLOC_N( char , sizeof(fftw_complex) * na->total ); 

  for(dim=0 ; dim<(na->rank) ; dim++){
    if(arg_v[dim]==FFT_YES){
      if(dim == 0){
        p = fftw_create_plan(na->shape[dim], dir, FFTW_ESTIMATE);
        if(p==NULL) 
            rb_raise(rb_eRuntimeError,"cannot allocate FFTW plan");

        idist   = na->shape[0];
        istride = 1;
        howmany = total / na->shape[0];

        fftw(p, howmany,
           (fftw_complex *)ptr    , istride, idist,
           (fftw_complex *)new_ptr, istride, idist);

        fftw_destroy_plan(p);
      }else if(dim == rank - 1){
        p = fftw_create_plan(na->shape[dim], dir, FFTW_ESTIMATE);
        if(p==NULL)
            rb_raise(rb_eRuntimeError,"cannot allocate FFTW plan");

        idist   = 1;
        istride = 
        howmany = total / na->shape[rank - 1];

        fftw(p, howmany,
           (fftw_complex *)ptr    , istride, idist,
           (fftw_complex *)new_ptr, istride, idist);
        fftw_destroy_plan(p);

      }else{
        idist   = 1;
        istride = 1; /* initialize */
        for(m=0 ; m<dim ; m++)
          istride = istride * na->shape[m];
        howmany = istride;
        length = istride * na->shape[dim];

        p = fftw_create_plan(na->shape[dim], dir, FFTW_ESTIMATE);
        if(p==NULL) 
            rb_raise(rb_eRuntimeError,"cannot allocate FFTW plan");

        for(m=0 ; m<total ; m += length){
          fftw(p, howmany,
             (fftw_complex *)&ptr[m]    , istride, idist,
             (fftw_complex *)&new_ptr[m], istride, idist);
        }
        fftw_destroy_plan(p);
      }
    }
  }

  ret_na = na_alloc_struct(NA_FFTW_COMP, rank, na->shape);
  ret_na->ptr = new_ptr;

  return ret_na;
}

NA *
fft_SCOMPLEX(int arg_c, int *arg_v, NA *na, fftw_direction dir)
{
  int n;
  NA *dcmp_na;
  dcomplex *dcmp_ptr;
  scomplex *scmp_ptr;

  dcmp_na = na_alloc_struct(NA_DCOMPLEX, na->rank, na->shape);
  dcmp_ptr = ALLOC_N(dcomplex ,na->total);
  scmp_ptr = (scomplex *)na->ptr;

  for(n=0 ; n<(na->total) ; n++){
    dcmp_ptr[n].r = scmp_ptr[n].r;
    dcmp_ptr[n].i = scmp_ptr[n].i;
  }

  dcmp_na->ptr = (char *)dcmp_ptr;
  return fft_DCOMPLEX(arg_c, arg_v, dcmp_na, dir);
}

static void 
check_argv(int argc, int *argv, int rank)
{
  int n;

  if(argc>rank) 
    rb_raise(rb_eArgError, "too many arguments");

  if(argc==1){
    if( (*argv)>(rank-1) )
      rb_raise(rb_eArgError, "invalid value");
  }else if(argc>1){
    for(n=0 ; n<argc ; n++){
      if( argv[n]>(rank-1) || argv[n]<0 )
    rb_raise(rb_eArgError, "invalid value");
    }
  }
  return;
}

void 
Init_fft(void)
{
  rb_define_method(cNArray, "fft"     ,fft_fftw    , -1);
  rb_define_method(cNArray, "ffti"    ,ffti_fftw   , -1);
/*rb_define_method(cNArray, "rfft"    ,rfft_fftw   , -1); */
/*rb_define_method(cNArray, "rffti"   ,rffti_fftw  , -1); */
}

6.4.2 na_fft.h

　以下に,na_fft.h全文を示す.　

#include <ruby.h>

enum NArray_Types {
  NA_NONE,
  NA_BYTE,      /* 1 */
  NA_SINT,      /* 2 */
  NA_LINT,      /* 3 */
  NA_SFLOAT,    /* 4 */
  NA_DFLOAT,    /* 5 */
  NA_SCOMPLEX,  /* 6 */
  NA_DCOMPLEX,  /* 7 */
  NA_ROBJ,      /* 8 */
  NA_NTYPES     /* 9 */
};

struct NARRAY {
  int    rank;    /* # of dimension */
  int    total;   /* # of total element */
  int    type;    /* data type */
  int   *shape;
  char  *ptr;     /* pointer to data */
  VALUE  ref;     /* NArray object wrapping this structure */
};

typedef struct { float  r,i; } scomplex;
typedef struct { double r,i; } dcomplex;
typedef struct NARRAY NA;

/* fft_fftw.c */
VALUE cNArray;
VALUE fft_fftw(int argc, VALUE *argv, VALUE self);
VALUE ffti_fftw(int argc, VALUE *argv, VALUE self);
void Init_fft(void);

/* fft_narray.c */
struct NARRAY *na_alloc_struct(int type, int rank, int *shape);
void na_free(struct NARRAY *ary);

/* copy from narray.h */

#define GetNArray(obj,var)  Data_Get_Struct(obj, struct NARRAY, var)

6.5 現状,および今後の予定

　本ライブラリは,未だ開発の途上であり,さらなる拡張の余地を残している.現在のバージョンでは,FFTWライブラリの複素フーリエ正変換および逆変換に対応しているのみで,実数正変換,および逆変換には対応しておらず,ISPACKについては完全に未実装である.

　今後の予定としては,以下のようになっている.

ISPACKライブラリを利用できるようにし,すべてのメソッドを利用可能にする.
完全に利用できるようになったあかつきには,NArray パッケージに patch を送り本ライブラリを融合する.
細かなバグフィックスとパッケージの維持を心がける.

　次回リリースでは,ISPACKによるサイン変換,コサイン変換を組みこむこと,NArrayライブラリに対するパッチの形式でリリースすることの二点を目指して開発を続けて行きたい.

6. 開発技術文書

6.1 開発環境,動作確認

6.2 必要な環境

6.3 ソースの構成

6.4 ソースコード

6.4.1 fft_fftw.c

6.4.2 na_fft.h

6.5 現状,および今後の予定