定義式通りのDFTの実装

クラスクラスタの具体クラスのまず一番最初、式通りのFourier変換を実装する。クラスクラスタの他のクラスのほとんどは抽象クラスではなくてこのクラスのサブクラスになる。結局やることは皆同じだし。

定義通りのDFT

まず、Nが小さな素数の場合用に定義通りのDFTを書いておく。これは以前書いたものと実質的に同じである。

#import "ASDFT1D.h"
@interface ASIntactDFT1D : ASDFT1D {
    unsigned int    N;
    float complex   *wcoef;
    float complex   *buf;    
}
- (unsigned int)length;
//  Private methods
- (void)allocWcoefForDirection:(int)dir;
@end

「intact」というのは素のままとか穢されていないとか言う意味の形容詞。抽象クラスを親クラスにしてメソッドを上書きする。

インスタンスメソッドとして長さと回転因子の配列とバッファを持つ。インスタンス変数に大文字を使うのはスタイル違反だけど、今回は特別。ノーテーションになるべく合わせた方がいい。バッファはさらにこのサブクラスを作ったときにも使い回す予定。

この実装、まずinit...メソッドは

- (id)initWithNumberOfPoints:(unsigned int)np andDirection:(int)dir
{
    self = [super init];
    N = np;
    if (N > 0)
        [self allocWcoefForDirection:dir];
    buf = (float complex *)malloc(N * sizeof(float complex));
    return self;
}

という簡単なもの。長さに0が渡されると0の割り算が発生するのでチェックしている。

allocWcoefForDirection:というメソッドは回転因子の配列を確保して値を詰める作業をしているだけ。ただし回転因子の計算は精度に影響を与えるので注意が必要である。

DFTの計算ではひとつの要素の値を得るのに足し算を$N$回することになるが、ここで桁落ちする可能性がある。DFTの応用の多くは要素ごとの絶対精度が必要なわけではなく、要素間の比の精度が確保できれば十分な場合が多い（DFT後の、最大の要素と最小の要素の両方に同じ有効桁数が必要になることは少ない）が、回転因子に丸め誤差が含まれているとそれがここで累積していく。回転因子の計算では、float幅であればdouble幅で計算して最後に丸めるとか、割り算は最後にするとかで計算精度が違ってくる場合がある。

Fourier変換の計算は

- (void)performForComplexArray:(float complex *)dat withStride:(int)stride
{
    int n, k;
    float complex   *bp = buf;
    float complex   *dap = dat;
    for (k = 0 ; k < N ; k ++) {
        *bp = *dap;
        dap += stride;
        bp ++;
    }
    for (k = 0 ; k < N ; k ++) {
        float complex   tmp = 0;
        bp = buf;
        for (n = 0 ; n < N ; n ++)
            tmp += (*bp) * wcoef[(k * n) % N];
        *dat = tmp;
        dat += stride;  bp ++;
    }
}

バッファにデータをコピーして定義式通りの計算をする。簡単。