您当前的位置:五五电子网电子知识单片机-工控设备DSP/FPGA技术利用低功耗微控制器开发FFT应用 正文
利用低功耗微控制器开发FFT应用

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

点击数:7397 次   录入时间:03-04 12:00:14   整理:http://www.55dianzi.com   DSP/FPGA技术

       清单1. 两种ADC采样伪码算法。第二种算法避免了第一种的问题——采样间隔抖动。

// ALGORITHM 1: INCONSISTENT SAMPLING Frequency - BAD!
// sample[] is an array of 16-bit va riables
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的位反转的展开循环。



www.55dianzi.com 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; ix_n_re[i] = magnLUT[x_n_re[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次。在从频谱的实部和虚部取得索引号前,需首先将它们转换为绝对值。因此,符号位为零。



www.55dianzi.com        从式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方波

            FFT Graph软件显示的由低功耗μC计算出的频谱

                             图3. FFT Graph软件显示的由 低功耗 μC计算出的频谱。

       接下来干什么?

       有兴趣的读者还可以花费大量的时间来继续优化和重新配置该FFT应用。尽管在本文中我们选择了radix-2算法,还有很多其他算法可以显著降低加法和乘法运算量。很多本文所未提及的优化可以提升FFT的速度。例如,作为纯实数的输入采样,其虚部总为零,频谱中只有前半部分有实际意义。利用这一点,第一级和最后一级FFT的执行速度可进一步优化,但需要付出更多的程序空间。

       总之,本文所讨论的算法对于低功耗μC上的FFT应用而言,提供了一个很好的出发点。如果想了解更多信息和具体实现的细节,请查阅我们为本应用所提供的、带有详细注释的固件信息。



上一页  [1] [2] 


本文关键字:控制器  开发  DSP/FPGA技术单片机-工控设备 - DSP/FPGA技术