DSP

TI DSP FFT算法的应用

2019-07-13 10:35发布

这篇文章是应一些找我讨论DSP的同学所写,贴在这里大家一起学习。
曾有不少论坛上的同学(包括DSP算法讨论群里的一些同学)问过我关于TI的FFT库的使用,这里我将我使用过的一些经验说一下。
TI的这个FFT库在计算速度、计算精度以及数据存储等方面是做了不少优化的,比如数据存储,若作N点的FFT,供查表用的旋转因子必须有N/2点的正弦值与N/2点的余弦值,这个库将其压缩成3N/4点的正弦值,因此就节省了N/4点的存储空间;另外计算N点实数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)有如下一句: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 ipcb[N+2];
因为AD结果寄存器是12位的,在数据左对齐的情况下直接左移15位即可:
ipcb[ConversionCount] = ((unsigned long)AdcRegs.ADCRESULT0)<<15;
另外我记得有份文档提到在某些存储器下是右对齐,此时则需左移19位,大家在使用时注意这个问题。
当采样完成后,按照FFT库的文档上的说明或者仿照文档附带的例程,进行相应函数调用来实现自己的FFT计算,比如按计算点数来定义各个变量数组,是否加窗,是否求解幅值平方值等。
在进行存储器分配时,文档上要求(128点实数FFT为例):
FFTipcb ALIGN(256) : { } > L0L1RAM PAGE 1
FFTmag > L0L1RAM PAGE 1
FFTtf > 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若因存储冲突肯定会造成查表值不准确,那计算FFT时肯定就不对了。有个论坛帖子作者说一定要放在origin = 0x008000开始位置,其实也不对,大家可自己去试验;后面我也会给出我的存储配置文件(即.cmd文件内容)。
无图无真相,下面给出计算实例。
假设有一信号包含两个谐波频率值,分别为413.0Hz(幅值设为1.00V)和287.0Hz (幅值设为0.400V),利用函数发生器产生这两路信号再混合,加上偏置后送入AD采样。设采样频率1024Hz,共采样2048点。图1的采样点均是右对齐的12位采样结果值。



图 1 AD采样得到的采样点图

利用TI的FFT库进行计算,查看mag数组,得到图2.



图 2 FFT计算结果(mag数组)

以1024Hz采样2048点,采样时间2s,对两个谐波频率可采样到整数倍周期;从另一个角度理解,此时最小频率分辨率为0.5Hz,413.0Hz与287.0Hz均是其整数倍数,故不会发生频谱展宽或混叠情况,计算得到的频谱图应该为两根尖峰线。从图2结果也能看出这一点。
查看mag数组,可知第一根尖峰线下标574,第二根尖峰线下标826,故真实频率值分别为:
574*0.5=287Hz , 826*0.5= 413Hz
若要计算幅值,按照输出的Q30格式除相应系数即可。注意最好另外定义浮点数组来做除运算,因为整型数据做除运算(或者右移位操作)会丢失小数位数据。
若要验证DSP的FFT计算结果,可将AD数据从CCS导入到MATLAB中做对比计算。有些同学不清楚如何导入导出,下面说一下步骤。
1. 点菜单栏file-data-save,选择保存类型Integer(若是其他进制还需在MATLAB中转换),点“确定”后,在“address”栏填入要保存数据的起始地址(填变量数组名或真实存储器地址皆可),在“length”栏内填入数据长度,“page”肯定选“data”页了;全部设好后点“OK”
2. 找到保存的数据文件,将后缀改成.txt,再用记事本打开(也可以不改后缀直接用记事本打开),删去第一行数据;
3. 打开MATLAB,点开“load data file…”,选中刚才的数据文件,然后按照提示一步步往下导入即可。最后不妨用变量temp来保存这些数据。
在MATLAB内运行如下代码:
%%%%%%%%%%%%%%%%%%
f_sample = 1024; %采样频率
N = 2048; %采样点N
y=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/length(fft_result);

%作频谱图
stem(f(1:N/2),mag(1:N/2),'b.');grid;
%%%%%%%%%%%%%%%%%%



可将MATLAB计算结果与CCS内的结果做些对比。

附:
存储器分配(.cmd文件内容),我用的是2808的板子,2812或其他的稍作改动即可。
MEMORY
{
PAGE 0 :
BEGIN : origin = 0x000000, length = 0x000002 
RAMM0 : origin = 0x000002, length = 0x0003FE
PRAMH0 : origin = 0x3FA000, length = 0x002000 
RESET : origin = 0x3FFFC0, length = 0x000002
BOOTROM : 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 = 0x001000
HL0SARAM : origin = 0x3F8000, 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 = 0
ramfuncs : > RAMM0 PAGE = 0 
.text : > 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 = 1

FFTtf :> DRAMH0 , PAGE = 1


FFTmag :> L0L1RAM, PAGE = 1
.const : > DRAMH0, PAGE = 1
.bss : > DRAMH0, PAGE = 1
.stack : > RAMM1, PAGE = 1
.sysmem : > RAMM1, PAGE = 1

.ebss : > DRAMH0, PAGE = 1
.econst : > PRAMH0, PAGE = 0 
.esysmem : > L0L1RAM, PAGE = 1

IQmath : > PRAMH0, PAGE = 0
IQmathTables : > BOOTROM, type = NOLOAD, PAGE = 0
}

另外,有些同学一开始运行FFT附带的例程,会遇上两个问题:
一是缺少函数文件,这个到TI官网下那个sprc083_SGEN包就好了,或找我我用邮箱发给大家也行;
二是提示编译不成功,找到上述配置文件的这一行
.reset : > RESET, PAGE = 0, TYPE = DSECT /* not used, */
在后面加上TYPE = DSECT /* not used, */。

这些是我凭印象写的,因好久没做这个FFT,可能会写出错误,欢迎大家提出来:-) 不少同学加我QQ问我关于TI的FFT的使用方法,因为他们总会碰上各种各样的问题——我也不知道为什么实际使用会有这么多问题,有时我也答不上来其中缘由——至少对我来说,我当时还用得挺顺利的,没有这么多问题啊。
   因此这里我将我在实验项目中写过的FFT模块代码贴上来,与大家一起分享一起学习。其实这段代码并没有什么深度或难解的地方,就是不断地调用TI写好的函数而已。这里没有用到Acquisition Modules也没有用到加窗,因此相关代码都注释了。
  
  
  //File name: FFT_module.c
  //Programmed by Sun Zy @ XJTU on 2009-12-03
  
  #include "stb.h"
  #include "fft.h" 
  #include "MyInclude.h"
  
  #define N 2048
  #pragma DATA_SECTION(ipcb, "FFTipcb");
  #pragma DATA_SECTION(mag, "FFTmag");
  
  RFFT32 fft=RFFT32_2048P_DEFAULTS; 
  long ipcb[N+2];
  long mag[N/2+1]; 
  int FFTEndFlag = 0;
  
  //Configure the FFT module
  void FFT_config()
  {
   /* Initialize acquisition module 
   acq.buffptr=ipcb;
   acq.tempptr=ipcb;
   acq.size=N;
   acq.count=N;
   acq.acqflag=1;*/
  
  /* Initialize FFT module */
   fft.ipcbptr=ipcb;
   fft.magptr=mag; 
  // fft.winptr=(long *)win;
   fft.init(&fft);
  }
  
  //Do the FFT computation
  void Compute_FFT()
  {
  RFFT32_brev(ipcb,ipcb,N);
   fft.calc(&fft);
   fft.split(&fft);
   fft.mag(&fft); 
  
   FFTEndFlag = 1; //Set end flag
  }