利用FFT成功实现拓扑识别（七）--压缩FFT的运行时间

FFT在嵌入式应用的困境

• ​​嵌入式系统具有体积小、可靠性高功能强、灵活方便等许多优点，其应用已深入到工业、农业、教育、国防、科研以及日常生活等各个领域，注定了其必然要处理一些现场的紧急、并发的任务，更广泛一点，也就是检测到一个信息到系统处理完返回给用户的这个反应时间必须满足任务的处理时间要求。这个时间可以划分为几个阶段：
  1.中断响应时间；
  2.RTOS进程调度响应时间；
  3.应用程序响应时间；
  简单说就是程序能及时的解决、处理一些比较急的事情，而不会出现“卡机”等情况。总之，往往需要在下面几项中寻求平衡

  （1）实时性要求高
  （2）系统成本限制下资源有限
  （3）处理器性能瓶颈（例如没有浮点指令等等）

实时系统（Real-time operating system,RTOS）的正确性不仅依赖系统计算的逻辑结果，还依赖于产生这个结果的时间。而FFT的运算，尤其是对频率分辨率要求较高的应用场合，采样点数就需要达到一定的数量级才能满足，这就引起FFT运算时间成级数上升，影响到系统的实时性，以本人做过的某个系统为例，实测FFT用时如下：
大多数嵌入式系统中，tick周期多为1ms左右，也就是说，一次FFT计算过程中，系统往往要发生多次的任务调度，进一步拖长了FFT时间，而对于依赖FFT结果进行决策的任务，实时性就大打折扣。怎么解决这个问题昵？
下面结合本人的实际经验，从三个方面来探讨：

数据的准备

• FFT计算的前提是采集够预定点数的原始数据，而这些数据往往是有AD采集得到的，不同的处理器（或扩展ADC）采集的数据，要适配FFT算法的数据格式要求，如何更高效、准确的准备这些数据并和FFT时序匹配，也能一定程度提高整体FFT运行效率，我们可以从以下几个方面着手改进：

（1）DMA中断和FFT的巧妙衔接。

• 以此项目为例：

  在初始化中配置好DMA目标区地址和长度，DMA控制器自动传输AD数据到FFT原始数据区，达到采样点长度后，DMA中断发生，处理函数如下：

void FFTDataDmaIrqProcess(void)
{
//printf(“–>>>>>>>>DMA数据块对应FFT数字滤波OriginalSignal原始数据区>>>>>>>>\n\n\n”);
if((READ_REG(DMA_STATUSTFR) & 0x1)) /* DMA complete Interrupt ch0 /
{
/clear Interrupt/
WRITE_REG(DMA_CLEARTFR, 0x1);
Cnt++;
sAddrT = sAddrT + sLength2*4 ; //DMA结束后指向下一个区
if(DdSpiDataReadySem)
{
//发送FFT数据准备好信号量，启动FFT数据滤波
SystemSemaphoreGiveFromIsr(FFTFilterDataReadySem);
}
}
}

（2）数据滤波和FFT的并发处理

前文介绍过数据滤波的实现和对FFT结果的影响，通过设立独立的滤波器的OriginalSignal缓冲区和SignalFiltered缓冲区，实现了FFT运算和数字滤波的独立并发运行：

void IIR6_Buterworth(float OriginalSignal,floatSignalFilter,int num)
{
float a[7] = {1, -3.986026373811, 8.102670852719, -9.80794698246, 7.589104117838, -3.496678118687, 0.8216600080371};
float b[7] = { 0.0001074738548191, 0,-0.0003224215644572, 0, 0.0003224215644572, 0,-0.0001074738548191};

// 6阶滤波
SignalFilter[0]=b[0]*OriginalSignal[0];
SignalFilter[1]=b[0]*OriginalSignal[1] + b[1]*OriginalSignal[0] - a[1]*SignalFilter[0];
SignalFilter[2]=b[0]*OriginalSignal[2] + b[1]*OriginalSignal[1] + b[2]*OriginalSignal[0]- a[1]*SignalFilter[1]- a[2]*SignalFilter[0];
SignalFilter[3]=b[0]*OriginalSignal[3] + b[1]*OriginalSignal[2] + b[2]*OriginalSignal[1] + b[3]*OriginalSignal[0]- a[1]*SignalFilter[2]- a[2]*SignalFilter[2] - a[3]*SignalFilter[0] ;
SignalFilter[4]=b[0]*OriginalSignal[4] + b[1]*OriginalSignal[3] + b[2]*OriginalSignal[2] + b[3]*OriginalSignal[1]+ b[4]*OriginalSignal[0]
               - a[1]*SignalFilter[3]- a[2]*SignalFilter[2] - a[3]*SignalFilter[1] - a[4]*SignalFilter[0];
SignalFilter[5]=b[0]*OriginalSignal[5] + b[1]*OriginalSignal[4] + b[2]*OriginalSignal[3] + b[3]*OriginalSignal[2]+ b[4]*OriginalSignal[1]+ b[5]*OriginalSignal[1]
               - a[1]*SignalFilter[4]- a[2]*SignalFilter[3] - a[3]*SignalFilter[2] - a[4]*SignalFilter[1]- a[5]*SignalFilter[0];

for (int i=6;i<num;i++)
{
    SignalFilter[i]=b[0]*OriginalSignal[i] + b[1]*OriginalSignal[i-1] + b[2]*OriginalSignal[i-2]+ b[3]*OriginalSignal[i-3] + b[4]*OriginalSignal[i-4]+ b[5]*OriginalSignal[i-5]+ b[6]*OriginalSignal[i-6]
                   - a[1]*SignalFilter[i-1]- a[2]*SignalFilter[i-2]- a[3]*SignalFilter[i-3]- a[4]*SignalFilter[i-4]- a[5]*SignalFilter[i-5]- a[6]*SignalFilter[i-6];
}
//发送FFT数据准备好信号量，启动FFT
			SystemSemaphoreGiveFromIsr(FFTDataReadySem);

}

算法的改进

• 在FFT算法中，不同的数据结构、库函数的使用上，用时差异非常大，可以优化的措施有：

（1）对sin、cos函数的优化。

         在嵌入式设备中，浮点数的计算相对来说是十分缓慢的，sin与cos函数的计算更是缓慢。而在FFT计算过程中，以前文所述算法中，对sin和cos函数使用中：

             theta = -PI * i / n;
    		   a = cos(theta);
	          b = sin(theta);

可以看出：
（1）计算的点数是有限的，和采样点数相关
（2）数据循环使用
这就给算法提速提供了空间，对于sin与cos相关的大量浮点计算，我们可以在系统初始化的时候进行一次性的计算，建立对应的数据表格：

float SinTab[FFT_DOt/2],CosTab[FFT_DOt/2];

 用查表法代替实时的计算，用时更少，

（2）对乘法计算的优化。

我们知道，计算机计算移位运算，右移一位相当于除2，右移n位

• x >> n

相当于除以2的n次方。同理，左移一位相当于原数据乘以2，左移n位：

• x<<n;

相当于乘以以2的n次方。FFT算法中的乘2操作好多，可以通过此方法大大压缩计算时间。

任务的优化

嵌入式系统中，不仅仅有后台的AD，FFT任务，还有人机接口、通信、容错等等不同需求的任务需要调度。任务定义大同小异，大致如下：
SYS_TASK_DEF(NAME, TASKID, Proc, QUEUE_LENGTH, STACK_SIZE, PRIORITY, AUTOSTART)
各个参数对应：
NAME 任务名称
TASKID,任务ID号
Proc,任务的的处理进程
QUEUE_LENGTH, 任务的堆序列长度
STACK_SIZE,任务的栈大小
PRIORITY 任务的任务优先级
AUTOSTART 任务的启动方式

通过对各个task的优先级设定，可以影响系统调度的时间片长短安排，从而实现对紧急任务的及时处理，相关的文章很多，不在赘述。

采用上述方法优化过后，可以将FFT的运算时间缩短一、二个重量级，大大提高了系统的实时性性能。