FFT中的Rader位反转算法

最近在思考fft中的bit-reverse（位反转）算法，给定一个数列，通过反转每项数字的二进制表示来创建一个新的数列，这个新的数列可以用于FFT中的蝶形算法中。

例如，从1到15：

0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111

反转二进制位就是：

1000 0100 1100 0010 1010 0110 1110 0001 1001 0101 1101 0011 1011 0111 1111

 

一种常见的位反转算法是Gold Rader bit reversal algorithm。对于这个算法，搜到的资料要么一笔带过，要么解释不对，所以我决定把我对这个算法的思路写下来，希望对同样思考这个算法的童鞋有所帮助。

Rader算法伪代码如下：

for i = 0 ... n − 2 do
k = n/2
if i < j then
swap g(i) and g(j)
end if
while k ≤ j do
j ⇐ j − k
k ⇐ k/2
end while
j ⇐ j + k
end for

代码里面，通过遍历i，找到i的反转j，并在需要交换i和j的项的值（交换是为了输出数列，这里我们并不关心，我们只关心j是怎么来的）。

假设n=16，i和j在遍历过程中的变化如下：

i：0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15（算法中15被优化掉了，我们这里为了便于理解也写上）

j：0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15

i二进制：0000 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111

j二进制：0000 1000 0100 1100 0010 1010 0110 1110 0001 1001 0101 1101 0011 1011 0111 1111

对于i，太好理解了，正常的加法从0每次加1，加到15。

而j（即i的按位反转），是怎么得到的呢？其实也是“加法”，也是从0开始每次加1，加到“15”，不过不是正常的加法罢了，j的加法规则（加1）如下：

从高位加1，并向低位进位，比如：

  111011

+1         

------------

   000111

=======================

   001011

 +1

-------------

   101011

观察上面两个例子，这个算法用语言描述起来就是：

a，从最高位往低位查找，碰到第一个0，这个0左边（高位）的1都变成0（即加1后依次向低位进位，进到低位出现0为止，进位的位置变成0）。

b，a中碰到第一个0所在的那个位置变成1.（完成向低位进位）

用符合人脑习惯的算法描述就是：

//  n是数字个数，一般是2的正整数次方
unsigned int Mask = n;
//  从最高位往最低位数查找，碰到第一个0，这个0左边（高位）的1都变成0，即依次向低位进位，进到低位出现0为止。
while (j & (Mask >>= 1))
//  把高位连续出现的1都变成0
    j&= ~Mask;
//  碰到第一个0所在的那个位置变成1
j |= Mask;

去除位运算的算法描述就是前面Gold Rader算法的伪代码形式了。