利用低功耗微控制器开发FFT应用

摘要:今天的低功耗微控制器(μc)也开始集成原先只存在于大型微处理器、asic和dsp中的外设功能,使我们有可能以很低的功耗实现复杂的算术运算。本文讨论一种快速傅立叶变换(fft)应用,并在一个含有单周期硬件乘法器的低功耗μc上实现该应用。
这个fft应用实时计算一路输入电压(图1中的vin)的频谱。为完成该任务,用一片模数转换器(adc)对vin进行采样,获得的采样传送给µc。然后,µc对这些采样执行256点fft运算,获得输入电压的频谱。为便于检测,µc将计算出的频谱数据传送给pc,由pc实时显示出来。
图1. 利用fft应用计算输入电压的频谱。
该fft应用的固件针对maxq2000系列中的一款16位、低功耗µc用c语言编写。有兴趣的读者可以下载(zip,2.4kb)该项目的固件和电路原理图。 背景知识为确定输入信号采样的频谱,我们需要对这些输入采样进行离散傅立叶变换(dft)。dft的定义如下:
其中n是采样的数量,x(k)是频谱,x(n)是一组输入采样。利用欧拉等式展开求和符,并分离输入采样和频谱的实部和虚部,得到以下等式:
式2和3中,求和符中第二项的消失是由于输入采样全部为实数。假定我们有n个采样,直接计算式2和3需要2n²次乘法和2n(n - 1)次加法。这样,我们的256点输入采样dft将需要进行131,072次乘法和130,560次加法运算。我们还是将注意力转向fft吧!
有多种fft算法可供使用。本应用采用普通的radix-2算法,继续将dft分解为两个更小的dft。为此,n必须是2的指数。这种radix-2 fft算法的步骤可归纳如图2所示的蝶型运算。观察这些蝶型运算我们可以发现,radix-2算法仅需(n / 2)log2(n)次乘法和nlog2(n)次加法。图2中用到的参数wn就是通常所谓的“旋转因子”,可以在执行算法前预先计算出来。
图2. 利用蝶型运算实现n = 8的fft。
在图2中,fft的输入显示为一种特殊的排列顺序,这种序列是对原始序列索引号的二进制位反转后得到的。因此,当我们对n = 8个采样执行radix-2 fft算法时,需要将输入数据的原始序列:
0 (000b), 1 (001b), 2 (010b), 3 (011b), 4 (100b), 5(101b), 6(110b), 7(111b)
重新排列为:
0 (000b), 4, (100b), 2 (010b), 6 (110b), 1 (001b), 5 (101), 3 (011), 7 (111)
fft输出则以正确的顺序排列。图2还说明,每个单独的蝶型运算所得的结果,是下一级fft运算所需的唯一数据。由于运算过程可“即位”进行,新值可替代旧值,这样,计算n个采样的fft只需要2n个变量(因为每个数据都包括实部和虚部两部分)。
fft完成后,结果为复数形式。式4和5将结果转换为极坐标方式后表示为:
有关dsp的文献中可以找到很多优化方法,可使上述dft/fft算法更小或更快。其中最重要的一种优化方法(可能也是最容易实现的)源于这样一个事实,那就是作为一个实数信号,其dft幅度是相关于x(n / 2)对称的,因此:
编写fft代码绝非易事。低功耗µc的一些局限又进一步使该任务复杂化。
存储器:我们所选的µc有2kb的ram。已经知道该算法需要用到2n个16位变量来存储fft数据,这样,我们的µc可以执行n最高为512的fft。然而,固件的其他部分也要用到一些ram。因此,在此项目中,我们限制n于256。若采用16位变量来表示每个值的实部和虚部,fft数据总共需要1024字节的ram。
速度:低功耗µc尽管具有高mips/ma性能,仍然需要一些优化手段来使运行fft的指令数尽可能少。好在本应用所用的c编译器(iar的embedded workbench for maxq,见www.iar.com)可提供多种级别的优化和设置。高效地使用硬件乘法器可使代码优化到可以接受的水平。
无浮点能力:所选的µc不具备浮点能力(低功耗产品一般都不具备浮点能力)。因此,所有运算都必须采用定点算法。为了表示小数,固件采用带符号的q8.7表示法。这样,在固件中假定: 第0位至第6位代表小数部分 第7位至第14位代表整数部分 第15位代表符号位(二的补码) 这样的安排对于加法和减法没有影响,但在做乘法时必须注意将数据按照q8.7格式对齐。
所选的数据表示法还要适应fft算法可能遇到的最大数值,同时又要提供足够的精度。例如,我们的adc可提供带符号的8位采样,以二的补码表示。如果输入为最大幅度(对于带符号8位采样为127)的直流电压,则其能谱全部包含于x(0)中,用q8.7表示为32512。这个数值能够由单个带符号的16位数据表示。 固件以下部分讨论在低功耗µc上执行radix-2 fft的固件实现。信号采样由adc读出后被存储在x_n_re数组中。这个数组代表x(n)的实部。虚部存储在x_n_im数组中,在开始运行fft前初始化为零。完成fft后,计算结果取代原始采样数据,被存储在x_n_re和x_n_im中。 获取采样fft算法假定采样是以固定的取样频率获得的。在为fft获取采样时如果不加小心将会产生一些问题。例如,采样间隔的抖动就会给fft结果引入误差,应尽力减小之。
adc采样循环中的判决语句会造成采样间隔的抖动。例如,我们的系统从adc读取带符号的8位采样,并将其存储在一组16位变量中。在下面的程序清单1中给出了两种伪码算法,执行这种adc读取-存储功能。算法1给出的方法会造成采样间隔的抖动,因为负采样比正采样需要更多的时间来读取并存储。
清单1. 两种adc采样伪码算法。第二种算法避免了第一种的问题――采样间隔抖动。 // algorithm 1: inconsistent sampling frequency - bad! // sample[] is an array of 16-bit variables for i = 0 to (n-1) begin doadcsampleconversion() // instruct adc to sample vin sample[i] = read8bitsamplefromadc() // read 8-bit sample from adc if (sample[i] & 0x0080) // if the 8-bit sample was negative sample[i] = sample[i] + 0xff00 // make the 16-bit word negative end // algorithm 2: fixed sampling frequency - good! // sample[] is an array of 16-bit variables for i = 0 to (n-1) begin doadcsampleconversion() // instruct adc to sample vin sample[i] = read8bitsamplefromadc() // read 8-bit sample from adc end for i = 0 to (n-1) begin if (sample[i] & 0x0080) // if the 8-bit sample was negative sample[i] = sample[i] + 0xff00 // make the 16-bit word negative end 三角函数表本fft算法通过查表(lut)而非计算得到正弦或余弦函数值。程序清单2给出了对于正弦和余弦lut的申明。实际固件的注释中包含了自动生成这些lut的源代码,可由程序调用。两个lut均含有n / 2分量,因为旋转因子的索引号变化范围为从0至n / 2 - 1 (见图2)。
清单2. 正弦和余弦函数lut。 const int coslut[n/2] = {+128,+127,+127, ... ,-127,-127,-127}; const int sinlut[n/2] = {+0 ,+3 , +6, ... ,+9 , +6, +3}; 这些lut中的数组被声明为const,强制编译器将它们存储于代码空间而非数据空间。由于lut数值须采用q8.7表示法,它们由正弦和余弦的实际值乘以27后得到。 位反转位反转排序(n已知)可在运行时通过计算、查表或直接利用展开循环编写。所有这些方法都需要在源代码的尺寸和运行速度间进行折衷。本fft应用利用展开循环进行位反转,其源代码较长,但运行速度快。程序清单3显示了该展开循环的实现。本应用固件的注释中包含了用于程序自动生成展开循环的源代码。
清单3. 用于实现n = 256的位反转的展开循环。 i=x_n_re[ 1]; x_n_re[ 1]=x_n_re[128]; x_n_re[128]=i; i=x_n_re[ 2]; x_n_re[ 2]=x_n_re[ 64]; x_n_re[ 64]=i; i=x_n_re[ 3]; x_n_re[ 3]=x_n_re[192]; x_n_re[192]=i; i=x_n_re[ 4]; x_n_re[ 4]=x_n_re[ 32]; x_n_re[ 32]=i; ... i=x_n_re[207]; x_n_re[207]=x_n_re[243]; x_n_re[243]=i; i=x_n_re[215]; x_n_re[215]=x_n_re[235]; x_n_re[235]=i; i=x_n_re[223]; x_n_re[223]=x_n_re[251]; x_n_re[251]=i; i=x_n_re[239]; x_n_re[239]=x_n_re[247]; x_n_re[247]=i; radix-2 fft算法采样按照位反转方式重新排序后就可进行fft运算了。本radix-2 fft应用的固件通过三个主循环执行图2所示的蝶型运算。外循环计数log2(n)级fft运算。内循环执行每一级的蝶型运算。
fft算法的核心部分是执行蝶型运算的一小块代码。程序清单4给出了这一块代码,遗憾的是,它是本应用中唯一“不可移植”的固件。宏mul_1和mul_2利用µc的硬件乘法器执行单指令周期乘法运算。这些宏的内容专用于maxq2000,可在实际固件中全部看到。
清单4. 用c编写的蝶型运算。 /* (1) macro mul_1(a,b,c): c=a*b (result in q8.7)*/ /* (2) macro mul_2(a,c) : c=a*last_b (result in q8.7)*/ mul_1(coslut[tf],x_n_re[b],resultmulrecos); mul_2(sinlut[tf],resultmulresin); mul_1(coslut[tf],x_n_im[b],resultmulimcos); mul_2(sinlut[tf],resultmulimsin); x_n_re[b] = x_n_re[a]-resultmulrecos+resultmulimsin; x_n_im[b] = x_n_im[a]-resultmulresin-resultmulimcos; x_n_re[a] = x_n_re[a]+resultmulrecos-resultmulimsin; x_n_im[a] = x_n_im[a]+resultmulresin+resultmulimcos; 复数的极坐标转换为了便于确定vin频谱的幅度,我们须要将复数形式的x(k)转换为极坐标形式。实现该转换的固件示于程序清单5。幅度值取代了原始的fft结果,因为固件不再需要这些数据。
清单5. fft结果从复数形式转换为极坐标形式。 const unsigned char magnlut[16][16] = { {0x00,0x10,0x20, ... ,0xd0,0xe0,0xf0}, {0x10,0x16,0x23, ... ,0xd0,0xe0,0xf0}, ... {0xe0,0xe0,0xe2, ... ,0xff,0xff,0xff}, {0xf0,0xf0,0xf2, ... ,0xff,0xff,0xff} }; ... ... /* compute x_n_re=abs(x_n_re) and x_n_im=abs(x_n_im) */ ... ... x_n_re[0] = magnlut[x_n_re[0]>>11][0]; for(i=1; i>11][x_n_im[i]>>11]; x_n_re[n_div_2] = magnlut[x_n_re[n_div_2]>>11][0]; 频谱幅度并非根据式4计算得到,而是通过一个二维lut查表得到。第一索引为频谱实部的高4位(msb),第二索引为频谱虚部的高4位。为得到这些数据,可将带符号的16位数据右移11次。在从频谱的实部和虚部取得索引号前,需首先将它们转换为绝对值。因此,符号位为零。
从式6我们已经知道,频谱的幅度是关于x(n / 2)对称的,因此我们只需将前(n / 2) + 1个频谱数据转换为极坐标形式。还有,我们可以看到,对于实数输入采样,x(0)和x(n / 2)的虚部总为零。因此这两条谱线的幅度被单独计算。本项目实际固件的注释中包含了用于自动生成该lut的源代码,可由程序调用来计算x(k)的幅度。 hamming或hann窗此项目固件还包括了对输入采样加hamming或hann窗的lut (q8.7格式)。加窗函数可有效降低对时域采样x(n)的舍入操作所引起的频谱泄漏。hamming和hann窗函数分别如式7和8所示。
程序清单6给出了实现这些函数的代码。同样,本项目实际固件的注释中包含了用于自动生成这些lut的源代码,可由程序调用来实现这些窗函数。
清单6. 用来实现hamming和hann窗函数的lut。 const char hamminglut[n] = {+10, +10, +10, ... ,+10, +10, +10}; const char hannlut[n] = { +0, +0, +0, ... , +0, +0, +0}; ... ... for(i=0; i<256; i++) { #ifdef windowing_hamming mul_1(x_n_re[i],hamminglut[i],x_n_re[i]); // x(n)*=hamming(n); #endif #ifdef windowing_hann mul_1(x_n_re[i],hannlut[i]),x_n_re[i]); // x(n)*=hann(n); #endif } 测试结果为了测试该fft应用的性能,固件将x(k)幅度通过µc的uart端口上传给pc。专门编写的fft graph软件(随该项目固件一起提供)用于从pc串口读取这些幅值,并以图形方式实时显示频谱。图3显示了µc以200ksps采样四种不同输入信号并处理后,由fft graph所显示出来的结果: 4.3v直流信号 50khz正弦信号 70khz正弦信号 6.25khz方波
图3. fft graph软件显示的由低功耗µc计算出的频谱。 接下来干什么?有兴趣的读者还可以花费大量的时间来继续优化和重新配置该fft应用。尽管在本文中我们选择了radix-2算法,还有很多其他算法可以显著降低加法和乘法运算量。很多本文所未提及的优化可以提升fft的速度。例如,作为纯实数的输入采样,其虚部总为零,频谱中只有前半部分有实际意义。利用这一点,第一级和最后一级fft的执行速度可进一步优化,但需要付出更多的程序空间。
总之,本文所讨论的算法对于低功耗µc上的fft应用而言,提供了一个很好的出发点。如果想了解更多信息和具体实现的细节,请查阅我们为本应用所提供的、带有详细注释的固件信息。

大数据正在推动更大的容量
abb变频器历史故障怎么查
纸张在线缺陷检测仪的原理、参数及功能
i2c时序图的详细讲解
SpringBoot部署打包成jar和war有什么不同呢?
利用低功耗微控制器开发FFT应用
王毅与荷兰谈光刻机出口问题
新一轮的换机潮是由什么引起的
科锐推出业界首项LED模组驱动兼容计划
如何使用蓝牙将Android手机中的照片和视频副本发送到树莓派
分析师称苹果新款iPhone将采用触摸感应外壳
一图了解华为MateX 1.75万元到底值不值
马斯克的星链互联网已实现初步组网功能
龙源大丰海上风电项目,在江苏省已投运海上风电项目中独占鳌头!
直接使用惯性测量单元 (IMU)
飞睿智能人体存在感应雷达方案
英特尔至强可扩展处理器有助于医疗保健维护数据隐私和合规性
基于AMD嵌入式锐龙的GPD Win Max便携掌机
电动自行车的霍尔线是什么意思,它和电机有什么关系?
车载扬尘监测站:可以移动的气象站