【Stockham FFT 入門】
このホームページはFFT(高速フーリエ変換)のアルゴリズムの１つである Stockham のアルゴリズムについて解説するホームページです。 Cooley-Tukey のアルゴリズムは理解できたけど、Stockham のアルゴリズムは、 何故あの 漸化式 が出てくるのか分からないという人に向けて書きました。 みなさんは Stockham のアルゴリズムが実は Cooley-Tukey のアルゴリズムの単純な変形であることをご存知でしょうか？ 天下り式に与えられるあの漸化式も Cooley-Tukey のバタフライ演算から簡単に得られます。このホームページを見れば、 あなたもきっと Stockham のアルゴリズムに納得がいくことでしょう。
まず、Cooley-Tukey のアルゴリズムは理解していることを前提に話を進めます。 Cooley-Tukey のアルゴリズムって何？そもそも FFT って何？という人は こちらのページ をご覧ください。
それでは、Cooley-Tukey のアルゴリズムを実装したコードを示しておきましょう。 Cooley-Tukey アルゴリズムを再帰呼び出しを使って実装すると以下のようになります。
#include <complex>
#include <cmath>
#include <utility>
typedef std::complex<double> complex_t;
void butterfly(int n, int q, complex_t* x) // バタフライ演算
// n : 系列長
// q : ブロックの開始位置
// x : バタフライ演算する系列(入出力)
{
const int m = n/2;
const double theta0 = 2*M_PI/n;
if (n > 1) {
for (int p = 0; p < m; p++) {
const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
const complex_t a = x[q + p + 0];
const complex_t b = x[q + p + m];
x[q + p + 0] = a + b;
x[q + p + m] = (a - b) * wp;
}
butterfly(n/2, q + 0, x);
butterfly(n/2, q + m, x);
}
}
void bit_reverse(int N, complex_t* x) // ビット反転並べ替え
// N : 系列長
// x : ビット反転並べ替えする系列(入出力)
{
for (int i = 0, j = 1; j < N-1; j++) {
for (int k = N >> 1; k > (i ^= k); k >>= 1);
if (i < j) std::swap(x[i], x[j]); // x[i]とx[j]を交換する
}
}
void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
butterfly(N, 0, x);
bit_reverse(N, x);
for (int k = 0; k < N; k++) x[k] /= N;
}
void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
for (int p = 0; p < N; p++) x[p] = conj(x[p]);
butterfly(N, 0, x);
bit_reverse(N, x);
for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}
係数 \(\frac{1}{N}\) はフーリエ変換の側に付けてあります(\(N\) は系列長)。 普通の流儀では逆フーリエ変換の方に付けるようです。しかし、 物理的な意味を考えるなら \(\frac{1}{N}\) はフーリエ変換の方に付けるのが自然です。 普通の流儀になれている人は読み替えてください。
Cooley-Tukey アルゴリズムの特徴は、リスト１
の中の bit_reverse() 関数です。
このアルゴリズムは Cooley-Tukey
のアルゴリズムの中でも周波数間引き(DIF)と言われるアルゴリズムですが、
この場合変換の最後に出力の順序を自然なフーリエ変換の順序にするために、
ビット反転並べ替えという並べ替えを行ってやる必要があります。
ここが Stockham のアルゴリズムとの違いで、
Stockham のアルゴリズムではビット反転並べ替えを行わなくとも自然な順序でフーリエ変換の結果が出力されます。
すなわち、bit_reverse() を省略できるようになるわけです。
しかし、出力を自然な順序にしたいだけなら、
bit_reverse() を用いなくとも、
Cooley-Tukey のアルゴリズムを使ってできないことはありません。
以下のようにしてやれば出力は自然な順序で並びます。
#include <complex>
#include <cmath>
typedef std::complex<double> complex_t;
void fft0(int n, int q, complex_t* x, complex_t* y)
// n : 系列長
// q : ブロック開始位置
// x : フーリエ変換する系列(入出力)
// y : 並べ替え用作業配列
{
const int m = n/2;
const double theta0 = 2*M_PI/n;
if (n > 1) {
for (int p = 0; p < m; p++) { // バタフライ演算
const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
const complex_t a = x[q + p + 0];
const complex_t b = x[q + p + m];
y[q + p + 0] = a + b;
y[q + p + m] = (a - b) * wp;
}
fft0(n/2, q + 0, y, x);
fft0(n/2, q + m, y, x);
for (int p = 0; p < m; p++) { // 偶数成分と奇数成分の合成
x[q + 2*p + 0] = y[q + p + 0]; // 偶数成分
x[q + 2*p + 1] = y[q + p + m]; // 奇数成分
}
}
}
void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
complex_t* y = new complex_t[N]; // 並べ替え用作業配列の確保
fft0(N, 0, x, y);
delete[] y;
for (int k = 0; k < N; k++) x[k] /= N;
}
void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
for (int p = 0; p < N; p++) x[p] = conj(x[p]);
complex_t* y = new complex_t[N]; // 並べ替え用作業配列の確保
fft0(N, 0, x, y);
delete[] y;
for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}
サイズが \(N = 2^L\) のフーリエ変換を \(F_N\) とし(\(L\) は 0 以上の整数)、 入力系列を \(x_p~(p = 0,1,\ldots,N-1)\)、フーリエ変換した結果を \(X_k~(k = 0,1,\ldots,N-1)\) とすると、フーリエ変換は \(X_k = F_N(x_p) = \frac{1}{N}\sum_{p=0}^{N-1}x_p W_N^{kp},~ W_N = \exp\left(-j\frac{2\pi}{N}\right),~j = \sqrt{-1}\) と書けます。今 \(N \ge 2\) と仮定し、これを FFT するために再帰で分解してやると、 周波数間引きでは以下のような関係が成り立つことはご存知ですね。
このように、フーリエ変換を再帰で分解してやると再帰の各段階で \(X_{2k}\) と \(X_{2k+1}\) という偶数成分と奇数成分が得られます。 ビット反転並べ替えが起こるのは、 この偶数成分と奇数成分をさらに偶数成分の偶数成分、 偶数成分の奇数成分というふうに分割していくからです。少し考えれば、 最終的にはビット反転順になることが分かるでしょう。
ならば出力を自然な順序にするには、各再帰の最後で偶数成分と奇数成分を合成して \(X_k\) にしてやれば良いことになります。それを行ったのが リスト２ です。しかし、このコードは遅いです。 そこで Stockham のアルゴリズムの出番ということになるのです。
Stockham のアルゴリズムを得るのに、先の自然な順序で出力が並ぶ Cooley-Tukey のアルゴリズムから無駄を省くことから始めてみましょう。 まず気がつくのは配列の読み書きが２ヶ所に分かれていることです。 一つはバタフライ演算を行っている所、 もう一つは偶数成分と奇数成分の合成を行っている所です。 バタフライ演算の所でせっかく配列から読み出しているのですから、 そのまま書き戻すのではなく、 同時に偶数成分と奇数成分の合成も行ってやったら良いのではないでしょうか。
しかし、バタフライ演算の所で偶数成分と奇数成分の合成を行うためには、 再帰呼び出ししているフーリエ変換をまたいで、 合成の操作を移動してやる必要があります。 当然フーリエ変換がそのままではうまくいきません。では、どうしたら良いか。 それは、偶数位置だけを取り出して変換するフーリエ変換と奇数位置だけを取り出して変換するフーリエ変換に変更してやれば良いのです。
再帰の段数が進むと偶数位置の偶数位置、 奇数位置の偶数位置などをフーリエ変換してやる必要が出てきます。 それを実現するには結局のところ配列を \(s\) 間隔 (\(s = 2^h~(h = 0,1,\ldots,L - 1)\)) でアクセスして変換するフーリエ変換が必要になります。 この \(s\) をストライドと呼ぶことにして リスト２ を書き換えてやると以下のようになります。
#include <complex>
#include <cmath>
typedef std::complex<double> complex_t;
void fft1(int n, int s, int q, complex_t* x, complex_t* y);
void fft0(int n, int s, int q, complex_t* x, complex_t* y)
// n : 系列長
// s : ストライド
// q : 偶数位置奇数位置の切り替え引数
// x : フーリエ変換する系列(入出力)
// y : 作業用配列
{
const int m = n/2;
const double theta0 = 2*M_PI/n;
if (n == 1) {}
else {
for (int p = 0; p < m; p++) { // バタフライして偶数成分と奇数成分を合成
const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
const complex_t a = x[q + s*(p + 0)];
const complex_t b = x[q + s*(p + m)];
y[q + s*(2*p + 0)] = a + b;
y[q + s*(2*p + 1)] = (a - b) * wp;
}
fft1(n/2, 2*s, q + 0, y, x); // 偶数位置のフーリエ変換(y が入力 x が出力)
fft1(n/2, 2*s, q + s, y, x); // 奇数位置のフーリエ変換(y が入力 x が出力)
}
}
void fft1(int n, int s, int q, complex_t* x, complex_t* y)
// n : 系列長
// s : ストライド
// q : 偶数位置奇数位置の切り替え引数
// x : フーリエ変換する系列(入力)
// y : 作業用配列兼出力配列(出力)
{
const int m = n/2;
const double theta0 = 2*M_PI/n;
if (n == 1) { y[q] = x[q]; }
else {
for (int p = 0; p < m; p++) { // バタフライして偶数成分と奇数成分を合成
const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
const complex_t a = x[q + s*(p + 0)];
const complex_t b = x[q + s*(p + m)];
y[q + s*(2*p + 0)] = a + b;
y[q + s*(2*p + 1)] = (a - b) * wp;
}
fft0(n/2, 2*s, q + 0, y, x); // 偶数位置のフーリエ変換(y が入力で同時に出力)
fft0(n/2, 2*s, q + s, y, x); // 奇数位置のフーリエ変換(y が入力で同時に出力)
}
}
void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
complex_t* y = new complex_t[N]; // 並べ替え用作業配列の確保
fft0(N, 1, 0, x, y);
delete[] y;
for (int k = 0; k < N; k++) x[k] /= N;
}
void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
for (int p = 0; p < N; p++) x[p] = conj(x[p]);
complex_t* y = new complex_t[N]; // 並べ替え用作業配列の確保
fft0(N, 1, 0, x, y);
delete[] y;
for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}
いきなり相互再帰が登場していますが、
これは配列の書き換えを効率的に行うためです。Stockham
のアルゴリズムでは並べ替えのための作業領域 y が必要です。
並べ替えた結果を一旦作業領域に保存します。保存した結果をすぐに
x に書き戻してやれば相互再帰にならずに済みますが、
それではバタフライ演算と偶数成分奇数成分の合成を一緒にして、
配列の読み書きを減らした意味がなくなります。x と y
を交互に書き換えれば余分な読み書きは最小限にできます。
そのため相互再帰になっています。
簡単に相互再帰の fft0, fft1
について説明しておきましょう。fft0(n, s, q, x, y) はサイズ
n のフーリエ変換で、入力 x
を受け取りストライド s でフーリエ変換した結果を再び
x に書き戻します。y は作業領域として使われます。
q は偶数位置と奇数位置を切り替えるために使われます。
fft1(n, s, q, x, y) はサイズ n のフーリエ変換で、
入力 x を受け取りストライド s
でフーリエ変換した結果を y に書き込みます。
q は偶数位置と奇数位置を切り替えるために使われます。
出力変数が y になっていることに注意してください。
これで x と y を交互に書き換え出来るようになります。
さて、今までの書き換えは理解できたでしょうか。 実はこれで Stockham のアルゴリズムになっています。 めったにお目にかかれないと言うか、 このサイト以外では見られないかもしれない再帰版の Stockham のアルゴリズムです。 試しに \(n=2^{L-h}\),\(~m=2^{-1}n\),\(~s=2^h\),\(~x_h(q,p)=x_h[q+sp]\) とおいてバタフライ演算を 漸化式 に書き換えてみてください。 以下のようになるはずです。\(x_h[~]\) は計算 \(h\) 段目の配列です。 入力系列 \(x_0[~]\) から始めて \(x_L[~]\) を求めるとフーリエ変換出来ます(係数 \(\frac{1}{N}\) をかける必要があります)。
fft0(n,s,q,x,y) あるいは fft1(n,s,q,x,y)
の呼ばれ方に注目して、s=1,2,4,8,... の時、
q の取る値の範囲を考えると、
以下のようになることに注意してください。
再帰０段目 s=1:q=0
再帰１段目 s=2:q=0, 1
再帰２段目 s=4:q=0, 2, 1, 3
再帰３段目 s=8:q=0,4,2,6,1,5,3,7
...
え？何？Stockham のアルゴリズムにならない？う〜ん、そうかもしれませんね。 私もこの変形を最初に思いついたときはそう思いました。「やべぇ！俺、 フーリエ変換なんていう基本的な分野で新しいアルゴリズムを思いついちゃったよ。 俺って天才？うひ！」とそれはもう興奮したもんです。いやマジで。しかし、 よくよく調べてみるとこれは周波数間引き(DIF)の Stockham のアルゴリズムというやつだったんですね。残念！このアルゴリズムが Stockham のアルゴリズムにならないと思った人はたぶん、 時間間引き(DIT)の Stockham のアルゴリズムを想定していたんでしょう。 ベースとなる Cooley-Tukey のアルゴリズムが DIT の場合、 DIT 版の Stockham になります。実際に書いてみると以下のようになります。
#include <complex>
#include <cmath>
typedef std::complex<double> complex_t;
void fft1(int n, int s, int q, complex_t* x, complex_t* y);
void fft0(int n, int s, int q, complex_t* x, complex_t* y)
// n : 系列長
// s : ストライド
// q : 偶数位置奇数位置の切り替え引数
// x : フーリエ変換する系列(入出力)
// y : 作業用配列
{
const int m = n/2;
const double theta0 = 2*M_PI/n;
if (n == 1) {}
else {
fft1(n/2, 2*s, q + 0, y, x); // 偶数位置のフーリエ変換(x が入力 y が出力)
fft1(n/2, 2*s, q + s, y, x); // 奇数位置のフーリエ変換(x が入力 y が出力)
for (int p = 0; p < m; p++) {
const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
const complex_t a = y[q + s*(2*p + 0)];
const complex_t b = y[q + s*(2*p + 1)] * wp;
x[q + s*(p + 0)] = a + b;
x[q + s*(p + m)] = a - b;
}
}
}
void fft1(int n, int s, int q, complex_t* x, complex_t* y)
// n : 系列長
// s : ストライド
// q : 偶数位置奇数位置の切り替え引数
// x : 作業用配列兼出力配列(出力)
// y : フーリエ変換する系列(入力)
{
const int m = n/2;
const double theta0 = 2*M_PI/n;
if (n == 1) { x[q] = y[q]; }
else {
fft0(n/2, 2*s, q + 0, y, x); // 偶数位置のフーリエ変換(y が入力で同時に出力)
fft0(n/2, 2*s, q + s, y, x); // 奇数位置のフーリエ変換(y が入力で同時に出力)
for (int p = 0; p < m; p++) {
const complex_t wp = complex_t(cos(p*theta0), -sin(p*theta0));
const complex_t a = y[q + s*(2*p + 0)];
const complex_t b = y[q + s*(2*p + 1)] * wp;
x[q + s*(p + 0)] = a + b;
x[q + s*(p + m)] = a - b;
}
}
}
void fft(int N, complex_t* x) // フーリエ変換
// N : 系列長
// x : フーリエ変換する系列(入出力)
{
complex_t* y = new complex_t[N]; // 並べ替え用作業配列の確保
fft0(N, 1, 0, x, y);
delete[] y;
for (int k = 0; k < N; k++) x[k] /= N;
}
void ifft(int N, complex_t* x) // 逆フーリエ変換
// N : 系列長
// x : 逆フーリエ変換する系列(入出力)
{
for (int p = 0; p < N; p++) x[p] = conj(x[p]);
complex_t* y = new complex_t[N]; // 並べ替え用作業配列の確保
fft0(N, 1, 0, x, y);
delete[] y;
for (int k = 0; k < N; k++) x[k] = conj(x[k]);
}
どうでしょう。漸化式 は Stockham のアルゴリズムになりましたね。おっと、 時間間引きの場合 \(m=2^h\),\(~t=2^{L-h}\),\(~s=2^{-1}t\),\(~x_h(q,p)=x_h[q+tp]\) とおくのを忘れないでください。
何？まだ Stockham のアルゴリズムにならない？そうですね。実はその後の調査で、 ここで示した方法とは違った Stockham のアルゴリズムがあるらしいことが分かりました。 申し訳ありませんが私も詳しいことは分かっていません。 もう一つの Stockham のアルゴリズムを想定していた人にはごめんなさい。 とりあえず、Stockham のアルゴリズムの種類に関して こちら でまとめてあります。
さて、これで一応 Stockham のアルゴリズムが求まりました。 これでめでたしめでたしとなればいいんですが、実はそんなに甘くありません。 このアルゴリズムは Stockham のアルゴリズムが Cooley-Tukey のアルゴリズムからの単純な変形だと説明するのには有効ですが、 依然として遅いままです。 ようするにこのままでは Stockham のアルゴリズムを用いる意味がないのです。 次章、 セクション2 では、さらなる書き換えを行って高速なStockham のアルゴリズムを導きます。