引き続きFFT

Program

高速フーリエ変換。
Cooley-Tukeyタイプの基数４の周波数間引きFFTの場合の流れの処理としては

Radix4のFFT処理
Radix2のFFT処理
ビットリバース処理

んで、Radix2のFFT処理については桁数Nが4^kの場合にはいらない。
つまり桁数が4*4とか4*4*4の形にできる場合ってことで
逆に4*4*2とか4*4*4*2とかだと最後にRadix2を行う。

これをプログラムに落とす場合、
周波数間引きの場合は先にRadix4処理を行うので
そのたびにN/=4としていくことで最後N==0なのかN==2で
Radix2処理を行うかどうかが簡単に判定できる。

イメージで書くとこんな感じで実に簡単。



m=N;

//Radix4処理

while((m>>=2)>=4)//m÷4が4以上の場合

{

}
//Radix2処理

if(m==2)

{

}

それに対して同じCooley-TukeyタイプRadix4の時間間引きFFTの場合は
周波数間引きの処理の流れとまったく逆になる。つまり

ビットリバース処理
Radix2のFFT処理
Radix4のFFT処理

この場合、Radix2処理を行うかどうかの判定を先にしなければならないため
Nを4で割りつづけて最後が0なのか2なのかはあらかじめ計算しなければならない。
たとえば4*4は(2^2)*(2^2)=2^(2+2)=2^4だし、4*4*2は(2^2)*(2^2)*2=2^(2+2+1)=2^5
という流れを考えれば、4で割りつづけて0か2かというのは
結局Nが2^nとおける場合にnが奇数なのか偶数なのかというのを調べればいい。
これをシンプルに表現すればまずN=2^mという形になおしてあげて(mを求めて)
mが偶数か奇数か、つまりmの最下位ビットが1か0かを調べればいい。



m=0;//N=2^mのmを求める

while(N>1)

{

    N>>=1;//N÷2

    ++m;

}

//mが偶数？奇数？

m&=1;

ところでこのN=2^mという数字について考えるともう少しシンプルにできる
2^mという数字は2進数のビット表現で見たときに
最下位ビット(右側)からずっと0が続いて最後に1がくる数字。
たとえば2^1=10とか2^3=1000とか2^4=10000とか2^6=1000000とか。
で、1がくる場所は右側から0番目,1番目,2番目・・と数えてm番目の位置になる。
結局mが偶数か奇数かは1の場所が偶数番目か奇数番目にあるかどうかで判定できる。
今調べたいのはNが4で割り切れずRadix2を処理する場合、つまりmが奇数の場合な訳だから
1010というフラグとNでANDしたーげればいい。
2^0を無視して考えれば1010は16進数Aなので



m=(N&0xAAAAAAAA)?1:0;

と簡単な形で求めることができる。