博文
遗忘了的，FFT的一个小角落
||
快速傅立叶变换（FFT）中最常见的基2算法，有着蝶形的结构。基2的变换序数有着倒序的特点。所说的倒序，就是序数的2进制表达，变换后的序数恰好是其表达的中心镜像对称。如16点的基2变换 ，序数“11”二进制1011，那么变换后的序数在11这个位置就是“13”二进制1101。
蝶状图的序数的算法方式：一级分为奇数偶数两组，偶数组在上方，奇数组在下方。二级把原先分好的偶数组除2，再分奇数组与偶数组，原先的奇数组减一后除2，再分奇数组，偶数组。依次类推，直到分组无法再分。应用FFT许多年了，却没认真考虑过为什么它的序数最后是二进制的倒序。
以16点2基为例
一级：
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
二级：
0 2 4 6 8 10 12 14 1 3 5 7 9 11 13 15
三级：
0 4 8 12 2 6 10 14 1 5 9 13 3 7 11 15
四级：
0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15
一级与四级二进制表达分别：
0： 0000 0： 0000
1： 0001 8： 1000
2： 0010 4： 0100
3： 0011 12：1100
4： 0100 2： 0010
5： 0101 10：1010
6： 0110 6： 0110
7： 0111 14：1110
8： 1000 1： 0001
9： 1001 9： 1001
10：1010 5： 0101
11：1011 13：1101
12：1100 3： 0011
13：1101 11：1011
14：1110 7： 0111
15：1111 15：1111
上述的序数的算法结构，用MATLAB语言描述大约如下：
function h=ap(h,n) % n: 2^n
for i=0:n-1
T=2^(n-i);
for j=0:(2^i-1)
h((1+j*T):(j*T+2^(n-i)))=ac(h,1+j*T,n-i);
end
end
end
function tmp=ac(h,start,n)
tmp=zeros(1,2^n);
for i=1:2^(n-1)
tmp(i)=h(start+2*(i-1));
tmp(i+2^(n-1))=h(start+2*i-1);
end
end
>> h=0:15;
>> ap(h,4)
ans =
0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15
倒序方式的MATLAB程序：>> N=4;
>> p(1:2^N,N:-1:1)=dec2bin(0:(2^N-1));
>> re=bin2dec(char(p));
re =
0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15
昨儿，莫名想到这个倒序序数的事（也许是 11.11二进制形式和蝴蝶的各种联想吧）。之后动笔排了好些次，也就得到为什么是倒序的原因了。
C++程序描述倒序算法如下：
int sequ( int *re, unsigned int n )
{
unsigned int i;
int x = 0 ;
do
{
for( i = 0;i < n; i++ )
{
re[ x ] |= ( ( x >> i ) & 0x01 ) << ( n - i - 1 );
}
}while( ++x < ( 1 << n ) );
return 0;
}
当一个序数x，按奇偶分组（x<16）：
最后一级，相对于频域的序数的位置是红还是黑，取决于(x>>3)&1，也就是变换后位置(H[x]>>0)&1，第三级取决于(x>>2)&1，也就是变换后位置(H[x]>>1)&1；第二级取决于(x>>1)&1，也就是变换后的位置(H[x]>>2 )&1；第一级是红是黑取决于(x>>0)&1，也就是变换后的位置(H[x]>>3)&1。同样的，其它2基变换的序数也满足这样的规律，由0到2^n-1的顺序序数依次奇偶分组后的序数必然是倒序。
以此文悼念那些业已失去的，浑浑噩噩的日子。
https://blog.sciencenet.cn/blog-567861-631966.html
上一篇：关于整数拆分的几个Matlab函数
下一篇：求正弦函数值
