[C言語]高速フーリエ変換の実装(周波数間引き型FFT)[コード付き]

どうもこんにちは〜　シュモクザメです。

前に離散フーリエ変換をC言語で実装する記事を書いたのですが、今回は離散高速フーリエ変換を実装していきます！
高速フーリエ変換自体の仕組みというよりは、コードを組む方に重きを置いているので注意してください！

追記(2023/5/1)
C++でもFFTを作ってみました！
ぜひ参考にしてください！

高速フーリエ変換 (FFT)とは？
一旦全体のコード
高速フーリエ変換の流れ
コード実装
終わり
最後にもう一度全体のコード

高速フーリエ変換 (FFT)とは？

高速フーリエ変換とはその名の通り、
普通のフーリエ変換より高速で計算できるフーリエ変換です。

んー…つまり計算式が簡単ってこと？

いえ、そういう意味ではなくて計算量が普通のフーリエ変換より少ないという意味です。
具体的には、普通のフーリエ変換はO(n^2)に対し、高速フーリエ変換はO(n*log n)で計算できます。

ただデメリットとしては、

変換する信号の数が2の累乗でなければならない
コードが難しくなる

という点があげられます。
特に変換する信号の数が2のn乗でなければならない点は重要なので要注意です！

またこの記事の以下では高速フーリエ変換のことをFFT(Fast Fourier Transform)と略して表記するのでお願いします。

一旦全体のコード

とりあえず全体のコードを貼っておきます

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

#define swap(a,b) do{int t=a; a=b; b=t;}while(0) //入れ替えを関数く形式マクロで表現
#define NAME 64//ファイル名の長さの上限です。64に特に意味はありません

typedef struct complex {//複素数を構造体で実現
    double r;
    double i;
}complex;

void cre_w(int N,complex w[]){//1のn乗根を三角関数を用いて生成
    for(int i=0; i<N; i++){
        w[i].r = cos(2*i*M_PI/N);
        w[i].i = sin((-2*i*M_PI)/N);
    }
}

void cal(complex a, complex b, complex *y){//複素数同士の積の計算
    y->r = a.r*b.r - a.i*b.i; //yはポインタで渡して計算結果をmainのに戻るようにする
    y->i = a.r*b.i + a.i*b.r;
}

int ff(int i){//渡されたiが2の何乗なのかをintで返す
    int r=0;
    for(; i>1;){
        i /= 2;
        r++;
    }
    return r;
}

int power2(int i){//２のi乗を返す
    int r=1;
    for(int l=0; l<i; l++)
        r *= 2;
    return r;
}

void change2(int N, int k, int bit[ff(N)]){//kを2進数に変換して配列に格納
    for(int j=0; j<k; j++){
        for(int l=0; l<ff(N); l++){
            bit[l]++;
            if(bit[l]<2)//筆算の容量で位上げを実装
                break;
            else
                bit[l]=0;
        }
    }
}

int change10(int N, int bit[ff(N)]){//2進数の配列を10進数に戻す
    int r=0;
    for(int i=0; i<ff(N); i++){
        r += bit[i] * power2(i);
    }
    return r;
}

void bit_revers(int N, int bit[ff(N)]){//2進数の配列をビットリバースする 配列の長さはff(N)
    
    for(int i=0; i<(int)(ff(N)/2); i++)
        swap(bit[i], bit[ff(N)-1-i]);
}

int rev(int N, int k){//10進数の数字を渡すとその数字をビットリバースした数字を返す
    int bit[ff(N)];
    
    for(int i=0; i<ff(N); i++)
        bit[i] = 0;
    
    change2(N,k,bit);
    bit_revers(N, bit);
    
    return change10(N, bit);
}

void clear_com(int N, complex a[N]){//complexの中身をすべて0にする
    for(int i=0; i<N; i++)
        a[i].r = a[i].i = 0;
}

void FFT(int N, complex x[N], complex X[N]){//FFT xには最初は信号を渡す、Xには求める変換が返ってくる
    complex w[N];
    clear_com(N, w);
    cre_w(N, w);//N個の1のN乗根を生成
    
    complex y[N];
    clear_com(N, y);//今回のループの計算結果を格納し、次のループに渡す用の配列
    
    complex X1[N/2];
    complex X2[N/2];//次のループに渡す用の配列
    clear_com(N/2, X1);
    clear_com(N/2, X2);
    
    if(N != 2){//Nが2ではないとき、つまりまだループを繰り返す場合
        int i;
        for(i=0; i<N/2; i++){
            y[i].r = x[i].r + x[N/2+i].r;
            y[i].i = x[i].i + x[N/2+i].i;//バタフライ演算の上の部分
        }
        
        FFT(N/2, y, X1);//上半分の計算結果を次のFFTループに渡す
        
        complex k;
        k.r=0; k.i=0;
        
        for(i=0; i<N/2; i++){
            k.r = x[i].r - x[N/2+i].r;
            k.i = x[i].i - x[N/2+i].i;
            cal(w[i], k, &y[i]);//バタフライ演算の下の部分を計算しyに格納
            
        }
         FFT(N/2, y, X2);//下部分も計算結果を次のFFTループに渡す
         
        for(int j=0; j<N/2; j++){//返ってきたX1,X2をXに順に格納する
            X[j] = X1[j];
            X[N/2+j] = X2[j];
        }
        
    }else{//N=2のとき、つまりループの末端の場合
        
        X[0].r = x[0].r + x[1].r;
        X[0].i = x[0].i + x[1].i;
        X[1].r = x[0].r - x[1].r;
        X[1].i = x[0].i - x[1].i;//バタフライ演算をして１つ上のループにXとして返す
      
    }
}

void SCAN(int N, complex x[N],char fn[NAME]){ //読み込み用関数
    FILE *fp;
    int i;
    if((fp=fopen(fn,"r"))==NULL){
        printf("ファイルを読み込めません\n");
    }else{
        for(i=0; i<N; i++){
            fscanf(fp,"%lf",&x[i].r);
        }
    }
    fclose(fp);
}

void WRITE(int N, complex x[N],char fn[NAME]){ //書き込み用関数
    FILE *fp;
    int i;
    
    if((fp=fopen(fn,"w"))==NULL)
        printf("ファイルに書き込めません\n");
    else{
        for(i=0; i<N; i++){
            fprintf(fp,"%lf %lf\n",x[i].r,x[i].i);
        }
    }
    fclose(fp);
}

int main(int argc, const char * argv[]) {
  
    int N=256;//ここはフーリ変換にかける信号長の長さを自分で設定する
    char fn[NAME];
    complex x[N];//信号を格納
    complex X[N];//関数FFTで最終的に求める変換を格納する配列
    complex X_t[N];//ビットリバースした数字に対応させる配列
    
    clear_com(N, x);
    clear_com(N, X);
    clear_com(N, X_t);
    
    printf("どのファイルから読み込みますか?");
    scanf("%s",fn);
    SCAN(N, x, fn);
    
    FFT(N, x, X);
    printf("hhhh\n");
     printf("%f %fa\n%f %f\n",X[2].r,X[2].i,X[3].r,X[3].i);
    for(int i=0; i<N; i++){//X_tにXの結果を入れる
        X_t[rev(N, i)] = X[i];
    }
    printf("hhhh\n");
    printf("%f %fa\n%f %f\n",X_t[2].r,X_t[2].i,X_t[3].r,X_t[3].i);
    printf("FFTが完了しました\n");
    
    printf("どのファイルに書き込みますか?");
    scanf("%s",fn);
    WRITE(N, X_t, fn);
    
    return 0;
}

最初の方に関数多すぎだし
途中のFFTっていう関数ややこしすぎるなあ

まあかなりややこしいっすよね。
一つずつ説明していくので頑張っていきましょう。

高速フーリエ変換の流れ

最初に大雑把に流れを確認します。

回転因子を求める
ビットリバースを考える
バタフライ演算

という３点に分けられます。
回転因子とビットリバースはバタフライ演算で用いるもので、主にバタフライ演算がこの高速フーリエ変換の山場です。

以下では実際のコードを実装しながら上の三点を簡単に解説していきます。

コード実装

下準備

まず今回のコードを組むにあたって必要な関数や構造体を先に紹介していきます。

まずは複素数関連
構造体と積計算の関数です。

typedef struct complex {//複素数を構造体で実現
    double r;
    double i;
}complex;

void cal(complex a, complex b, complex *y){//複素数aとbの積の計算 結果は複素数yに入れる
    y->r = a.r*b.r - a.i*b.i; //yはポインタで渡して計算結果を戻るようにする
    y->i = a.r*b.i + a.i*b.r;
}

詳しく知りたい方はこちらの記事をどうぞ↓

次は進数変換関連

int ff(int i){//渡されたiが2の何乗なのかをintで返す
    int r=0;
    for(; i>1;){
        i /= 2;
        r++;
    }
    return r;
}

int power2(int i){//２のi乗を返す
    int r=1;
    for(int l=0; l<i; l++)
        r *= 2;
    return r;
}

void change2(int N, int k, int bit[ff(N)]){//kを2進数に変換して配列に格納
    for(int j=0; j<k; j++){
        for(int l=0; l<ff(N); l++){//信号長Nのff(N)の長さの2進数を作る
            bit[l]++;
            if(bit[l]<2)//筆算の容量で位上げを実装
                break;
            else
                bit[l]=0;
        }
    }
}

int change10(int N, int bit[ff(N)]){//2進数の配列を10進数に戻す
    int r=0;
    for(int i=0; i<ff(N); i++){
        r += bit[i] * power2(i);
    }
    return r;
}

power2は後に使います。
2進数の状態ではbitという配列に各桁ずつ値を格納してますが、bit自体の長さ (つまり考えている2進数の桁数)は信号長Nの因数2の数です(コード内の関数ff()にNを渡したもの)

次はファイル操作関連

void SCAN(int N, complex x[N],char fn[NAME]){ //読み込み用関数
    FILE *fp;
    int i;
    if((fp=fopen(fn,"r"))==NULL){
        printf("ファイルを読み込めません\n");
    }else{
        for(i=0; i<N; i++){
            fscanf(fp,"%lf",&x[i].r);
        }
    }
    fclose(fp);
}

void WRITE(int N, complex x[N],char fn[NAME]){ //書き込み用関数
    FILE *fp;
    int i;
    
    if((fp=fopen(fn,"w"))==NULL)
        printf("ファイルに書き込めません\n");
    else{
        for(i=0; i<N; i++){
            fprintf(fp,"%lf %lf\n",x[i].r,x[i].i);
        }
    }
    fclose(fp);
}

以上の関数等をこれからバキバキに使っていきます！

回転因子の導出

まず回転因子とは以下の複素数です

分割数？いまいちイメージが湧かんな…

以下の画像だとよくわかると思います

これは複素数平面上の単位円を4分割している点を表しています。
一般化するとこのような感じですね。

こんな感じでオイラーの公式を用いて実部と虚部に分解して一般化できます。

この一般化の式を元に関数を作ります。

void cre_w(int N,complex w[]){//回転因子を配列に格納する
    for(int i=0; i<N; i++){
        w[i].r = cos(2*i*M_PI/N);
        w[i].i = sin((-2*i*M_PI)/N);
    }
}

見たまんまですね。
下準備で作った複素数の構造体の配列にに実部と虚部に分けて代入していきます。

ビットリバース

ビットリバースとはある二進数をその並びの真ん中を中心に入れ替えることです

今回は以下のように連続した整数の並びを入れ替えるためにビットリバースを用います

以上のコードがこちら

#define swap(a,b) do{int t=a; a=b; b=t;}while(0) //入れ替えを関数く形式マクロで表現

void bit_revers(int N, int bit[ff(N)]){//2進数の配列をビットリバースする 配列の長さはff(N)
    
    for(int i=0; i<(int)(ff(N)/2); i++)
        swap(bit[i], bit[ff(N)-1-i]);
}

int rev(int N, int k){//10進数の数字を渡すとその数字をビットリバースした数字を返す
    int bit[ff(N)];
    
    for(int i=0; i<ff(N); i++)
        bit[i] = 0;
    
    change2(N,k,bit);
    bit_revers(N, bit);
    
    return change10(N, bit);
}

関数bit_reversでは2進数が入った配列bitの中身を先頭と末尾から順に一つずつ関数形式マクロswapで入れ替えて行っています。
入れ替える回数は配列bitの長さff(N0)を2で割ったものです。これでビットリバースが完了しましたね。

関数revの流れは単純に

2進数変換後のビットの長さを関数ff()で調べて
受け取った10進数の整数kを関数change2()でその長さの2進数にする
その2進数を関数bit_reversでビットリバースする
ビットリバースしたものをchange10()で10進数に戻したものを返り値とする

ここで紹介した関数ff()、change2、change10は上の下準備の項目でコードを紹介してるので確認しといてください。

バタフライ演算

今回の山場のバタフライ演算です。

とりあえず一般化したバタフライ演算のアルゴリズムどうぞ
N個の信号を処理するときのFFTをスケールNのFFTとしているのでご注意ください

この演算の画像のを文字に起こしてみると、

L個の信号の上半分{0番~(N/2)-1番}と下半分{(N/2)番~L-1番}に分ける。
手順１で並べた信号の、i番目と(N/2)+i番目を足したものを(i=0,1,2…(N/2)-1)の順で上から並べ、
手順１で並べた信号の、i番目から(N/2)+i番目を引いたものを(i=0,1,2…(N/2)-1)の順で手順２で並べたものの下に並べる。
手順2で並べた(N/2)個の信号にスケール(N/2)のFFTを施す。(Nを半分にして手順1に戻る)
手順3で並べた(N/2)個の信号の上から順に回転因子w(N,i)を掛けて、スケール(N/2)FFTを施す。(Nを半分にして手順１に戻る)

ここのアルゴリズムは絶対に理解しましょう。

手順4,5の通り、FFT実行中に半分のスケールのFFTを実行します。
その新たなFFTの中でもさらに半分のFFTを呼び出して…
つまり再帰的な計算ということですね。
といった感じで以下のようなイメージがわくと思います。

どんどん再帰していって、最後はどのタイミングで終了するの？

イメージ図にある通り新しく実行するFFTのスケール(N/2スケール)が2の時は、そのFFTを最後の再帰とします。まあ2より小さい場合はバタフライ演算できませんからね。

一旦イメージ図の一番右上のスケール2FFTを考えましょう。
スケール2FFTでは新しいFFTを実行しないので計算の値が出て処理を終えます。そうすると一個前のスケール4FFTに戻ります。
この時点でこのスケール4FFTの手順４が終わったことになるので手順５つまり下半分のスケール2FFTを開始します。

手順5の処理(スケール2FFT)も計算の値が出て処理が次のFFTに行かないので、一個前のN=4のFFTに戻ります。
この時点で手順５が終わったことになるので、つまりこのスケール4FFTが終わったことになります。

となると今度はさらにひとつ前のスケール8のFFTに戻るのですが、この時点でスケール8FFTの手順４が終わったことになるので今度は手順５、つまり下半分のスケール4FFTが始まって….

といった感じで最初はどんどん再帰していきスケール2FFTまで行ったら処理を一つ前のFFTに戻して…という流れをすべてのスケール2FFT(N/2回分)が終わるまでやります。

つまりこの演算の全体のイメージとしては
スケール2FFTにたどり着くまで再帰して、たどり着いたらそこでバタフライ演算を実行して値を二つ出して、を繰り返してすべてのスケール２FFTの値を回収する
ですかね

FFTのコードです。

void FFT(int N, complex x[N], complex X[N]){//FFT xには最初は信号を渡す、Xには求める変換が返ってくる
    complex w[N];
    clear_com(N, w);
    cre_w(N, w);//N個の1のN乗根を生成
    
    complex y[N];
    clear_com(N, y);//今回のループの計算結果を格納し、次のループに渡す用の配列
    
    complex X1[N/2];
    complex X2[N/2];//次のループに渡す用の配列
    clear_com(N/2, X1);
    clear_com(N/2, X2);
    
    if(N != 2){//Nが2ではないとき、つまりまだループを繰り返す場合
        int i;
        for(i=0; i<N/2; i++){
            y[i].r = x[i].r + x[N/2+i].r;
            y[i].i = x[i].i + x[N/2+i].i;//バタフライ演算の上の部分
        }
        
        FFT(N/2, y, X1);//上半分の計算結果を次のFFTループに渡す
        
        complex k;
        k.r=0; k.i=0;
        
        for(i=0; i<N/2; i++){
            k.r = x[i].r - x[N/2+i].r;
            k.i = x[i].i - x[N/2+i].i;
            cal(w[i], k, &y[i]);//バタフライ演算の下の部分を計算しyに格納
            
        }
         FFT(N/2, y, X2);//下部分も計算結果を次のFFTループに渡す
         
        for(int j=0; j<N/2; j++){//返ってきたX1,X2をXに順に格納する
            X[j] = X1[j];
            X[N/2+j] = X2[j];
        }
        
    }else{//N=2のとき、つまりループの末端の場合
        
        X[0].r = x[0].r + x[1].r;
        X[0].i = x[0].i + x[1].i;
        X[1].r = x[0].r - x[1].r;
        X[1].i = x[0].i - x[1].i;//バタフライ演算をして１つ上のループにXとして返す
    }
}

特に説明する点は

if文でスケール2か判断している
処理した信号を次のFFTに渡す仮引数x
スケール2のFFTで処理した信号を回収する仮引数X
上半分のFFTで回収する信号はX1、下半分のFFTで回収する信号はX2として、それをX1X2の順で２倍の配列に格納して、一つ前のFFTに回収させている

ですかね。
まあ解説の部分を見れば今何をやっているかはなんとなくわかると思います。

仕上げ

実はまだ終わりじゃないんですね。
最後にバタフライ演算で回収した値の配列をビットリバースで並び変える必要があります。

まあでもビットリバースもすでに作ってあるのでここは簡単ですね。
main関数内で以下の処理をして配列の要素の順番を入れ替えた配列を作ります

for(int i=0; i<N; i++){//X_tにXの結果を入れる
        X_t[rev(N, i)] = X[i];
    }

終わり

あ〜　疲れた。
再帰的な関数を説明するところがマジで大変でした..

まあでもそれなりに有用なコードが書けたと思います。
みなさんぜひ使ってみてください！

最後にもう一度全体のコード

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

#define swap(a,b) do{int t=a; a=b; b=t;}while(0) //入れ替えを関数く形式マクロで表現
#define NAME 64//ファイル名の長さの上限です。64に特に意味はありません

typedef struct complex {//複素数を構造体で実現
    double r;
    double i;
}complex;

void cre_w(int N,complex w[]){//1のn乗根を三角関数を用いて生成
    for(int i=0; i<N; i++){
        w[i].r = cos(2*i*M_PI/N);
        w[i].i = sin((-2*i*M_PI)/N);
    }
}

void cal(complex a, complex b, complex *y){//複素数同士の積の計算
    y->r = a.r*b.r - a.i*b.i; //yはポインタで渡して計算結果をmainのに戻るようにする
    y->i = a.r*b.i + a.i*b.r;
}

int ff(int i){//渡されたiが2の何乗なのかをintで返す
    int r=0;
    for(; i>1;){
        i /= 2;
        r++;
    }
    return r;
}

int power2(int i){//２のi乗を返す
    int r=1;
    for(int l=0; l<i; l++)
        r *= 2;
    return r;
}

void change2(int N, int k, int bit[ff(N)]){//kを2進数に変換して配列に格納
    for(int j=0; j<k; j++){
        for(int l=0; l<ff(N); l++){
            bit[l]++;
            if(bit[l]<2)//筆算の容量で位上げを実装
                break;
            else
                bit[l]=0;
        }
    }
}

int change10(int N, int bit[ff(N)]){//2進数の配列を10進数に戻す
    int r=0;
    for(int i=0; i<ff(N); i++){
        r += bit[i] * power2(i);
    }
    return r;
}

void bit_revers(int N, int bit[ff(N)]){//2進数の配列をビットリバースする 配列の長さはff(N)
    
    for(int i=0; i<(int)(ff(N)/2); i++)
        swap(bit[i], bit[ff(N)-1-i]);
}

int rev(int N, int k){//10進数の数字を渡すとその数字をビットリバースした数字を返す
    int bit[ff(N)];
    
    for(int i=0; i<ff(N); i++)
        bit[i] = 0;
    
    change2(N,k,bit);
    bit_revers(N, bit);
    
    return change10(N, bit);
}

void clear_com(int N, complex a[N]){//complexの中身をすべて0にする
    for(int i=0; i<N; i++)
        a[i].r = a[i].i = 0;
}

void FFT(int N, complex x[N], complex X[N]){//FFT xには最初は信号を渡す、Xには求める変換が返ってくる
    complex w[N];
    clear_com(N, w);
    cre_w(N, w);//N個の1のN乗根を生成
    
    complex y[N];
    clear_com(N, y);//今回のループの計算結果を格納し、次のループに渡す用の配列
    
    complex X1[N/2];
    complex X2[N/2];//次のループに渡す用の配列
    clear_com(N/2, X1);
    clear_com(N/2, X2);
    
    if(N != 2){//Nが2ではないとき、つまりまだループを繰り返す場合
        int i;
        for(i=0; i<N/2; i++){
            y[i].r = x[i].r + x[N/2+i].r;
            y[i].i = x[i].i + x[N/2+i].i;//バタフライ演算の上の部分
        }
        
        FFT(N/2, y, X1);//上半分の計算結果を次のFFTループに渡す
        
        complex k;
        k.r=0; k.i=0;
        
        for(i=0; i<N/2; i++){
            k.r = x[i].r - x[N/2+i].r;
            k.i = x[i].i - x[N/2+i].i;
            cal(w[i], k, &y[i]);//バタフライ演算の下の部分を計算しyに格納
            
        }
         FFT(N/2, y, X2);//下部分も計算結果を次のFFTループに渡す
         
        for(int j=0; j<N/2; j++){//返ってきたX1,X2をXに順に格納する
            X[j] = X1[j];
            X[N/2+j] = X2[j];
        }
        
    }else{//N=2のとき、つまりループの末端の場合
        
        X[0].r = x[0].r + x[1].r;
        X[0].i = x[0].i + x[1].i;
        X[1].r = x[0].r - x[1].r;
        X[1].i = x[0].i - x[1].i;//バタフライ演算をして１つ上のループにXとして返す
      
    }
}

void SCAN(int N, complex x[N],char fn[NAME]){ //読み込み用関数
    FILE *fp;
    int i;
    if((fp=fopen(fn,"r"))==NULL){
        printf("ファイルを読み込めません\n");
    }else{
        for(i=0; i<N; i++){
            fscanf(fp,"%lf",&x[i].r);
        }
    }
    fclose(fp);
}

void WRITE(int N, complex x[N],char fn[NAME]){ //書き込み用関数
    FILE *fp;
    int i;
    
    if((fp=fopen(fn,"w"))==NULL)
        printf("ファイルに書き込めません\n");
    else{
        for(i=0; i<N; i++){
            fprintf(fp,"%lf %lf\n",x[i].r,x[i].i);
        }
    }
    fclose(fp);
}

int main(int argc, const char * argv[]) {
  
    int N=256;//ここはフーリ変換にかける信号長の長さを自分で設定する
    char fn[NAME];
    complex x[N];//信号を格納
    complex X[N];//関数FFTで最終的に求める変換を格納する配列
    complex X_t[N];//ビットリバースした数字に対応させる配列
    
    clear_com(N, x);
    clear_com(N, X);
    clear_com(N, X_t);
    
    printf("どのファイルから読み込みますか?");
    scanf("%s",fn);
    SCAN(N, x, fn);
    
    FFT(N, x, X);
    printf("hhhh\n");
     printf("%f %fa\n%f %f\n",X[2].r,X[2].i,X[3].r,X[3].i);
    for(int i=0; i<N; i++){//X_tにXの結果を入れる
        X_t[rev(N, i)] = X[i];
    }
    printf("hhhh\n");
    printf("%f %fa\n%f %f\n",X_t[2].r,X_t[2].i,X_t[3].r,X_t[3].i);
    printf("FFTが完了しました\n");
    
    printf("どのファイルに書き込みますか?");
    scanf("%s",fn);
    WRITE(N, X_t, fn);
    
    return 0;
}