fftReal.cpp

00001 //------------------------------------------------------------------------------
00002 //  データが実数の場合の FFT class
00003 //      FFT     複素共役になる部分は計算しない
00004 //      IFFT    複素共役の部分を除いた値から計算する
00005 //
00006 //  2020/12/12, Copyright (c) 2020 MIKAMI, Naoki
00007 //------------------------------------------------------------------------------
00008 
00009 #include "fftReal.hpp"
00010 
00011 namespace Mikami
00012 {
00013     // コンストラクタ
00014     FftReal::FftReal(int16_t n)
00015             : N_FFT_(n), N_INV_(1.0f/n), wTable_(n/2), bTable_(n), u_(n)
00016     {
00017         // __clz(): リーディング 0 の数を取得する命令に対応
00018         uint32_t shifted = n << (__clz(n) + 1);
00019         MBED_ASSERT(shifted == 0);  // n が 2 のべき乗であることをチェック
00020 
00021         // 回転子の表を作成
00022         Complex arg = Complex(0, -6.283185f/N_FFT_);
00023         for (int k=0; k<N_FFT_/2; k++)
00024             wTable_[k] = exp(arg*(float)k);
00025 
00026         // ビット逆順のための表を作成
00027         uint16_t nShift = __clz(n) + 1;
00028         for (int k=0; k<n; k++)
00029             // __rbit(k): ビットの並びを逆にする命令に対応
00030             bTable_[k] = __rbit(k) >> nShift;
00031     }
00032 
00033     // FFT の実行
00034     //      有効な部分 y[0] ～ y[N_FFT_/2]
00035     //      それ以外の部分は複素共役になるので，計算をしていない
00036     void FftReal::Execute(const float x[], Complex y[])
00037     {
00038         for (int n=0; n<N_FFT_; n++) u_[n] = x[n];
00039 
00040         // 最終ステージを除いた部分
00041         ExcludeLastStage();
00042 
00043         // 最終ステージ
00044         y[0] = u_[0] + u_[1];
00045         y[N_FFT_/2] = u_[0] - u_[1];
00046         for (int k=2; k<N_FFT_; k+=2) u_[k] = u_[k] + u_[k+1];
00047 
00048         // ビット逆順の並べ替え
00049         for (int k=1; k<N_FFT_/2; k++) y[k] = u_[bTable_[k]];
00050     }
00051 
00052     // IFFT の実行
00053     //      このクラスで計算された FFT の結果の IFFT を計算する
00054     //      実行結果    x[0] ～ x[N_FFT_-1]
00055     void FftReal::ExecuteIfft(const Complex y[], float x[])
00056     {
00057         int half = N_FFT_/2;
00058 
00059         for (int n=0; n<=half; n++) u_[n] = y[n];
00060         for (int n=half+1; n<N_FFT_; n++)
00061             u_[n] = conj(y[N_FFT_-n]);      // 後半は複素共役になっているものとする
00062 
00063         // 最終ステージを除いた部分
00064         ExcludeLastStage();
00065 
00066         // 最終ステージとビット逆順の並べ替え処理
00067         x[0] = N_INV_*(u_[0].real() + u_[1].real());
00068         x[half] = N_INV_*(u_[0].real() - u_[1].real());
00069 
00070         for (int n=2; n<N_FFT_; n+=2)
00071         {
00072             float un  = u_[n].real();
00073             float un1 = u_[n+1].real();
00074             x[Index(n)]   = N_INV_*(un + un1);
00075             x[Index(n+1)] = N_INV_*(un - un1);
00076         }
00077     }
00078 
00079     // 最終ステージを除いた処理
00080     void FftReal::ExcludeLastStage()
00081     {
00082         uint16_t nHalf = N_FFT_/2;
00083         for (int stg=1; stg<N_FFT_/2; stg*=2)
00084         {
00085             uint16_t nHalf2 = nHalf*2;
00086             for (int kp=0; kp<N_FFT_; kp+=nHalf2)
00087             {
00088                 uint16_t kx = 0;
00089                 for (int k=kp; k<kp+nHalf; k++)
00090                 {
00091                     // バタフライ演算
00092                     Complex uTmp = u_[k+nHalf];
00093                     u_[k+nHalf] = (u_[k] - uTmp)*wTable_[kx];
00094                     u_[k] = u_[k] + uTmp;
00095                     kx = kx + stg;
00096                 }
00097             }
00098             nHalf = nHalf/2;
00099         }
00100     }
00101 }