二次电缆寻线仪及保护序列化测试仪FFT算法优化研究

史春旻　沈心怡　张玮　俞家融　王天序

摘 要：针对快速傅里叶变换（Fast Fourier Transform，FFT）算法在单片机中运行耗时较长的问题，开展了耗时原因分析。针对现有耗时较长的三角函数计算、2的N次方计算，优化为查表法计算；平方和的根号计算，优化为运行速度更快的简化计算。该优化方法简单实用，可运用于其他需要计算优化的场景。最后对该优化方案进行了比较，研究结果表明，该优化方法显著提升了FFT的计算效率，具有很强的工程实用性。

关键词：二次电缆寻线仪；保护序列化测试仪；FFT算法优化

中图分类号：TP33    文献标志码：A    文章编号：1671-0797（2023）13-0059-04

DOI：10.19514/j.cnki.cn32-1628/tm.2023.13.015

0    引言

FFT（Fast Fourier Transform）算法是离散傅里叶变换的快速算法，由伟大的数学家高斯在计算小行星轨道时首次提出，高斯去世后，后人在其手稿中整理并发表在《高斯全集》第三卷，但当时并没有FFT算法的使用场景，FFT算法就这样被遗忘在历史的长河中。20世纪60年代，由于检测地下核试验振动波等应用需求，美国投入了大量资源来研究傅里叶变换的计算方法。新的FFT算法由美国数学家和计算机科学家John W. Tukey于1965年提出[1]，它是一种用于计算复杂信号的离散傅里叶变换的快速方法，主要优点是计算速度快，而且可以使用计算机高效地实现。FFT算法是一种数字信号处理（DSP）技术，用于将一个信号从时域转换为频域，它的主要作用是从信号中提取出频率特性，从而用于信号的处理、检测、分析和改善等。

在声波二次电缆寻线仪及保护序列化测试仪设备研发过程中，需要在单片机中运行FFT算法程序，为了优化和提升FFT算法程序在单片机中的计算速度，本文将研究国内电力系统继电保护和测控装置中普遍使用的FFT算法，并提出改进方法，在计算精度和计算速度之间取得平衡，为今后涉及信号处理的相关产品研发和算法优化提供一定的优化思路。

1    分裂基FFT算法介绍

本文采用了国内继电保护装置厂家在中低压保护装置中普遍使用的分裂基FFT算法。

分裂基算法（Split-Radix FFT Algorithm，SRA）由Duhamel和Hollman于1984年提出，最初是针对规模为2的幂的DFT问题的计算优化。Split-Radix FFT將N点DFT分解为一个N/2点DFT和两个N/4点DFT。分裂基的基本思路是对偶序列输出使用基2算法，对奇序列输出使用基4算法，将基2和基4分解组合在一起。在FFT计算中用较大的基可以进一步减少复数乘法，N/4点基4算法比N/2点基2算法更有效。Split-Radix FFT是众多FFT算法中乘法和加法次数最少的算法之一，其运算结构好，运算程序短，在国内继电保护装置中被广泛采用，研究表明其乘法复杂度接近于理论最小值[1]。尽管该库的计算速度很快，但该库在RP2040 Arm-Cortex M0单片机中的执行速度依然有一定的延时，为了能在二次电缆寻线仪及保护序列化测试仪项目的信号处理中有更好的表现，本文对分裂基FFT算法进行了一定程度的优化。

2    FFT算法介绍

2.1    傅里叶变换原理

傅里叶变换，将被测试信号分解为不同强度和频率的正弦波，把所有的正弦波全部叠加起来则是原始信号。原始信号可能看起来完全不像正弦波，这就是傅里叶变换的神奇之处。

如果想知道原始信号中某个频率的信号的强度，只需要将这个信号乘以这个正弦波，再计算曲线和坐标轴的面积之和。

比如某个信号f（t），包含以下分量：

f（t）=a1×sin t+a2×sin（2t）+a3×sin（3t）

如果要求sin t在f（t）的分量，需要在等式两边同时乘以sin x。

2.3    快速傅里叶变换计算耗时的原因分析

检查FFT运算的源码，经过统计可知，在计算8点FFT时，经过了一次2的N次方运算、24次循环计算，每次循环包含4次sin和cos运算；最后求信号幅值和相位时，需要完成4次平方和再开根号运算、4次arctan函数运算。因为乘法、除法运算量几乎不可压缩，本文将重点关注FFT运算过程中数学函数调用计算量的优化。

3    FFT算法优化方案

3.1    使用查表法实现sin和cos函数优化

为了加快单片机三角函数的计算速度，需要避免使用math.h头文件中的数学函数。本文选择使用查表法计算sin和cos，在8位单片机中使用8位的sin和cos精度，在32位单片机中使用32位精度。

图1以8位精度为例，构造一个sin 0°到sin 90°的数组。为了降低单片机的内存使用，这里统一使用8位内存整型值数组保存sin值。实际项目中，根据单片机的内存大小，在内存充足的情况下，既可以使用16位整数数组，也可以使用32位浮点数组保存sin的值，使用浮点数组时，fast_sin函数不用再除以255，可减少一次除法运算。图1中的每个值，除以255就可以得到对应的sin值。

根据图1构造的数组，编写sin函数，所有超出0°～90°的值，根据sin函数的周期性和对称性，做相应的平移和符号变换，再从数组中获取对应的值。查表法sin函数实现方案如图2所示。

同样的方法可以构造查表法cos函数，经过测试，该方法相比传统的数学函数库，能将sin和cos的计算时间优化到原有时间的3%，极大地提升了正弦函数和余弦函数的计算效率。

3.2    使用查表法求解2的N次方

FFT求解过程中，需要求解2的N次方。根据需要求解信号的采样长度，提前构造2的N次方数组，能大幅度提升该计算的速度。图3为最大12阶2的N次方构造数组示例（实际项目中需根據情况确定数组的长度）。

经测试，使用查表法求解2的N次方，能有效节约FFT运算时间1～2 μs。

3.3    使用近似法快速计算平方和根

在FFT运算中，需要求解复数的模√Re2+Im2，即平方和的平方根求解，此运算是比较耗时的运算，使用近似求解平方根的方法，能有效提升计算速度。图4是快速平方和根号函数的C语言代码实现，在fastRSS（Root of Squared Sum）中，当两数的比例大于4倍的时候，直接返回最大数，两数比例小于4倍的时候，使用近似求解。近似法快速计算平方和根函数实现方案如图4所示。

该实现方案和精确解的误差是多少呢？表1为Im固定为10，Re=1至Re=10时，fastRSS函数结果和精确解的误差分析。

从表1的结果可以看出，使用fastRSS函数，最大误差没有超过2%，具有很高的精度，同时避免使用数学函数库中的根号运算，有效提升了FFT算法计算速度。

4    FFT优化后运算时间和精度分析

树莓派RP2040单片机价格低廉，资源丰富。本文通过在RP2040单片机上运行优化前后的FFT代码，进行时间和精度分析。通过构造两个电压和电流信号，并将采样值在50 Hz周期中进行32点采样，将采样保存在测试数组中，再进行FFT运算，分析FFT运行的时间。表2是离散后的32点采样值。

使用表2数据，分别代入优化前后的程序，分析两者数据的差异。表3为优化前和优化后RP2040单片机中FFT程序输出结果。优化前程序计算时间为3 639 μs，优化后的执行时间为2 124 μs，优化效果显著，节约单片机运算时间约42%。

从表3的结果可知，最大误差为50 Hz电流分量结果，优化前20.5 A，优化后20.32 A，误差为0.88%。优化后的FFT算法具有非常高的精度。

5    结论和建议

本文对FFT计算源码过程中的函数调用进行了研究，针对部分数学函数提出了优化方案，在一定程度上提升了单片机中FFT算法的运算效率，平均效率提升约42%。改进措施简单实用，对单片机中其他计算量较大的运用场景也有一定的借鉴作用。

单片机中FFT算法还有其他的改进方法，比如在ARM-cortex M4处理器中调用SIMD指令加速乘法及三角函数运算[4]。本文仅提出单一的函数简化计算方法，未来还能在现有的基础之上，进一步进行算法优化，提高单片机中FFT算法运行速率，降低硬件设计成本，提升芯片工作效率。

[参考文献]

[1] 陈暾.基于ARM_V8平台的FFT实现与优化研究[D].长沙：湖南师范大学，2018.

[2] 王琳，胡耀，王世元.从采样的角度谈信号与系统中的傅里叶变换[J].安徽师范大学学报（自然科学版），2022，45（4）：332-337.

[3] 操庐宁.基于X86-64计算平台的FFT实现与并行优化[D].长春：吉林大学，2019.

[4] 姚建宇，张祎维，张广婷，等.基于SIMD的三角函数高性能实现与优化[J].计算机科学，2021，48（12）：29-35.

收稿日期：2023-03-14

作者简介：史春旻（1981—），男，江苏东台人，高级工程师，研究方向：变电站二次系统运维。