FFTにおけるビット逆転

なぜビットが逆転するのかよくわかりません。
バタフライ演算を用いたシグナルフローを見れば、逆転しているのがよくわかるのですが、数式だけ見ると、よくわからなくなります。
ビット逆転を数式で証明できる方、おりましたらよろしくお願いいたします。

No.1 回答者： keyguy 回答日時： 2004/07/23 16:30

FFTは並べ替え方によって2通り有るが

最初に奇数項と偶数項をそれぞれまとめる方法について言うと

奇数項だけを取り出して数が半分のFFTをつくり
偶数項だけを取り出して数が半分のFFTを作る
この奇数項と偶数項を順に並べたものの順番（０からふる）のビット表示は元の順番の1位桁と最上位桁を交換したものである

半分になったFFTはさらに同じようにしてその半分の半分のFFT（元の1/4のFFT)になる
以下2個（1個でもよい）のFFTになるまで繰り返す

以上の事からビット交換の仕組みが分かるだろう

こんなのは黒板で説明すればすぐに理解できるし説明するほうもほとんど負担が無い

もう1つの方法について補足したまえ

この回答への補足

説明不足ですみません。
時間間引きでも周波数間引きでもどっちでもいいのですが、
基数２のＦＦＴを考えた場合、フーリエ変換値をＸとすると、

Ｘ(ｋ2，ｋ1，ｋ0)＝ΣΣΣｘ0(ｎ2，ｎ1，ｎ0) Ｗ

となります(Ｗは回転子です。次数は省略します)。
ここで、ｋ（基本周波数の整数倍）・ｎは、２進数で考えて、
ｋ＝４ｋ2＋２ｋ1＋ｋ0　　（ｋ2，ｋ1，ｋ0＝０，１）
ｎ＝４ｎ2＋２ｎ1＋ｎ0　　（ｎ2，ｎ1，ｎ0＝０，１）
であり、
1つ目のΣはｎ0＝０～１までの和（外側）
2つ目のΣはｎ1＝０～１までの和（真ん中）
3つ目のΣはｎ2＝０～１までの和（内側）
です。
この式の３つのΣを内側、真ん中、外側で分けて考えたとき、

ｘ1＝(ｋ0，ｎ1，ｎ0)＝Σｘ0(ｎ2，ｎ1，ｎ0) Ｗ
ｘ2＝(ｋ0，ｋ1，ｎ0)＝Σｘ1(ｋ0，ｎ1，ｎ0) Ｗ
ｘ3＝(ｋ0，ｋ1，ｋ2)＝Σｘ2(ｋ0，ｋ1，ｎ0) Ｗ

となります。よって、

Ｘ(ｋ2，ｋ1，ｋ0)＝ｘ3＝(ｋ0，ｋ1，ｋ2)

となり、ｋの順序が逆になっていることが分かります。これがまさにビット逆転なのですが、これらの式展開は、シグナルフローを見て考えれば理解できるのですが、式だけ眺めていると何故かわからなくなってしまいます（ビット逆転を仮定して考えると理解できるが、最終的にビット逆転されることを考慮に入れないと、わからないのです）。
分かりづらかったらすみません。

つまり、初めの式を分割した式のところで、ｋとｎが混合されるところがなぜか分からないのです。
ご指導の方、よろしくお願いいたします。

補足日時：2004/07/23 20:38

No.2 回答者： keyguy 回答日時： 2004/07/23 21:52

用は出来上がりのバタフライ図を見ているから分からなくなるだけで

１つのFFTが２つのFFTで構成されている家庭を分析すれば分かるはずです
あとはその単純な繰り返しです

入力側から分割していく方法がNo.1で出力側から分割していく方法が補足です
（時間、周波数と呼ぶよりもこの方が直接的でしょう）

2つを理解するのならば簡単でせう
2つへの分解は一般性があるのでそれで十分
後は単純な繰り返し

この回答へのお礼

うーん。。。ちょっとわからないですね。サンプリング点数をもっと少なくして、もう少し考えてみます。
ご回答（というかアドバイスですね。）ありがとうございました。

お礼日時：2004/07/23 22:04

No.3 回答者： Rossana 回答日時： 2004/09/28 01:01

ANo.1の補足回答

2進数のビットで見て行きましょう．
X(k2, k1, k0)＝ΣΣΣx0(n2, n1, n0)exp[-j(2π/8)kn]
＝ΣΣΣx0(n2, n1, n0)exp[-j(2π/8)k(4n2＋2n1＋n0)]
＝ΣΣΣx0(n2, n1, n0)exp[-j(2π/8)k(4n2)]exp[-j(2π/8)k(2n1)]exp[-j(2π/8)kn0]
＝ΣΣΣx0(n2, n1, n0)exp[-jπkn2]exp[-j(π/2)kn1]exp[-j(π/4)kn0]

まず，n2の段について計算
x1(k0, n1, n0)＝Σx0(n2, n1, n0)exp[-jπn2(4k2＋2k1＋k0)]
＝Σx0(n2, n1, n0)exp[-jπn2k0]
＝x0(n2, n1, 0)＋x0(n2, n1, 1)exp[-jπk0]

次に，n1の段について計算
x2(k0, k1, n0)＝Σx1(k0, n1, n0)exp[-j(π/2)n1(4k2＋2k1＋k0)]
＝Σx1(k0, n1, n0)exp[-j(π/2)n1(2k1＋k0)]
＝x1(k0, 0, n0)＋x1(k0, 1, n0)exp[-jπk1]exp[j(π/2)k0]

次に，n0の段について計算
x3(k0, k1, k2)＝Σx2(k0, k1, n0)exp[-j(π/4)n0(4k2＋2k1＋k0)]
＝Σx2(k0, k1, n0)exp[-j(π/4)n0(4k2＋2k1＋k0)]
＝x2(k0, k1, 0)＋x2(k0, k1, n0)exp[-jπk2]exp[j(π/4)(2k1＋k0)]

よって，
X(k2, k1, k0)＝x3(k0, k1, k2)
とXとx3の配列がビット逆順の関係になっています．

この回答へのお礼

久しぶりに見たら回答がありうれしく思います。
ありがとうございます。

お礼日時：2005/02/04 09:22

No.4 回答者： Rossana 回答日時： 2004/09/28 01:07

すいません．ちょっとミスがありました．訂正お願いします．

●訂正箇所

まず，n2の段について計算
x1(k0, n1, n0)＝Σx0(n2, n1, n0)exp[-jπn2(4k2＋2k1＋k0)]
＝Σx0(n2, n1, n0)exp[-jπn2k0]
＝x0(0, n1, n0)＋x0(1, n1, n0)exp[-jπk0]

次に，n0の段について計算
x3(k0, k1, k2)＝Σx2(k0, k1, n0)exp[-j(π/4)n0(4k2＋2k1＋k0)]
＝Σx2(k0, k1, n0)exp[-j(π/4)n0(4k2＋2k1＋k0)]
＝x2(k0, k1, 0)＋x2(k0, k1, 1)exp[-jπk2]exp[j(π/4)(2k1＋k0)]

まだどこかミスがあるかもしれません．あったら指摘下さい．

この回答へのお礼

ありがとうございます。
このように展開する方法は解ります。間違いもないと思います。

ただ、理論的に説明するとどうなるのでしょうか？

わからないのは、n2の段について計算において、出力x1(k0, n1, n0)に、なぜ「k0」が入ってくるのかがよく解らないのです。
確かにここを「k0」ではなく「k1」にしてしまうと、計算できなくなってしまうので、なんとなくは解るのですが、論理的に説明できないのです。

よろしくお願いいたします。

お礼日時：2005/02/04 09:21

No.5 回答者： Rossana 回答日時： 2005/02/22 01:34

>n2の段について計算において、出力x1(k0, n1, n0)に、なぜ「k0」が入ってくるのかがよく解らないのです。

x1(k0, n1, n0)＝x0(0, n1, n0)＋x0(1, n1, n0)exp[-jπk0]
とk0の関数になっているからというのでは説明不足でしょうか？どういうふうに分からないかもう少し具体的に教えて頂けますか？私も真に完全に理解できているわけではないのですが。実際にプログラムを作ってみると分かるかもしれませんよ。

この回答へのお礼

こんにちは。度々ありがとうございます。

そうですね。k0の関数になっているからというのは、今の僕の中の結論としてあるんです（一応…）。

ただ、何かしっくりこないんですよね(^^;)
これは論理的なのでしょうか？

回転子にk0が現れるからといって、「出力の最上位桁にk0が来る理由はない」「他のパラメーターではダメなのか？」と思ってしまいます。

考えすぎでしょうか？

お礼日時：2005/02/22 13:10

No.6 回答者： Rossana 回答日時： 2005/02/22 17:45

>回転子にk0が現れるからといって、「出力の最上位桁にk0が来る理由はない」「他のパラメーターではダメなのか？」と思ってしまいます。

>考えすぎでしょうか？

考えすぎではなくごもっともな意見だと思います．天下り的に多くの本などでは書かれていると思いますが，私も天下り的に素直に納得できないほうです．x1はk0,n1,n0の関数なのですから当然
x1(k0, n1, n0)ではなくx1(n1, k0, n0)とかk0が一番最上位でなくてもいいじゃないかというご質問でよろしいでしょうか？

私もそこまで完全によく分かっていないので，うまく説明できませんが，プログラムを作成したときに例えば，
x[9]という配列があるとします．これは_D(10進表示)，_B(2進表示)として
9_D=(1001)_B
ですから，x(9)=x(1,0,0,1)とでも表示できますでしょうか．
例えばC言語のプログラムではfor(i=0; i<N; i++)などのようにインクリメントして配列を繰り上げていきますが，x(1,0,0,1)というのをx(9)というのに対応させれば自然ですが，位置をずらして例えばx(9)=x(0,0,1,1)という対応をさせてしまうとプログラムに新たな規則を作ったりしなければならず，難解なプログラムになってしまうと考えられます．
以上，単純に説明しましたが，素人の考えなので，誤りもあるかと思います．その際は識者の皆様のご指摘を頂きたく思います．
何か分からないことがあったら指摘下さい．共に考えていきましょう．

この回答へのお礼

回答ありがとうございます。

先も申しましたが、疑問なのは「パラメーターの位置」と「パラメーターの種類」です。
今回のRossanaさんのご回答は、「位置」に関することですよね。確かに桁の位置を変えると、x(9)＝x(1,0,0,1)もx(9)＝x(0,0,1,1)もあり得ますよね。
ただ、x(9)＝x(1,0,0,1)に対して、x(9)≠x(0,0,1,1)になってしまいますよね。
これは、回転子の周期性・対称性を考慮に入れるからだと僕は解釈しています（しっくりきませんが…苦笑）。例えば８点FFTの場合、1段目のバタフライ演算において、バタフライ演算の対象となる
2点は、４つ飛びの値になります（x[0]とx[4]、x[1]とx[5]、x[2]とx[6]、x[3]とx[7]）。すなわちこれは、最上位桁に依存しています（最上位桁の重みは、2^2=4であるから）。2段目のバタフライ演算の対象は、２つ飛びの値になります（x1[0]とx1[2]、x1[1]とx1[3]、x1[4]とx1[6]、x1[5]とx1[7]）。すなわちこれは、最上位桁より１つ下の桁に依存しています。同様に3段目のバタフライ演算に対象は、１つ飛びの値となり最下位桁に依存します。このように、周期性・対称性を考慮に入れると納得はできるのですが、もし、「周期性・対称性」という概念がなかったら、……。となってしまいます。（涙）

「パラメーターの種類」に関して、僕は消去法（？）で考えてみたんですけど、（出力の最上位桁に）k0が来ないと絶対におかしくなるんですよ。例えば、k0ではなくk1,n0,n1を入れてたのですが、おかしくなります。結論として、k0が出力の最上位桁に来るしかないんですよねぇ…。

もっと論理的に説明できないでしょうか…。

お礼日時：2005/02/23 08:30

No.7 回答者： Rossana 回答日時： 2005/02/24 23:15

論理的にと言われると難しいですね…．

私のキャパシティではご期待に添えられないかもしれません（悲）．
この質問はもう半年前くらいなので，この議論を見ている方はいないのでは（悲）と思います．
もう一度，新たにご質問されるのもいいと思います．
以下は，理解の助けになるでしょうか．
☆x1(k0, n1, n0)=x0(n2, n1, 0)＋x0(n2, n1, 1)exp[-jπk0]について
x1(2n1+n0)=x1(0, n1, n0)＝x0(0, n1, n0)+x0(1, n1, n0) … x1(0),x1(1),x1(2),x1(3)に対応
x1(4+2n1+n0)=x1(1, n1, n0)＝x0(0, n1, n0)-x0(1, n1, n0) … x1(4),x1(5),x1(6),x1(7)に対応

☆x2(k0, k1, n0)=x1(k0, 0, n0)＋x1(k0, 1, n0)exp[-jπk1]exp[j(π/2)k0]について
x2(n0)=x2(0, 0, n0)=x1(0, 0, n0)＋x1(0, 1, n0) … x2(0),x2(1)に対応
x2(2+n0)=x2(0, 1, n0)=x1(0, 0, n0)-x1(0, 1, n0) … x2(2),x2(3)に対応
x2(4+n0)=x2(1, 0, n0)=x1(1, 0, n0)＋jx1(1, 1, n0) … x2(4),x2(5)に対応
x2(4+2+n0)=x2(1, 1, n0)=x1(1, 0, n0)-jx1(1, 1, n0) … x2(6),x2(7)に対応

☆x3(k0, k1, k2)=x2(k0, k1, 0)＋x2(k0, k1, 1)exp[-jπk2]exp[j(π/4)(2k1＋k0)]について
x3(0)=x3(0, 0, 0)=x2(0, 0, 0)＋x2(0, 0, 1) … x3(0)に対応
x3(1)=x3(0, 0, 1)=x2(0, 0, 0)-x2(0, 0, 1) … x3(1)に対応
x3(2)=x3(0, 1, 0)=x2(0, 1, 0)＋jx2(0, 1, 1) … x3(2)に対応
x3(2+1)=x3(0, 1, 1)=x2(0, 1, 0)-jx2(0, 1, 1) … x3(3)に対応
x3(4)=x3(1, 0, 0)=x2(1, 0, 0)+exp[j(π/4)]x2(1, 0, 1) … x3(4)に対応
x3(4+1)=x3(1, 0, 1)=x2(1, 0, 0)-exp[j(π/4)]x2(1, 0, 1) … x3(5)に対応
x3(4+2)=x3(1, 1, 0)=x2(1, 1, 0)+exp[j(3π/4)]x2(1, 1, 1) … x3(6)に対応
x3(4+2+1)=x3(1, 1, 1)=x2(1, 1, 0)-exp[j(3π/4)]x2(1, 1, 1) … x3(7)に対応

これを見る限り例えば，
x1(*,n1,n0)⇔x0(*,n1,n0)という対応
x2(*,*,n0)⇔x1(*,*,n0)という対応
x3(*,*,*)⇔x2(*,*,*)という対応
となっています．

うまく説明できませんね．なんか狐に摘まれた感じがします．
natsumitokaの命名された「消去法」なるもので考えるのはいい作戦ですね．私は単純に「種類」の方はx1(k1,n1,n0)などは
右辺にk1なるパラメータを含んでいないから単純にありえないと考えてしまったのですが，これも考えるとさらに分からなくなりますね．
「位置」の方に関して
x0(0,n1,n0)±x0(1,n1,n0)
との対応としては他にx1(n0,k0,n1)など3!-1=5通りありますが．

この回答へのお礼

またまたありがとうございます．
う～ん…，わかんないですねぇ…．
本当にすみません．
僕の力不足ですね…．

今のところ「消去法（？）」が一番答えに近いですかねぇ…．
また質問しなおしてみます．ありがとうございました．

お礼日時：2005/03/01 22:28