2の基数のFFTの実装

具体的なサブクラスの実装の続き。今日はいちばん一般的なCooly-Tukey型のFFTの実装。20年前も似たようなことをやっていたような．．．．

2の基数のFFT

2の基数のFFTはASIntactDFT1Dのサブクラスにする。

#import "ASIntactDFT1D.h"
@interface ASRadix2FFT1D : ASIntactDFT1D {
    unsigned int    lengLog2;
    unsigned int    *bitrevTable;
}
- (id)initWithLog2OfPoints:(unsigned int)log2OfLength andDirection:(int)dir;
//  Private methods
- (void)setBitReverseTable;
@end

専用の初期化メソッド（単に2のベキの値をもらうだけ）とビット逆順のテーブルを作るメソッドを追加する。

その初期化メソッドは

- (id)initWithLog2OfPoints:(unsigned int)log2OfLength andDirection:(int)dir
{
    lengLog2 = log2OfLength;
    self = [super initWithNumberOfPoints:(1 << lengLog2) andDirection:dir];
    [self setBitReverseTable];
    return self;
}

- (id)initWithNumberOfPoints:(unsigned int)np andDirection:(int)dir
{
    if (! isPowerOf2(np)) {
        self = [super initWithNumberOfPoints:0 andDirection:dir];
        bitrevTable = nil;
        [self release];
        return [[ASDFT1D alloc] initWithNumberOfPoints:np
                                          andDirection:dir];
    }
    return [[[self class] alloc] initWithLog2OfPoints:log2OfN(np)
                                         andDirection:dir];
}

で、いちおう親クラスの初期化メソッドも上書きしておく。そこではNが2のベキでなかった場合は抽象クラスを呼んで自分を破棄する。このクラスを直接呼ばないかぎりはここへは来ないはず。

回転因子の配列は親クラスであるASIntactDFT1Dの初期化のときに計算されるのでそれをそのまま使う。

具体的な計算は

- (void)performForComplexArray:(float complex *)dat withStride:(int)stride
{
    fftstep(dat, lengLog2, stride, wcoef, 1);
    int k;
    for (k = 0 ; k < (1 << lengLog2) ; k ++) {
        int br = bitrevTable[k];
        if (k < br) {
            float complex   tmp = dat[br];
            dat[br] = dat[k];
            dat[k] = tmp;
        }
    }
}

で、C関数を呼んで、ビット逆順の並び替えをやるだけ。本体のC関数は

static void    fftstep(float complex *fin, int lengLog2, int stride, float complex *wcoef, int wstride)
{
    int nbar = (1 << (lengLog2 - 1)) * stride;
    int n;
    for (n = 0 ; n < nbar ; n ++) {
        int nind = n * stride;
        float complex   temp0 = fin[nind];
        float complex   temp1 = fin[nbar + nind];
        fin[nind] = temp0 + temp1;
        fin[nbar + nind] = (temp0 - temp1) * wcoef[n * wstride];
    }
    if (lengLog2 > 1) {
        fftstep(fin, lengLog2 - 1, stride, wcoef, wstride * 2);
        fftstep(fin + nbar, lengLog2 - 1, stride, wcoef, wstride * 2);
    }
}

でアルゴリズム通りに再帰呼び出しになってる。これが一番簡単。