Application Report
ZHCA414 – Jan. 2012
基于TMS320C64x+DSP的FFT实现
冯华亮
DSP System
ABSTRACT
This application note introduces the implementation of FFT (Fast Fourier Transform) on TMS320C64x+ DSP family, the performance is also discussed.
摘要
本文介绍基于TI TMS320C64x+ DSP 的FFT（快速傅立叶变换）的实现，并讨论相关性能。
1
ZHCA414
Contents
1
快速傅立叶变换介绍 ....................................................................................................................... 3 1.1 基2 FFT ................................................................................................................................. 3 1.2 基4 FFT ................................................................................................................................. 6 1.3 混合基4和基2 FFT ............................................................................................................. 10 1.4 混合基FFT的常规C语言实现 ............................................................................................. 13 1.5 快速傅立叶逆变换(IFFT) ....................................................................................................... 17 2 TMS320C64x+ DSP简介及FFT相关特性 ................................................................................... 17 3 基于TMS320C64x+的FFT实现 .................................................................................................. 21
3.1 C64x+ DSP库里的FFT函数 ................................................................................................ 21 3.2 对旋转因子访问的优化 .......................................................................................................... 22 3.3 蝶形运算的优化 .................................................................................................................... 25 3.4 Cache冲突的优化 ................................................................................................................ 28 4 执行速度指标 ............................................................................................................................... 29 5 数据缩放和精度 ............................................................................................................................ 33 References ......................................................................................................................................... 35 Figures
Figure 1. 基2，N=8 点FFT .............................................................................................................. 6 Figure 2. 基4，N=16点FFT ............................................................................................................ 9 Figure 3. N=32混合基4和基2 FFT ............................................................................................... 12 Figure 4. TMS320C64x+ DSP内核框图 ......................................................................................... 18 Figure 5. C64x+ DSP上1024点基2 FFT比特反转 ....................................................................... 20 Figure 6. C64x+ DSP上1024点基4 FFT比特反转 ....................................................................... 21 Figure 7. FFT函数命名规则............................................................................................................ 22 Figure 8. 优化实现中的基本蝶形运算 .............................................................................................. 26 Figure 9. FFT函数执行的指令周期数 ............................................................................................. 30 Figure 10. IFFT函数执行的指令周期数 ............................................................................................ 32
Tables
基2，8点FFT序号比特反正 ............................................................................................ 6 基4，16点FFT序号比特反转 .......................................................................................... 9 旋转因子组和FFT级之间的关系 ..................................................................................... 12 32点混合基FFT序号比特反转 ........................................................................................ 12 C64x+ DSP上用于FFT的指令 ....................................................................................... 18 C64x+ DSP库里的FFT函数 ........................................................................................... 21 可连续访问的旋转因子表 ................................................................................................. 23 FFT函数执行的指令周期数 ............................................................................................. 29 IFFT函数执行的指令周期数 ............................................................................................ 31 Cache “touch”对DSP_fft16x16性能的影响 .................................................................. 32 数据动态缩放的DSP_fft16x16的性能 ............................................................................ 34
Table 1. Table 2. Table 3. Table 4. Table 5. Table 6. Table 7. Table 8. Table 9. Table 10. Table 11.
2 基于TMS320C64x+DSP的FFT实现
ZHCA414
1 快速傅立叶变换介绍
序列x(n), n = 0...N ?1的离散傅立叶变换 (Discrete Fourier Transform, DFT) X(k), k = 0...N ?1 可定义为:
knX(k)=DFTN(x(n))=∑x(n)WN,n=0N?1k=0,....,N?1
kn其中，WN=e?j(2π/N)kn是旋转因子。
反离散傅立叶(Inverse DFT, IDFT)可定义为: 其中
?kn=ej(2π/N)kn WN1x(n)=N∑X(k)Wk=0N?1?knN,n=0,....,N?1
它们的输入输出序列都是复数。如果我们定义一个基本的运算单位为两个复数相乘再相加，则直接根据以上公式进行DFT计算需要N2次基本运算。当N比较大时，N2会变得非常大。这使得直接的DFT计算方法对很多应用来说都不现实。
1.1 基2 FFT
然而，如果N是2的整数次方，一种快速傅立叶变换（Fast Fourier Transform, FFT）算法可显著的减少运算量。基2 FFT算法利用了旋转因子的以下周期性特性:
0=e?j(2π/N)0=e?j(0)=cos(0)+jsin(0)=1 WNkN/2WN=e?j(2π/N)kN/2=e?j(πk)=(cos(?π)+jsin(?π))k=(?1)k kn2kn=e?j(2π/N)2kn=e?j(2π/(N/2))kn=WNWN/2
利用这些特性，可以把N点DFT分解为以下两个N/2点DFT:
基于TMS320C64x+DSP的FFT实现 3
ZHCA414
X(k)=∑x(n)W
n=0N/2?1N?1
nkN
N/2?1
=
∑x(n)W
n=0
nkN
+
n=N/2
∑x(n)W
N?1
nkN
=
∑[x(n)W
n=0
nkN(n+N/2)k
+x(n+N/2)WN]
N/2?1
===
∑[x(n)W
n=0
nkNkN/2nk
+x(n+N/2)WNWN]
N/2?1n=0
∑[x(n)W
n=0
nkNnk
+(?1)kx(n+N/2)WN]
N/2?1
∑[x(n)+(?1)
k
nk
x(n+N/2)]WN
让我们把输出序列X(k), k= 0,…, N-1 分解成连个序列：
偶数序列: X(2r), r= 0,…,N/2-1
N/2?1X(2r)==
∑[x(n)+(?1)n=02r2nrx(n+N/2)]WNN/2?1n=0∑[x(n)+x(n+N/2)]W∑[x(n)+x(n+N/2)]Wn=02nrN
nrN/2N/2?1==DFTN/2(x(n)+x(n+N/2))
奇数序列: X(2r+1), r= 0,…,N/2-1
N/2?1X(2r+1)=N/2?1∑[x(n)+(?1)n=0(2r+1)n*(2r+1)x(n+N/2)]WN=
∑[x(n)?x(n+N/2)]WWnNn=02nrN
N/2?1=∑{[x(n)?x(n+N/2)]Wn=0nNnr}WN/2n)=DFTN/2((x(n)?x(n+N/2))WN按以上方法反复的分解N点DFT成2/N点DFT直到N=2，使得N点傅立叶变换的运算复杂度由
原来的N2降到Nlog2N，这是非常显著的改进。这个算法也被称为频域抽取(Decimate-In-Frequency, DIF)FFT算法。图1演示了N=8点FFT的分解运算过程。上面的两个箭头（每个箭头上都有一个数字“1”）表示两个数相加，下面的两个箭头（一个箭头上数字“1”，一个箭头
nn
上数字“-1”）表示两个数相减。后面的WN表示加减的结果再与旋转因子WN相乘。
4 基于TMS320C64x+DSP的FFT实现
ZHCA414
X(0)X(1)4 pointsEven11DFTSequenceX(2)X(3)X(4)1W08-1W1X(5)84 pointsOddX(6)W28DFTSequenceX(7)W38X(0)11112 pointsEven of evenDFTSequenceX(1)X(2)1W004(W8)-12 pointsOdd of evenSequenceX(3)W1(W248)DFTX(4)1W08-12 pointsEven of oddW1X(5)8DFTSequenceX(6)W28W04(W08)2 pointsOdd of oddX(7)W38W1(W248)DFTSequenceX(0)111111X(0)X(1)1-1X(4)X(2)1W04(W08)-1X(2)X(3)W124(W8)X(6)X(4)1W08-1X(1)W1X(5)8X(5)X(6)W208W4(W08)X(3)X(7)W38W14(W28)X(7) 基于TMS320C64x+DSP的FFT实现 5
ZHCA414
Figure 1. 基2，N=8 点FFT
从图1可以看出，输出数据的顺序被打乱了。这种乱序其实有规律，我们把顺序的序号用二进制数列在表1中的左边，把乱序的序号用二进制数列在表1中的右边。从表中我们可以看出，乱序序号的二进制码可由顺序序号的二进制码反转得到，所以，这种规律被叫做比特反转。
Table 1.
Normal order of
index n 0 1 2 3 4 5 6 7
基2，8点FFT序号比特反正
Reversed bits of
index n 000 100 010 110 001 101 011 111
Bit-reversed of order
index n 0 4 2 6 1 5 3 7
Binary bits of index n 000 001 010 011 100 101 110 111
1.2 基4 FFT
如果FFT点数N是4的整数次方，采用基4FFT算法可以进一步减少运算量。基4 FFT算法算法利用了旋转因子的以下周期性特性::
kN/4WN=e?j(2π/N)kN/4=e?j(πk/2)=(cos(?π/2)+jsin(?π/2))k=(?j)k 2kN/4WN=e?j(2π/N)2kN/4=e?j(πk)=(cos(?π)+jsin(?π))k=(?1)k
3kN/4WN=e?j(2π/N)3kN/4=e?j(3πk/2)=(cos(?3π/2)+jsin(?3π/2))k=(j)k nk4nkWN=e?j(2π/N)4kn=e?j(2π/(N/4))kn=WN/4
利用这些特性，可以把N点DFT分解为以下4个N/4点DFT:
X(k)=∑x(n)Wn=0N/4?1N?1nkNN/4?1=∑x(n)Wn=0nkN2N/4?1+n=N/4∑x(n)WnkN3N/4?1+n=2N/4∑x(n)WnkN+n=3N/4∑x(n)WN?1nkN=
∑[x(n)Wn=0n=0nkN(n+N/4)k(n+2N/4)k(n+3N/4)k+x(n+N/4)WN]+x(n+2N/4)WN+x(n+3N/4)WNN/4?1
kN/4N==∑[x(n)+x(n+N/4)Wkn=02kN/43kN/4nk+x(n+2N/4)WN+x(n+3N/4)WN]WNN/4?1∑[x(n)+(?j)nkx(n+N/4)+(?1)kx(n+2N/4)+(j)kx(n+3N/4)]WN让我们把输出序列X(k), k= 0,…, N-1 分解成四个序列, X(4r), X(4r+1), X(4r+2), X(4r+3), r= 0,…,N/4-1
6 基于TMS320C64x+DSP的FFT实现
ZHCA414
N/4?1X(4r)==
∑[x(n)+(?j)n=04r4nrx(n+N/4)+(?1)4rx(n+2N/4)+(j)4rx(n+3N/4)]WNN/4?1n=0∑[x(n)+x(n+N/4)+x(n+2N/4)+x(n+3N/4)]WnrN/4=DFTN/4[x(n)+x(n+N/4)+x(n+2N/4)+x(n+3N/4)]=DFTN/4([x(n)+x(n+2N/4)]+[x(n+N/4)+x(n+3N/4)])=DFTN/4(A1(n)+B1(n))
其中,
A1(n)=x(n)+x(n+2N/4); B1(n)=x(n+N/4)+x(n+3N/4)
N/4?1X(4r+1)=N/4?1∑[x(n)+(?j)n=0(4r+1)(4r+1)nx(n+N/4)+(?1)(4r+1)x(n+2N/4)+(j)(4r+1)x(n+3N/4)]WN==∑[x(n)?jx(n+N/4)?x(n+2N/4)+jx(n+3N/4)]Wn=0nN4nrWNN/4?1n=0∑[x(n)?jx(n+N/4)?x(n+2N/4)+jx(n+3N/4)]WnNnrWN/4n)=DFTN/4([x(n)?jx(n+N/4)?x(n+2N/4)+jx(n+3N/4)]WNn)=DFTN/4({[x(n)?x(n+2N/4)]?j[x(n+N/4)?x(n+3N/4)]}WNn)=DFTN/4([A2(n)?jB2(n)]WN 其中,
A2(n)=x(n)?x(n+2N/4); B2(n)=x(n+N/4)?x(n+3N/4)
基于TMS320C64x+DSP的FFT实现 7
ZHCA414
N/4?1X(4r+2)=N/4?1∑[x(n)+(?j)n=0(4r+2)(4r+2)nx(n+N/4)+(?1)(4r+2)x(n+2N/4)+(j)(4r+2)x(n+3N/4)]WN==∑[x(n)?x(n+N/4)+x(n+2N/4)?x(n+3N/4)]Wn=02nN4nrWNN/4?1n=0∑[x(n)?x(n+N/4)+x(n+2N/4)?x(n+3N/4)]W2nNnrWN/42n)=DFTN/4([x(n)?x(n+N/4)+x(n+2N/4)?x(n+3N/4)]WN2n)=DFTN/4({[x(n)+x(n+2N/4)]?[x(n+N/4)+x(n+3N/4)]}WN2n)=DFTN/4([A1(n)?B1(n)]WN 其中,
A1(n)=x(n)+x(n+2N/4); B1(n)=x(n+N/4)+x(n+3N/4)
N/4?1X(4r+3)=N/4?1∑[x(n)+(?j)n=0(4r+3)(4r+3)nx(n+N/4)+(?1)(4r+3)x(n+2N/4)+(j)(4r+3)x(n+3N/4)]WN==∑[x(n)+jx(n+N/4)?x(n+2N/4)?jx(n+3N/4)]Wn=03nN4nrWNN/4?1n=0∑[x(n)+jx(n+N/4)?x(n+2N/4)?jx(n+3N/4)]W3nNnrWN/43n)=DFTN/4([x(n)+jx(n+N/4)?x(n+2N/4)?jx(n+3N/4)]WN3n)=DFTN/4({[x(n)?x(n+2N/4)]+j[x(n+N/4)?x(n+3N/4)]}WN3n)=DFTN/4([A2(n)+jB2(n)]WN其中,
A2(n)=x(n)?x(n+2N/4); B2(n)=x(n+N/4)?x(n+3N/4)
按以上方法反复的分解N点DFT成N/4点DFT直到N=4，使得N点傅立叶变换的运算复杂度由
原来的N2降到Nlog4N。图2演示了N=16点FFT的分解运算过程。
8 基于TMS320C64x+DSP的FFT实现
ZHCA414
X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)X(8)X(9)X(10)X(11)X(12)X(13)X(14)X(15)1j-1-j1-11-11-j-1j0W1611111-j-1j1111X(0)X(4)1-11-1X(8)X(12)X(1)X(5)1j-1-j1W162W161-j-1j11113W161j-1-j1-11-1X(9)X(13)X(2)X(6)0W162W164W161-j-1j1111W6161j-1-j1-11-1X(10)X(14)X(3)X(7)0W163W166W161-j-1j11119W161j-1-j1-11-1X(11)X(15)
Figure 2. 基4，N=16点FFT
从图2可以看出，输出数据的顺序也被打乱了。这种乱序其实有规律，我们把顺序的序号用二进制数列在表1中的左边，把乱序的序号用二进制数列在表1中的右边。从表中我们可以看出，乱序序号的二进制码可由顺序序号的二进制码以2比特为单位反转得到。
Table 2.
Normal order of
index n 0 1
基4，16点FFT序号比特反转
Reversed bits of
index n 00 00 01 00
Bit-reversed of order
index n 0 4
Binary bits of index n 00 00 00 01
基于TMS320C64x+DSP的FFT实现 9
ZHCA414
2 3 4 5 6 7 8 9 10 11 12 13 14 15
00 10 00 11 01 00 01 01 01 10 01 11 10 00 10 01 10 10 10 11 11 00 11 01 11 10 11 11
10 00 11 00 00 01 01 01 10 01 11 01 00 10 01 10 10 10 11 10 00 11 01 11 10 11 11 11
8 12 1 5 9 13 2 6 10 14 3 7 11 15
1.3 混合基4和基2 FFT
如果N= 4M*2, 傅立叶变换可被分解成log4N= M个基4 FFT和最后一级的一个基2 FFT. 图3演示了N=32点FFT的分解运算过程。
10 基于TMS320C64x+DSP的FFT实现
ZHCA414
X(0)X(0)X(1)X(16)X(2)W0(W0832)X(4)X(3)W1(W4832)X(20)W0(W0832)X(4)X(8)X(5)W2W84(32)X(24)X(6)W0(W0832)X(12)X(7)W3(W12832)X(28)X(8)W032X(1)X(9)W132X(17)X(10)W232W0W08(32)X(5)X(11)W31432W8(W32)X(21)X(12)W432W0(W0832)X(9)28X(13)W532W4(W32)X(25)0X(14)W632W(W0832)X(13)X(15)W7W312328(W32)X(29)X(16)W032X(2)W2X(17)32X(18)X(18)W4032W(W0832)X(6)W6X(19)32W148(W32)X(22)X(20)W80032W8(W32)X(10)W10X(21)32W2(W8432)X(26)X(22)W1232W008(W32)X(14)X(23)W1432W3(W12832)X(30)X(24)W032X(3)X(25)W332X(19)X(26)W6W0(W032832)X(7)X(27)W932W148(W32)X(23)X(28)W1232W0(W0832)X(11)X(29)W1532W284(W32)X(27)X(30)W1832W0(W0832)X(15)X(31)W21332W(W12832)X(31) 基于TMS320C64x+DSP的FFT实现
11
ZHCA414
Figure 3. N=32混合基4和基2 FFT
值得注意的是，在第一级中，每个蝶形运算的旋转因子都不一样；在第二级，四组蝶形运算使用相同的旋转因子。下一级每一组中的不同蝶运算仍使用不同的旋转因子。重复以上规则，可得到下面的表：
Table 3.
级数 1 2 … log4N
旋转因子组和FFT级之间的关系
蝶形组的个数 1 4 … N/4
蝶的总数数 N/4 N/4 … N/4
每组中蝶的个数 N/4 N/16 … 1
FFT的实现代码利用以上关系来减少旋转因子的访问，即每一级中相同的旋转因子只读取一次，然后给所有使用这个相同旋转因子的蝶形运算用。
由图3可看出，数据输出的顺序被打乱了。表4可帮我们理解这种乱序的规则。它将正常序号二进制编码以2比特为单位反转。因为总比特个数不是偶数，最后的一个比特单独反转。
Table 4.
Normal order of
index n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
32点混合基FFT序号比特反转
Reversed bits of
index n
0 00 00 1 00 00 0 01 00 1 01 00 0 10 00 1 10 00 0 11 00 1 11 00 0 00 01 1 00 01 0 01 01 1 01 01 0 10 01 1 10 01 0 11 01 1 11 01 0 00 10 1 00 10 0 01 10 1 01 10
Binary bits of index n
00 00 0 00 00 1 00 01 0 00 01 1 00 10 0 00 10 1 00 11 0 00 11 1 01 00 0 01 00 1 01 01 0 01 01 1 01 10 0 01 10 1 01 11 0 01 11 1 10 00 0 10 00 1 10 01 0 10 01 1 Bit-reversed of order
index n 0 16 4 20 8 24 12 28 1 17 5 21 9 25 13 29 2 18 6 22
12 基于TMS320C64x+DSP的FFT实现
ZHCA414
20 21 22 23 24 25 26 27 28 29 30 31
10 10 0 10 10 1 10 11 0 10 11 1 11 00 0 11 00 1 11 01 0 11 01 1 11 10 0 11 10 1 11 11 0 11 11 1
0 10 10 1 10 10 0 11 10 1 11 10 0 00 11 1 00 11 0 01 11 1 01 11 0 10 11 1 10 11 0 11 11 1 11 11
10 26 14 30 3 19 7 23 11 27 15 31
1.4 混合基FFT的常规C语言实现
下面的C代码是混合基FFT的常规C语言实现，它用FFT算法的主要发明人Cooley和Turkey命名。它要求正常顺序的输入x[]，产生正常顺序的输出y[]。输入，输出和旋转因子wn[]都是复数，实部和虚部存在数组中相邻的单元，实部在偶数单元，虚部在奇数单元。所有数据数据都是16位的Q.15格式的小数。
基于TMS320C64x+DSP的FFT实现 13
ZHCA414
void CooleyTukeyFft16x16(int N, short wn[], short x[], short y[]) { int n1, n2, ie, ia1, ia2, ia3; int i0, i1, i2, i3, j, k;
short co1, co2, co3, si1, si2, si3; short xt0, yt0, xt1, yt1, xt2, yt2; short xh0, xh1, xh20, xh21, xl0, xl1,xl20,xl21;
n2 = N; ie = 1;
/*Performs log4(N)-1 stages of radix4 FFT */ for (k = N; k > 4; k >>= 2) { n1 = n2; n2 >>= 2; ia1 = 0;
/*loop through butterflies with different twiddle factors j * i0 = N/4 */ for (j = 0; j < n2; j++)
{ ia2 = ia1 + ia1; ia3 = ia2 + ia1;
co1 = wn[2 * ia1 ]; si1 = wn[2 * ia1 + 1]; co2 = wn[2 * ia2 ]; si2 = wn[2 * ia2 + 1]; co3 = wn[2 * ia3 ]; si3 = wn[2 * ia3 + 1]; ia1 = ia1 + ie;
/*loop through the groups with same twiddle factors*/ for (i0 = j; i0< N; i0 += n1) { i1 = i0 + n2; i2 = i1 + n2; i3 = i2 + n2;
xh0 = x[2 * i0 ] + x[2 * i2 ]; //calculate real(A1) xh1 = x[2 * i0 + 1] + x[2 * i2 + 1]; //calculate imag(A1) xl0 = x[2 * i0 ] - x[2 * i2 ]; //calculate real(A2) xl1 = x[2 * i0 + 1] - x[2 * i2 + 1]; //calculate imag(A2)
xh20 = x[2 * i1 ] + x[2 * i3 ]; //calculate real(B1) xh21 = x[2 * i1 + 1] + x[2 * i3 + 1]; //calculate imag(B1) xl20 = x[2 * i1 ] - x[2 * i3 ]; //calculate real(B2) xl21 = x[2 * i1 + 1] - x[2 * i3 + 1]; //calculate imag(B2)
x[2 * i0 ] = (xh0 + xh20); //calculate real(A1+B1) x[2 * i0 + 1] = (xh1 + xh21); //calculate imag(A1+B1)
xt0 = xh0 - xh20; //calculate real(A1-B1) yt0 = xh1 - xh21; //calculate imag(A1-B1) xt1 = xl0 + xl21; //calculate real(A2-jB2) yt1 = xl1 - xl20; //calculate imag(A2-jB2) yt2 = xl1 + xl20; //calculate imag(A2+jB2) xt2 = xl0 - xl21; //calculate real(A2+jB2)
x[2 * i1 ] = (xt1 * co1 + yt1 * si1) >> 15;
14 基于TMS320C64x+DSP的FFT实现
ZHCA414
x[2 * i1 + 1] = (yt1 * co1 - xt1 * si1) >> 15; x[2 * i2 ] = (xt0 * co2 + yt0 * si2) >> 15; x[2 * i2 + 1] = (yt0 * co2 - xt0 * si2) >> 15; x[2 * i3 ] = (xt2 * co3 + yt2 * si3) >> 15; x[2 * i3 + 1] = (yt2 * co3 - xt2 * si3) >> 15; }
}
ie <<= 2; }
/*performs either a radix2 or radix4 transform on the */
/*last stage depending on \ /*then this last stage is also a radix4 transform, otherwise it is a */ /*radix2 transform.*/ if(4==k) {
i0= 0;
for(j= 0; j< N/4; j++) {
i1= i0+ 1; i2= i1+ 1; i3= i2+ 1;
xh0 = x[2 * i0 ] + x[2 * i2 ]; //calculate real(A1) xh1 = x[2 * i0 + 1] + x[2 * i2 + 1]; //calculate imag(A1) xl0 = x[2 * i0 ] - x[2 * i2 ]; //calculate real(A2) xl1 = x[2 * i0 + 1] - x[2 * i2 + 1]; //calculate imag(A2)
xh20 = x[2 * i1 ] + x[2 * i3 ]; //calculate real(B1) xh21 = x[2 * i1 + 1] + x[2 * i3 + 1]; //calculate imag(B1) xl20 = x[2 * i1 ] - x[2 * i3 ]; //calculate real(B2) xl21 = x[2 * i1 + 1] - x[2 * i3 + 1]; //calculate imag(B2)
x[2 * i0 ] = xh0 + xh20; //calculate real(A1+B1) x[2 * i0 + 1] = xh1 + xh21; //calculate imag(A1+B1) x[2 * i2 ] = xh0 - xh20; //calculate real(A1-B1) x[2 * i2 + 1] = xh1 - xh21; //calculate imag(A1-B1) x[2 * i1 ] = xl0 + xl21; //calculate real(A2-jB2) x[2 * i1 + 1] = xl1 - xl20; //calculate imag(A2-jB2) x[2 * i3 + 1] = xl1 + xl20; //calculate imag(A2+jB2) x[2 * i3 ] = xl0 - xl21; //calculate real(A2+jB2) i0+= 4; } }
else if(2==k) {
i0= 0;
for(j= 0; j< N/2; j++) {
i1= i0+ 1;
xh0 = x[2 * i0 ] + x[2 * i1 ]; xh1 = x[2 * i0 + 1] + x[2 * i1 + 1]; xl0 = x[2 * i0 ] - x[2 * i1 ]; xl1 = x[2 * i0 + 1] - x[2 * i1 + 1];
x[2*i0]= xh0; x[2*i0+ 1]= xh1;
基于TMS320C64x+DSP的FFT实现 15
ZHCA414
x[2*i1]= xl0; x[2*i1+ 1]= xl1;
i0+= 2; } }
//reverse the output order for (k = 0; k < N; k++) {
j= bitReverse(k);
y[2*j]= x[2*k];
y[2*j+1]= x[2*k+1]; }
} 前面的基4 FFT由三个主循环实现。最外层的”k”代表FFT的级，第二层的”j”代表不同的旋转因子，最内的”i0”代表使用相同旋转因子的不同的蝶。最后一级可能是基4或基2，放在最后单独实现
根据Table 3的关系我们可以得出以下规律:
a) 最内层循环次数从1到N/4变化，每级扩大4倍。 b) 第二层循环次数从N/4到1变化，每级缩小1/4。 因此，规律a)和b)正好相反。
c) 如果最内层循环和第二层循环合并成一层循环，则循环次数固定为N/4。在DSP上的优化实现正是这样做的。
FFT的旋转因子按如下公式产生：
nWN=e?j2π*n/N=cos(?2π*n/N)+jsin(?2π*n/N)=cos(2π*n/N)?jsin(2π*n/N)
生成旋转因子wn[]的代码如下：
void gen_wn_16(short *w, int N) {
int n;
float M = 32767;
float PI = 3.14159265358979323846;
for (n = 0; n < N; n++) {
w[2*n] = M * cos(2.0 * PI * n / N); w[2*n+1] = M * sin (2.0 * PI * n / N); } }
请注意，保存在旋转因子表中的旋转因子的虚部是sine而不是-sine，“sine”前面的负号在FFT复数乘法运算时被补进去。
上面FFT实现中，每个蝶形运算所用的三个旋转因子在存储空间上是不连续的，这会降低访问的效率。在DSP的优化实现中，旋转因子表被重新组织，从而可以连续访问。
16 基于TMS320C64x+DSP的FFT实现
ZHCA414
1.5 快速傅立叶逆变换(IFFT)
逆向DFT也可以类似地按前面介绍的方法来实现。IFFT和FFT实现上的主要区别是： 1, 输出结果除N
kn?kn2, 旋转因子WN变成WN
另外，IFFT也可以用FFT函数来计算：
1x(n)=NN?1k=0∑X(k)W?knN1N?1*1kn*=[∑X(k)WN]={DFT[X*(k)]}* Nk=0N实现方法是，先把输入数据取复共轭，执行FFT，再对结果取复共轭并除N。但，如果输入和输
出数据的复共轭计算不方便的话，就需要用专门的代码来实现IFFT。 本文档仅专注于FFT的优化实现。 相同的方法可应用到IFFT的实现上。
2 TMS320C64x+ DSP简介及FFT相关特性
本节简介TMS320C64x+ DSP，探讨DSP内为FFT实现而做的增强功能。
TMS320C64x+ DSP是Texas Instruments最新的DSP系列。它运行速度可达到1GHz，有8个运算单元，因此每个周期可执行8条指令，每秒钟最高可执行8G个16bitsx16bits的乘累加运算。下图是C64x+ DSP核对框图，更详细的信息请参阅”TMS320C64x C64x+ DSP CPU and Instruction Set Reference Guide (spru732)”.
基于TMS320C64x+DSP的FFT实现 17
ZHCA414
256 bits program bus C64x+ Core Instruction Fetch SPLOOP Buffer Control Registers In-Circuit Emulation Data Path 2 B Register File B0 –B31 D1 D2 + + M2 S2 x x x x + + + + Interrupt Control Instruction Dispatch Instruction Decode Data Path 1 A Register File A0 – A31 L1 + + + + S1 + + + + M1 L2 + + + + x x x x x x + + x x 64-Bit Data Bus 64-Bit Data Bus 32KB L1D Cache 256-Bit Data Bus L2 SRAM Figure 4. TMS320C64x+ DSP内核框图 C64x+有两组64-bit数据总线，可直接访问32KB一级数据存储器(L1D)。L1D中的0～32KB可以被配置为cache，其它部分可被用作普通RAM。一级存储器的速度和内核一样，二级(L2)存储器的速度是DSP核的一半。
为了充分利用64-bit数据总线，算法的数据访问要被优化成每次访问64-bit(double word，双字)；为了充分利用L1D cache，数据应该尽量的被连续访问。
C64x+ DSP针对FFT做了专门的增强，下表总结了主要的对FFT运算有帮助的指令：
Table 5.
FFT运算 蝶形运算中的加法和减法运算
C64x+ 指令 ADD2
C64x+ DSP上用于FFT的指令
说明
两个16-Bit整数加。用于16 bit复数加： (a+jb) + (c+jd)= (a+c)+ j(b+d)
SUB2
两个16 Bit整数减。用于16 bit复数减：
18 基于TMS320C64x+DSP的FFT实现
ZHCA414
(a+jb) - (c+jd)= (a-c)+ j(b-d)
ADDSUB2
对相同输入数据的并行的ADD2和SUB2运算。用于16 bit complex加和减：
(a+jb) + (c+jd)= (a+c)+ j(b+d) (a+jb) - (c+jd)= (a-c)+ j(b-d)
ADDSUB AVG2
并行的ADD和SUB运算，用于32-bit FFT。
两个16-Bit整数取平均。用于蝶形运算第一个加法输出的结果除2：
(a+jb) + (c+jd)= (a+c)/2+ j(b+d)/2
加减运算的结果乘旋转因子
CMPYR
16-bit复数乘法并对结果四舍五入和饱和: (a+jb) * (c+jd)= (ac-bd)>>16+ j(ad+bc)>>16
MPY2IR SMPY32
输出数据顺序比特反转
BITR DEAL SHFL NORM
两个16-Bit x 32-Bit运算，结果右移16位并四舍五入。用于32-bit数据， 16-bit旋转因子的FFT 32-Bit x 32-Bit乘法。用于32-bit FFT.
32-bit word比特顺序反转。
DEAL对32-bit数据做比特交织。SHFL是DEAL的逆运算。
确定数据的最高有效位。也就是计算log2(N).
关于这些指令的细节，请参阅”TMS320C64x C64x+ DSP CPU and Instruction Set Reference Guide (SPRU732)”。
基2FFT的比特反转可由以下代码实现：
power2= 30- _norm(N); //power2= log2(N) j= _bitr(k)>>(32-power2);
这里的k是原始序号，j是重排序的序号，N是FFT点数。
上面代码中的“_bitr”和“_norm”被称作指令BITR和NORM的intrinsic，这是在DSP C语言程序中指定使用某条DSP指令的方法。DSP编译器并不把它们当作函数调用，而是把它们编译成相应的指令。几乎每条DSP指令都有相应的intrinsic，这使得C语言程序的优化更便捷。
基于TMS320C64x+DSP的FFT实现 19
ZHCA414
20下图演示了在C64x+ DSP上1024点基2FFT比特反转的实现，图中每个字母代表一个比特。
310ab cd ef gh ij_bitr()310ji hg fe dc ba_bitr()>>(32-power2)310ji hg fe dc ba
Figure 5. C64x+ DSP上1024点基2 FFT比特反转
基4 FFT的比特反转可用以下代码实现：
power2= 30- _norm(N); //power2= log2(N)
j= _shfl(_rotl(_deal(_bitr(k)>>(32-power2)),16))
下图演示了在C64x+ DSP上1024点基4 FFT比特反转的实现。
310ab cd ef gh ij_bitr()310ji hg fe dc ba_bitr()>>(32-power2)310ji hg fe dc ba_deal()310j h f d b i g e c a_rotl(,16)310 i g e c a j h f d b_shfl()310ij gh ef cd ab
基于TMS320C64x+DSP的FFT实现
ZHCA414
Figure 6. C64x+ DSP上1024点基4 FFT比特反转
混合基FFT的比特反转可用以下代码实现：
power2= 30- _norm(N); //power2= log2(N) j= _shfl(_rotl(_deal(_bitr(k)>>(32-power2)),16)) if(last stage is radix-2) //swap the last single bit j= ((k&1)<<(power2-1))|(j&((1< 3 基于TMS320C64x+的FFT实现 TI提供的C64x+ DSP算法库中包含多种FFT的优化实现。所以这些库函数的可在C语言中调用。通过调用这些优化的函数，可以获得比标准C语言程序快得多的执行速度。TI同时也提供这些库函数的原代码，以方便客户在此基础上修改以得满足特殊要求的FFT。 3.1 C64x+ DSP库里的FFT函数 下表总结了C64x+ DSP库里的FFT函数： Table 6. 函数 void DSP_fft16x16( short *w, int nx, short *x, short *y) void DSP_fft16x16_imre( short *w, int nx, short *x, short *y) void DSP_fft16x16r( int nx, short *x, short *w, short *y, int radix, int offset, int n_max) void DSP_fft16x32( short *w, int nx, int *x, int *y) void DSP_fft32x32( int *w, int nx, int *x, int *y) void DSP_fft32x32s( int *w, int nx, int *x, int *y) C64x+ DSP库里的FFT函数 说明 16 bits 复输入输出数据，复数按先实部/后虚部交织存放 除最后一级外，每级的结果右移1并四舍五入 16 bits 复输入输出数据，复数按先虚部/后实部交织存放 除最后一级外，每级的结果右移1并四舍五入 16 bits 复输入输出数据，复数按先实部/后虚部交织存放 除最后一级外，每级的结果右移1并四舍五入 针对比较小的Cache优化 32 bits 复输入输出数据，复数按先实部/后虚部交织存放 32 bits 复输入输出数据，复数按先实部/后虚部交织存放 32 bits旋转因子 32 bits 复输入输出数据，复数按先实部/后虚部交织存放 32 bits旋转因子 除最后一级外，每级的结果右移1 void DSP_ifft16x16( 16 bits 复输入输出数据，复数按先实部/后虚部交织存放 基于TMS320C64x+DSP的FFT实现 21 ZHCA414 short *w, int nx, short *x, short *y) Void DSP_ifft16x16_imre( short *w, int nx, short *x, short *y) void DSP_ifft16x32( short *w, int nx, int *x, int *y) void DSP_ifft32x32( int *w, int nx, int *x, int *y) 除最后一级外，每级的结果右移1并四舍五入 16 bits 复输入输出数据，复数按先虚部/后实部交织存放 除最后一级外，每级的结果右移1并四舍五入 32 bits 复输入输出数据，复数按先实部/后虚部交织存放 32 bits 复输入输出数据，复数按先实部/后虚部交织存放 32 bits旋转因子 Figure 7 说明的FFT函数的命名规则： DSP_ fft 16 x 16 ( )前缀fft或ifft旋转因子精度 16 bit 32 bit变型: _imre: 复数按先虚部/后实部交织存放 r:针对比较小的Cache优化输入输出数据精度 16 bit 32 bit s:结果右移1并四舍五入 Figure 7. FFT函数命名规则 所以函数都按混合基FFT实现。除最后一级外，都是基4运算，最后一级可能是基2或基4。FFT点数要求是2或4的整数次方，且16 ≤ nx ≤ 65536。 所以数据都是Q.15或Q.31的复小数，实部虚部交织存放。 输出数据是经过排序后的数据。FFT蝶形运算可以直接在输入数组x[]上执行，但比特反转的结果必须存放的另一数组，所以最终输出结果要求存放的另一个数组y[]。 关于DSP库的详细信息，请参阅”TMS320C64x+ DSP Library Programmer's Reference (SPRUEB8)”. 下面章节介绍DSP库里面的FFT函数实现的优化细节。 3.2 对旋转因子访问的优化 在第一部分介绍的普通FFT实现里，对旋转因子表的访问不是连续的。为了利用cache的特点，旋转因子表被按如下所示的方式重排，以使它能被连续访问： 22 基于TMS320C64x+DSP的FFT实现 ZHCA414 Table 7. 可连续访问的旋转因子表 第一级的3N/4个旋转因子 第二级的3N/16个旋转因子 第三级的3N/64个旋转因子 …… 0 对每个基4蝶，第一个旋转因子都是WN=1，它不需要被存储，只有其它三个需要存储。所以， 第一级只用3N/4个旋转因子。在下一级，运算被分成四个小FFT，每个小FFT使用相同的旋转因子，所以第二级只用3N/16个旋转因子。以此类推，总共需要的旋转因子的个数为： 3N/4+ 3N/16+ 3N/64+ … = 3N(1/4+ 1/16+ 1/64+ …) <= 3N*(1/4)/(1-(1/4))= N 所以，旋转因子表的大小也是N。 下面的代码产生用于DSP_fft16x16函数的优化的旋转因子表: 基于TMS320C64x+DSP的FFT实现 23 ZHCA414 int gen_twiddle_fft16x16(short *w, int n) { int i, j, k; float M = 32767; float PI = 3.14159265358979323846; /*loops through the stages, loop count is log4(N)-1, Twiddle factor in last stage is always 1 , is not necessary to generate*/ for (j = 1, k = 0; j < n >> 2; j = j << 2) { /*generate twiddle factors for two butterflies in one loop, W[0-1, 4-5, 8-9] for the first butterfly, the others for the second butterfly, loop count= (N/4)/(2j)*/ for (i = 0; i < n >> 2; i += j << 1) { /*twiddle factor for the second output of first butterfly*/ w[k + 1] = M * cos(2.0 * PI * (i ) / n); w[k + 0] = M * sin(2.0 * PI * (i ) / n); /*twiddle factor for the second output of second butterfly*/ w[k + 3] = M * cos(2.0 * PI * (i + j) / n); w[k + 2] = M * sin(2.0 * PI * (i + j) / n); /*twiddle factor for the third output of first butterfly*/ w[k + 5] = -M * cos(4.0 * PI * (i ) / n); w[k + 4] = -M * sin(4.0 * PI * (i ) / n); /*twiddle factor for the third output of second butterfly*/ w[k + 7] = -M * cos(4.0 * PI * (i + j) / n); w[k + 6] = -M * sin(4.0 * PI * (i + j) / n); /*twiddle factor for the fourth output of first butterfly*/ w[k + 9] = M * cos(6.0 * PI * (i ) / n); w[k + 8] = M * sin(6.0 * PI * (i ) / n); /*twiddle factor for the fourth output of second butterfly*/ w[k + 11] = M * cos(6.0 * PI * (i + j) / n); w[k + 10] = M * sin(6.0 * PI * (i + j) / n); k += 12; } } return k; } 请注意，上面的代码里，为了FFT乘法的方便，有些旋转因子的符号被调整。FFT实现在做复数乘法运算时都巧妙的补偿了这些调整的符号。 其它FFT函数需要的旋转因子表可能略有不同，在DSP库里面，每个FFT函数都提供了一个对应的旋转因子生成函数。 通常，旋转因子生成函数在软件初始化阶段调用，生成的旋转因子表被存到数组里。对于不同的FFT点数，要生成不同的旋转因子表。FFT函数则在运行的过程中反复调用，旋转因子表也被反复使用。 24 基于TMS320C64x+DSP的FFT实现 ZHCA414 3.3 蝶形运算的优化 C64x+ DSP library里FFT的实现利用了Table 3中的旋转因子组和FFT级之间的关系，标准C代码里的两个内层循环\和 \被合成一个内层循环。 为了进一步利用C64x+ DSP的8个并行功能单元和两条64bit数据总线，内层循环被展开一倍，每次循环处理两个蝶形，下面的伪代码描述了基本的处理过程： for(loop through the stages) for(i=0; i< N; i+=8) { //every loop calculate two butterflies, 8 points ...... } 我们以DSP_fft16x16的优化实现为例，来说明实现方法。下图描述了每次循环内两个蝶形的运算过程。 X[0] + jX[1] X[2] + jX[3] 1111X[h2]+ jX[h2+1]X[h2+2]+ jX[h2+3]1-j-1jw2[0] + j w2[1] w2[2]+ j w2[3]X[l1] + j X[l1+1] 1-11-1w2[4] + j w2[5]w2[6] + j w2[7]X[l1+2] + j X[l1+3] X[l2] + j X[l2+1] X[l2+2] + j X[l2+3] 1j-1-jw2[8] + j w2[9]w2[10] + j w2[11] 基于TMS320C64x+DSP的FFT实现 25 ZHCA414 Figure 8. 优化实现中的基本蝶形运算 上图中，输入数组定义为： Int16 x[2*N]; 请注意，每个复数的实部和虚部存储在相邻的两个16bit单元，实部在偶数单元，虚部在奇数单元。如上图所示，第一个输入复数为x[0]+jx[1]. 两个蝶形的8个输入点可以用4个64bit (double word) load指令读取： x_3210 = _amemd8(&x[0]); x_h2_32_x_h2_10 = _amemd8(&x[h2]); x_l1_32_x_l1_10 = _amemd8(&x[l1]); x_l2_32_x_l2_10 = _amemd8(&x[l2]); _amemd8 intrinsic 代表double word load指令. _amemd8前面的\说明读取的8字节是按8字节对齐(align)的。 旋转因子的定义如下： Int16 w2[2*N]; 两个蝶形运算用到的6个旋转因子用64bit (double word) load指令读取: co11si11_co10si10 = _amemd8_const(w2); co21si21_co20si20 = _amemd8_const(&w2[4]); co31si31_co30si30 = _amemd8_const(&w2[8]); _amemd8_const intrinsic也执行double word wide load指令, 它的后缀“const”代表被访问的数据不会被修改。 请注意，蝶形运算的第一个输出用到旋转因子始终是1，所以不必存取和计算。 由于两个蝶的运算相同，我们以第一个蝶为例。蝶形内的运算包括加法，减法和乘j。具体运算分解如下： First output= x + x_h2 + x_l1 + x_l2= (x + x_l1) + (x_h2 + x_l2) Third output= x - x_h2 + x_l1 - x_l2= (x + x_l1) - (x_h2 + x_l2) Second output= x - j*x_h2 - x_l1 + j*x_l2= (x - x_l1) - j*(x_h2 - x_l2) Fourth output= x + j*x_h2 - x_l1 - j*x_l2= (x - x_l1) + j*(x_h2 - x_l2) 下面的代码计算 (x + x_l1) 和 (x - x_l1)， (x_h2 + x_l2) 和 (x_h2 - x_l2) 26 基于TMS320C64x+DSP的FFT实现 ZHCA414 xh1_0_xh0_0_xl1_0_xl0_0 = _addsub2((_lo(x_3210)), (_lo(x_l1_32_x_l1_10))); xh1_0_xh0_0 = (_hill(xh1_0_xh0_0_xl1_0_xl0_0)); xl1_0_xl0_0 = (_loll(xh1_0_xh0_0_xl1_0_xl0_0)); xh21_0_xh20_0_xl21_0_xl20_0 = _addsub2((_lo(x_h2_32_x_h2_10)), (_lo(x_l2_32_x_l2_10))); xh21_0_xh20_0 = (_hill(xh21_0_xh20_0_xl21_0_xl20_0)); xl21_0_xl20_0 = (_loll(xh21_0_xh20_0_xl21_0_xl20_0)); 上面的 _lo 和 _loll intrinsic取64 bit数据的低32 bit。每64 bits输入数据中，低32 bit是第一个蝶的数据。类似地，_hi 和 _hill intrinsic取64 bit数据的高32 bit。 那么，蝶形运算第一个输出的计算如下： x_1o_x_0o = _avg2(xh21_0_xh20_0, xh1_0_xh0_0); 上面的_avg2把两个数高低16位分别相加，然后右移一位 (>>1)并四舍五入。 蝶形运算第三个输出的计算如下： myt0_0_mxt0_0 = _sub2(xh21_0_xh20_0, xh1_0_xh0_0); xl1_1_0 = _cmpyr((_lo(co21si21_co20si20)), myt0_0_mxt0_0); _cmpyr把两个16bits Q.15的复数 (蝶形运算的第三个输出和其旋转因子)相乘，得到32bits Q.30的数据，加上0x8000 (四舍五入)，然后取高16位，得到16bit Q.14数据，这实际等效于把Q.15格式的数据右移一位2 (>>1)并做四舍五入。 蝶形运算的第二个和第四个输出的计算因为包含j*(x_h2 - x_l2)而有一点复杂。它可以用复数乘法指令实现，但效率较低。 对于一个16比特复数(a+jb)，实部a存储在低16bit，虚部b存储在高16bit，j*(a+jb)= -b+ja，所以，乘j可以由高低16bit互换得到b+ja，然后与旋转因子相乘。注意，b前面没有负号，而前面我们提到，生成旋转因子时某些旋转因子数值前面的符号也做了特殊的变化，这就正好补偿了b前面的负号。 蝶形运算的第二个和第四个输出的计算如下： xl0_0_xl1_0 = _rotl(xl1_0_xl0_0, 16); //swap the real/imaginary part xt1_0_yt2_0_xt2_0_yt1_0 = _addsub2(xl0_0_xl1_0, xl21_0_xl20_0); xt1_0_yt2_0 = (_hill(xt1_0_yt2_0_xt2_0_yt1_0)); xt2_0_yt1_0 = (_loll(xt1_0_yt2_0_xt2_0_yt1_0)); yt1_0_xt1_0_yt2_0_xt2_0 = _dpackx2(xt1_0_yt2_0, xt2_0_yt1_0); yt1_0_xt1_0 = (_hill(yt1_0_xt1_0_yt2_0_xt2_0)); yt2_0_xt2_0 = (_loll(yt1_0_xt1_0_yt2_0_xt2_0)); xh2_1_0 = _cmpyr((_lo(co11si11_co10si10)), yt1_0_xt1_0); xl2_1_0 = _cmpyr((_lo(co31si31_co30si30)), yt2_0_xt2_0); 请在DSPLib的\\dsplib_v210\\src\\DSP_fft16x16目录查看完整的源代码来研究进一步的细节。其它FFT函数的实现代码也可以从相应的目录找到。 基于TMS320C64x+DSP的FFT实现 27 ZHCA414 3.4 Cache冲突的优化 由于FFT运算过程中有三个数组需同时访问，在基于Cache访问的系统，如C64x+，这可能会导致Cache冲突，从而明显的降低效率。 按上面介绍的实现方法，FFT是按下面的方式逐级计算的： Process N data in stage 1 Process N data in stage 2 …… Process N data in last stage 由于每级运算都会用到大部分数据，如果所有用到的数据合起来大于Cache的容量，则每级处理都会出现Cache冲突。例如，2048点16x16 FFT需要2K*4B=8KB输入数据，8KB旋转因子，和8KB输出数据。如果DSP的L1D Cache大小为32KB，则刚好能满足2048点FFT，处理大于2048点的FFT则会引起L1D Cache的冲突。 我们前面讨论过，对于N点基4 FFT，第一级之后，它被分解成4个N/4点FFT。所以，FFT可按以下方式计算： Perform first stage of N point FFT Perform N/4 point FFT Perform N/4 point FFT Perform N/4 point FFT Perform N/4 point FFT 在第一步，所有的数据都被用到；但在后面，只有N/4的数据被用到。所以，Cache冲突可能只发生在第一步。这样就减轻了Cache冲突的影响。 DSP_fft16x16r函数就是基于以上考虑来实现的。它的定义是： void DSP_fft16x16r(int nx, short * restrict x, short * restrict w, short * restrict y, int radix, int offset, int nmax) 这个函数执行在nx点数据上执行分解后的部分的FFT运算，分解前的主FFT的点数是nmax。如我们前面讨论的，FFT是逐级分解的，本函数一直将FFT分解到点数为radix。如果要函数执行的最后一级，则radix为4或2，nx点个输出会按偏移offset存储到输出数组。例如： DSP_fft16x16r (8192, &x[0], &w[0], y, 4, 0, 8192); 上面的函数调用执行完整的8192点FFT直到radix=4。基于32KB的Cache，每一级都会发生Cache冲突。 同样的功能可以分解成以下函数调用： DSP_fft16x16r(8192, &x[0], &w[0] , y, 2048, 0, 8192); DSP_fft16x16r(2048, &x[0], &w[2*6144], y, 4, 0, 8192); DSP_fft16x16r(2048, &x[2*2048],&w[2*6144], y, 4, 1*2048, 8192); DSP_fft16x16r(2048, &x[4*2048],&w[2*6144], y, 4, 2*2048, 8192); DSP_fft16x16r(2048, &x[6*2048],&w[2*6144], y, 4, 3*2048, 8192); 第一个函数将8192点FFT分解到2048点，也就是，处理基4 FFT的第一级。其它的四个函数执行4个完整的2048点FFT（直到radix=4）并且根据偏移offset将结果存到输出数组。请注意，如前面3.2节讨论的，第二级运算用到的旋转因子从旋转因子表的3*N/4位置开始。 28 基于TMS320C64x+DSP的FFT实现 ZHCA414 基于32KB Cache，Cache冲突只在执行第一个函数时发生，后4个函数执行时不会发生Cache冲突。 按上面的分解方法，点数为N，最后一级为基”rad”的FFT，可以被分解为： DSP_fft16x16r(N, &x[0], &w[0], y, N/4, 0, N); DSP_fft16x16r(N/4,&x[0], &w[2*(3*N/4)], y, rad, 0, N); DSP_fft16x16r(N/4,&x[2*(N/4)], &w[2*(3*N/4)], y, rad, N/4, N); DSP_fft16x16r(N/4,&x[2*(N/2)], &w[2*(3*N/4)], y, rad, N/2, N); DSP_fft16x16r(N/4,&x[2*(3*N/4)], &w[2*(3*N/4)], y, rad, 3*N/4, N); 4 执行速度指标 文档”TMS320C64x+ DSP Library Programmer's Reference (SPRUEB8)” 提供了每个FFT函数理论上的执行速度指标，它没有考虑实际系统中cache miss的开销。本节提供在DSP硬件上的实测数据供参考。 Table 8 和 Figure 9 是在TMS320TCI6484 EVM板上测试的各种FFT执行指令周期数。L1D cache 大小设为 32KB。每个测试之前Cache都被清空。 Table 8. FFT函数执行的指令周期数 256 11501 1595 1601 2071 2362 2882 3307 3204 512 26529 3559 3565 4613 4950 6276 7259 7106 1024 52937 6883 6889 8897 9242 12316 14283 13988 2048 120141 15459 15465 19779 20170 27362 31679 31090 4096 240149 30267 30273 38427 38838 61100 74681 72582 8192 617637 79140 79142 96511 94490 291540 305807 352544 16384 2626845 339383 339385 403297 317907 597642 625685 716664 Point 64 128 5788 837 843 1060 1344 1499 1759 1698 CooleyTukey 2557 fft16x16 DSP_fft16x16 383 DSP_fft16x16_imre 389 full DSP_fft16x16r 478 decomposed 722 DSP_fft16x16r DSP_fft16x32 775 DSP_fft32x32 831 DSP_fft32x32s 763 基于TMS320C64x+DSP的FFT实现 29 ZHCA414 FFT cycles benchmark400003500030000CooleyTukeyFft16x16 DSP_fft16x16 DSP_fft16x16_imre full DSP_fft16x16r decomposed DSP_fft16x16rDSP_fft16x32 DSP_fft32x32 DSP_fft32x32s cycles2500020000150001000050000025651276810241280153617922048pointsFFT cycles benchmark800000700000600000CooleyTukeyFft16x16 DSP_fft16x16 DSP_fft16x16_imre full DSP_fft16x16r decomposed DSP_fft16x16rDSP_fft16x32 DSP_fft32x32 DSP_fft32x32s cycles50000040000030000020000010000000204840966144points819210240122881433616384 Figure 9. FFT函数执行的指令周期数 在上面的表和图里，“full DSP_fft16x16r”是一次执行完FFT运算；“decomposed DSP_fft16x16r” 是按前面3.4节介绍的Cache优化方法，将N点FFT分解为4个N/4点FFT。 从上面的数据可以看出，标准C实现的CooleyTukey FFT比DSP Library里优化的FFT实现消耗的指令周期多很多。对于少于4096点的FFT，消耗的指令周期和点数几乎成正比；而当点数大于4096点时，指令周期数增长加快，这是由于Cache冲突造成的。可以看出“decomposed DSP_fft16x16r”在点数较大时消耗的指令周期少于” full DSP_fft16x16r”。 30 基于TMS320C64x+DSP的FFT实现 ZHCA414 Table 9 和 Figure 10是在TMS320TCI6484 EVM板上测试的各种IFFT执行指令周期数。 Table 9. Point DSP_ifft16x16 DSP_ifft16x16_imre DSP_ifft16x32 DSP_ifft32x32 64 402 439 683 783 128 856 900 1461 1740 IFFT函数执行的指令周期数 256 512 3575 3615 6424 7238 1024 6899 6941 12680 14262 2048 15475 15521 28412 31656 4096 30283 30331 60450 74648 8192 79151 79207 323618 304706 16384 339397 339449 668534 623634 1611 1647 2872 3288 IFFT cycles benchmark400003500030000DSP_ifft16x16 DSP_ifft16x16_imreDSP_ifft16x32 DSP_ifft32x32 cycles25000200001500010000500000256512768points10241280153617922048 IFFT cycles benchmark800000700000600000DSP_ifft16x16 DSP_ifft16x16_imreDSP_ifft16x32 DSP_ifft32x32 cycles50000040000030000020000010000000204840966144points819210240122881433616384 基于TMS320C64x+DSP的FFT实现 31 ZHCA414 Figure 10. IFFT函数执行的指令周期数 用TI的指令周期消耗分析工具，我们发现上面的数据中大概有1/3浪费在L1D cache miss上。主要由访问输入数组，旋转因子和输出数组引起，它发生在这些数据从L2搬到L1D时。 由于相同的旋转因子在多次调用FFT函数时被反复使用，所以，一个好的解决旋转因子L1D miss的方法是，把L1D的一部分当作普通RAM使用，直接把旋转因子表放到L1D RAM 中。 为了解决输入输出数据引起的L1D miss，我们也可以在L1D中为它们划分专门的空间。由于输入输出数据是不断变化的，可以在DSP执行当前的FFT时，同时用DMA将下一组数据搬到L1D中。但需要注意的是，在L1D中分配固定的数据块会减少L1D Cache的大小，这可能会降低其它需要较大L1D Cache的算法或函数的性能。所以，这需要在系统设计上做折中。 另外一个方法可减少Cache miss的影响，它把需要用到的数据在访问之前”touch”到Cache中，关于“touch”方法的更多信息，请参阅”TMS320C64x+ DSP Cache User's Guide (SPRU862)”。以上所以测试都没有用”touch”方法。 Take DSP_fft16x16 as example, Table 10 shows the impact of touching different data before perform DSP_fft16x16. Table 10. Cache “touch”对DSP_fft16x16性能的影响 64 128 256 512 1024 2048 4096 8192 16384 No touch 413 Touch twiddle factor 412 Touch twiddle factor and input data Touch twiddle factor and input, output data 857 834 793 1591 1500 1389 3555 3334 3087 6875 6392 5873 15447 14442 13379 11997 22.33% 30279 28258 26115 31795 13.75% 79160 79132 81416 339391 344064 348701 411 431 768 1272 2785 5204 Best improvement 0.48% 10.39% 20.05% 21.66% 24.31% 86208 353363 0.04% N/A 上表中，第二行数据是下面两个函数总共消耗的指令周期数： touch(w16, nx*4); DSP_fft16x16(w16, nx, x16, y16); 上面的代码在调用DSP_fft16x16之前先把旋转因子”touch”到Cache中。 上表中，第三行数据是下面三个函数总共消耗的指令周期数： touch(w16, nx*4); touch(x16, nx*4); DSP_fft16x16(w16, nx, x16, y16); 上面的代码在调用DSP_fft16x16之前先把旋转因子和输入数据”touch”到Cache中。 上表中，第四行数据是下面四个函数总共消耗的指令周期数： 32 基于TMS320C64x+DSP的FFT实现 ZHCA414 touch(w16, nx*4); touch(x16, nx*4); touch(y16, nx*4); DSP_fft16x16(w16, nx, x16, y16); 上面的代码在调用DSP_fft16x16之前先把旋转因子，输入数据和输出数据”touch”到Cache中。 最后一行表示最好情况和没有”touch”的情况相比性能改进的百分比。最后的性能改进可达到24%。 通常，在FFT函数调用之前”touch”数据可减少FFT函数的执行周期数。但调用”touch”函数会消耗额外的指令周期。所以，只有当FFT函数上节省的周期数大于”touch”函数消耗的周期数才能带来整体性能的改进。 通常，被”touch”的数组越大，节省的周期数会越多。但总共被”touch”的数据不能超过cache的大小。这就是上表中点数大于4096的情况，Cache冲突会在”touch”时发生，性能不但不能改善，反而有可能恶化，所以”touch”不应该在这种情况下使用。 在实际的应用中，数据结构可能比这里的测试更复杂，设计者应该尝试不同的数据”touch”的组合来得到最好的结果。另外，在实际应用中，FFT函数的输入可能是它之前的另一个函数的输出，这时，FFT的输入数据可能已经在Cache之中，那么就没有必要再去”touch”这些数据。如果FFT的输出会被另一个函数使用，那么执行FFT之前”touch” FFT的输出数据数组可能更有意义。 5 数据缩放和精度 FFT计算中的乘法和加法可能会导致溢出。为了解决乘法溢出，所有的数据都被表示为小数，并按小数相乘，因为两个小数相乘的结果总是小于或等于1。 第二种溢出来自于加法。根据前面讲到的基4 FFT算法，每个蝶形输出由四个数相加得到，它最多带来2bit的进位。对于基2运算(混合基最后一级)，每级可能产生1bit进位。解决这个问题最直接的方法是每级的输入预先右移2bit或1bit。 每级都对输入数据右移2bit或1bit会消耗额外的指令周期。一个简化的方法是一次性对输入数据右移足够的bit，以保证后面所有运算都不会溢出。N点FFT整个运算可能带来的进位总数为log2(N)，所以一次性把输入数据右移log2(N) bit就可防止FFT运算溢出。 然而，右移防止溢出是以损失精度为代价的。例如，1024点FFT，为防止溢出，输入数据必须被右移10bit。如果输入数据是16bit的，那么，只有6bit有效精度。 在DSP library里，大部分FFT函数在每级(最后一级除外)FFT运算后都对输出右移一位并四舍五入。这由使用带右移功能的加法和乘法指令实现，如AVG2 和 CMPYR, 这并不增加额外的指令周期。对于这种实现，要防止FFT运算溢出，只需要对输入数据右移: log2(N)- (#stages-1)= log2(N)- ceil[log4(nx)?1] FFT点数越多，移位所导致的精度损失越大。如果16bit FFT不能满足精度要求，可以用32bit的FFT。通常，32bit FFT消耗的指令周期是16bit FFT的两倍，占用的存储器也是两倍。 基于TMS320C64x+DSP的FFT实现 33 ZHCA414 另外一个解决的办法是检查每级输出的进位，然后确定右移几位来防止溢出。这会使精度损失最小，但它增加了实现的复杂度。TI DSP library不直接支持这种实现。由于TI DSP Library提供了所有函数的源代码，客户可以在这些源代码的基础上做一些修改来实现特殊功能。让我们以DSP_fft16x16为基础，实现一个叫DSP_fft16x16_scale的带动态缩放的FFT函数。 我们假设输入数据是16-bit Q.15格式的小数，表示为： S.ABC DEFG HIJK LMNO, 上面，每个字母代表一个bit，最高位，S，是符号位。A…O是数据位。为了防止基4 FFT加法溢出，我们需要两个保护bit，也就是，输入数据的bit 14 (A) 和 bit 13 (B)应该等于符号位(S)。如果它们不等于符号位，我们就需要右移来得到两个保护bit。 为了确定bit 13 和 14 是否和符号位相同，我们把输入数据先右移两位，得到如下新数据： S.SSA BCDE FGHI JKLM 把它和原始数据异或(XOR)，如果原始bit 14 或 13 与符号位不同则结果的bit 14 或 13就是1。 对每级的所有输出数据进行这种检测，代码如下： next_scale|= ((x0)^_shr2(x0,2))|; ((x1)^_shr2(x1,2))| ((x2)^_shr2(x2,2))| ((x3)^_shr2(x3,2)); 上面的 x0, x1, x2, x3是蝶形运算的四个输出。请注意，它们是32位的，低16位是复数的实部，高16位是复数的虚部。所以，我们用intrinsic _shr2()对高低16位同时右移。在每个数据上的检测结果被或(OR)在一起，因此，如果任何一个数据的bit14或13和符号位不同，则最终结果的bit14或13就为1。在每一级运算的最后，我们可以确定要右移多少位，代码如下： next_scale= next_scale|(next_scale>>16); //OR higher 16 bit with lower 16 bit if(next_scale&0x4000) //if bit 14 is 1, we should scale 2 bit scale= 2; else if(next_scale&0x2000) //if bit 13 is 1, we should scale 1 bit scale= 1; else scale= 0; 在每个数据被下一级使用前，我们对它右移“scale”，这样就可以保证下一级运算不溢出，代码如下： x_3210 = _amemd8(&x[0]); x_3210_l= _shr2(_lo(x_3210),scale); x_3210_h= _shr2(_hi(x_3210),scale); 这个实现使得精度损失最小化，但它增加了指令周期数。下表比较了DSP_fft16x16_scale 和 DSP_fft16x16的指令周期数。这里两种情况都使用Cache “touch”。 Table 11. 数据动态缩放的DSP_fft16x16的性能 64 128 256 512 1024 2048 4096 8192 DSP_fft16x16 407 744 1250 2763 5182 11964 31753 86186 34 基于TMS320C64x+DSP的FFT实现 ZHCA414 DSP_fft16x16_scale 813 1266 2662 4840 11522 21886 56056 131958 从上面的数据可看出，DSP_fft16x16_scale消耗的周期数大概是DSP_fft16x16的两倍。 References 1. TMS320C64x+ DSP Library Programmer's Reference (SPRUEB8) 2. TMS320C64x C64x+ DSP CPU and Instruction Set Reference Guide (SPRU732) 3. TMS320C64x+ DSP Cache User's Guide (SPRU862) 4. Extended-Precision Complex Radix-2 FFT_IFFT Implemented on TMS320C62x (SPRA696) 5. Auto-scaling Radix-4 FFT for TMS320C6000 DSP (SPRA654) 6. Understanding Digital Signal Processing, Richard G. Lyons, 2004 基于TMS320C64x+DSP的FFT实现 35 重要声明 德州仪器(TI)及其下属子公司有权在不事先通知的情况下,随时对所提供的产品和服务进行更正、修改、增强、改进或其它更改，并有权随时中止提供任何产品和服务。客户在下订单前应获取最新的相关信息,并验证这些信息是否完整且是最新的。所有产品的销售都遵循在订单确认时所提供的TI销售条款与条件。 TI保证其所销售的硬件产品的性能符合TI标准保修的适用规范。仅在TI保证的范围内,且TI认为有必要时才会使用测试或其它质量控制技术。除非政府做出了硬性规定,否则没有必要对每种产品的所有参数进行测试。 TI对应用帮助或客户产品设计不承担任何义务。客户应对其使用TI组件的产品和应用自行负责。为尽量减小与客户产品和应用相关的风险，客户应提供充分的设计与操作安全措施。 TI不对任何TI专利权、版权、屏蔽作品权或其它与使用了TI产品或服务的组合设备、机器、流程相关的TI知识产权中授予的直接或隐含权限作出任何保证或解释。TI所发布的与第三方产品或服务有关的信息，不能构成从TI获得使用这些产品或服务的许可、授权、或认可。使用此类信息可能需要获得第三方的专利权或其它知识产权方面的许可，或是TI的专利权或其它知识产权方面的许可。对于TI的产品手册或数据表，仅在没有对内容进行任何篡改且带有相关授权、条件、限制和声明的情况下才允许进行复制。在复制信息的过程中对内容的篡改属于非法的、欺诈性商业行为。TI对此类篡改过的文件不承担任何责任。 在转售TI产品或服务时，如果存在对产品或服务参数的虚假陈述，则会失去相关TI产品或服务的明示或暗示授权，且这是非法的、欺诈性商业行为。TI对此类虚假陈述不承担任何责任。 TI产品未获得用于关键的安全应用中的授权，例如生命支持应用（在该类应用中一旦TI产品故障将预计造成重大的人员伤亡），除非各方官员已经达成了专门管控此类使用的协议。购买者的购买行为即表示，他们具备有关其应用安全以及规章衍生所需的所有专业 技术和知识，并且认可和同意，尽管任何应用相关信息或支持仍可能由TI提供，但他们将独力负责满足在关键安全应用中使用其产品及TI产品所需的所有法律、法规和安全相关要求。此外，购买者必须全额赔偿因在此类关键安全应用中使用TI产品而对TI及其代表造成的损失。TI产品并非设计或专门用于军事/航空应用，以及环境方面的产品，除非TI特别注明该产品属于“军用”或“增强型塑料”产品。只有TI指定的军用产品才满足军用规格。购买者认可并同意，对TI未指定军用的产品进行军事方面的应用，风险由购买者单独承担，并且独力负责在此类相关使用中满足所有法律和法规要求。 TI产品并非设计或专门用于汽车应用以及环境方面的产品，除非TI特别注明该产品符合ISO/TS16949要求。购买者认可并同意，如果他们在汽车应用中使用任何未被指定的产品，TI对未能满足应用所需要求不承担任何责任。可访问以下URL地址以获取有关其它TI产品和应用解决方案的信息： 产品 数字音频放大器和线性器件数据转换器DLP?产品 DSP-数字信号处理器时钟和计时器接口逻辑电源管理微控制器(MCU)RFID系统 OMAP机动性处理器无线连通性 www.ti.com.cn/audiowww.ti.com.cn/amplifierswww.ti.com.cn/dataconverterswww.dlp.comwww.ti.com.cn/dsp www.ti.com.cn/clockandtimerswww.ti.com.cn/interfacewww.ti.com.cn/logicwww.ti.com.cn/power www.ti.com.cn/microcontrollerswww.ti.com.cn/rfidsyswww.ti.com/omap www.ti.com.cn/wirelessconnectivity德州仪器在线技术支持社区 www.deyisupport.com IMPORTANTNOTICE通信与电信计算机及周边消费电子能源工业应用医疗电子安防应用汽车电子视频和影像 应用 www.ti.com.cn/telecomwww.ti.com.cn/computerwww.ti.com/consumer-appswww.ti.com/energywww.ti.com.cn/industrialwww.ti.com.cn/medicalwww.ti.com.cn/securitywww.ti.com.cn/automotivewww.ti.com.cn/video 邮寄地址：上海市浦东新区世纪大道1568号，中建大厦32楼邮政编码：200122 Copyright?2012德州仪器半导体技术（上海）有限公司