二次電纜尋線儀及保護(hù)序列化測(cè)試儀FFT算法優(yōu)化研究

史春旻　沈心怡　張瑋　俞家融　王天序

摘 要：針對(duì)快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）算法在單片機(jī)中運(yùn)行耗時(shí)較長(zhǎng)的問(wèn)題，開(kāi)展了耗時(shí)原因分析。針對(duì)現(xiàn)有耗時(shí)較長(zhǎng)的三角函數(shù)計(jì)算、2的N次方計(jì)算，優(yōu)化為查表法計(jì)算；平方和的根號(hào)計(jì)算，優(yōu)化為運(yùn)行速度更快的簡(jiǎn)化計(jì)算。該優(yōu)化方法簡(jiǎn)單實(shí)用，可運(yùn)用于其他需要計(jì)算優(yōu)化的場(chǎng)景。最后對(duì)該優(yōu)化方案進(jìn)行了比較，研究結(jié)果表明，該優(yōu)化方法顯著提升了FFT的計(jì)算效率，具有很強(qiáng)的工程實(shí)用性。

關(guān)鍵詞：二次電纜尋線儀；保護(hù)序列化測(cè)試儀；FFT算法優(yōu)化

中圖分類號(hào)：TP33? ? 文獻(xiàn)標(biāo)志碼：A? ? 文章編號(hào)：1671-0797（2023）13-0059-04

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

0? ? 引言

FFT（Fast Fourier Transform）算法是離散傅里葉變換的快速算法，由偉大的數(shù)學(xué)家高斯在計(jì)算小行星軌道時(shí)首次提出，高斯去世后，后人在其手稿中整理并發(fā)表在《高斯全集》第三卷，但當(dāng)時(shí)并沒(méi)有FFT算法的使用場(chǎng)景，F(xiàn)FT算法就這樣被遺忘在歷史的長(zhǎng)河中。20世紀(jì)60年代，由于檢測(cè)地下核試驗(yàn)振動(dòng)波等應(yīng)用需求，美國(guó)投入了大量資源來(lái)研究傅里葉變換的計(jì)算方法。新的FFT算法由美國(guó)數(shù)學(xué)家和計(jì)算機(jī)科學(xué)家John W. Tukey于1965年提出[1]，它是一種用于計(jì)算復(fù)雜信號(hào)的離散傅里葉變換的快速方法，主要優(yōu)點(diǎn)是計(jì)算速度快，而且可以使用計(jì)算機(jī)高效地實(shí)現(xiàn)。FFT算法是一種數(shù)字信號(hào)處理（DSP）技術(shù)，用于將一個(gè)信號(hào)從時(shí)域轉(zhuǎn)換為頻域，它的主要作用是從信號(hào)中提取出頻率特性，從而用于信號(hào)的處理、檢測(cè)、分析和改善等。

在聲波二次電纜尋線儀及保護(hù)序列化測(cè)試儀設(shè)備研發(fā)過(guò)程中，需要在單片機(jī)中運(yùn)行FFT算法程序，為了優(yōu)化和提升FFT算法程序在單片機(jī)中的計(jì)算速度，本文將研究國(guó)內(nèi)電力系統(tǒng)繼電保護(hù)和測(cè)控裝置中普遍使用的FFT算法，并提出改進(jìn)方法，在計(jì)算精度和計(jì)算速度之間取得平衡，為今后涉及信號(hào)處理的相關(guān)產(chǎn)品研發(fā)和算法優(yōu)化提供一定的優(yōu)化思路。

1? ? 分裂基FFT算法介紹

本文采用了國(guó)內(nèi)繼電保護(hù)裝置廠家在中低壓保護(hù)裝置中普遍使用的分裂基FFT算法。

分裂基算法（Split-Radix FFT Algorithm，SRA）由Duhamel和Hollman于1984年提出，最初是針對(duì)規(guī)模為2的冪的DFT問(wèn)題的計(jì)算優(yōu)化。Split-Radix FFT將N點(diǎn)DFT分解為一個(gè)N/2點(diǎn)DFT和兩個(gè)N/4點(diǎn)DFT。分裂基的基本思路是對(duì)偶序列輸出使用基2算法，對(duì)奇序列輸出使用基4算法，將基2和基4分解組合在一起。在FFT計(jì)算中用較大的基可以進(jìn)一步減少?gòu)?fù)數(shù)乘法，N/4點(diǎn)基4算法比N/2點(diǎn)基2算法更有效。Split-Radix FFT是眾多FFT算法中乘法和加法次數(shù)最少的算法之一，其運(yùn)算結(jié)構(gòu)好，運(yùn)算程序短，在國(guó)內(nèi)繼電保護(hù)裝置中被廣泛采用，研究表明其乘法復(fù)雜度接近于理論最小值[1]。盡管該庫(kù)的計(jì)算速度很快，但該庫(kù)在RP2040 Arm-Cortex M0單片機(jī)中的執(zhí)行速度依然有一定的延時(shí)，為了能在二次電纜尋線儀及保護(hù)序列化測(cè)試儀項(xiàng)目的信號(hào)處理中有更好的表現(xiàn)，本文對(duì)分裂基FFT算法進(jìn)行了一定程度的優(yōu)化。

2? ? FFT算法介紹

2.1? ? 傅里葉變換原理

傅里葉變換，將被測(cè)試信號(hào)分解為不同強(qiáng)度和頻率的正弦波，把所有的正弦波全部疊加起來(lái)則是原始信號(hào)。原始信號(hào)可能看起來(lái)完全不像正弦波，這就是傅里葉變換的神奇之處。

如果想知道原始信號(hào)中某個(gè)頻率的信號(hào)的強(qiáng)度，只需要將這個(gè)信號(hào)乘以這個(gè)正弦波，再計(jì)算曲線和坐標(biāo)軸的面積之和。

比如某個(gè)信號(hào)f（t），包含以下分量：

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

如果要求sin t在f（t）的分量，需要在等式兩邊同時(shí)乘以sin x。

2.3? ? 快速傅里葉變換計(jì)算耗時(shí)的原因分析

檢查FFT運(yùn)算的源碼，經(jīng)過(guò)統(tǒng)計(jì)可知，在計(jì)算8點(diǎn)FFT時(shí)，經(jīng)過(guò)了一次2的N次方運(yùn)算、24次循環(huán)計(jì)算，每次循環(huán)包含4次sin和cos運(yùn)算；最后求信號(hào)幅值和相位時(shí)，需要完成4次平方和再開(kāi)根號(hào)運(yùn)算、4次arctan函數(shù)運(yùn)算。因?yàn)槌朔ā⒊ㄟ\(yùn)算量幾乎不可壓縮，本文將重點(diǎn)關(guān)注FFT運(yùn)算過(guò)程中數(shù)學(xué)函數(shù)調(diào)用計(jì)算量的優(yōu)化。

3? ? FFT算法優(yōu)化方案

3.1? ? 使用查表法實(shí)現(xiàn)sin和cos函數(shù)優(yōu)化

為了加快單片機(jī)三角函數(shù)的計(jì)算速度，需要避免使用math.h頭文件中的數(shù)學(xué)函數(shù)。本文選擇使用查表法計(jì)算sin和cos，在8位單片機(jī)中使用8位的sin和cos精度，在32位單片機(jī)中使用32位精度。

圖1以8位精度為例，構(gòu)造一個(gè)sin 0°到sin 90°的數(shù)組。為了降低單片機(jī)的內(nèi)存使用，這里統(tǒng)一使用8位內(nèi)存整型值數(shù)組保存sin值。實(shí)際項(xiàng)目中，根據(jù)單片機(jī)的內(nèi)存大小，在內(nèi)存充足的情況下，既可以使用16位整數(shù)數(shù)組，也可以使用32位浮點(diǎn)數(shù)組保存sin的值，使用浮點(diǎn)數(shù)組時(shí)，fast_sin函數(shù)不用再除以255，可減少一次除法運(yùn)算。圖1中的每個(gè)值，除以255就可以得到對(duì)應(yīng)的sin值。

根據(jù)圖1構(gòu)造的數(shù)組，編寫sin函數(shù)，所有超出0°～90°的值，根據(jù)sin函數(shù)的周期性和對(duì)稱性，做相應(yīng)的平移和符號(hào)變換，再?gòu)臄?shù)組中獲取對(duì)應(yīng)的值。查表法sin函數(shù)實(shí)現(xiàn)方案如圖2所示。

同樣的方法可以構(gòu)造查表法cos函數(shù)，經(jīng)過(guò)測(cè)試，該方法相比傳統(tǒng)的數(shù)學(xué)函數(shù)庫(kù)，能將sin和cos的計(jì)算時(shí)間優(yōu)化到原有時(shí)間的3%，極大地提升了正弦函數(shù)和余弦函數(shù)的計(jì)算效率。

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

FFT求解過(guò)程中，需要求解2的N次方。根據(jù)需要求解信號(hào)的采樣長(zhǎng)度，提前構(gòu)造2的N次方數(shù)組，能大幅度提升該計(jì)算的速度。圖3為最大12階2的N次方構(gòu)造數(shù)組示例（實(shí)際項(xiàng)目中需根據(jù)情況確定數(shù)組的長(zhǎng)度）。

經(jīng)測(cè)試，使用查表法求解2的N次方，能有效節(jié)約FFT運(yùn)算時(shí)間1～2 μs。

3.3? ? 使用近似法快速計(jì)算平方和根

在FFT運(yùn)算中，需要求解復(fù)數(shù)的模√Re2+Im2，即平方和的平方根求解，此運(yùn)算是比較耗時(shí)的運(yùn)算，使用近似求解平方根的方法，能有效提升計(jì)算速度。圖4是快速平方和根號(hào)函數(shù)的C語(yǔ)言代碼實(shí)現(xiàn)，在fastRSS（Root of Squared Sum）中，當(dāng)兩數(shù)的比例大于4倍的時(shí)候，直接返回最大數(shù)，兩數(shù)比例小于4倍的時(shí)候，使用近似求解。近似法快速計(jì)算平方和根函數(shù)實(shí)現(xiàn)方案如圖4所示。

該實(shí)現(xiàn)方案和精確解的誤差是多少呢？表1為Im固定為10，Re=1至Re=10時(shí)，fastRSS函數(shù)結(jié)果和精確解的誤差分析。

從表1的結(jié)果可以看出，使用fastRSS函數(shù)，最大誤差沒(méi)有超過(guò)2%，具有很高的精度，同時(shí)避免使用數(shù)學(xué)函數(shù)庫(kù)中的根號(hào)運(yùn)算，有效提升了FFT算法計(jì)算速度。

4? ? FFT優(yōu)化后運(yùn)算時(shí)間和精度分析

樹莓派RP2040單片機(jī)價(jià)格低廉，資源豐富。本文通過(guò)在RP2040單片機(jī)上運(yùn)行優(yōu)化前后的FFT代碼，進(jìn)行時(shí)間和精度分析。通過(guò)構(gòu)造兩個(gè)電壓和電流信號(hào)，并將采樣值在50 Hz周期中進(jìn)行32點(diǎn)采樣，將采樣保存在測(cè)試數(shù)組中，再進(jìn)行FFT運(yùn)算，分析FFT運(yùn)行的時(shí)間。表2是離散后的32點(diǎn)采樣值。

使用表2數(shù)據(jù)，分別代入優(yōu)化前后的程序，分析兩者數(shù)據(jù)的差異。表3為優(yōu)化前和優(yōu)化后RP2040單片機(jī)中FFT程序輸出結(jié)果。優(yōu)化前程序計(jì)算時(shí)間為3 639 μs，優(yōu)化后的執(zhí)行時(shí)間為2 124 μs，優(yōu)化效果顯著，節(jié)約單片機(jī)運(yùn)算時(shí)間約42%。

從表3的結(jié)果可知，最大誤差為50 Hz電流分量結(jié)果，優(yōu)化前20.5 A，優(yōu)化后20.32 A，誤差為0.88%。優(yōu)化后的FFT算法具有非常高的精度。

5? ? 結(jié)論和建議

本文對(duì)FFT計(jì)算源碼過(guò)程中的函數(shù)調(diào)用進(jìn)行了研究，針對(duì)部分?jǐn)?shù)學(xué)函數(shù)提出了優(yōu)化方案，在一定程度上提升了單片機(jī)中FFT算法的運(yùn)算效率，平均效率提升約42%。改進(jìn)措施簡(jiǎn)單實(shí)用，對(duì)單片機(jī)中其他計(jì)算量較大的運(yùn)用場(chǎng)景也有一定的借鑒作用。

單片機(jī)中FFT算法還有其他的改進(jìn)方法，比如在ARM-cortex M4處理器中調(diào)用SIMD指令加速乘法及三角函數(shù)運(yùn)算[4]。本文僅提出單一的函數(shù)簡(jiǎn)化計(jì)算方法，未來(lái)還能在現(xiàn)有的基礎(chǔ)之上，進(jìn)一步進(jìn)行算法優(yōu)化，提高單片機(jī)中FFT算法運(yùn)行速率，降低硬件設(shè)計(jì)成本，提升芯片工作效率。

[參考文獻(xiàn)]

[1] 陳暾.基于ARM_V8平臺(tái)的FFT實(shí)現(xiàn)與優(yōu)化研究[D].長(zhǎng)沙：湖南師范大學(xué)，2018.

[2] 王琳，胡耀，王世元.從采樣的角度談信號(hào)與系統(tǒng)中的傅里葉變換[J].安徽師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2022，45（4）：332-337.

[3] 操?gòu)]寧.基于X86-64計(jì)算平臺(tái)的FFT實(shí)現(xiàn)與并行優(yōu)化[D].長(zhǎng)春：吉林大學(xué)，2019.

[4] 姚建宇，張祎維，張廣婷，等.基于SIMD的三角函數(shù)高性能實(shí)現(xiàn)與優(yōu)化[J].計(jì)算機(jī)科學(xué)，2021，48（12）：29-35.

收稿日期：2023-03-14

作者簡(jiǎn)介：史春旻（1981—），男，江蘇東臺(tái)人，高級(jí)工程師，研究方向：變電站二次系統(tǒng)運(yùn)維。