混合基数のFFTの実装

任意点数に対応するFFTをここ数週間考えてきた。その実装もこれで一段落。今日は、まえに整理した混合基数のFFTのアルゴリズムをObjective-Cで実装する。

混合基数のFFT

混合基数はちょっと面倒。
ほかと同じようにASIntactDFT1Dのサブクラスにする。

#import "ASIntactDFT1D.h"
@interface ASMixedRadixFFT1D : ASIntactDFT1D {
    unsigned int    radixQ;
    id              successiveRadix;
    id              dftOrChirpz;
}
@end

インスタンス変数のradixQは式-37のQである。

ほかの変数は、まずsuccessiveRadixについては初期化メソッドを見ればわかる。

- (id)initWithNumberOfPoints:(unsigned int)np andDirection:(int)dir
{
    self = [super initWithNumberOfPoints:np andDirection:dir];
    radixQ = smallestFactorInteger(np);
    if (radixQ != np)
        successiveRadix = [[[self class] alloc]
                           initWithNumberOfPoints:(np / radixQ) 
                                     andDirection:dir];
    else
        successiveRadix = nil;

    if ([[super class] chirpZDesirableFor:radixQ])
        dftOrChirpz = [[ASChirpzFFT1D alloc] initWithNumberOfPoints:radixQ
                                                       andDirection:dir];
    else {
        if (successiveRadix)
            dftOrChirpz = [[ASIntactDFT1D alloc] initWithNumberOfPoints:radixQ
                                                           andDirection:dir];
        else
            dftOrChirpz = nil;
    }
    return self;
}

なにをしてるかと言うと、与えられた長さの一番小さな素因数をひとつとりだしてそれをradixQとしている。smallestFactorInteger()はCの関数で引数が素数ならそれ自身を、そうでなければ小さい素数から順に割って一番最初に割り切れた素数を返すと言うもの。

長さが素数でないなら、自分と同じクラスのインスタンスを作って式-37のRの値を渡してインスタンス変数のsuccessiveRadixにセットしている。これが入れ子になって素因数分解されることになる。入れ子の最後はsuccessiveRadixがnilにセットされて、素因数の一番最後であることを表す。

dftOrChirpz変数はFFTの計算を見ればわかる。

- (void)performForComplexArray:(float complex *)dat withStride:(int)stride
{
    if (successiveRadix == nil) {
        if (dftOrChirpz == nil)
            [super performForComplexArray:dat withStride:stride];
        else
            [dftOrChirpz performForComplexArray:dat withStride:stride];
    }
    else {
        unsigned int    radixR = (N / radixQ);
        float complex   *bp;
        int             q, s, p;
        
        for (q = 0 ; q < radixQ ; q ++) {
            [successiveRadix performForComplexArray:(dat + q * stride)
                                         withStride:stride * radixQ];
            for (s = 0 ; s < radixR ; s ++) {
                int     idx = s * radixQ + q;
                buf[idx]  = dat[idx * stride] * wcoef[(s * q) % N];
            }
        }
        bp = buf;
        for (s = 0 ; s < radixR ; s ++) {
            [dftOrChirpz performForComplexArray:(buf + s * radixQ) withStride:1];
            for (p = 0 ; p < radixQ ; p ++)
                dat[(p * radixR + s) * stride] = *(bp ++);
        }
    }
}

前半のif文の中身は最後の素因数だった場合で、これはそのままDFTを計算する。そのとき、素因数がN₀より大きい場合にはChirp z変換のFFTを使っている。そうでない場合は親クラスの素のDFTを使う。

後半のelseの中身は、素因数に分解したその途中の場合で、式-42をそのまま展開している。わかりにくいけど並び替えをしながら式-42を2次元のFFTとして計算している。dftOrChirpzは初期化メソッドの中でQがN₀より大きい場合はChirp z変換、小さい場合は素のDFTのインスタンスが設定されている。その判定は親の抽象クラスのchirpZDesirableFor:メソッドでチェックしている。

アルゴリズムの使い分けが全部抽象クラスの初期化メソッドに集中できれば良かったが、そのようにうまくできなかった。後から見て「なにこれ？」とか言いそう。

これは初期化メソッドを

- (id)initWithNumberOfPoints:(unsigned int)np andDirection:(int)dir
{
    self = [super initWithNumberOfPoints:np andDirection:dir];
    radixQ = smallestFactorInteger(n);
    if (radixQ != np)
        successiveRadix = [[[self class] alloc] initWithNumberOfPoints:(np / radixQ) 
                                                          andDirection:dir];
    else
        successiveRadix = nil;
    dftOrChirpz = [[ASDFT1D alloc] initWithNumberOfPoints:radixQ
                                             andDirection:dir];
    return self;
}

としてしまえば、抽象クラスに選ばせるようなコードにすることができるけど、そうすると自分でできる場合（自分自身もASIntactDFT1Dのサブクラスだから）があるのに、わざわざインスタンスを作ってそれにやらせることになる。それはさすがに無駄な気がした。

このへんは効率よりもわかりやすさ、美しさという点から考え直す必要があるな。