在实现FFT计算的时候,第一步要做的就是实现倒位序的实现,这里有一种算法,叫做雷德(Rader)算法。
下面假如使用A[I]存的是顺序位序,而B[J]存的是倒位序。I<J的时候需要变序,I>J的时候就不用,不然就白忙活了。
例如 N = 8 的时候,
倒位序 顺序 二进制表示 倒位序 顺序
0 0 000 000
4 1 100 001
2 2 010 010
6 3 110 011
1 4 001 100
5 5 101 101
3 6 011 110
7 7 111 111
由上面的表可以看出,按自然顺序排列的二进制数,其下面一个数总是比其上面一个数大1,即下面一个数是上面一个数在最低位加1并向高位进位而得到的。而倒位序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位进位而得到。
I、J都是从0开始,若已知某个倒位序J,要求下一个倒位序数,则应先判断J的最高位是否为0,这可与k=N/2相比较,因为N/2总是等于100..的。如果k>J,则J的最高位为0,只要把该位变为1(J与k=N/2相加即可),就得到下一个倒位序数;如果K<=J,则J的最高位为1,可将最高位变为0(J与k=N/2相减即可)。然后还需判断次高位,这可与k=N4相比较,若次高位为0,则需将它变为1(加N4即可)其他位不变,既得到下一个倒位序数;若次高位是1,则需将它也变为0。然后再判断下一位。。。。
下面用C代码实现雷德算法。
#include <stdio.h>
double x[16] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
double y[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
int N = 16;
void main()
{
int i,j,k;
double TempR,TempI;
int c = 0,m = 0;
for(j=0,i=0;i<N-1;i++) //这里实现了奇偶前后分开排序
{
printf("i = %d ",i);
printf("j = %d ",j);
printf("k = %d ",k);
if(i<j) //如果i<j,即进行变址
{
TempR = x[j];
TempI = y[j];
x[j] = x[i];
y[j] = y[i];
x[i] = TempR;
y[i] = TempI;
}
k = N/2; //求j的下一个倒位序
while(k<(j+1)) //如果k<=j,表示j的最高位为1
{
j = j-k; //把最高位变成0
k = k/2; //k/2,比较次高位,依次类推,逐个比较,直到某个位为0
}
j = j+k; //把0改为1
for(m=0;m<N/2;m++)
{
for(c=0;c<2;c++) //一行打印两个数
printf(" %10.7f+J %10.7f",x[2*m+c],y[2*m+c]);
printf(" ");
}
printf(" ");
}
}