【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%的全部濾掉。不過為了后面可能需要的積分計算准確性更高,保留着也是可以的。

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

        (未完待續)

后續還有如何均值化,方波化等等,最終得到非常貼近實際的步數結果


注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2021 ITdaan.com