FFTのひとつの例

＜FFTのひとつの例＞

AD変換して求めた入力信号の周波数特性を求めるため、２のｎ乗のポイント数を計算する高速フーリエ変換をC++で記述したサンプルプログラムを参考に下記に示す。
FFTを計算するにあたり、ハ二ングの窓を掛けて、滑らかに両端の値をゼロにして、入力信号を見かけ上有限の長さに押し込めている。

このサンプルを実際に動作させるためには、おそらく、
ご使用されるプログラミング環境に合わせて修正・変更・追加などが必要です。
万一、このサンプルを動作させる場合は、あなたの責任においておこなってください。

#include <math.h>

//--------------------------------------------------------
TCFFT::TCFFT(int nzyou)
// データのポイント数が2のnzyou乗である
// FFTを計算するたの準備を行う。
//
//　　nzyou=8の時は　計算するポイント数は 256
//    nzyou=9    　　　　　　　　　　　　 512
//    nzyou=10   　　　　　　　　　　　　 1024
//    nzyou=11   　　　　　　　　　　　　 2048
{


    value_nzyou= nzyou;
    // データのポイント数をvalue_pointsに格納する。
    value_points=f_zyou(value_nzyou);

    // 準備1　
    // 繰り返して使用するものを、
    // あらかじめ、計算して値を求めておく。
    f_mkdat();
    f_bitrev();

    // 準備2　窓の準備。窓はFFTで必修ではないが、実用上、使用する。
    f_set_hanig();

    // 作業領域を確保
    works=(double *)malloc(sizeof(double) * value_points * 2);
    datas=(double *)malloc(sizeof(double) * value_points * 2);

}
//--------------------------------------------------------
void TCFFT::f_fft()
{
// 入出力は共通で下記のような複素数とする。
// data[2*i ] 実数部
// data[2*i+1]　虚数部　＊実数の場合はゼロを設定する。
// iは、0,1,...,value_points-1 の整数である。

   int n,nh,i,i1,i2;
   double *ws;

   // ポイント数
   nh=value_points;

   // sin,cos のテーブルをセット
   ws=mk;

   // 計算
   for(i1=0;i1<value_nzyou;++i1)
   {
         nh= nh/2;
         n= f_zyou(i1);

         for(i2=0;i2<n;++i2)
         {
            f_fftsub(nh,2*nh*i2,n,datas,ws);
         }
   }

   // 入れ替え

   for(i=0;i<value_points;++i)
   {
     i1= table[i];
     works[2*i1] = datas[2*i];
     works[2*i1+1]= datas[2*i+1];
   }
   for(i=0;i<value_points;++i)
   {
     datas[2*i]=works[2*i];
     datas[2*i+1]=works[2*i+1];
   }


}
//--------------------------------------------------------
int TCFFT::f_zyou( int nzyou)
// ２のn乗を出力する
{
   return( 1 << nzyou);
}
//---------------------------------------------------------
void TCFFT::f_mkdat( )
// sin,cosの値を事前に計算しておく
{
   int i;
   double x;
   double *p;

   mk = (double *)malloc( sizeof(double) * value_points *2 );
   p=mk;

   for(i=0;i<value_points;++i)
   {
       x= M_PI * 2.0 *(double)i / (double)value_points;
       *p= cos(x);
       ++p;
       *p= -sin(x);
       ++p;
   }
}
//---------------------------------------------------------
void TCFFT::f_bitrev()
//   bit-reversed ビット　リバースド　入れ替え順番の設定
{
   int i,j,x1,x;
   int out;
   int *p;

   table=(int *)malloc( sizeof(int) * value_points);
   p=table;

   for(i=0;i<value_points;++i)
   {
     out=0;
     x1=i;
     for(j=value_nzyou;j>0;--j)
     {
        x = x1 & 0x0001;
        out = out + x*f_zyou(j-1);
        x1>>=1;
     }
     *p=out;
     ++p;
   }
}
//---------------------------------------------------------
void TCFFT::f_multc( double *x, double *ws, double *z)
// 複素数の掛け算
{
   z[0]     = x[0] * ws[0] - x[1] * ws[1];
   z[1]     = x[0] * ws[1] + x[1] * ws[0];
}
//---------------------------------------------------------
void TCFFT::f_fftsub(int nh, int offset,int nstep, double datas[], double ws[])
// 小さい分割部分の計算
{
   double a[4];
   int i;
   int x1,y1;
   int x2,y2;

   for(i=0;i<nh;++i)
   {
      x1=2*(i+offset);
      y1=x1+1;
      x2=2*(i+nh+offset);
      y2=x2+1;

      a[0]= datas[x1];
      a[1]= datas[y1];
      datas[x1] = a[0] + datas[x2];
      datas[y1]= a[1] + datas[y2];
      a[0]= a[0] - datas[x2];
      a[1]= a[1] - datas[y2];
      f_multc( a, (ws+2*i*nstep), (datas + x2));
   }
}
//--------------------------------------------------------
//
// 以下は、AD変換して求めた入力信号のデータに窓を掛けて
// 値を複素数に拡張する。FFTを求める前の準備用。
//
//--------------------------------------------------------
void TCFFT::f_dohan(short *wave, int mode)
//   入力信号の16ビット整数データ(wave)にハニングの窓を掛けて、
//   複素数(datas)に変換する。
//
//   mode が HANNING 時はハニングの窓を掛ける。
//   mode が HANNING 以外ならば、窓を掛けない。
{
   int i;

   for(i=0;i<value_points;++i)
   {
      if( mode == HANNING)
      {
        datas[2*i] = (double)wave[i] * hng_data[i];   // 実数部
      }
      else
      {
        datas[2*i] = (double)wave[i];
      }
      datas[2*i+1] = 0.0;                   // 虚数部
   }
}
//--------------------------------------------------------
void TCFFT::f_set_hanig()
// ハニングの窓を計算する
{
    int i,z;
    double y,y2;
    double *p;

    hng_data=(double *)malloc( sizeof(double) * value_points);
    p=hng_data;

    y2=(double)value_points;
    for(i=0;i<value_points;++i)
    {
      y = (double)i/y2 ;
      y = y * 2.0 * M_PI;
      *p= 0.5 - 0.5 * cos(y);
      ++p;
    }

}

//--------------------------------------------------------
double *TCFFT::return_datas_pointer()
{
   return( datas);
}
//--------------------------------------------------------
int TCFFT::return_value_points()
{
   return(value_points);
}
//--------------------------------------------------------
TCFFT::~TCFFT()
{
     free(mk);
     free(table);
     free(hng_data);
     free(works);
     free(datas);
}
//--------------------------------------------------------
//
//
//--------------------------------------------------------

#define HANNING 1

class TCFFT
{
public:
        void f_dohan(short *wave, int mode);
        void f_fft();

        //　クラスTCFFTの内部の値を返す関数
        // FFTの計算結果が格納されるdatasを返す。
        double *return_datas_pointer();
        // FFT計算のためのデータポイント数を返す。
        int return_value_points();

        double *datas;

        TCFFT(int nzyou);
      virtual ~TCFFT();

private:
        int f_zyou(int nzyou);
        void f_multc( double *x, double *ws, double *z);
        void f_fftsub(int nh, int offset,int nstep, double datas[], double ws[]);
        void f_mkdat();
        void f_bitrev();
        void f_set_hanig();

        double *mk;
        int    *table;
        double *hng_data;
        double *works;

        int    value_nzyou;
        int    value_points;
};

No.2c 2006年1月12日