【Delphi】如何在三軸加速器的頻譜分析中使用FFT(快速傅里葉變換)算法

       關於傅里葉變換的做用，網上說的太過學術化，且都在說原理，以及如何編碼實現，可能不少人有個模糊印象，在人工智能，圖像識別，運動分析，機器學習等中，頻譜分析成爲了必備的手段，可將離散信號量轉換爲數字信息進行歸類分析。

今天這裏講的不是如何實現，而是如何使用傅里葉變換

       但頻譜分析中，涉及到的信號處理知識對大部分軟件開發的人來講，太過於晦澀難懂，傅里葉變換，拉普拉斯，卷積，模相，實數，虛數，複數，三角函數等等，已經能讓軟件工程師望而卻步，形成懂知識的人沒法開發，懂開發的人沒法分析，而同時具有兩種技能的人，除了專業的研究生博士，剩下的少之又少。

       站在軟件開發人員的角度，可否拋開這些晦澀難懂概念，進行純軟件方面的信號處理？這個問題可能沒有答案，最好可以經過實踐來證實。並且如果拋開這些概念，讓那些有興趣的開發人員學習到信號處理，頻譜分析，先無論可不可能，光想一想，其做用也是有的，至少「高深」的東西可以被用上了。

三軸加速器步數計算分析

       下面就以三軸加速器（X，Y，Z三個加速器）的實際應用開發爲例，來說一講（限於我的水平，可能比較粗略）。

       時下運動APP也有很多，很多APP開發人員就會接觸到陀螺儀，加速器等，拋開計步的準確性不講，咱們這裏來看看這些計步分析的原理，如何利用這些傳感器進行行爲分析，首先須要獲取傳感器的數據（如何獲取就不在這裏講了，有Android，iOS，獨立藍牙/WIFI設備等），將獲取的數據輸出到平面座標圖表中：

      

       如上圖，即爲時域座標曲線圖，X,Y,Z三種顏色分別表明三軸加速器的三個方向加速度隨時間的變化曲線。座標的橫軸爲時間軸，縱軸爲加速度值（該值由加速度傳感器的電壓值轉換而來，具體由硬件開發者決定）。

       咱們的採集頻率爲50Hz，即橫軸（時間軸）每個點爲50分之1秒，能夠看到前4秒（0-200）沒有太明顯走動，4秒後到12秒走動速度較慢，12秒走路稍微加快。在這裏不得不感嘆人腦簡直是宇宙超級無敵，光看圖表就可以直接知道行爲，沒法想象將來誰可以讓人工智能超越人腦。

        在軟件開發中，咱們要如何經過代碼分析獲得上面「人腦的分析結果「？ 要直接上代碼不是個人做風，百度能夠搜到一篇騰訊的文章《利用三軸加速器的計步測算方法》，此文章的分析方法能夠說大體可以理解，但對於讀者可能幫助不是特別大，由於涉及到的計算方法卻語焉不詳，好比，在沒有進行濾波前，上圖的X軸波峯波谷並不平滑，能夠看出每一個大波峯上又有2個小波谷，若是以此來計算，得出來的步數偏差很是大，如何解決偏差？

        可能一些作過信號處理的人會說，作一個「均值滾動濾波+低通濾波+高通濾波」就差很少了，我只能說對於複雜的「週期性」信號分析是可行的，對於三軸加速器的步數計算來講，實際中並不存在什麼週期性信號，人的上坡下坡走路非直線，時快時慢等動做都會加大分析難度，每每理論提及來容易，實現起來難，會發現實際狀況和理論不太同樣。

        先不說這麼多了，拿到信號原始數據，立刻進行傅里葉變換，Delphi可使用FPC一個開源的庫AlgLib，最新版好像不開源了且分爲免費版和收費版，但codetyphon收錄了之前的開源版本能夠用。

        alglib中關於快速傅里葉變換是在u_fft.pp文件中，delphi只把後綴pp改爲pas就可使用了。使用複數FFTC1D和實數FFTR1D函數進行變換，這裏僅作單軸分析，因此使用實數FFTR1D函數就行了，複數傅里葉變化應用在三維波動曲線的分析中，這裏不用就不關心了（針對運動曲線來說，至於無線電信號的複合器或者多維矩陣數據分析的狀況就多了），2個函數中的N即爲採集的數據量，固然實際中咱們有時候採集了幾十秒，而N的最佳範圍是在128-2048範圍（既充分又計算量不高且作代碼運算方便），因此咱們每10秒（50Hz即500個數據，最好是直接擴展到512個。2的n次方）進行一次傅里葉變換分析便可。

題外話：

        專業人士看到這裏可能會特別說明FFT僅僅是離散傅里葉變換(DFT)，因爲計算機是由時間片定時執行一次的，因此採集的數據也是每隔一段時間採集一次，固然間隔時間越短越貼近連續採集，總的來說，就是非連續的，離散的數據，計算機也只能作離散傅里葉變換分析了，電子電路設備光學設備等才能作連續傅里葉變換。

         至於爲什麼使用實數傅里葉變換，由於咱們採集的數據是一個一個採集的（三軸其實是分開獨立來分析的），因此採集的數據只是一維（時間維度不計在內），使用實數傅里葉變換再適合不過了。

FFTR1D函數的調用方法

        好比採集到的數據爲[3.4,  4.1, 6.7, 3.1...],則以下代碼（僅參考用）：

       

var
  A : TComplex1DArray; 
  R : TReal1DArray;
  I, N : Integer;
begin
   N := 50*10 // 10秒（建議使用512，2的N次方）
   SetLength(R, N);
   for I:=0 to N-1 do
   begin     // 時域的複合幅度值
      R[I]:=Data[nIdx];//按順序填上採集的數據[3.4,  4.1, 6.7, 3.1...]
   end;
   FFTR1D(R, N, A);  // 使用實數傅里葉變換方法（速度理論上比複數快一倍）
   ////////////////////////// 
   // 到這裏已經完成傅里葉變換，如何使用變換後的數據？？
end;

 

複數數組A即爲傅里葉變換結果，可是這個結果有啥用？這裏先引用百科上的一段話來幫助理解後面如何利用傅里葉變換結果的代碼：

假設FFT以後某點n用複數a+bi表示，
那麼這個複數的模就是An=根號a*a+b*b，
相位就是Pn=atan2(b,a)。
根據以上的結果，就能夠計算出n點（n≠1，且n<=N/2）對應的信號的表達式爲：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。
對於n=1點的信號，是直流份量，幅度即爲A1/N。

因爲FFT結果的對稱性，一般咱們只使用前半部分的結果，即小於採樣頻率一半的結果。

       通過變換後的複數數組A，便是變換結果，但並非分析的數據，而是一箇中間數據，由該數據能夠獲得模值和相位，根據傅里葉變換原理，模值計算以下：

// 對A數組每一個點的複數求模，即實部和虛部平方和再開根號。 
nAValue := Math.Power(Abs(A[I].X*A[I].X)+Abs(A[I].Y*A[I].Y), 0.5); 

　　至於相位則是根據公式Pn=atan2(b,a)計算，可是頻譜分析比較有用的是模值（主要是濾波，這裏區別於物理硬件上的濾波-採用共振/折射等方法），因此其餘無論。

       前面講到使用10秒的數據進行計算，這樣獲得的傅里葉變換結果複數也有500個，因爲傅里葉的對稱性，咱們只取0-250的結果來計算模值便可，由於251-500是對稱的結果。

       計算模值後再輸出平面座標圖表（頻域圖），以下：

  

固然幅度也有做用，幅度的計算說明以下（來自百度百科）：

假設原始信號中某個頻率的峯值（幅度）爲V，那麼FFT的結果的每一個點（除了第一個點直流份量以外）的模值就是V的N/2倍。而第一個點就是直流份量，它的模值就是直流份量的N倍。

計算幅度 = 模值*2/N（第一個點= 模值/N），後再輸出平面座標圖表（頻域圖），以下：

 

如上圖，網絡上各類傅里葉變換文章中的標誌性的結果圖就出來了（請忽略圖中的時間刻度與文章的10秒500個數據等相關性，由於我用的圖是隨便一段數據來的。）

看到上面的傅里葉變換結果圖，會不會以爲悲劇纔開始發生。沒錯，結果是出來了，但咱們該如何分析呢？

如何分析傅里葉變換結果

    首先傅里葉變換後的結果要看得懂，直白的講，變換後的結果反應了各個頻率的模值，模值越高表明該頻率產生的「做用」越大，也就表明越切合實際存在該頻率的信號在起做用，而圖中的縱座標就是模值（或幅度），橫座標就是頻率，但頻率是整數的，實際頻率與前面的採集頻率和數據量有關，關係以下：

       若採樣頻率是50Hz，採集50個數聽說明咱們可以分辨出1Hz的頻率，能夠簡單認爲分辨率爲1，而500個點表明頻率的分辨率是10，也說明這個傅里葉變換結果可以分辨出0.1Hz的頻率，在這裏表明每個橫座標的整數點對應的頻率刻度爲0.1Hz。

       因此看到上圖15-20範圍內的X加速器（藍色線）的模值很高，表明1.5-2Hz的頻率最切合實際，說明在這一部分的分析結果人的走路速度是1秒鐘1.5-2步（走得比較慢），而咱們的業餘半馬運動員步頻在180步/分鐘左右，專業的190步/分鐘左右。

       繼續話題，實際中咱們開發的軟件是須要自動化處理，不可能跟人腦同樣高級，一看圖就知道結果，固然若是靠人腦，傅里葉變換都不須要就能夠分析了，因此別想那麼多了，在人工智能沒法代替人腦的時代，軟件的信號處理，到這裏纔剛剛開始。

       如何進行下一步分析呢，咱們能夠取一個範圍，如1-5Hz內的模值數據，其餘的濾波處理掉，同時根據採樣頻率的二分之一爲帶寬(奈奎斯特頻率)之間的關係，50Hz的採樣頻率必須把高於25Hz的頻段濾掉（幸虧人行走最高5Hz，非要跑到6Hz以上的瘋子就無論了吧），獲得簡化後的模值複數數組。

       另外，小幅度波動也能夠濾掉，至於幅度值多少能夠去掉，能夠自行統計判斷，可能不一樣加速度出來的值也不同，簡單的按照最大幅度的百分比算，低於10%的所有濾掉。不過爲了後面可能須要的積分計算準確性更高，保留着也是能夠的。

利用傅里葉反變換進行濾波

        （未完待續）

後續還有如何均值化，方波化等等，最終獲得很是貼近實際的步數結果