在MSP430F1611上实现周期图谱分析及校正

2 防止溢出，保证精度

FFT中的蝶形运算如图3所示。设蝶形输入为：X1(k)，实部为X1(k)r，虚部为X1(k)i；X2(k)，实部为X2(k)r，虚部为 X2(k)i。设蝶形输出为：X(k)，实部为X(k)，，虚部为X(k)i；X(k+N／2)，实部为X(k+N／2)r，虚部为X(k+N／2)i。则有：





由式(1)和式(2)可以看出，蝶形输出的实部和虚部是由3个数相加得到的，因此数据可能会放大3倍。如果计算过程中的数据始终使用定点数表示，随着级数的增加，就会发生溢出。例如，使用16位定点数表示，其最高位(从左数第1位)为符号位，其表示的数据范围为-32 768～32 767。如果采样得到的数据最大值为4 096，经过两级计算后蝶形最大输出就可能为4096×3×3=36 864，超出了16位定点数的表示范围。

下面介绍保证数据计算精度的方法。

为了提高计算速度，系统中使用定点数法运算FFT，且使用Q13表示数据。蝶形运算中，其蝶形输出的数据的实部和虚部都使用3次加法运算，即每级蝶形运算都可能使数据扩大3倍，因此，蝶形输出的实部和虚部结果都需要右移2位(缩小4倍)以防止溢出。但随着计算级数的增加，移位将会使数据变得越来越小。例如，128点FFT，总共需要7级运算，数据最终将移位2×7=14位(缩小47=16 384倍)，因此当信号幅值不够大时，经过多级运算可能会无法分辨出主信号频率。

设FFT运算结果的主信号频率点的对应实部为r，虚部为i，其幅值为A(ADC的量化值)，参与运算的数据点数为N，由FFT功率谱计算的性质可得：



设经过定点FFT运算，也就是运算过程中有移位，则该主信号频率点的模为K，即：



联立式(3)和式(4)，得



由于功率谱估计是找出功率谱中的最大值，确定主信号的频率，根据经验，使用定点数运算FFT，当实部和虚部的模的平方K2为2时，就无法由功率谱分辨出主信号频率。由式(5)可得：



因此，当K2为2，N为128时，A=128×1．414=180．992=181，即当信号的幅值为18l／4 096×2．538=112 mV，就分辨不出主信号频率。考虑K2为2的极限情况，当A为724，N为512时，即给定信号幅值为724／4 096×2．538=449 mV时，就分辨不出主信号频率。

为了防止计算结果经过多次移位后，数据太小无法分辨主信号，系统针对定点FFT运算采取如下处理：由于FFT定点运算中，一般情况下，为了处理方便，每级蝶形运算中乘法结果都限制在-1～1范围内，即乘法运算的结果始终为小数(只有经过加法运算，数据才有可能超出-1～1范围)，因此，通过判断蝶形输出的结果，决定是否移位。当发现超出-1～1范围，就将本级的所有蝶形运算的输出结果右移2位，没有超出就不进行移位。

3 内存分配

由式(3)可知，功率谱估算时需要另外开辟一段内存空间存储功率谱结果。例如，当进行2048点基于FFT的功率谱分析时，需用1024个浮点数存放功率谱计算结果，这将占有很大一段内存。但实际运算中，每个频率点功率，只与其FFT运算结果中的对应频率点的实部、虚部有关，而与其他频率点无关。因此功率谱运算中，可采取以下步骤将存放实部的空间存放功率谱：

①实部、虚部数据平方计算。由于MSP430F1611内部集成了硬件乘法器，因此可将乘法器的第一操作数寄存器(OP1)、第二操作数寄存器(OP2) 写入相同的数据实现平方运算。

②平方结果移位。平方结果需要右移13位，使用Q13表示，同时使用16位的临时变量将平方结果保存。

③功率谱计算结果保存。实部平方结果、虚部平方结果相加后再存人原来的实部单元。

经过上述步骤后，就可将原来存放实部、虚部数据的内存单元再次利用。

定点FFT运算过程中，还可将用来存放采集数据的内存空间，再次用作存放FFT运算过程中的实部数据，另外再开辟同等大小的内存空间，存放虚部数据。例如，对于RAM空间为10 KB的MSP430F16ll来说，使用16位定点数运算FFT，最多能够运算2 048点。因为实部、虚部结果都需4 096 KB，故共需8．192 KB，正好小于10KB；而运算4 096点FFT时，共需16．384 KB，超出10 KB。

4 程序实现

算法实现时使用如下方法简化了程序运算过程：

①C程序调用汇编FFT程序，同时为了处理方便将功率谱运算过程也用C语言实现。为了使汇编程序中使用的内存空间与C程序中的内存空间地址不发生冲突，汇编程序中所需的变量都在C文件中定义。

②由于实部、虚部都使用C语音数组来存储，当计算点数很多时，数组将很大。例如，当运算2 048点FFT时，就需定义两个长度为2 048的整形数组，这两段数组不能用堆栈局部空间存储，只能用全局数组，由于C语言规定全局变量默认初始化为0，MSP430的IAR编译环境，进入 main函数之前的cstart函数中就用cstar_inh_zero函数对全局变量进行初始化，由于定义的数组太长，初始化需要很长时间，导致程序还没有进入main函数，看门狗就已经复位。因此定义全局数组时，加上_no_init关键字。例如，定义一个数据长度为2 048的不需要初始化的整型数组，使用语句no_init int fft[2048]。

③为节约RAM内存空间，将旋转因子对应的正余弦表制作成表格存放在ROM空间中，而蝶形运算时可通过查表得到旋转因子后再进行FFT运算。

④由于FFT运算过程中的旋转因子是通过左移14位取整得到的，因此蝶形运算过程中需要乘法运算结果右移14位。MSP430F1611单片机乘法器的结果寄存器，由高16位乘法结果寄存器RESFII、低16位乘法结果寄存器RESLO组成，右移14位操作就是保留高位结果寄存器所有内容和低位结果寄存器中的高两位，因此RESHI、RESLO一起向左移2位，然后保留高位结果寄存器作为乘法结果就可实现右移14位过程。