Chirp z変換によるFFTの実装

今回面白いと思ったChirp z変換のFFT。その実装はあんがい当たり前になってしまった。まあそういうもんだろ。

Chirp z変換を使ったFFT

Chirp z変換のFFTも同じようにする。やっぱり定義通りのDFTを親クラスにする。

#import "ASIntactDFT1D.h"
#import "ASRadix2FFT1D.h"

@interface ASChirpzFFT1D : ASIntactDFT1D {
    unsigned int    trueN;
    unsigned int    fftN;
    ASRadix2FFT1D   *fft;
    ASRadix2FFT1D   *fftInv;
    float complex   *pArray;
    float complex   *qArray;
    float complex   *qConj;
}
...

インスタンス変数として畳み込みのための2の基数のFFTを正逆2個持つ。また、式-29のp(n)とq(n)、それからq^†(n)（q(n)の共役）の配列を持たせる。

初期化の実装は

- (id)initWithNumberOfPoints:(unsigned int)np andDirection:(int)dir
{
    self = [super initWithNumberOfPoints:np * 2 andDirection:dir];
    //  note the number of points to be doubled
    trueN = np; //  note that trueN != N
    fftN = 1 << log2OfN(2 * np - 1);
    [self allocArrays];
    fft = [[ASRadix2FFT1D alloc] initWithNumberOfPoints:fftN
                                           andDirection:dir];
    fftInv = [[ASRadix2FFT1D alloc] initWithNumberOfPoints:fftN
                                              andDirection:-dir];
    [self allocArrays];
    return self;
}

というもので、未来の自分に向けてコメントがあるように、親クラスを倍の長さで初期化している。これは式-29のq(n)が2倍の長さ用の回転因子になるためである。ほんとの長さはわざわざtrueNというインスタンス変数に保持しておく。

log2OfNはC関数で

unsigned int    log2OfN(unsigned int n)
{
    unsigned int p = 0;
    while ((1 << p) < n)
        p ++;
    return p;
}

というようなもの。与えられたnより小さくない2のベキの数のlog₂の値を返す。

allocArraysメソッドは

- (void)allocArrays
{
    pArray = (float complex *)malloc(fftN * sizeof(float complex));
    qArray = (float complex *)malloc(fftN * sizeof(float complex));
    qConj = (float complex *)malloc(fftN * sizeof(float complex));
    float complex   *qa = qArray;
    float complex   *qc = qConj;
    int k;
    
    for (k = 0 ; k <= fftN / 2 ; k ++) {
        *qa = wcoef[(k * k) % N];
        *qc = conjf(*qa);
        qc ++;  qa ++;
    }
    float complex   *qp = qa - 2;
    for (; k < fftN ; k ++) {
        *qa = *qp;
        *qc = conjf(*qa);
        qp --;  qc ++;  qa ++;
    }
    [fft performForComplexArray:qConj withStride:1];
    //  qConj has already been transformed here
}

というようなもので式-29をそのまま計算している。親クラスで回転因子が計算できているのでそれを使う。ここにも未来の自分用にコメントがあってq†(n)用の配列はここでFFTしてしまっている。式-30を見ればわかるようにq†(n)はFourier変換先のものしか使わないからこうできる。

FFTは

- (void)performForComplexArray:(float complex *)dat withStride:(int)stride
{
    int k;
    [self setPArrayWithComplexArray:dat withStride:stride];
    [fft performForComplexArray:pArray withStride:1];
    for (k = 0 ; k < fftN ; k ++)
        pArray[k] *= qConj[k];
    [fftInv performForComplexArray:pArray withStride:1];
    for (k = 0 ; k < trueN ; k ++) {
        *dat = pArray[k] * qArray[k] / fftN;
        dat += stride;
    }
}

で式-30そのまま。最後に2の基数のFFTの長さで割っているのは正規化係数の調整のため。

setPArrayWithComplexArray:withStride:メソッドは式-29のp(n)を計算している。中身は要素ごとにかけ算するだけなので見るまでもない。