FFT（高速フーリエ変換）振動解析のSTM32マイコン用サンプルコード

1.概略

　STマイクロエレクトロニクスでは、STM32マイコン用にさまざまなサンプルコードを無償で提供していますが、FFT（高速フーリエ変換）振動解析に関しても複数のサンプルコードを提供しています。　現在リリースされているFFT振動解析のサンプルコードを下表にまとめました。いずれも、MotionSPというSTのミドルウェアを使用しています。
　この記事では、VibrationAnalysisという、下表の黄色の枠で囲ったサンプルコードの内容について、簡単に紹介したいと思います。異なる２つのSoftwareパッケージに含まれていますが、基本的には同じ内容です。このサンプルコードは、下表に示したサンプルコードの中でも、比較的シンプルでありながら、実装内容も充実しているため、この記事で取りあげることにしました。

＊ 対応評価ボード：STEVAL-BFA001V2Bについて
このボードは、広帯域・低ノイズの３軸振動センサ（加速度センサ）：IIS3DWB¹を搭載し、産業用機器間に使用される通信規格IO-Linkに対応した、産業機器向けのモニタリング用評価ボードです。少し用途が特化したボードですので、汎用的な評価ボードへの移植に関する記事も参考にしてみてください。

*表中の各ソフトウェアパッケージの詳細は、以下のリンク先を参照ください。
・X-CUBE-MEMS1: https://www.st.com/ja/embedded-software/x-cube-mems1.html
・FP-IND-PREDMNT1: https://www.st.com/ja/embedded-software/fp-ind-predmnt1.html
・STSW-PROTEUS: https://www.st.com/ja/embedded-software/stsw-proteus.html

MotionSPとは？
　振動解析を通じて各種機械・装置などの故障の予兆を検出する予知保全向けとして、STが提供している信号処理のソフトウェア・ライブラリです。具体的には、センサのデータからFFT演算を行い、スペクトルピーク値などの周波数領域の特徴量の算出や、速度RMS値や加速度のピーク値などの時間領域の特徴量の算出も行います。詳細は、下のリンク先を参照ください。
MotionSP Middleware Implementation（PDF英文資料）

2.MotionSPを使用したサンプルコード：VibrationAnalysis

2.1 動作概要

2.1.1 ボード上での動作

　まず最初に、VibrationAnalysisのサンプルコードをインストールしたボードが、どのような動作をするのかを簡単に紹介します。下図はこの記事に関連するハードウェアのみを、ブロック図にしたものです。

[接続] サンプルコード動作時のボードとPCとの接続は、STLINK²のUSB-シリアル変換の機能を使用したVirtual COM接続です。データの送受信はテキストデータなので、PC側は、TeraTermなどの端末ソフトが使用できます。

[動作] ボードをSTLINK経由でPCに接続した後、端末ソフト上で、次の動作内容を繰り返します。

① デフォルト値に設定済みの振動解析用パラメータがリスト表示され、ユーザからの入力待ちで停止。
② ここで、パラメータの変更がキーボード入力で可能。
③ パラメータ変更終了後、サンプルコードの振動解析が開始され、設定したデータ取得時間経過後、その結果が端末ソフトにテキスト表示。

2.1.2 FFT振動解析処理の大まかな流れ

ソースコードの中身に見ていく前に、FFT振動解析の大まかなステップに触れておきます。

1. 加速度センサから加速度値のデータをFFTサイズ分取得
2. 取得したデータに対して窓関数を乗算し、FFT演算
3. 設定したデータ取得時間の間、1.と2.を繰り返す
4. データ取得時間経過後、FFT演算結果の平均値を算出し、正規化とピーク値の算出を行う

2.2 サンプルコードの構成

　ここからは、サンプルコードの中身を見ていきます。サンプルコード全体の大まかな流れとしては、下図のように、周辺機能などの初期化コードの後、mainループに入り、加速度センサからの割り込み処理を伴いながら、振動解析の最上位層の関数（ MotionSP_Manager関数 ） が動作します。

ソースコードの構成は下表のようになっています。 MotionSP_MainManager( )などの振動解析の上位層の関数は、主にMotionSP_Manager.cの中に含まれ、それらの関数から、振動解析処理のコアとなる、ミドルウェア：MotionSPのAPIをコールする形で構成されています。

以降、振動解析処理の初期化とメインの処理内容について、順次説明していきます。

2.3 プラットフォームの初期化

　プラットフォームの初期化では、上のフローチャート図のように、UARTやLED用GPIOの初期化以外に、MotionSPに関する初期化（MotionSP_Init()）を行っています。その内容としては、大まかに次の２つに分けられます。

• 使用するセンサに関する基本的な設定（主にI/Fの確立）
• 振動解析用パラメータをヘッダファイルにあるデフォルト値に設定。加速度センサのパラメータとMotionSP用のパラメータの２つに対して以下のように設定。

  • AcceleroParamsデフォルト設定値：

  • MotionSP_Parametersデフォルト設定値：

2.4 振動解析処理

　振動解析処理の最上位の関数：MotionSP_MainManager() のフローチャート図を、以下に示しました。

＊MotionSPを使用した振動解析動作状態（動作中：１／停止中：０）を示すグローバル変数が、MotionSP_Runningです。

この動作フローの内容を見ていく上で、次の３つのステップ：
(1) 設定 　⇒　(2) 解析 　⇒　(3) 出力
に分けて順次説明していきます。

2.4.1 設定

ここで説明する対象の関数は次の２つとなります。
・MotionSP_VibrationAnalysisInit()
・MotionSP_VibrationAnalysisStart()

2.4.1.1 MotionSP振動解析初期化 :[ MotionSP_VibrationAnalysisInit() ]

①振動解析パラメータの変更： [ MotionSP_VibrationAnalysisSetParam() ]
　「プラットフォームの初期化」時に、デフォルト値に設定した振動解析パラメータ：
・AcceleroParams
・MotionSP_Parameters
のそれぞれに対して、TeraTerm等の端末ソフトからの指示に従って、設定値の変更を行います。

②振動解析パラメータの設定： [ MotionSP_VibrationAnalysisVariableInit() ]

振動解析パラメータ（MotionSP_Parameters）の設定値に基づいて、MotionSPで必要となる変数を初期化

• FFT解析結果格納変数（AccMagResults）の初期化
• 窓関数データ（Filter_Params）の生成 (MotionSP_SetWindFiltArray())
• 時間領域解析データ結果格納変数（TimeDomainData）の初期化
• MotionSPミドルウェア内で使用する変数の初期化
• ARM CMSIS FFT演算関数の初期化（arm_rfft_fast_init_f32()）

③加速度センサの設定とODRの計測 : [ MotionSP_AcceleroConfig() ]

• AcceleroParamsの振動解析パラメータに基づいて、加速度センサを設定（フルスケール値, ODRの設定など）
• FFT演算時には実際のデータの取得周波数（ODR）を把握することが必要なため、この関数の中で、加速度センサを取得するタイミングの周波数を、加速度センサからの割り込み（INT1）回数で計測しています。この時に使用する割り込みハンドラはMotionSP_DataReady_IRQ_Rtn()です。

[割り込み処理] 加速度センサには２つの割り込み用の出力端子(INT1, INT2)があります。両方ともこのサンプルコードでは使用しており、以下のように割り込み内容を割り当てています。（2.4.2 ①加速度センサ内のFIFOデータの読み出しの図を参照）

• INT1 : 加速度センサ内の出力データの準備ができたことを知らせるDataReady
• INT2 : 加速度センサ内にあるFIFO内のデータ量が設定した閾値まで達したことを知らせるFifoFull

2.4.1.2 センサのFIFOと循環バッファの設定:[ MotionSP_VibrationAnalysisStart() ]

この関数にて、割込みが有効になり、振動解析動作を開始

• 循環バッファ（Motion_Data.AccCircBuff）の設定
  • メモリ確保とゼロリセット
  • 循環バッファ構造体変数(MotionSP_Data.AccCircBuff)とFIFO割り込みのフラグの初期化
• 加速度センサのFIFOに関する設定
  • FIFOからの読み出し用バッファのメモリ確保
  • FIFOの容量制限をウォータマークのレベルに設定
  • センサの出力をFIFOモードに設定
  • INT2の割込みを、FIFO-FULLに割り当て、有効にする
  • BDR（Batch Data Rate:FIFOへの書き込みデータレート）を設定
  • データ取得時間計測用に、Systickタイマーのタイムスタンプの取得を開始

　すでに、似たような名前のパラメータや変数が複数登場しましたので、少し整理するために、振動解析の関数との関係をイメージ図にしてみました。

＊AcceleroParamsの変数：MotionSP_Manager.cの中で定義
＊それ以外の上図の変数：MotionSP.cの中で定義

2.4.2 解析

　VibrationAnalysisのメインとなる関数：MotionSP_MainManager() の説明に移っていきます。ここでは、加速度センサ内のFiFoからデータを読み出し、時間領域と周波数領域の解析処理を行っていますが、この記事では、周波数領域の解析処理に絞って説明します。 また、実際のサンプルコードでは、X軸,Y軸,Z軸それぞれについて振動解析処理を行っていますが、説明を簡単にするため、X軸のみを例に説明していきます。

周波数領域の解析処理の大まかな流れとしては、次の通りです。

①加速度センサ内のFIFOデータの読み出し(MotionSP_FifoFull_IRQ_Rtn())
②循環バッファのデータ量をチェック(MotionSP_FrequencyDomainAnalysis())

③FFT演算処理（MotionSP_fftExecution()）

では、それぞれに関して少し詳しく見ていきます。

---

①加速度センサ内のFIFOデータの読み出し： [ MotionSP_FifoFull_IRQ_Rtn() ]

　加速度センサは、Fifo内のデータ量が所定の閾値に達したことを知らせるため、割込み端子（INT2）を使って、マイコンに対して割り込みを要求します。
　その割り込みハンドラ(MotionSP_FifoFull_IRQ_Rtn())では、下図のように、読み出したセンサデータを加速度の単位に変換し、循環バッファに格納後、読み出し完了のフラグ（MotionSP_Data.FifoEmpty）をセットします。

②循環バッファのデータ量をチェック：[ MotionSP_FrequencyDomainAnalysis() ]

　上記の割り込みハンドラにてセンサからのデータ読み出しが完了すると、そのフラグを検出し、まずは、FFT演算に必要となるデータ量が、循環バッファに蓄積されているかを判断する必要があります。
　まず、FFT演算で必要となるデータ取得タイミングのイメージとしては、下図のようになります。循環バッファのデータを設定に基づいて、オーバーラップさせながら、順次データを取得していきます。

データの取り出しタイミングは、次の２つ循環バッファのインデックス値（変数） を比較して判断しています。

• accCircBuffIndexTmp ：必要な格納データ量を示す目標のインデックス値
• accCircBuffIndexForFftTmp：現在の格納データ量を示すインデックス値
  　
  具体的には、目標のインデックス値は、FFT演算に必要なデータ量を示し、
  （現在のインデックス値）>=（目標のインデックス値）であれば、FFT演算に必要なデータ量があり、データを取り出すタイミングになります。 このタイミングで、FFT演算処理関数内（MotionSP_fftExecution())に移行します。

循環バッファ関連の変数：
・循環バッファのオーバーフローを示すフラグ：MotionSP_Data.AccCircBuff.Ovf
・次の目標インデックスまでの差分：accCircBuffIndexForFftDelta
・循環バッファのインデックス値としての変数：次の４つあります。

　比較用のインデックスが用意されていますが、これは、目標値が循環バッファのサイズを超えた際にもインデックスの大小比較で判断できるようにするためです。

＜例：FFTサイズ：4096、オーバーラップ：75％＞
（１回目） 解析の最初は、FFTサイズ分のデータを取り出す必要があるため、
（循環バッファに格納したデータ量）>=（FFTサイズ） の条件で、FFTサイズ分のデータを循環バッファから取り出します。そのため、目標インデックス値（accCircBuffIndexForFftTmp）には、FFTサイズー１（例：4095）の値が格納されています。

（２回目以降） オーバーラップしない新規データ分を補充して、FFTサイズとするため、（新たに格納するデータ量）>= (オーバーラップしないデータ量） の条件で、FFTサイズ分のデータを循環バッファから取り出します。そのため、目標インデックス値（accCircBuffIndexForFftTmp）には、前回の目標インデックス値に対して、オーバーラップしないデータ量（=accCircBuffIndexForFftDelta=4096 x (1-0.75) =1024）が加算された値（4095 + 1024 = 5119）となります。

③FFT演算処理：[ MotionSP_fftExecution() ]

　FFTサイズ分データを取得する度に、FFT演算処理を行います。この（データ取得）＋（FFT演算処理）の繰り返しは、設定したデータ取得時間の間行われ、設定時間を経過した際は、演算結果の平均化等を行い、出力データ格納変数に格納します。ここでは、それらの処理を行うMotionSP_fftExecution()の中を見ていきます。

（１）データ取得時間内の処理：

　以下の一連の処理が、１回の加速度データのFFT演算処理になりますが、精度を上げるために設定したデータ取得時間の間、この一連の処理を繰り返し実行しています。

1.循環バッファからのデータ取得(MotionSP_fftInBuild())
2.窓関数を乗算(motionSP_fftUseWindow())
3.FFT演算(MotionSP_fftCalc())
4.FFT演算結果の平均化前処理(MotionSP_fftAverageCalcTime())

下図は、一連の処理の簡単なイメージ図です。

図中のプロット図に関しては、配列のインデックス値をX軸にしてグラフにしたイメージ図です。

　「１.循環バッファからのデータ取得」に関しては、memcpyを使ったデータのコピーとなりますが、それ以外は、ほぼCMSIS-DSPライブラリを使用した演算処理です。
　データ取得時間の計測には、Systickタイマーの割り込み回数をカウントして計測しています。MotionSP_FrequencyDomainAnalysis（）の中でそのカウント値のチェックをしています。

（２）データ取得時間経過後の処理：

ここからは、FFT演算出力データに対する処理となります。主な内容は以下の通りです。

5.FFT演算結果の平均化(MotionSP_fftAverageCalcTime())
6.FFT演算結果のスケール補正(MotionSP_fftAdapting())
7.FFT演算結果のピーク値算出(MotionSP_fftPeakFinding())

• 平均化：
  「（１）データ取得時間内の処理」で各インデックス毎に累積してあったFFT出力データの合計値をそのデータ数（FFT演算回数）で除算して平均値を求めます。
  AccMagResults.Array.X[i] = AccMagResults.Array.X[i] / AccMagResults.MagSize

• スケール補正：
  窓関数によるFFT出力振幅の減衰を補正するため、窓関数の補正値を乗算し、更にFFT出力振幅がFFTサイズに依存しように除算することで、FFT出力振幅の正規化をしています。
  AccMagResults.Array.X[i] = (AccMagResults.Array.X[i] x WSF) / AccMagResults.MagSize

• ピーク値算出：
  最後にFFT演算処理で得られた結果（FFTサイズ分の配列データ）から、arm_max_f32()関数を用いて、最大値を算出し、AccMagResults.Max.X.valueに格納。

2.4.3 出力

解析処理終了後、MotionSP_VibrationAnalysisSendResult() 関数を使用して、振動解析結果の内容を、UART経由でTeraTerm等の端末ソフトにテキスト出力します。出力内容は、FFTの振幅スペクトルグラフのプロットデータを出力します。

• 縦軸のデータ：FFTの振幅値として、上記のAccMagResults.Array.Xの配列データをそのまま出力
• 横軸のデータ：次式より、各インデックスに対応する周波数値を算出して出力
  インデックス値　x　周波数分解能（AccMagResults.BinFreqStep）

＊周波数分解能は、FFT演算結果で得られる周波数範囲（ センサデータ出力レート（ODR）の半分 ）をFFTサイズで除算した値であり、（１）設定のMotionSP_VibrationAnalysisInit（） ー　MotionSP_AcceleroConfig() 関数内で以下のように定義されています。

AccMagResults.BinFreqStep = (AcceleroODR.Frequency / 2) / AccMagResults.MagSize;

　
PC端末ソフト（TeraTerm）上に表示される内容は、下図のようになります。

VibrationAnalysisサンプルコードの内容に関する説明は、これで以上となります。

　今回のサンプルコードでは、振動解析結果をテキストデータとして表示するだけですが、他のサンプルコードでは、解析結果をPCやスマホに送信して、GUIやアプリ上に表示することも可能な内容となっています。＜関連記事＞に、その動作内容の一部を紹介しています。

＜関連記事＞：

• 「IoT向け小型センサ･ノード開発キット「STWIN」が予知保全に使える！？ 　スマホ編（基本）」

• 「STM32F4用の振動解析FFTサンプルコードをSTM32U5マイコン用にポーティング」

1. IIS3DWB製品情報：https://www.st.com/ja/mems-and-sensors/iis3dwb.html ↩

2. STLINKは、STM32マイコン用のデバッグツールですが、いくつか種類があります。一部を除き、殆どの製品で、USB-シリアル変換アダプタ（Virtual COM）の機能が搭載されています。詳細は以下のリンク先の資料を参照ください。https://www.st.com/resource/en/product_presentation/stlink-debugging-and-programming-tools-overview.pdf ↩