C言語でIFFTを実装したい

実現したいこと

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

前提

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

該当のソースコード

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

C
1#include <stdio.h>
2#include <stdlib.h>
3#include <math.h>
4#include <complex.h>
5
6void fft(double complex* x, int n) {
7
8    if (n == 1) {
9
10        return;
11    }
12
13    double complex* xe = malloc(n / 2 * sizeof(double complex));
14    double complex* xo = malloc(n / 2 * sizeof(double complex));
15
16    for (int i = 0; i < n / 2; i++) {
17
18        xe[i] = x[2 * i];
19        xo[i] = x[2 * i + 1];
20    }
21
22    fft(xe, n / 2);
23    fft(xo, n / 2);
24
25    for (int i = 0; i < n / 2; i++) {
26
27        double complex w = cexp(-2 * M_PI * I * i / n);
28        double complex t = w * xo[i];
29
30        x[i] = xe[i] + t;
31        x[i + n / 2] = xe[i] - t;
32    }
33
34    free(xe);
35    free(xo);
36
37    return;
38}
39

試したこと

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

C
1#include <stdio.h>
2#include <stdlib.h>
3#include <math.h>
4#include <complex.h>
5
6void ifft(double complex* x, int n) {
7
8    if (n == 1) {
9
10        return;
11    }
12
13    double complex* xe = malloc(n / 2 * sizeof(double complex));
14    double complex* xo = malloc(n / 2 * sizeof(double complex));
15
16    for (int i = 0; i < n / 2; i++) {
17
18        double complex w = cexp(-2 * M_PI * I * i / n);
19
20        xe[i] = x[i] + x[i + n / 2];
21        xo[i] = (x[i] - x[i + n / 2]) * w;
22    }
23
24    fft(xe, n / 2);
25    fft(xo, n / 2);
26
27    for (int i = 0; i < n / 2; i++) {
28
29        x[i] = xe[i];
30        x[i + n / 2] = xo[i];
31    }
32
33    free(xe);
34    free(xo);
35
36    return;
37}
38

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

Debian 11.5
gcc 10.2.1 (Debian 10.2.1-6)

ps_aux_grep

2023/03/22 06:15

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

iFQ7Vj

2023/03/22 06:35 編集

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

ps_aux_grep

2023/03/22 06:48

はい，であれば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

2023/03/22 07:22 編集

すれ違いのため削除

行動規範の内容に同意します

回答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倍にして完成です．

C
1#include <stdio.h>
2#include <stdlib.h>
3#include <math.h>
4#include <complex.h>
5
6void fft(double complex* x, int n) {
7
8    if (n == 1) {
9
10        return;
11    }
12
13    double complex* xe = malloc(n / 2 * sizeof(double complex));
14    double complex* xo = malloc(n / 2 * sizeof(double complex));
15
16    for (int i = 0; i < n / 2; i++) {
17
18        xe[i] = x[2 * i];
19        xo[i] = x[2 * i + 1];
20    }
21
22    fft(xe, n / 2);
23    fft(xo, n / 2);
24
25    for (int i = 0; i < n / 2; i++) {
26
27        double complex w = cexp(-2 * M_PI * I * i / n);
28        double complex t = w * xo[i];
29
30        x[i] = xe[i] + t;
31        x[i + n / 2] = xe[i] - t;
32    }
33
34    free(xe);
35    free(xo);
36
37    return;
38}
39
40void _ifft(double complex* x, int n) {
41    if (n == 1) {
42
43        return;
44    }
45
46    double complex* xe = malloc(n / 2 * sizeof(double complex));
47    double complex* xo = malloc(n / 2 * sizeof(double complex));
48
49    for (int i = 0; i < n / 2; i++) {
50
51        xe[i] = x[2 * i];
52        xo[i] = x[2 * i + 1];
53    }
54
55    _ifft(xe, n / 2);
56    _ifft(xo, n / 2);
57
58    for (int i = 0; i < n / 2; i++) {
59
60        double complex w = cexp(2 * M_PI * I * i / n);
61        double complex t = w * xo[i];
62
63        x[i] = xe[i] + t;
64        x[i + n / 2] = xe[i] - t;
65    }
66
67    free(xe);
68    free(xo);
69
70    return;
71}
72
73void ifft(double complex* x, int n) {
74    _ifft(x, n);
75    for (int i = 0; i < n; i++) x[i] /= n;
76    return ;
77}
78
79int main() {
80    double complex x[] = {1, 0, 1, 0, 1, 0, 1, 0};
81    printf("Input: ");
82    for (int i = 0; i < 8; i++) {
83        printf("(%f + %fi) ", creal(x[i]), cimag(x[i]));
84    }
85    printf("\n");
86
87    fft(x, 8);
88    
89    printf("FFT: ");
90    for (int i = 0; i < 8; i++) {
91        printf("(%f + %fi) ", creal(x[i]), cimag(x[i]));
92    }
93    printf("\n");
94    
95    ifft(x, 8);
96    
97    printf("IFFT: ");
98    for (int i = 0; i < 8; i++) {
99        printf("(%f + %fi) ", creal(x[i]), cimag(x[i]));
100    }
101    return 0;
102}

shell:stdout
1Input: (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) 
2FFT: (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) 
3IFFT: (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)

C:短く書く方
1void _fft(double complex* x, int n, int inv) {
2
3    if (n == 1) {
4
5        return;
6    }
7
8    double complex* xe = malloc(n / 2 * sizeof(double complex));
9    double complex* xo = malloc(n / 2 * sizeof(double complex));
10
11    for (int i = 0; i < n / 2; i++) {
12
13        xe[i] = x[2 * i];
14        xo[i] = x[2 * i + 1];
15    }
16
17    _fft(xe, n / 2, inv);
18    _fft(xo, n / 2, inv);
19
20    for (int i = 0; i < n / 2; i++) {
21
22        double complex w = cexp((inv ? 1 : -1) * 2 * M_PI * I * i / n);
23        double complex t = w * xo[i];
24
25        x[i] = xe[i] + t;
26        x[i + n / 2] = xe[i] - t;
27    }
28
29    free(xe);
30    free(xo);
31
32    return;
33}
34
35void fft(double complex* x, int n) {
36    _fft(x, n, 0);
37    return ;
38}
39
40void ifft(double complex* x, int n) {
41    _fft(x, n, 1);
42    for (int i = 0; i < n; i++) x[i] /= n;
43    return ;
44}

投稿2023/03/22 06:52

編集2023/03/22 08:19

ps_aux_grep

総合スコア1579

iFQ7Vj

2023/03/22 07:35 編集

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

iFQ7Vj

2023/03/22 08:08

ご教授頂いたところ恐れ多いですが、少し問題があったので書かせて頂きます。 ifft内でifftを呼び出すはずですが、fftが呼ばれています ifft内のnで割るところですが、再帰的に呼び出されるので2で割るようにしなければなりません {1, 0, 1, 0, 1, 0, 1, 0}では期待した結果が得られますが、その他のケースではうまく動きません