基于FFT的频谱分析方法可以从含有噪声的信号中提取有用的信息，在仪器仪表的数据处理中具有重要的应用价值。为了保证频谱分析的精度，往往进行多点的FFT运算，例如，1024点、2048点等，这样运算量大、所占内存也大，只有采用DSP(数字信号处理器)才能实现实时的处理。目前，在工业现场普遍使用的两线制、低功耗自动化仪表，由于仪表本身消耗电流必须控制在4 mA之内，所以无法采用DSP等运算能力强的芯片，只能采用低功耗单片机;而低功耗单片机的运算速度和内存容量都很有限，所以，至今未见用其进行多点数FFT运算的报道。为了能够用低功耗单片机实时做FFT运算，以提高自动化仪表信息处理的能力，我们用汇编语言编制FFT程序，在程序中用定点数运算(以下简称定点FFT)，采取措施防止数据溢出，保证计算精度，合理分配内存。测试结果表明，我们编制的程序在MSP430F、1611单片机上，完成一次2048点的基于FFT的频谱分析和校正只需要500 ms，精度也达到要求，可以用于以低功耗单片机为核心的仪表中，实时完成信号处理任务。
1 定点运算
1.1 数据表示
在MSP430中使用C语言实现FFT运算，其乘法和加法运算都是默认使用浮点实现的。于MSP430属于定点单片机，因此浮点运算必须由大量的定点指令模拟，这将耗费大量的时间。因此我们针对MSP430的特点，使用汇编语言编制FFT程序，在程序中用定点数运算，并将数据统一使用16位定点数表示。16位定点数中最高位(左边的第1位)作为符号位，剩下的15位用于存放数值。数据格式如图1所示。
1.2 数据定标
定点单片机参与数值运算的数都是16位的整型数，但是运算过程中的数不一定都是整数。那么，定点计算过程中如何处理小数呢?这其中的关键就是由程序员来确定一个数的小数点处于16位中的哪一位。这就是数的定标。
通过设定小数点在16位数中的不同位置，就可以表示不同大小和不同精度的小数了。数的定标有Q表示法和S表示法两种。表1列出了一个16位数的16种Q表示、S表示及它们所能表示的十进制数值范围。
从表1中可以看出，同样一个16位数，若小数点设定的位置不同，它所表示的数也就不同。例如，十六进制数2000H=8192，用Q0表示;十六进制数2000H=O.25，用Q15表示;但对于定点运算来说，处理方法是完全相同的。下面简要介绍如何使用定点数乘法运算模拟浮点数乘法。
设浮点乘法运算的表达式为：float x，y，z;z=xy。假设经过统计(这里“统计”的意思是所有计算中数据范围都在定标范围内)后x的定标值为Qx，y的定标值为Qy，乘积z的定标值为Qz，则z=xy;zq×2-Qz=xq×yq×2-(Qx+Qy);zq=(xqyq)2Qz-(Qx+Qy)。所以，定点表示的乘法为：
1.3 FFT计算过程中的数据定标
为了在以MSP43F1611为处理器的仪表系统上进行基于FFT的功率谱估计，必须先由MSP430F1611的ADC进行采样，而ADC采样得到的数据需要经过定标后才能进行定点计算。定标过程为：ADC的采样电压范围为0～2.5 V，因此，采样过程实际上就是将信号电压除以2.5进行归一化，使得采样得到的数据范围为O～1 V，此时数据就可用Q15表示，即将ADC的12位采样结果寄存器中的数据右移4位保存起来，维持12位精度，转换为Q15定点数表示。由于FFT运算过程中，蝶形输出相对蝶形输入数据被放大了3倍，因此蝶形输出数据范围为一3～+3 V。此时数据如果仍然使用Q15表示，就会发生溢出，故改用Q13表示数据，即将12位ADC数据右移1位。实际上经过处理后，ADC数据精度没有变化，但使用Q13表示数据比用Q12表示数据，其蝶形输出的数据精度高。这是由于定点计算时需要对蝶形输出右移以防止溢出，而使用Q13表示数据比使用Q12表示数据少右移了1位，因此多了1位有效数据。FFT运算过程中使用Q13表示数据，就使得加法乘法运算都可以直接使用定点指令实现，减少了很多判断处理，提高了运算速度。使用Q13表示数据，即最高位(左边的第1位)是符号位，剩下的15位表示数据的大小。表示数据大小的15位中，高2位(左边的第2位和第3位)用来表示数据中的整数部分，在计算中作为保护位;最低的13位(右边的13位)表示数据中的小数部分，如果经过某次蝶形单元运算后，最大值正好被放大3倍，此时数据就由13位扩大到15位，保证数据不会增大到16位，冲走符号位，发生溢出。运算完成后将FFT计算过程中的这一级所有结果都右移2位，就能够使得这一级计算结果的最大值仍然可用13位表示，同时也可将这一级所有蝶形运算输出的数据同时缩小，保证下级计算。表示数值大小的15位数据的数据格式如图2所示。
1.4 旋转因子数据定标
FFT运算过程使用定点计算，且使用有符号乘法，必须始终保留1位作为符号位;而旋转因子范围为-1～1，因此可定标为Q14，转换为16位定点数。其转换过程为：根据参与FFT运算的数据点数计算出旋转因子的正余弦表，然后将正余弦表乘以16384，即左移14位，最后四舍五人取整。如果使用Q15表示数据，即需要左移15位，那么正余弦表中最大值1，经过上述处理后成为-1，发生溢出。
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，当实部和虚部的模的平方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位过程。
5 算法测试
为了验证算法的实时性和正确性，通过信号发生器产生标准信号送入所研制的基于MSP430f1611为核心的处理系统，对算法进行了全面的测试。
(1)测试算法运行时间
测试对2 048点数据进行功率谱估计所需要的总时间。预先设置MSP430F1611单片机的P5.4引脚为普通I/O，且为输出方式，接着，循环执行FFT运算和功率谱估计程序，且每次开始FFT计算和功率谱估算前，将P5.4输出电平翻转，因此P5.4输出的相邻两次翻转电平的时间间隔就是FFT运算和功率谱估计的总时间。通过数字存储示波器观测P5.4引脚输出的信号波形，如图4所示。
每次高低电平翻转的时间间隔约为500 ms，即对2 048点数据进行功率谱估计总共需要500 ms。
(2)测试算法计算精度
由于FFT运算的最大误差发生在非整周期采样时，所以，选择这些最大误差点来进行测试。由于泄漏误差，信号基频表示为
f0=(k+d)fs/N (7)
式(7)中，k为整数，d为小数(定义d为泄漏误差系数)。由于泄漏误差不超过频率分辨率的二分之一，所以|d|≤O.5。当d=O时，即为整周期采样情况;当d=O.5时，就是最大非整周期采样的地方。因为所研制的基于FFT的频谱分析方法将应用于数字涡街流量计，在此，针对气体40口径频率范围为69～1 380 Hz，设定采样频率为3 717.472 199 Hz，数据点数为2 048，选择不同的k值得到不同的频率信号，由信号发生器产生幅值为60 mV的这些信号，送入两线制涡街流量计信号处理系统低通滤波器前端，然后经过预放大电路和低通滤波电路后，送入MSP430F1611进行频率估计。
京公网安备 11010802033920号
