FFT在低功率微程序控制器中的应用

0 引言
    在以前，外围设备是更大的微处理器如特定用途集成电路和DSP中的内存。现在，低功率微处理器包括了外围设备，这样就有机会在低功耗的情况下进行复杂的运算。本文介绍了在低功率的微处理器中执行快速傅里叶变换(FFT)，其中微处理器包括一周期的硬件乘法器。这个应用可以实时计算输入电压的频谱。
    为了完成此任务，一个模数转换器对输入信号进行采样然后传输到微处理器。微处理器再对样本信号进行256点的FFT，这样就获得输入电压的频谱。为了测试其有效性，微处理器计算频谱的幅值然后实时地传输给示波器。

1 背景
    为了确定输入样本信号的频谱信息，需要计算输入样本的离散傅里叶变换(DFT)。离散傅里叶变换定义为：
   
式中：N是样本点数；X(k)是频谱，与x(n)代表输入样本。利用欧拉方程的一致性将这个求和公式中的输入样本与频谱分离为它们的实部与虚部，可得以下方程式：
   
    因为输入样本只是考虑实部。式(2)与式(3)中的求和公式的第二项消失了，假设有N个样本，直接计算式(2)、式(3)需要2N2次乘法及2N(N-1)次加法。因此256点输入样本的DFT将要求131072次乘法和130560次加法。
    已经出现了很多种FFT算法。普通的以基为2的算法连续将DFT分解成2个更小的DFT。为了使其变成可能，N必须分解为2的整数幂。转化为以2为基的FFT的步骤见图1的蝶形计算。从图1的蝶形计算中可观察到，获得基为2的FFT算法的解只需要(N／2)log2N次乘法与Nlog2N次加法。在图l中的值WH通常认为是旋转因子且能够在执行FFT前计算得到。

    在图1中，FFT的输入具有特殊的形式。它是具有位倒置下标的原始顺序。因此，当计算基为2的N=8的FFT时，输入数据的记录顺序要求为O(000b)，1(001b)，…，0(000b)，4(100b)，…。
    FFT是以正确的顺序作为输出。图1同样揭示了单一的蝶形计算的结果只是FFT的下一阶段的输入。因为计算是在适当的位置中完成的，旧值可以代替新获得值且在计算N点的FFT只是需要2N个变量样本(需要2N个变量是因为每一个变量值都有一个实部与虚部)。
    当完成FFT时，结果是以复数为记法的。式(4)和式(5)将复数表示形式转变为以极坐标表示：
   
    在DSP的文章里介绍了很多关于DFT／FFT的优化方法，使其计算速度更快且需要的计算量更少，其中一个比较重要的优化方法(也可能是最容易执行的)。
    从观察DFT中可以获得，因为具有N点的实值信号的DFT是以X(N／2)为对称的，因此有：
    [p]
2 执行要点
    写代码实现DFT不是一件容易的事，因为用低功率的微处理器实现DFT算法的实际情况是相当复杂的。例如，这些微处理器通常：
    (1)有限的内存。选择的微处理器只有2 KB的RAM。从上面叙述可知实现FFT至少需要2N×16 B变量。微处理器不能执行样本点数N大于512的FFT。这是不现实的，因为别的固件同样需要一些字节的RAM。因此在实际执行的过程中，通常将样本点数局限在256点。使用16 B的变量表示每一个值的实部与虚部，这种情况下对于FFT的数据需要1024B的RAM。
    (2)有限的速度。尽管低功率的微处理器具有高达每秒百万条指令的速度，仍然需要一些优化方法来减少在执行FFT过程中所用到的指令。所幸的是在应用过程中。C编译器包括很多优化的级别设置。小心使用芯片的硬件乘法同样可以使得代码优化到一个可以接受的水平。
    (3)没有浮点数功能。所选择的微处理器特别是那些低功率的微处理器没有浮点数功能。因此所有的计算都需要定点算法。为了表示分数，固件将使用有符号的Q8．7标记。因此固件将假设：O～6 B表示每一数字的分数部分；7～14 B表示每一数字的整数部分；第15字节是这个数字的符号位。
    这种形式对于加法和减法是没有影响的，但是对于乘法则必须注意，使所有数据排成Q8．7的形式。例如对于Q8．7的乘法如下：

    为了获得比较精确的FFT结果，Q8．7排列形式的一致性同样适用于具有比较大样本点数的FFT。例如，模／数转换器以实部和虚部互补的形式提供8位的符号数。如果输入的是直流电压(+127为有符号的8位样本数)，从X(0)中将会完全获得其频谱，以Q8．7标记等于32512。这个值很适合于用16位的符号数表示。

3 固件
    下面介绍计算基为2的FFT所需的固件。当从模／数转换器中读取样本数后，存储在数组x_n_re中。这个数组表示x(n)的实部。在执行FFT前，虚部的值初始化为零，存储在数组x_n_im中。当执行完FFT时，频域的幅值将代替原来的样本值，且存储在x_n_re和x_n_im中。
3．1 采集样本
    FFT算法假设以固定采样率来采集样本的。尽管这是在本文考虑范围之外，但是如果认真对待采集样本的代码同样会产生问题。例如，不稳定的采样率将会产生错误的FFT结果，所以应该尽量使该情况最小化。对模／数转换器采样的原码每一次循环以及输出结果命令都有可能对采样率产生不稳定性。例如，系统从摸／数转换器中读取8位的有符号数，然后存储在16位的数组变量中。
    下面列出了关于从模／数转换器中读取及存储数据的2个伪码算法。第1个记为算法l，将会引起采样率的不稳定。因为负数样本比正的样本需要更多的时间来读取及存储。中断同样不能保证采样代码的健全。
    模／数转换器采样(ADC)的伪码：
    算法1：不一致的采样率。
   
3．2 三角法来查寻表格
    FFT利用查寻表的方法(LUTs)来代替直接计算正弦与余弦的值。下面叙述中给出了对正弦与余弦的LUTs的声明。固件中的声明包括在程序中自动产生这些表格的原始代码。LUTs中的正弦与余弦都具有N／2样本，因为旋转因子的下标从0～(N／2)-1变化。
    正弦与余弦函数的LUTs：

    包括这些LUTs的声明为常量，强迫编译器将这些数据存储在码区而不是数据区。由于微处理器中的RAM的有限性这样做是很重要的。由于LUTs的值必须以Q8．7方式排列，因此与正弦和余弦对应的值应该乘以27。
3．3 位倒置
    位倒置(N是已知的)可以在运行中计算，利用1个查寻表格标记或者直接用一个打开的环来写。每1种方法有其各自的大小与执行速度的平衡。本文利用开环直接写的方法来执行位倒置。实际的固件由原码来自动产生这个开环。 [p]
    利用开环来获得位倒置(其中N=256)：

3．4 基为2的FFT算法
    在对样本执行位倒置后，即可以计算FFT了。利用蝶形方法(见图1)计算基为2的FFT的固件包括3个主要的循环。在循环之外包括log2N的FFT计算阶段。在每一阶段循环内部执行单独的蝶形运算。FFT算法的核心是执行每一蝶形运算的块码。不幸的是，这一块码的计算是不轻便的。MUL_1与MUL_2利用微处理器的硬件乘法来执行乘法。
3．5 复数转化为极坐标表示
    为了计算输入信号的幅值，必须将复数X(k)转换为极坐标来表示。幅值将代替在固件中不再需要的FFT中的原来的值。
    式(4)中决定了二维的LUT其幅值而不是其计算。第1个值是频谱实部4个最重要的位，而第2个值是频谱虚部4个最重要的位。为了获得这些最重要的位，16位有符号数右移11次。频谱的实部与虚部都可以被用作下标时，它们被其绝对值代替了(因此符号位将是零)。
    从式(6)中，可以知道频谱是以X(N／2)对称的，只有前N／2+1个幅值被转化为极坐标表示。同样，对于输入样本为实数的频谱虚部中的X(0)与X(N／2)通常为零。因此这两个幅值通常分别是单独计算的。固件中用来自动计算X(k)的包括原码。
3．6 海明窗或汉宁窗
    为了实现此任务的固件包括LUTs(Q8．7形式)将海明窗或汉宁窗应用到输入样本中。加窗对于防止泄漏是有用的。加窗可以在时域上对输入样本截短。海明窗的方程如方式(8)所示，而汉宁窗则如式(9)所示。
   
    同样，这些实际固件的注释包括自动为这些窗函数产生LUTs的原码。

4 测试结果
    为了对FFT结果测试，利用微处理器中的通用异步收发报机端口。固件将幅值X(k)上传到PC上。这个程序包括FFT图表，即利用一个窗口应用程序从PC中的串行端口中读取这些幅值，实时地将计算得到的幅值利用图表画出来。利用微处理器对4个不同输入电压信号进行200Kb／s采样，并在图2中输出其FFT图像。

5 结语
    当然你可以利用无限的时间来优化及计算FFT。但是本文采取了基为2的FFT，这种方法可以很大限度地减少计算FFT所需要的加法与乘法。由于篇幅问题，有很多提高执行FFT速度的优化方法并没有在本文中给出。例如，对于实数的输入样本信号，输入样本的虚部通常为零，而且只有开始的一半频谱是重要的。利用这个信息，第一和最后一个FFT阶段就可以更快速的执行(但可能将需要编写更多的代码)。对于低功率的微处理器，在本文中所提到的关于FFT的算法是一个好的开始。