TI库实现2048点或更大点数FFT

1、在DSP运算中，经常需要把输入时域信号在频域进行处理之后，再还原为时域信号，这样就需要进行FFT和IFFT运算：x(n) -> FFT -> X(f) -> 频域处理 -> Y(f) -> IFFT -> y(n)而一般的DSP芯片只支持整数运算，也就是说只能进行定点小数计算。N点FFT计算出0 N-1，N个复数： 0，A，N/2，A*，A为(N/2-1)个复数，A*为A的共轭复数。FFT的公式为：            

2、;           N     X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N)、1 < =  k < =  N.                

3、;      n = 1IFFT的公式为：                          N     x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N)、1 < =  n <

4、=  N.                         k = 1假设我们对ADC转换器转换的数字信号进行FFT运算，若输入数据为16bit的短整型数，我们可以把它看作Q15的从-1到1之间的小数。根据FFT的公式我们可以知道，FFT变换之后的结果将超出这个范围。例如在matlab中输入fft(sin(1:8*0.5)，可以看到结果：2.859

5、7，-0.8019 - 3.0216i，0.4312 - 0.8301i，0.5638 - 0.3251i，0.5895，0.5638 + 0.3251i，0.4312 + 0.8301i，-0.8019 + 3.0216i实际上，FFT变换之后的数据的范围在-N到N之间，N为FFT的点数。为了正确地表达-N到N之间的数值，输出数据的Q值将变小，例如若N=1024，输入数据为Q15的话，那么输出数据则必须为Q5才能够确保结果不会溢出。这样的结果将丢失很多信息，以至于IFFT无法还原为原来的数据。如下图所示：Q15 -> 1024FFT -> Q5 -> 1024IFFT -&

6、gt; Q5这样，经过FFT和IFFT变换之后，数据从Q15变成了Q5，丢失了10bit的信息。为了使得定点FFT和IFFT之后能够还原为原来的数据，必须使用32bit定点FFT，输入数据虽然是16bit，在进行FFT的时候将将它转换为32bit运算，这样输入的数据为Q31，1024点FFT的输出数据为Q21，满足16bit的精度要求。这样整个计算流程变为：Q15 -> Q31 -> 1024FFT32bit -> Q21 -> 1024IFFT32bit -> Q21 -> Q15在TI的DSPLIB中的FFT和IFFT都提供了两种选择：SCALE和NOS

7、CALE。仍以1024点计算为例，它们的意思分别为：Q15 -> NOSCALE FFT -> Q15 (结果可能溢出)Q15 -> SCALE FFT    -> Q5 (结果不会溢出，但是精度降低)Q15 -> NOSCALE IFFT -> Q15 (结果可能溢出)Q15 -> SCALE IFFT    -> Q5 (结果不会溢出，但是精度降低)32位定点FFT，IFFT与此类似，不再重复。因此，用DSPLIB进行FFT和IFFT计算时，注意必须采用

8、32bit精度，而且FFT变换时采用SCALE，而IFFT变换时采用NOSCALE。关于TIFFT库的使用  2010-01-25 12:53 星期一这篇文章是应一些找我讨论DSP的同学所写，贴在这里大家一起学习。曾有不少论坛上的同学（包括DSP算法讨论群里的一些同学）问过我关于TI的FFT库的使用，这里我将我使用过的一些经验说一下。TI的这个FFT库在计算速度、计算精度以及数据存储等方面是做了不少优化的，比如数据存储，若作N点的FFT，供查表用的旋转因子必须有N/2点的正弦值与N/2点的余弦值，这个库将其压缩成3N/4点的正弦值，因此就节省了N/4点的存储空间；另外计算N

9、点实数FFT时，一般简单的做法是将N点实数的虚部全化为0来处理，而这个库则把N 点实数数据打包成N/2点复数数据来处理，在计算速度和存储空间都有很大改进。之前我浏览helloDSP论坛的帖子，有很多人发出疑问：为什么我计算出的mag值全为零？这样的帖子真不少见；TI的官方工程师论坛不少老外也发问：Why did I get all zeros ?事实上我想主要原因是输入数据格式不对。我认为使用这个库主要注意一下两点：1. 数据输入输出的Q格式；2. 存储空间分配。下面以32位实数FFT为例来说明。注意到文档的40页有如下说明：1. 在函数void calc(RFFT32_handle)有如下一

10、句：Note that the input and output data are in Q31 format.；2. 在函数void mag(FFT128R_handle)有如下一句：Note that the magnitude output is stored in Q30 format.因为28x系列DSP是定点处理器，而FFT计算涉及到不少浮点计算，TI使用Q格式来解决这个问题（Q格式说明可参考sprc087_IQmath）。事实上输入数据采用Q31格式能在避免计算溢出前提下获得最好的计算精度。对AD采样的数据进行FFT计算，定义计算缓冲区数组：long ipcbN+2;因为AD结果

11、寄存器是12位的，在数据左对齐的情况下直接左移15位即可：ipcbConversionCount = (unsigned long)AdcRegs.ADCRESULT0)<<15;另外我记得有份文档提到在某些存储器下是右对齐，此时则需左移19位，大家在使用时注意这个问题。当采样完成后，按照FFT库的文档上的说明或者仿照文档附带的例程，进行相应函数调用来实现自己的FFT计算，比如按计算点数来定义各个变量数组，是否加窗，是否求解幅值平方值等。在进行存储器分配时，文档上要求（128点实数FFT为例）：FFTipcb ALIGN(256) : > L0L1RAM PAGE 1FFTm

12、ag > L0L1RAM PAGE 1FFTtf > NVMEM PAGE 0 /* Non volatile memory */.econst >NVMEM PAGE 0 /* Non volatile memory */注意两点：FFT计算缓冲区FFTipcb需在page1上连续分配2N个位置（以ALIGN来指定），FFTtf（旋转因子）需放在Non volatile memory的page0内（事实上如何才能为Non volatile我也不清楚）。FFTtf位置这点我之前在调试时对计算FFT影响很大，因为twiddle factor若因存储冲突肯定会造成查表值不准确，那

13、计算FFT时肯定就不对了。有个论坛帖子作者说一定要放在origin = 0x008000开始位置，其实也不对，大家可自己去试验；后面我也会给出我的存储配置文件（即.cmd文件内容）。无图无真相，下面给出计算实例。假设有一信号包含两个谐波频率值，分别为413.0Hz（幅值设为1.00V）和287.0Hz （幅值设为0.400V），利用函数发生器产生这两路信号再混合，加上偏置后送入AD采样。设采样频率1024Hz，共采样2048点。图1的采样点均是右对齐的12位采样结果值。图 1 AD采样得到的采样点图利用TI的FFT库进行计算，查看mag数组，得到图2.以1024Hz采样2048点，采样时间2s

14、，对两个谐波频率可采样到整数倍周期；从另一个角度理解，此时最小频率分辨率为0.5Hz，413.0Hz与287.0Hz均是其整数倍数，故不会发生频谱展宽或混叠情况，计算得到的频谱图应该为两根尖峰线。从图2结果也能看出这一点。查看mag数组，可知第一根尖峰线下标574，第二根尖峰线下标826，故真实频率值分别为：574*0.5=287Hz , 826*0.5= 413Hz若要计算幅值，按照输出的Q30格式除相应系数即可。注意最好另外定义浮点数组来做除运算，因为整型数据做除运算（或者右移位操作）会丢失小数位数据。若要验证DSP的FFT计算结果，可将AD数据从CCS导入到MATLAB中做对比计算。有些

15、同学不清楚如何导入导出，下面说一下步骤。1.点菜单栏file-àdata-àsave，选择保存类型Integer(若是其他进制还需在MATLAB中转换)，点“确定”后，在“address”栏填入要保存数据的起始地址（填变量数组名或真实存储器地址皆可），在“length”栏内填入数据长度，“page”肯定选“data”页了；全部设好后点“OK”2.找到保存的数据文件，将后缀改成.txt，再用记事本打开（也可以不改后缀直接用记事本打开），删去第一行数据；3.打开MATLAB，点开“load data file”，选中刚才的数据文件，然后按照提示一步步往下导入即可。最后不妨用变量

16、temp来保存这些数据。在MATLAB内运行如下代码：%f_sample = 1024; %采样频率N = 2048; %采样点Ny=temp; %temp即为导入数据的变量名n = 0:N-1;t = n/f_sample;%做采样点图plot(t,y);figure;stem(t,y,'.');figure;%fft变换并作图fft_result = fft(y);mag = abs(fft_result)/(N/2);mag(1)=0; %为观察谐波分量，此处特意将直流分量置为0 %求解真实频率值f = (0:length(fft_result)-1)*f_sample/

17、length(fft_result);%作频谱图stem(f(1:N/2),mag(1:N/2),'b.');grid;%可将MATLAB计算结果与CCS内的结果做些对比。附：存储器分配（.cmd文件内容），我用的是2808的板子，2812或其他的稍作改动即可。MEMORYPAGE 0 :BEGIN : origin = 0x000000, length = 0x000002 RAMM0 : origin = 0x000002, length = 0x0003FEPRAMH0 : origin = 0x3FA000, length = 0x002000 RESET : orig

18、in = 0x3FFFC0, length = 0x000002BOOTROM : origin = 0x3FF000, length = 0x000FC0 TESARAM : origin = 0x008000, length = 0x002000 PAGE 1 : BOOT_RSVD : origin = 0x000400, length = 0x000080 RAMM1 : origin = 0x00A000, length = 0x001000 L0L1RAM : origin = 0x00B000, length = 0x001000HL0SARAM : origin = 0x3F8

19、000, length = 0x001010 DRAMH0 : origin = 0x3F9010, length = 0x000ff0 SECTIONS/* Setup for "boot to SARAM" mode: The codestart section (found in DSP28_CodeStartBranch.asm)re-directs execution to the start of user code. */codestart : > BEGIN, PAGE = 0ramfuncs : > RAMM0 PAGE = 0 .text :

20、 > PRAMH0, PAGE = 0.cinit : > RAMM0, PAGE = 0.pinit : > RAMM0, PAGE = 0.switch : > RAMM0, PAGE = 0.reset : > RESET, PAGE = 0, TYPE = DSECT /* not used, */FFTipcb ALIGN(4096): > HL0SARAM, PAGE = 1FFTtf :> DRAMH0 , PAGE = 1FFTmag :> L0L1RAM, PAGE = 1.const : > DRAMH0, PAGE =

21、 1.bss : > DRAMH0, PAGE = 1.stack : > RAMM1, PAGE = 1.sysmem : > RAMM1, PAGE = 1.ebss : > DRAMH0, PAGE = 1.econst : > PRAMH0, PAGE = 0 .esysmem : > L0L1RAM, PAGE = 1IQmath : > PRAMH0, PAGE = 0IQmathTables : > BOOTROM, type = NOLOAD, PAGE = 0另外，有些同学一开始运行FFT附带的例程，会遇上两个问题：一是缺少函数

22、文件，这个到TI官网下那个sprc083_SGEN包就好了，或找我我用邮箱发给大家也行；二是提示编译不成功，找到上述配置文件的这一行.reset : > RESET, PAGE = 0, TYPE = DSECT /* not used, */在后面加上TYPE = DSECT /* not used, */。这些是我凭印象写的，因好久没做这个FFT，可能会写出错误，欢迎大家提出来：-）有问题就解决问题之 TI库实现2048点或更大点数FFT 2010-05-18 13:20 星期二TI提供的FFT库（sprc081）使用起来很方便，只要直接调用即可；更重要的是其中做了好些优化，而且是用

23、汇编写成，运算速度肯定要比自己写的C代码更快好多。这个库最大能运算1024点的实数FFT，但实际中我们可能还需要更大点数的FFT，那如何去修改呢？其实修改的地方并不多，以前我曾花了一段时间浏览了一下那些头疼的汇编代码，发现模块化写得很好啊（佩服TI的工程师！），这样修改起来就快多了。1 2048点FFT直接在工程中加入原来的库文件fft.lib，并在代码中作如下修改：第一步：在有关FFT代码的头文件fft.h中增添如下一段代码：#define RFFT32_2048P_DEFAULTS (long *)NULL, (long *)NULL, 1024, 10, (long *)NULL, (l

24、ong *)NULL, 0, 0, 1, (void (*)(void *)CFFT32_init, (void (*)(void *)CFFT32_calc, (void (*)(void *)RFFT32_split, (void (*)(void *)RFFT32_mag, (void (*)(void *)RFFT32_win结构体RFFT32关键说明：1 第一个(long *)NULL,，是计算缓冲区（ipcb）的指针；2 第二个(long *)NULL,，是查表旋转因子（ft）的指针；3 1024指的是FFT的大小，注意这里是2048点实数FFT被打包成1024点的复数，这是这个库

25、的实数FFT运算的特点，这样节省一半运算时间和存储空间；4 10即1024是2的10次方，因为这里采用的是基2时间抽取算法；5 下面的两个(long *)NULL,分别是运算结果中存放幅值（mag）以及加窗表格（win）的指针；6 下面的两个0分别是计算结果中的峰值幅值与峰值频率，这里初始化0；7 下面的1指的是查表的步进值，因在这个sprc081文档的库已经有Q31格式（以及Q30格式）的1024点的正弦因子表，因此对于1024点复数FFT查表步进值设为1即可。8 下面其他的就是要调用的函数指针了，不再赘述。明白上述结构体的定义后，改写其他点数就很容易了，比如这个库最少能做到128点的实数F

26、FT，我们还可以相应改成64点或者更低点数的实数FFT，比如64点的结构体就是如下的一段：#define RFFT32_64P_DEFAULTS (long *)NULL, (long *)NULL, 32, 5, (long *)NULL, (long *)NULL, 0, 0, 32, (void (*)(void *)CFFT32_init, (void (*)(void *)CFFT32_calc, (void (*)(void *)RFFT32_split, (void (*)(void *)RFFT32_mag, (void (*)(void *)RFFT32_win第二步：修改源

27、程序代码里的变量定义：#include "fft.h"#define N 2048 /定义N的点数，2048点/分别在数据存储区为ipcb和mag分配足够的地址，具体见下面的连接器文件#pragma DATA_SECTION(ipcb, "FFTipcb"); #pragma DATA_SECTION(mag, "FFTmag");/定义变量RFFT32 fft=RFFT32_2048P_DEFAULTS; long ipcbN+2;long magN/2+1;/其他部分函数调用都是一样的，这里不再列出第三步：在连接器（.cmd）文件

28、里分配存储地址：这里给出我在2808上的分配方法：MEMORYPAGE 1 : L0L1RAM : origin = 0x00B000, length = 0x001000 HL0SARAM : origin = 0x3F8000, length = 0x001010 DRAMH0 : origin = 0x3F9010, length = 0x000ff0SECTIONS FFTipcb ALIGN(4096): > HL0SARAM, PAGE = 1 FFTtf :> DRAMH0 , PAGE = 1 FFTmag :> L0L1RAM, PAGE = 1好，到此为止

29、就修改好了，可以使用这个FFT库实现2048点实数FFT了。2 4096点或更大点数的FFT第一步：修改汇编代码并制成库文件 其实要实现4096点的FFT需要改动的地方也不多，主要因为原来的这个库只提供了1024点的正弦值，我们需要自己算出2048点的Q31（或Q30）格式的一周期正弦值，并将其代替原来汇编代码里的1024点的表格。汇编代码的其他地方不需要改动，不然可能产生不可预知的后果。在汇编代码文件sel_q.asm里查看设置：; Select Twiddle factor Q formatTF_QFMAT .set Q30 若是TF表设为Q30，则只需将计算出的正弦值移位成Q30格式即可

30、。另外需注意的是，2048点的FFT需要2048点的TF表，但这个库将其优化成只需2048×3/4 = 1536点，即1.5K大小的表格（见sprc081文档第10页，注意到正弦的第二部分与余弦的第一部分是相重合的，这样便节省了1/4周期的点数）。计算好正弦值并加入到汇编代码里，然后将所有的汇编代码导进CCS里做成一个新的库，到时直接添加到工程里调用。为了方便这里我已将其做成一个新的库命名为FFT_new.lib，大家直接用即可。第二步：修改头文件fft.h 结构体定义增添如下一段代码：#define RFFT32_4096P_DEFAULTS (long *)NULL, (long *)NUL