贝斯手蒲羽整容:可用于嵌入式计算的定点FFT算法

时域抽取法的定点FFT程序3.0 ！

简介

本示例程序演示了如何计算定点和浮点各个长度和位数的按时间抽取FFT算法。所有程序在VC6下编译通过。

最精确的FFT计算最好采用浮点，但浮点对于嵌入式设备计算量太大，有的朋友想采用定点，但网上资料不多，且长度位数要求不一，于是本人制作了这个版本。

此程序可以实现对长度为64，128，256，512，1024，2048的数据进行定点FFT（如果需要其他长度可自行修改）以及浮点版本。

根据处理器处理能力不同，对于数据类型的长度，提供了8位版本，16位版本和32位版本，可用于单片机等不同位数的嵌入式处理器。当然，三个版本主要是给定了不同的处理精度。

此程序采用了运算级间动态缩放舍去最低位，保证了输出不溢出，且提供了尽可能高的精度。具有自动适应输入值大小的能力。因此有较高的信噪比。

使用说明

FFT计算内容整合在timefft.h中，可将其加入到图形环境中或者控制台中，正如示例所示。

有两个值需要根据环境修改：

长度通过#define N 长度值

决定，长度值可取上文所述六种，主要是考虑嵌入式设备处理能力，一般来讲8位机用256点，这样可以保证下标在8bits以内，16位机可以考虑稍大的点数。

类型数据位数通过：

#define BIT 位数

决定，只能取8，16，32

1.对于图形界面演示示例，信号源和绘图程序在WaveAnalysisView.cpp中的void CWaveAnalysisView::OnDraw(CDC* pDC)函数中。

测试数据采用两个正弦波的叠加，即sin(2*3.14159*i/8)+sin(2*3.14159*i/64);

并将其扩大到最大值之间。

绘图的版本频域图像模值没有开根号，表示的是模值的平方。

如果时域直流量过大，则会引起频率的0点值过大，图像变得缩放的太小，如果为了绘图方便，可将其归零。

2.控制台下示例版本是将timefft.h文件独立出来，和一个控制台结合形成的。 其演示的是一个门函数引发的sinc函数。

定点算法原理和精度分析

由于定点运算位数为定值，因此不能像浮点一样进行各种乘法和加法处理，因为这两种运算会产生溢出。同时，对于旋转因子的小数模式也需要修整为可以运算的整数模式，这便是FFT定点的两个需要处理的难点。

本算法采用了运算量较大但效果较好的动态溢出检查。也就是说在计算每一级之前，对所有值进行遍历，如果发现有超过溢出允许输入值的变量，则对整个数组进行右移缩小处理。这样就保证了数值较小的输入不会被过度缩放，而数值很大的数据也能够通过FFT运算。由于对所有的数组值进行缩放，因此不会造成失真，但是幅度大小会有变化。因此可以理解为定点算法只求解相对值。在使用过程中，尽量不要叠加过多直流分量，因为这将导致频率0点处值很大，其它的值被缩放殆尽。影响精度。

下面详细介绍缩放方法：

一个乘加运算可以理解为：

X1+X2*WN

我们需要找出在什么情况下输出会最大，并在设定输出在此数据类型允许值下时，输入最大取多少。

由于X1，X2，WN都为复数，且abs(X2)=abs(X2*WN)，因此最坏的情况是，在模一定的情况下，原来X2的实部和虚部的数值经过旋转因子，被全部转移到了实部，而虚部为零。或者反之。

而由于当X2的实部和虚部相等时，它的模值最大，假设全部转移到实部，在加上X1的实部的值就构成了最大的输出。

我们另这个最大的输出等于类型最大允许值（8bit时为0x7f，16bit时为0x7fff），由于我们需要对所有输入做统一的限定值以便程序判断，所以只能有X1.real=X2.real

Sqrt(X2.real^2+X2.imag^2)+X1.real= Sqrt(2*X2.real^2)+X1.real=(sqrt(2)+1)X1.real

我们将此输出值限定在类型数据允许的最大整数上，则有

(sqrt(2)+1)X1.real

其中(sqrt(2)+1)为常数=2.4142

得到X1.real最大不能超过MAX/2.4142，虚部一样。我们定义此值为OVS

蝶形运算每算完一级，就要进行溢出判断，如果发现其中任何一个值大于此值，则需要对所有值缩放处理。

这时，我们面临怎样缩放的问题，如果采用溢出时所有值/2的方法，那么，它的输入端允许的最大信号幅值就是几乎每级计算完后都需要除以2才通过所有计算。此时

输入/2*2.4142=输出

如果只计算一级，那么输入不得超过OVS/1.207=52(以8bit为例)，

而每计算一级都要乘以1.207，成为底数大于1等比数列，此时如果计算1024点的FFT，则要保证输出不溢出，则：

输入*1.207^log(2,N)=输入*1.207^10

以八位为例max=127（0x7f），则输入不得超过19，可见，信噪比太低，不能应用。

继续观察运算我们发现，最多可能在同一级上溢出两个bit，我们动态的调整缩放的量级则会改善上述情况。也就是在溢出两位时除以4，溢出一位时除以2，经过计算机模拟，发现此时只要满足单极不溢出的条件（如8bit时为52）则总体就不会溢出，但是，付出的代价是两个bit的精度被处理掉了。因此整体的精度为7bit-2bit（8bit-符号位）=5bit，当然误差可能由于多次这样的操作而累积。

程序中，实现时在每次进入新的一级运算之前，将所有值遍历，如果超出2倍的ovs值则除以4，超出一倍除以2，否则不处理。这样，乘加运算在数据类型内部不会造成溢出了。

旋转因子的缩放：由于在计算旋转因子时，不能用小数，因此，将小数扩大到与位数相当，之后与x[l]相乘，结果再以相同比例缩小，

比如32位版本中：首先将其扩大

wn[k].real=(int)(cos(2*3.14159*k/N)*0x7fffffff);

在计算关于旋转因子的乘积项时乘完的结果再缩小。

XkWn.real=(((__int64)x[l].real*wn[wi].real)>>31)-(((__int64)x[l].imag*wn[wi].imag)>>31);

XkWn.imag=(((__int64)x[l].imag*wn[wi].real)>>31)+(((__int64)x[l].real*wn[wi].imag)>>31);

(__int64)是由于VC不支持32位乘法时自动切换到64为导致了数据丢失，因此做强制转换，其它类型不需要。

下载地址：UploadFiles/2009-2/231626252129.rar