1 前言DSP结构可以分为定点和浮点型两种。

其中,定点型DSP可以实现整数、小数和特定的指数运算,它具有运算速度快、占用资源少、成本低等特点;灵活地使用定点型DSP进行浮点运算能够提高运算的效率。

目前对定点DSP结构支持下的浮点需求也在不断增长,主要原因是:实现算法的代码往往是采用C/C++编写,如果其中有标准型的浮点数据处理,又必须采用定点DSP器件,那么就需要将浮点算法转换成定点格式进行运算。

同时,定点DSP结构下的浮点运算有很强的可行性,因为C语言和汇编语言分别具有可移植性强和运算效率高的特点,因此在定点DSP中结合C语言和汇编语言的混合编程技术将大大提高编程的灵活度,以及运算速度。

大多数DSP的开发工具只是在C语言的基础上支持标准的浮点运算,而定点DSP硬件一般都是面向定点的运算,不支持标准的浮点运算,缺乏硬件的支持极大地限制了浮点的应用,因而标准的浮点运算在实际定点DSP应用中并不多见。

C5509是一款16位定点DSP。

在本文中,对C5509输入FTSK信号,用C语言和汇编语言混合编程的方式对输入浮点型的FTSK信号进行相关运算,并输出浮点运算结果。

这里叶变换(FFT是一种高效实现离散傅里叶变换(DFT的快速算法,是数字信号处理中最为重要的工具之一,它在声学,语音,电信和信号处理等领域有着广泛的应用。

2 方案设计2.1 方案的提出DSP (数字信号处理器与一般的微处理器相比有很大的区别,它所特有的系统结构、指令集合、数据流程方式为解决复杂的数字信号处理问题提供了便利,本文选用TMS320C54X 作为DSP 处理芯片,通过对其编程来实现FFT 的浮点DSP实现。

DSP 应用系统设计的一般流程如下图所示:图2.1 DSP 系统流程图 2.2方案的论证旋转因子W N 有如下的特性。

对称性:/2k k N N N W W +=-周期性:kk NN NW W +=利用这些特性,既可以使DFT 中有些项合并,减少了乘法积项,又可以将长序列的DFT 分解成几个短序列的DFT 。

FFT 就是利用了旋转因子的对称性和周期性来减少运算量的。

FFT 的算法是将长序列的DFT 分解成短序列的DFT 。

例如:N 为偶数时,先将N 点的DFT 分解为两个N/2点的DFT ,使复数乘法减少一半:再将每个N/2点的DFT 分解成N/4点的DFT ,使复数乘又减少一半,继续进行分解可以大大减少计算量。

最小变换的点数称为基数,对于基数为2的FFT 算法,它的最小变换是2点DFT 。

一般而言,FFT 算法分为按时间抽取的FFT (DIT FFT和按频率抽取的FFT(DIF FFT 两大类。

D IF FFT 算法是在时域内将每一级输入序列依次按奇/偶分成2个短序列进行计算。

而DIF FFT 算法是在频域内将每一级输入序列依次奇/偶分成2个短序列进行计算。

两者的区别是旋转因子出现的位置不同,得算法是一样的。

在DIF FFT 算法中,旋转因子k N W 出现在输入端,而在DIF FFT算法中它出现在输入端。

假定序列x(n的点数N 是2的幂,按照DIF FFT 算法可将其分为偶序列和奇序列。

偶序列:1(0,(2,(4,(N-2,(2,0,1,/21x x x x x x r r N ==- 即奇序列:2(1,(3,(5,(1,(21,0,1,/21x x x x N x x r r N -=+=- 即则x(n的DFT 表示为11/21/212(2100/21/212212((((2(21(((2N N nknkNN n n N N rk r k N Nr r N N rk k rk N N Nr r X k x n Wx n W n n x r Wx r Wx r WWx r W--==--+==--===+=++=+∑∑∑∑∑∑为偶数为奇数由于22(2//2j N NN W eW π-⎡⎤===⎣⎦,则(3式可表示为/21/211/22/212(((((0,1,/21(3N N rkk rk N NN r r kN X k x r WWx r WX k W X k k N --===+=+=-∑∑式中, (1k X 和(2k X 分别为(1n x 和(2n x 的N/2的DFT 。

由于对称性,,2/KN N k N W W -=+则12(/2((k N X k N X k W X k +=-。

因此,N 点(k X 可分为两部分:前半部分:12/,1,0(((21-=+=N k k X W k X k X kN (4后半部分:12/,1,0((2/(21-=-=+N k k X W k X N k X k N (5从式(4和式(5可以看出,只要求出0~N/2-1区间(1k X 和(2k X 的值,就可求出0~N-1区间(k X 的N 点值。

以同样的方式进行抽取,可以求得N/4点的DFT ,重复抽取过程,就可以使N 点的DFT 用上组2点的 DFT 来计算,这样就可以大减少运算量。

基2 DIF FFT 的蝶形运算如图(a 所示。

设蝶形输入为(1p x m -和(1q x m -,输出为(p x m 和(q x m ,则有k N m m m W q x p x p x (((11--+= (6k N m m m W q x p x q x (((11---= (7在基数为2的FFT 中,设N=2M,共有M 级运算,每级有N/2个2点FFT 蝶形运算,因此,N 点FFT 总共有N N 2log 2/(个蝶形运算。

(1q x m - (p x m(1q x m - (q x m-1图(a 基2 DIF FFT 的蝶形运算例如:基数为2的FFT ,当N=8时,共需要3级,12个基2 DIT FFT 的蝶形运算。

其信号流程如图(b所示。

x(0x(0W N0x(4x(1-1W N0x(2x(2-1W N0W N2x(6x(3-1 -1W N0x(1x(4-1W N0W N1x(5 x(5 -1 -1W N0W N2x(3x(6-1 -1W N0W N2W N3x(7 x(7 -1 -1 -1图(b 8点基2DIF FFT蝶形运算从图(b可以看出,输入是经过比特反转的倒位序列,称为位码倒置,其排列顺序为7(4(xx,0(xxx。

输出是按自然顺序排列,其顺序为x xx,5(3(,,,,1(2(,6(0(xxx 。

,x,7(,1(,6(2.3 方案选择此相关运算的输人是浮点型数据,相关系数是小于1的单精度浮点型数。

对于定点DSP,由于不能直接进行浮点数的乘法运算,因此必须对输入数据进行类型转换。

首先,相关运算的输入数据是FTSK浮点数据。

在C语言中,单精度浮点数据是以IEEE754标准存储的32位数据,而C5509中C语言调用汇编语言,是通过寄存器AR0从C语言传递给汇编语言的是数据指针,这个指针是指向16位数据的,所以相关的输入32位浮点数要先转化为16位整型数据。

本文这样实现:C程序中先把浮点数据乘以10后(提高运算精度,强制类型转化为整型数据,然后把此16位数据的指针赋给调用汇编的入口参数,即通过寄存器AR0传递到汇编程序中。

然后,在汇编程序中,相关的系数是小于l 的小数;在DSP 中,汇编语言直接定义的格式是将其转换为16位二进制2的补码表示形式(例如0.8用8×32 768/lO 来表示。

从汇编程序入口进入的、经过强制类型转换的整型数据也是以16位二进制形式存储的,通过与16位的小数相乘得到的是32位数,存储在累加器A 中。

其中,前16位是运算结果的整数部分,后16位是小数部分。

由于从汇编语言程序返回C 程序的参数是16位的,故取运算结果的高16位(此前已经把输入数据乘以lO ,最大限度地提高了运算精度,这里直接取高16位。

把这16位数据返回C 程序,得到整型数据,再强制类型转化为单精度浮点型数据,再除以10,即得到了最后相关运算的结果。

经实际运算检验,通过这种方法在C5509里进行浮点运算,最终结果实现了很高的精度,而且通过调用汇编语言,极大地提高了运算的效率。

在IEEE754单精度浮点标准中,明确包含了符号位,第32位用作符号位。

尾数进行了归一化,以产生一个1.f 格式的数,f 是小数部分,占用分配的23位。

因为规格化的数最左一位总是所以不需要存储该位,在该格式中它是隐式的。

这样一个n 位的尾数实际上存放了一个n+l 位数。

为使尾数规格化,指数被适当增减,来跟踪规格化所需的左右移位数以及小数点。

最常用的是用8位指数表示0~255,即127(12(1.se f --⨯其中:s 是符号位,0为正数,1为负数;e 是指数位,无符号8位;f 是尾数的小数部分,23位。

例如:IEEE754格式下浮点正数00110001001111l000000001000000000的十进制表示为:符号位=0(因为它是一个正数尾数=23456151222222------++++++=1.48440551758指数=65222++十进制等值数=(981271.484405517582-⨯=92.7649207368110-⨯IEEE754格式下浮点负数1110000101010110000000000000000的十进制表示为: 符号位=1(因为它是一个负数尾数=135612222----++++指数=761222++十进制等值数=(1941271.68752--⨯=202.4903104499510-⨯3 硬件设计3.1 CS 开发环境CCS 提供了配置、建立、调试、跟踪和分析程序的工具,它便于实时、嵌入式信号处理程序的编制和测试,它能够加速开发进程,提高工作效率。

CCS 提供了基本的代码生成工具,它们具有一系列的调试、分析能力。

CCS 支持如下图1.1所示的开发周期的所有阶段。

图 3.1 CCS 开发周期阶段图3.2 SEED-DEC2812开发实验箱SEED-DECxxxx 系列嵌入式DSP 开发板本着模块化、总线型、开放式、系列化的设计思想,采用统一的系统结构、模块结构和机械结构,以多种典型DSP 处理器构成具有标准总线和相同物理尺寸的高性能嵌入式DSP 开发板。

SEED-DEC2812 嵌入式DSP 开发板原理框图如图1.2所示:图 3.2 DSP 开发板原理框图4程序代码DATE:09/15/2006AUTHOR:CHEN PENGVERSION:ORIGINAL 09/15/2006void fft_r4_real_512( short x[], short w[], short rs[], short points { int n1, n2, n3, ie;int iw1, iw2, iw3;int rw1, rw2, rw3;int rx0, rx1, rx2, rx3;int ix0, ix1, ix2, ix3;int iStage, iButterfly, counter;int xe_real, xe_imag, xo_real, xo_imag;long real0, real1, real2;long imag0, imag1, imag2;/* Exception handling *//* points must be 128, 512 or 2048 */if( points != 512 || points != 2048 || points != 128 {printf("#### Error: Points of input data must be 128, 512 or 2048\n";return;}n2 = (points >> 1; // n2 is used for butterfly countsie = 1; // ie is used for the steps for the index of twist factorscounter = 0;/* Complex FFT calculation *//** X(4r = FFT[A(n] A(n = x(n + x(n + N/2 + x(n + N/4 + x(n + 3N/4* X(4r + 2 = FFT[B(n] B(n = [x(n + x(n + N/2 - x(n + N/4 - x(n + 3N/4]*Wn(2n* X(4r + 1 = FFT[C(n] C(n = [x(n - x(n + N/2 - jx(n + N/4 + jx(n + 3N/4]*Wn(n * X(4r + 3 = FFT[D(n] D(n = [x(n - x(n + N/2 + jx(n + N/4 - jx(n + 3N/4]*Wn(3n */for( iStage = (points >> 1; iStage > 1; iStage >>= 2 {n1 = n2;n2 >>= 2;n3 = n2 << 1;rw1 = 0;iw1 = 1;for( iButterfly = 0; iButterfly < n2; iButterfly++ {rw2 = (rw1 + rw1 << 1;rw3 = (rw2 + rw1 << 1;iw2 = rw2 + 1;iw3 = rw3 + 1;for( rx0 = ( iButterfly<<1 ; ix0 < points; ix0 += n1 {rx1 = rx0 + n3; //real parts indexrx2 = rx2 + n3;rx3 = rx2 + n3;ix0 = rx0 + 1; //image parts indexix1 = rx1 + 1;ix2 = rx2 + 1;ix3 = rx3 + 1;/*************************\*even points calculations\*************************/real0 = x[rx1] + x[rx3];imag0 = x[ix1] + x[ix3];real1 = x[rx0] + x[rx2];imag1 = x[ix0] + x[ix2];real2 = x[rx0] - x[rx2];imag2 = x[ix0] - x[ix2];x[rx0] = (real0 + real1 >> rs[counter]; //A(n real partx[ix0] = (imag0 + imag1 >> rs[counter]; //A(n image partreal1 = real1 - real0;imag1 = imag1 - imag0;real0 = x[rx1] - x[rx3];imag0 = x[ix1] - x[ix3];x[rx1] = ((real1 * w[rw2] - (imag1 * w[iw2] >> rs[counter];//B(n real partx[ix1] = ((imag1 * w[rw2] + (real1 * w[iw2] >> rs[counter];//B(n image part /*************************\*odd points calculations\*************************/real1 = real2 + imag0;imag1 = imag2 + (-real0;real2 = real2 - imag0;imag2 = imag2 - (-real0;x[rx2] = ((real1 * w[rw1] - (imag1 * w[iw1] >> rs[counter];//C(n real partx[ix2] = ((imag1 * w[rw1] + (real1 * w[iw1] >> rs[counter];//C(n image part x[rx3] = ((real2 * w[rw3] - (imag2 * w[iw3] >> rs[counter];//D(n real partx[ix3] = ((imag2 * w[rw3] + (real2 * w[iw3] >> rs[counter];//D(n image part }//end of inner for(rw1 += ie;rw1 <<= 1;iw1 = rw1 + 1;}//end of outer for(counter++;ie <<= 2;}//end of outest for(/** Bit-reverse*///Add the code here/* Output separation *//** [X^(k + X^*(-k] / 2 = Xe(k* [X^(k - x^*(-k] / 2j = Xo(k* X(k = Xe(k + Xo(k * e(-j*2*pi*k/N * N = points / 2* X^*(-k = X^*(N - k*/rx1 = 0; //index for kix1 = 1;rx2 = points - 2; //index for N-kix2 = points - 1;rw1 = 0; //index for twist factor iw1 = 1; for( ; rx1 < points; {xe_real = ( x[rx1] + x[rx2] >> 1; //Xe(k real part xe_imag = ( x[ix1] + (-x[ix2] >> 1; //Xe(k image part xo_real = ( x[ix1] - (-x[ix2] >> 1; //Xo(k real partxo_imag = (-( x[rx1] - x[rx2] >> 1; //Xo(k image partx[rx1] = xe_real + (xo_real * w[rw1] - xo_imag * w[iw1]; //X(k real partx[ix1] = xe_imag + (xo_real * w[iw1] + xo_imag * w[rw1];//X(k image partrx1 += 2;rx2 -= 2;ix1 = rx1 + 1;ix2 = rx2 + 1;rw1++;iw1 = rw1 + 1;}return;}5系统仿真结果图5.1输出波形图(一图5.2输出波形图(二6设计总结本实验通过学习快速傅里叶变换(FFT的原理,然后在CCS平台下编程对其进行模拟仿真,对快速傅里叶变换(FFT有个一个较深刻的理解。