どうもー
今回は昔に書いたC言語で高速フーリエ変換を実装する記事の続きです！
↓以前の記事
C++で実装する上でのメリット
まずここから考えていきましょう。
CとC++の違いとして今回着目するのはVector型です！
vector型は動的配列とも呼ばれているC++のみ使える機能です。
普通の配列とは違い、簡単に要素数を後ろに足したり、要素数を取得したり、、などが標準で実装されています。（気になる人はぜひ調べて）
さらに今回着目したいのはvectorは関数の返り値として使える点です！
普通の配列は関数の返り値にすることはできないので、もし配列の中身を変更するような処理のコードを書く際は仮引数にその配列を入れるのが定石でしたね。
でもパッと関数だけを見た時に、返り値に求めているものが来た方が分かりやすいよね、、
て事でゴリゴリにvectorを使っていきます。
vectorの基本機能(push_backとかinsertとか)については解説しないのでわからないことがあったらググってください
FFTのアルゴリズム
高速フーリエ変換については前回の記事で説明したのでここでは簡単に。
普通のフーリエ変換(DFT)は計算量がO(n^2)なのに対して高速フーリエ変換はO(n*log(n))で実現できます。
なんでそうなるかというとDFTにおける同じ計算をまとめることで処理の数が減って、結果計算量の短縮につながっているんですね。
(↓DFTの記事)
ただ制約として対象のデータの長さが２の冪乗である必要があります。
今回の記事ではデータの長さが２の冪乗であるという前提で進めていきますので悪しからず。
必要な計算の流れとしては
- 回転因子の計算
- バタフライ演算
- ビットリバースによる並べ替え
です。以下のコードの実装ではこれらを簡単に説明しながらコードを解説していきます。
コードの実装
ここからコードの実装ですー
include
まずは必要な機能をincludeしていきましょう
#include <iostream>
#include <vector>
#include <complex>
#include <bitset>
#include "math.h"
using namespace std;とりあえず<iostream>は必要ですよね。
それに加えてvector型を使うために<vector>
フーリエ変換で得られるスペクトルや回転因子は複素数なので<complex>
ビットリバースでビット表現を用いるため<bitset>
三角関数などを用いるための”math.h”
をincludeしておきます。
あとstdを省略のためにusing namespace std;も書いています。
回転因子の導出
はいここからが関数です。
まず回転因子の定義から
て見せられてもよくわからないですよね。
てことでイメージ図です。
こんな感じで複素平面上の単位円をN分割したそれらを回転因子と呼びます。
この関数ではその回転因子の座標を全て取得して返すのが目的です。
てことでコードです！
vector<complex<double>> calcRotationFactor(int N){
vector<complex<double>> w;
for(int i=0; i<N; i++){
double real = cos(2*i*M_PI/N);
double img = sin((-2*i*M_PI)/N);
w.push_back(complex<double>(real,img));
}
return w;
}まず返り値の型ですが、vector<complex<double>>となっています！
？？？？？？？？？？
ってなりますよね。
まあじっくり見ていけば簡単です。
まず中のcomplex<double>はdouble型の複素数を格納する型です
最初にincludeした<complex>で定義されている複素数を表す型ですね。
詳しくはcomplex型で調べてください、、
そんで外枠のvector< >はcomplex<double>を要素とする動的配列を定義っていう意味です。
今回の記事では基本的にこの型を用いるのでしっかりイメージをつけときましょう。引数は先ほど説明したN分割のNを仮引数として渡します。
そして本文を見ていきましょう
vector<complex<double>> w;
for(int i=0; i<N; i++){
double real = cos(2*i*M_PI/N);
double img = sin((-2*i*M_PI)/N);
w.push_back(complex<double>(real,img));
}
return w;まず返り値を格納するようにvector<complex>> wを定義します。
そこからはfor文で回転因子の定義通りにwに座標を計算して格納していきます。
その時の格納ではpush_back()で格納していくのですが、このvectorの要素はcomplex型なのでcomplex<double>(real,img)というようにしています。
これで回転因子は終わりです！
ビットリバース
まずはビットリバースってなんぞやから
簡単に説明すると二進数表示した時の0と1(ビット)の並びを真ん中を中心に入れ替えるってことです。
こんな感じですね。注意したいのは0と1を反転させるのではなく、あくまで入れ替えるってことですね。
そんで入れ替えてなにをするかというと
このようにある数列に対して、全てを二進数表示した後にそれらをビットリバースする。その後また10進数表示に直すと新たな数列が得られますよね。
今回のFFTでは次に解説するバタフライ演算を終えたデータに対して、このビットリバースによる配置変換を行う必要があります。
てことでコードです
bitSwap
void bitSwap(bitset<20> *A, int a, int b){
//bit配列Aの位置aとbの値を交換する
bool tmp;
tmp = A->test(a);
A->set(a,A->test(b));
A->set(b,tmp);
}まずはビットの位置を入れ替える関数bitSwapから。
引数ではビットを入れ替える整数を指すポインタbitset<20> *A
入れ替える位置aとbを受け取ります。
？？？？？？？？？？
そうなりますよね、bitset<20> *Aとはなんぞやと
これはincludeした<bitset>に入っているビット形式で整数を保持する型bitset<>のポインタ型です。
詳しいことはbitsetについて調べて欲しいのですが、とりあえず“bitset<20> A”は整数値Aを20桁のビット表示にしたものと思ってもらって大丈夫です。
その20桁ってのはどこから来た数字なの？
痛いところを突かれましたね。
本当はAの大きさに合わせて桁表示を変えたいところなのですが、どうやらbitset<>で指定する桁数は明確な数字でないとerrorになってしまうとのこと、、
てことで泣く泣く20桁のビット表示にしてます。
（まあ20桁のビット表示でデータ長1,048,576まで対応できるので足りないってことはないと思います、、）
そして本文ではビットの入れ替えを行なっています。
//bit配列Aの位置aとbの値を交換する
bool tmp;
tmp = A->test(a);
A->set(a,A->test(b));
A->set(b,tmp);
ここでの入れ変えの手順はよくある一度tmpに値を預けて上書きしていく方法ですね。
今回はビットなのでtmpの型はboolで十分です。またAから値を取り出したり、値を上書きする際は<bitset>で定義されている.test()や.set()を用いています。
詳しくはbitsetで検索してください、、、、
bitReverse
次はビットリバースを行なって数列を取得する関数です！
void bitReverse(int N, int br[N]){
int bitsize = log2(N);
for(int i=0; i<N; i++){
bitset<20> A = i;
for(int j=0; j<(bitsize/2); j++)
bitSwap(&A, j, bitsize-1-j);
br[i]=(int)A.to_ulong();
}引数にはデータ長であるNとビットリバースによって得られた数列を格納する配列br[N]を受け取ります。
なんでbrがvectorではなく配列なのかというとbrは確定で長さが決まっているからですね。
本文ではまずlog2(N)によってビットにした際の長さを計算します。
int bitsize = log2(N);前提よりNは2の冪乗であるのでここで得られるbitsizeも必ず整数値になるはずです。
そして次は1からNまでをビットリバースした値を配列に格納していきます。
for(int i=0; i<N; i++){
bitset<20> A = i;
for(int j=0; j<(bitsize/2); j++)
bitSwap(&A, j, bitsize-1-j);
br[i]=(int)A.to_ulong();
}まず関数bitSwapの仮引数で受け取ったように対象の値を20桁のビット型としてAに格納します。
その後ビットリバースの定義通り中心までの距離が同じ値同士を入れ替えていきます。それがbitSwap(&A, j, bitsize-1-j);で表現されているわけですね。
全てのビットを入れ替えたら整数値に戻して配列に格納します。
そのときに用いるのが<bitset>で定義されている.to_ulong()ですね。これはbitset<>型をlong型に変換できる関数です。今回はこれをint型に変換してbrに格納しています。
これでビットリバースは終わりです！
バタフライ演算
はいー今回の最難関部ですね。
まずアルゴリズムから説明していきます。
説明のためにデータ長がNのデータに対するFFTをスケールNのFFTとします。
てことでデータ長Nのデータx0,x1,,,,,,xN-1に対するFFTは以下のようになります
文字に起こすと以下のような感じです
- N個のデータの上半分{0番~(N/2)-1番}と下半分{(N/2)番~N-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を半分にして手順１に戻る)
まあいわゆる再帰的な処理ですね。
こんなイメージです
簡単に処理の流れを説明すると、上の手順に書いてある通りどんどん再帰的に半分に分けていきます。そんで分け続けると以下のようにデータ長が2になるはずです。
データ長が2になったらその2つのデータを用いてバタフライ演算を行なって処理を終えます
一旦イメージ図の一番右上のスケール2FFTを考えましょう。
スケール2のFFTでは新しいFFTを実行しないので計算の値が出て処理を終えます。そうすると一個前のスケール4のFFTに戻ります。
この時点でこのスケール4のFFTの手順４が終わったことになるので手順５、つまり下半分のスケール2FFTを開始します。
手順5の処理(スケール2のFFT)も計算の値が出て処理が次のFFTに行かないので、一個前のスケール4のFFTに戻ります。
この時点で手順５が終わったことになるので、つまりこのスケール4のFFT全体が終わったことになります。
となると今度はさらにひとつ前のスケール8のFFTに戻るのですが、この時点でスケール8FFTの手順４が終わったことになるので今度は手順５、つまり下半分のスケール4FFTが始まって….
といった感じで最初はどんどん再帰していきスケール2FFTまで行ったら処理を一つ前のFFTに戻して…という流れをすべてのスケール2FFT(N/2回分)が終わるまでやります。
つまりこの演算の全体のイメージとしては
スケール2FFTにたどり着くまで再帰して、たどり着いたらそこでバタフライ演算を実行して値を二つ出して、を繰り返してすべてのスケール２FFTの値を回収する
ですかね
てことでコードです
vector<complex<double>> fftProcessing(vector<complex<double>> *data){
vector<complex<double>> tmp1,tmp2;//次に渡すよう
vector<complex<double>> collect1,collect2;//回収用
int datasize = (int)data->size();//データ長を取得
vector<complex<double>> w= calcRotationFactor(datasize);//回転因子を取得
if(datasize != 2){//データ長が２ではない時
//上半分
for(int i=0; i<(datasize/2); i++)
tmp1.push_back(data->at(i) + data->at(datasize/2+i));
collect1 = fftProcessing(&tmp1);
//下半分
for(int i=0; i<(datasize/2); i++)
tmp2.push_back((data->at(i) - data->at(datasize/2+i)) * w.at(i));//回転因子をかけたものをtmp2に格納
collect2 = fftProcessing(&tmp2);
//回収したスペクトルを結合
collect1.insert(collect1.end(),collect2.begin(),collect2.end());
}else{
collect1.push_back(data->at(0)+data->at(1));
collect1.push_back(data->at(0)-data->at(1));
}
return collect1;
}返り値、引数
まずは関数の返り値と引数を確認しましょう
vector<complex<double>> fftProcessing(vector<complex<double>> *data);返り値は上でも解説した型vector<complex<double>> です。これはフーリエ変換を計算した値を返すためですね。
そして仮引数dataもvector<complex<double>>のポインタ型です。これはFFTのアルゴリズムの中の半分にしたデータを次のFFTに渡すためですね。
次は本文です。
vector<complex<double>> tmp1,tmp2;//次に渡すよう
vector<complex<double>> collect1,collect2;//回収用
int datasize = (int)data->size();//データ長を取得
vector<complex<double>> w= calcRotationFactor(datasize);//回転因子を取得
とりあえず定義しまくってますね。
最初のvector<complex<double>> tmp1,tmp2は次のFFTに渡すためのデータを格納するためのvectorです。
次のvector<complex<double>> collect1,collect2は計算し終わったFFTを回収するためのvectorです。
次にdataから大きさをsize()で取得してint型datasizeに格納します。これがいわゆるこのFFTにおけるスケールNに当たりますね。dataをポインタ型で受け取っているのでdata.size()でなくdata->size()としている点に注意してください。
その次は同様にvector<complex<double>>型でwを定義して先ほど作成した回転因子を計算する関数によって回転因子を取得しています。その際早速datasizeを用いています。
再帰処理部
ここからがややこしいですね。
アルゴリズムの部分で説明したようにスケールがN=2かそうでないかで処理が変わってきます。まずはNが2ではない時、つまり再帰していく部分を考えましょう
if(datasize != 2){//データ長が２ではない時
//上半分
for(int i=0; i<(datasize/2); i++)
tmp1.push_back(data->at(i) + data->at(datasize/2+i));
collect1 = fftProcessing(&tmp1);
//下半分
for(int i=0; i<(datasize/2); i++)
tmp2.push_back((data->at(i) - data->at(datasize/2+i)) * w.at(i));//回転因子をかけたものをtmp2に格納
collect2 = fftProcessing(&tmp2);
//回収したスペクトルを結合
collect1.insert(collect1.end(),collect2.begin(),collect2.end());
}まずはアルゴリズムでいう手順4の部分ですね。コメントアウトの//上半分の部分です。
for文で0からdatasize/2まで回して、tmp1に計算式通りにpush_backしていきます。
for(int i=0; i<(datasize/2); i++)
tmp1.push_back(data->at(i) + data->at(datasize/2+i));それらが終わったらtmp1を次のFFTに渡します。つまり再帰呼び出しですね。
その際に返ってくる値は回収用のvectorであるcollect1に格納します。
collect1 = fftProcessing(&tmp1);これで手順4が終了ですね。次は手順５、つまりデータの下半分の処理も同様に行います。
//下半分
for(int i=0; i<(datasize/2); i++)
tmp2.push_back((data->at(i) - data->at(datasize/2+i)) * w.at(i));//回転因子をかけたものをtmp2に格納
collect2 = fftProcessing(&tmp2);基本的には先ほどと同じですね。違う点としては下半分の処理には回転因子が絡んでくるので忘れずにってところですかね。
そして手順4,5が終われば一つ上のスケールのFFTに戻ります。そのために回収した値を結合します。
collect1.insert(collect1.end(),collect2.begin(),collect2.end());これでcollect1は一つ上のスケールの大きさになっているはずです。
スケール2のFFT
次はN=2(datasize=2)の場合、つまり再起処理の終点の処理の部分です。
先ほどのif文の条件に引っ掛からなかった場合ですね。
else{
collect1.push_back(data->at(0)+data->at(1));
collect1.push_back(data->at(0)-data->at(1));
}ここは定義通りバタフライ演算を行なっていくだけです。
その際にcollect1にその値を格納しているのがポイントです。
return
最後はreturnです。
return collect1;どちらの場合でもcollect1に求めたデータを格納していたのでreturnもcollect1にします。
これでfftのバタフライ演算は終了です！
FFTの関数
次はこれまで作成してきた計算をまとめてFFTとして機能させるような関数です。
vector<complex<double>> fft(vector<double> signal){
int datasize = (int)signal.size();
vector<complex<double>> data(datasize);
vector<complex<double>> spectol(datasize);
//ビットリバースを計算
int index[datasize];
bitReverse(datasize, index);
//信号を複素数vectorに変換
for(int i=0; i<datasize; i++)
data.at(i)=(complex<double>(signal.at(i),0));
//fftの処理部
data = fftProcessing(&data);
//fftProcessingの計算結果をビットリバースして格納
for(int i=0; i<datasize; i++)
spectol.at(index[i])=data.at(i);
return spectol;
}返り値はvector<complex<double>> です。これはフーリエ変換を計算した値を返すためですね。
仮引数のsignalは注意して欲しいんですが、vector<double>です。フーリエ変換する前のデータは実部だけの１次元データなのでcomplexではないですからね。
最初はもろもろの定義ですね。
int datasize = (int)signal.size();
vector<complex<double>> data(datasize);
vector<complex<double>> spectol(datasize);dataが処理に渡す用のvectorでspectolがフーリエ変換の結果を格納するためのvectorです。
//ビットリバースを計算
int index[datasize];
bitReverse(datasize, index);そして先ほどのビットリバースの関数を用いて配列indexにビットリバース数列を格納します。
//信号を複素数vectorに変換
for(int i=0; i<datasize; i++)
data.at(i)=(complex<double>(signal.at(i),0));次は仮引数signalをdataに代入しています。
再帰処理の関数の仮引数の形がvector<complex<double>>なので合わせてあげるためですね。もちろん虚部は0になるようにしています。
//fftの処理部
data = fftProcessing(&data);そして処理部です。（たった一行なのが綺麗）
//fftProcessingの計算結果をビットリバースして格納
for(int i=0; i<datasize; i++)
spectol.at(index[i])=data.at(i);
return spectol;最後にbrに格納された数列に沿って計算結果を並べ替えて高速フーリエ変換の完成です！
まとめ(前回との比較)
やはりFFTの解説の記事は大変ですね、、
書いていて思ったのは
- Vectorを返り値にできるので見た目がすっきりする
- .size()で大きさを取得できるのでスケールNを再帰処理の際に渡す必要がない
- bitsetのおかげでビットリバースがスタイリッシュに
などなどC++で得られる恩恵がたくさんありますね。
また今度classを用いたFFTも実装して記事にしようかなと思います。
ではまたーー
全体のコード
#include <iostream>
#include <vector>
#include <complex>
#include <bitset>
#include "math.h"
using namespace std;
//回転因子の計算
void bitSwap(bitset<20> *A, int a, int b){
//bit配列Aの位置aとbの値を交換する
bool tmp;
tmp = A->test(a);
A->set(a,A->test(b));
A->set(b,tmp);
}
void bitReverse(int N, int br[N]){
int bitsize = log2(N);
for(int i=0; i<N; i++){
bitset<20> A = i;
for(int j=0; j<(bitsize/2); j++)
bitSwap(&A, j, bitsize-1-j);
br[i]=(int)A.to_ulong();
}
}
vector<complex<double>> calcRotationFactor(int N){
vector<complex<double>> w;
for(int i=0; i<N; i++){
double real = cos(2*i*M_PI/N);
double img = sin((-2*i*M_PI)/N);
w.push_back(complex<double>(real,img));
}
return w;
}
vector<complex<double>> fftProcessing(vector<complex<double>> *data){
vector<complex<double>> tmp1,tmp2;//次に渡すよう
vector<complex<double>> collect1,collect2;//回収用
int datasize = (int)data->size();//データ長を取得
vector<complex<double>> w= calcRotationFactor(datasize);//回転因子格納用
if(datasize != 2){//データ長が２ではない時
//上半分
for(int i=0; i<(datasize/2); i++)
tmp1.push_back(data->at(i) + data->at(datasize/2+i));
collect1 = fftProcessing(&tmp1);
//下半分
for(int i=0; i<(datasize/2); i++)
tmp2.push_back((data->at(i) - data->at(datasize/2+i)) * w.at(i));//回転因子をかけたものをtmp2に格納
collect2 = fftProcessing(&tmp2);
//回収したスペクトルを結合
collect1.insert(collect1.end(),collect2.begin(),collect2.end());
}else{
collect1.push_back(data->at(0)+data->at(1));
collect1.push_back(data->at(0)-data->at(1));
}
return collect1;
}
vector<complex<double>> fft(vector<double> signal){
int datasize = (int)signal.size();
vector<complex<double>> data(datasize);
vector<complex<double>> spectol(datasize);
//ビットリバースを計算
int index[datasize];
bitReverse(datasize, index);
//信号を複素数vectorに変換
for(int i=0; i<datasize; i++)
data.at(i)=(complex<double>(signal.at(i),0));
//fftの処理部
data = fftProcessing(&data);
//fftProcessingの計算結果をビットリバースして格納
for(int i=0; i<datasize; i++)
spectol.at(index[i])=data.at(i);
return spectol;
}
int main(int argc, const char * argv[]) {
vector<double> signal = {1,2,3,4,5,6,7,8};
vector<complex<double>> spectol;
spectol = fft(signal);
for(int i=0; i<spectol.size(); i++)
cout<<spectol.at(i)<<endl;
return 0;
}
コメント