質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.44%
C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

Q&A

解決済

1回答

1399閲覧

C言語でIFFTを実装したい

iFQ7Vj

総合スコア52

C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

0グッド

0クリップ

投稿2023/03/22 06:00

編集2023/03/22 06:39

実現したいこと

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)

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

PondVillege

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} だと理解しています
PondVillege

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 編集

すれ違いのため削除
guest

回答1

0

ベストアンサー

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
PondVillege

総合スコア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}では期待した結果が得られますが、その他のケースではうまく動きません
PondVillege

2023/03/22 08:15

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.44%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問