C言語でIFFTを実装したい

Question

実現したいこと

Cooley–Tukeyアルゴリズムを使用してIFFTを実装したい

前提

指数関数の符号が逆であるため、FFTの複素共役を取ると逆フーリエ変換を得られる事が分かります。
ただ、FFTと対称的にIFFTを実装する(再帰関数のみで実装する)方法が分からないのでご教授頂ければと思います。

該当のソースコード

一応、以下にCooley–Tukeyアルゴリズムを使用したFFTの実装を載せておきます

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

void fft(double complex* x, int n) {

    if (n == 1) {

        return;
    }

    double complex* xe = malloc(n / 2 * sizeof(double complex));
    double complex* xo = malloc(n / 2 * sizeof(double complex));

    for (int i = 0; i < n / 2; i++) {

        xe[i] = x[2 * i];
        xo[i] = x[2 * i + 1];
    }

    fft(xe, n / 2);
    fft(xo, n / 2);

    for (int i = 0; i < n / 2; i++) {

        double complex w = cexp(-2 * M_PI * I * i / n);
        double complex t = w * xo[i];

        x[i] = xe[i] + t;
        x[i + n / 2] = xe[i] - t;
    }

    free(xe);
    free(xo);

    return;
}

試したこと

IFFTは以下のように実装してみましたが、期待した結果が得られませんでした

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

void ifft(double complex* x, int n) {

    if (n == 1) {

        return;
    }

    double complex* xe = malloc(n / 2 * sizeof(double complex));
    double complex* xo = malloc(n / 2 * sizeof(double complex));

    for (int i = 0; i < n / 2; i++) {

        double complex w = cexp(-2 * M_PI * I * i / n);

        xe[i] = x[i] + x[i + n / 2];
        xo[i] = (x[i] - x[i + n / 2]) * w;
    }

    fft(xe, n / 2);
    fft(xo, n / 2);

    for (int i = 0; i < n / 2; i++) {

        x[i] = xe[i];
        x[i + n / 2] = xo[i];
    }

    free(xe);
    free(xo);

    return;
}

補足情報（FW/ツールのバージョンなど）

Debian 11.5
gcc 10.2.1 (Debian 10.2.1-6)

Comments

[ps_aux_grep]

複素共役をとることへのご自身の理解を併記してください

[iFQ7Vj]

指数関数の符号が逆ですから、
e^-ix = cos -x + i sin -x = cos x - i sin x = bar{e^ix}
だと理解しています

[ps_aux_grep]

はい，であればfft関数のあるコードの27行目
double complex w = cexp(-2 * M_PI * I * i / n);
にある回転因子wの肩の符号を逆転させて
double complex w = cexp(2 * M_PI * I * i / n);
にするだけでifft関数ができます．
ifft後は振幅がn倍されるので，これを1/nにして終わりです．

[iFQ7Vj]

すれ違いのため削除

Answer 1

fft関数のあるコードの27行目double complex w = cexp(-2 * M_PI * I * i / n);にある回転因子wの肩の符号を逆転させてdouble complex w = cexp(2 * M_PI * I * i / n);にするだけでifft関数ができます．

ifft後は振幅がn倍されるので，これを1/n倍にして完成です．

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

void fft(double complex* x, int n) {

    if (n == 1) {

        return;
    }

    double complex* xe = malloc(n / 2 * sizeof(double complex));
    double complex* xo = malloc(n / 2 * sizeof(double complex));

    for (int i = 0; i < n / 2; i++) {

        xe[i] = x[2 * i];
        xo[i] = x[2 * i + 1];
    }

    fft(xe, n / 2);
    fft(xo, n / 2);

    for (int i = 0; i < n / 2; i++) {

        double complex w = cexp(-2 * M_PI * I * i / n);
        double complex t = w * xo[i];

        x[i] = xe[i] + t;
        x[i + n / 2] = xe[i] - t;
    }

    free(xe);
    free(xo);

    return;
}

void _ifft(double complex* x, int n) {
    if (n == 1) {

        return;
    }

    double complex* xe = malloc(n / 2 * sizeof(double complex));
    double complex* xo = malloc(n / 2 * sizeof(double complex));

    for (int i = 0; i < n / 2; i++) {

        xe[i] = x[2 * i];
        xo[i] = x[2 * i + 1];
    }

    _ifft(xe, n / 2);
    _ifft(xo, n / 2);

    for (int i = 0; i < n / 2; i++) {

        double complex w = cexp(2 * M_PI * I * i / n);
        double complex t = w * xo[i];

        x[i] = xe[i] + t;
        x[i + n / 2] = xe[i] - t;
    }

    free(xe);
    free(xo);

    return;
}

void ifft(double complex* x, int n) {
    _ifft(x, n);
    for (int i = 0; i < n; i++) x[i] /= n;
    return ;
}

int main() {
    double complex x[] = {1, 0, 1, 0, 1, 0, 1, 0};
    printf("Input: ");
    for (int i = 0; i < 8; i++) {
        printf("(%f + %fi) ", creal(x[i]), cimag(x[i]));
    }
    printf("\n");

    fft(x, 8);
    
    printf("FFT: ");
    for (int i = 0; i < 8; i++) {
        printf("(%f + %fi) ", creal(x[i]), cimag(x[i]));
    }
    printf("\n");
    
    ifft(x, 8);
    
    printf("IFFT: ");
    for (int i = 0; i < 8; i++) {
        printf("(%f + %fi) ", creal(x[i]), cimag(x[i]));
    }
    return 0;
}

Input: (1.000000 + 0.000000i) (0.000000 + 0.000000i) (1.000000 + 0.000000i) (0.000000 + 0.000000i) (1.000000 + 0.000000i) (0.000000 + 0.000000i) (1.000000 + 0.000000i) (0.000000 + 0.000000i) 
FFT: (4.000000 + 0.000000i) (0.000000 + 0.000000i) (0.000000 + 0.000000i) (0.000000 + 0.000000i) (4.000000 + 0.000000i) (0.000000 + 0.000000i) (0.000000 + 0.000000i) (0.000000 + 0.000000i) 
IFFT: (1.000000 + 0.000000i) (0.000000 + 0.000000i) (1.000000 + 0.000000i) (0.000000 + 0.000000i) (1.000000 + 0.000000i) (0.000000 + 0.000000i) (1.000000 + 0.000000i) (0.000000 + 0.000000i) 

void _fft(double complex* x, int n, int inv) {

    if (n == 1) {

        return;
    }

    double complex* xe = malloc(n / 2 * sizeof(double complex));
    double complex* xo = malloc(n / 2 * sizeof(double complex));

    for (int i = 0; i < n / 2; i++) {

        xe[i] = x[2 * i];
        xo[i] = x[2 * i + 1];
    }

    _fft(xe, n / 2, inv);
    _fft(xo, n / 2, inv);

    for (int i = 0; i < n / 2; i++) {

        double complex w = cexp((inv ? 1 : -1) * 2 * M_PI * I * i / n);
        double complex t = w * xo[i];

        x[i] = xe[i] + t;
        x[i + n / 2] = xe[i] - t;
    }

    free(xe);
    free(xo);

    return;
}

void fft(double complex* x, int n) {
    _fft(x, n, 0);
    return ;
}

void ifft(double complex* x, int n) {
    _fft(x, n, 1);
    for (int i = 0; i < n; i++) x[i] /= n;
    return ;
}

Comments

[iFQ7Vj]

丁寧なご回答ありがとうございました。
IDFT, IFFTについて理解できました。

[iFQ7Vj]

ご教授頂いたところ恐れ多いですが、少し問題があったので書かせて頂きます。

ifft内でifftを呼び出すはずですが、fftが呼ばれています
ifft内のnで割るところですが、再帰的に呼び出されるので2で割るようにしなければなりません

{1, 0, 1, 0, 1, 0, 1, 0}では期待した結果が得られますが、その他のケースではうまく動きません

[ps_aux_grep]

すみません，修正しました．
普通，fft関数の引数にifftかどうかを示すフラグを付けて符号反転の操作し，ifft関数の末尾でnで除算するのが一般的です．