不韋 呂
FFT アナライザ このプログラムの説明は，CQ出版社「トランジスタ技術」の2021年10月号から開始された連載記事「STM32マイコンではじめるPC計測」の中にあります．このプログラムといっしょに使うPC側のプログラムについても同誌を参照してください．
Dependencies: Array_Matrix mbed SerialTxRxIntr DSP_ADDA UIT_FFT_Real Window
Revision 0:e5fc70976c00, committed 2021-09-09
- Comitter:
- MikamiUitOpen
- Date:
- Thu Sep 09 08:52:33 2021 +0000
- Child:
- 1:d9dbfbe95c8d
- Commit message:
- 1
Changed in this revision
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Array_Matrix.lib Thu Sep 09 08:52:33 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/Array_Matrix/#d3aa1ddb57e1
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DSP_ADDA.lib Thu Sep 09 08:52:33 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/DSP_ADDA/#a1dcee67c67e
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/DoubleBuffer.hpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,58 @@
+//--------------------------------------------------------
+// ダブル･バッファの template クラス
+// バッファに２次元配列（Matrix クラス）を使用
+//
+// 2021/05/22, Copyright (c) 2021 MIKAMI, Naoki
+//--------------------------------------------------------
+
+#ifndef DOUBLE_BUFFER_HPP
+#define DOUBLE_BUFFER_HPP
+
+#include "Matrix.hpp"
+using namespace Mikami;
+
+template<class T> class DoubleBuffer
+{
+public:
+ // コンストラクタ
+ explicit DoubleBuffer(int size, T initialValue = 0)
+ : N_(size), buf_(2, size, initialValue), ping_(0), pong_(1),
+ index_(0), full_(false) {}
+
+ // データを格納
+ void Store(T data) { buf_[ping_][index_++] = data; }
+
+ // 出力バッファからデータの取り出し
+ T Get(int n) const { return buf_[pong_][n]; }
+
+ // バッファが満杯でバッファを切り替える
+ void IsFullSwitch()
+ {
+ if (index_ < N_) return;
+
+ ping_ ^= 0x1; // バッファ切換えのため
+ pong_ ^= 0x1; // バッファ切換えのため
+ index_ = 0;
+ full_ = true;
+ }
+
+ // バッファが満杯で，true を返す
+ bool IsFull()
+ {
+ bool temp = full_;
+ if (full_) full_ = false;
+ return temp;
+ }
+
+private:
+ const int N_; // バッファのサイズ
+ Matrix<T> buf_; // バッファ
+ int ping_, pong_; // バッファ切替用
+ int index_; // 入力データのカウンタ
+ bool full_; // 満杯の場合 true
+
+ // コピー･コンストラクタおよび代入演算子の禁止のため
+ DoubleBuffer(const DoubleBuffer&);
+ DoubleBuffer& operator=(const DoubleBuffer&);
+};
+#endif // DOUBLE_BUFFER_HPP
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/IIR_Filter/Biquad.hpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,62 @@
+//--------------------------------------------------------------
+// 縦続形 IIR フィルタの構成要素として使う 2 次の IIR フィルタ
+// b0 は 1 と仮定している
+//
+// 2020/11/04, Copyright (c) 2020 MIKAMI, Naoki
+//--------------------------------------------------------------
+
+#include "mbed.h"
+
+#ifndef IIR_BIQUAD_HPP
+#define IIR_BIQUAD_HPP
+
+class Biquad
+{
+public:
+ // フィルタの係数をまとめて扱うための構造体
+ struct Coefs { float a1, a2, b1, b2; };
+
+ // デフォルト･コンストラクタ
+ // 係数は構造体 Ceofs で与える
+ Biquad(const Coefs ck = (Coefs){0, 0, 0, 0})
+ : a1_(ck.a1), a2_(ck.a2), b1_(ck.b1), b2_(ck.b2),
+ un1_(0), un2_(0) {}
+
+ // 係数を個別に与えるコンストラクタ
+ Biquad(float a1, float a2, float b1, float b2)
+ : a1_(a1), a2_(a2), b1_(b1), b2_(b2), un1_(0), un2_(0) {}
+
+ virtual ~Biquad() {}
+
+ // 2 次のフィルタを実行する
+ float Execute(float xn)
+ {
+ float un = xn + a1_*un1_ + a2_*un2_;
+ float yn = un + b1_*un1_ + b2_*un2_;
+
+ un2_ = un1_;
+ un1_ = un;
+
+ return yn;
+ }
+
+ // 係数を設定する
+ void SetCoefs(const Coefs ck)
+ {
+ a1_ = ck.a1;
+ a2_ = ck.a2;
+ b1_ = ck.b1;
+ b2_ = ck.b2;
+ }
+
+ // 内部変数（遅延器）のクリア
+ void Clear() { un1_ = un2_ = 0; }
+
+private:
+ float a1_, a2_, b1_, b2_; // フィルタの係数
+ float un1_, un2_; // 遅延器
+
+ // コピー･コンストラクタ禁止
+ Biquad(const Biquad&);
+};
+#endif // IIR_BIQUAD_HPP
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/IIR_Filter/Coefs_IIR_LP.hpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,47 @@
+//-----------------------------------------------------
+// 縦続形 IIR フィルタの次数と係数の定義
+//
+// 2021/06/01, Copyright (c) 2021 MIKAMI, Naoki
+//-----------------------------------------------------
+
+#include "Biquad.hpp"
+
+// FFT アナライザで使うフィルタ
+// 標本化周波数が 102.4 kHz の場合，5.24 kHz 以上で
+// 少なくとも 60 dB 減衰させる LPF
+
+// 低域通過フィルタ
+// 連立チェビシェフ特性
+// 次数 ： 10 次
+// 標本化周波数：102.4000 kHz
+// 遮断周波数 ： 5.0000 kHz
+// 通過域のリップル： 0.50 dB
+// 阻止域の減衰量 ：60.00 dB
+
+const int ORDER1_ = 10;
+const Biquad::Coefs CK1_[] = {
+ { 1.820963E+00f, -8.355312E-01f, -9.111580E-01f, 1.0f}, // 1段目
+ { 1.853721E+00f, -9.012503E-01f, -1.783341E+00f, 1.0f}, // 2段目
+ { 1.881627E+00f, -9.570190E-01f, -1.871513E+00f, 1.0f}, // 3段目
+ { 1.895648E+00f, -9.843304E-01f, -1.892515E+00f, 1.0f}, // 4段目
+ { 1.902771E+00f, -9.961831E-01f, -1.898312E+00f, 1.0f}}; // 5段目
+const float G01_ = 1.232249E-03f; // 利得定数
+
+// 白色雑音発生器で使うフィルタ
+// 標本化周波数が 100 kHz の場合，
+// 100/(4π) ≒8.15 kHz 以上で
+// 少なくとも 60 dB 減衰させる LPF
+
+// 低域通過フィルタ
+// 連立チェビシェフ特性
+// 次数 ： 6 次
+// 標本化周波数：102.4000 kHz
+// 遮断周波数 ： 5.0000 kHz
+// 通過域のリップル： 0.50 dB
+// 阻止域の減衰量 ：60.00 dB
+const int ORDER2_ = 6;
+const Biquad::Coefs CK2_[] = {
+ { 1.794327E+00f, -8.134116E-01f, -6.466692E-01f, 1.0f}, // 1段目
+ { 1.830789E+00f, -8.923580E-01f, -1.704269E+00f, 1.0f}, // 2段目
+ { 1.874859E+00f, -9.684694E-01f, -1.812154E+00f, 1.0f}}; // 3段目
+const float G02_ = 1.381274E-03f; // 利得定数
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/IIR_Filter/IirCascade.hpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,57 @@
+//---------------------------------------------------
+// 縦続形 IIR フィルタ
+//
+// 2020/11/04, Copyright (c) 2020 MIKAMI, Naoki
+//---------------------------------------------------
+
+#include "Biquad.hpp"
+#include "Array.hpp" // Array クラスが定義されている
+using namespace Mikami;
+
+#ifndef IIR_CASCADE_HPP
+#define IIR_CASCADE_HPP
+
+class IirCascade
+{
+public:
+ // コンストラクタ
+ IirCascade(int order, const Biquad::Coefs ck[], float g0)
+ : order_(order), hn_((order+1)/2)
+ { SetCoefs(order, ck, g0); }
+
+ // コンストラクタ
+ IirCascade(int order, const Biquad hk[], float g0)
+ : order_(order), hn_((order+1)/2, hk), g0_(g0) {}
+
+ virtual ~IirCascade() {}
+
+ // フィルタ処理を実行する
+ float Execute(float xn)
+ {
+ float yn = g0_*xn;
+ for (int k=0; k<(order_+1)/2; k++) yn = hn_[k].Execute(yn);
+ return yn;
+ }
+
+ // 係数の設定
+ void SetCoefs(int order, const Biquad::Coefs ck[], float g0)
+ {
+ if (order_ != order)
+ {
+ order_ = order;
+ hn_.SetSize((order+1)/2);
+ }
+ g0_ = g0;
+ for (int k=0; k<(order+1)/2; k++) hn_[k].SetCoefs(ck[k]);
+ }
+
+ // 内部変数（遅延器）のクリア
+ void Clear()
+ { for (int k=0; k<(order_+1)/2; k++) hn_[k].Clear(); }
+
+private:
+ int order_; // 次数
+ Array<Biquad> hn_; // Biquad クラスのオブジェクトの配列
+ float g0_; // 利得定数
+};
+#endif // IIR_CASCADE_HPP
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/MSeq16.hpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,42 @@
+//---------------------------------------------------------
+// M 系列信号発生器（N = 16）
+//
+// 2020/10/17, Copyright (c) 2020 MIKAMI, Naoki
+//---------------------------------------------------------
+
+#include "mbed.h"
+
+#ifndef MSEQ16_HPP
+#define MSEQ16_HPP
+
+namespace Mikami
+{
+ class MSeq16
+ {
+ public:
+ MSeq16() : reg_(1) {}
+
+ // 戻り値： 1 => 1, 0 => -1
+ int Execute()
+ {
+ if ((reg_ & B_M_) == B_M_)
+ {
+ reg_ = ((reg_ ^ XOR_) << 1) | 1; // 1 の場合の処理
+ return 1;
+ }
+ else
+ {
+ reg_ = reg_ << 1; // 0 の場合の処理
+ return -1;
+ }
+ }
+ private:
+ static const uint16_t XOR_ = (1 << (2-1))
+ | (1 << (3-1))
+ | (1 << (5-1)); // XOR の位置に対応する定数
+ static const uint16_t B_M_ = 1 << (16-1); // 16 段目に相当するビットを調べる
+
+ uint16_t reg_;
+ };
+}
+#endif // MSEQ16_HPP
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/MyFFT_Analyzer/FFT_Analyzer.cpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,30 @@
+//-------------------------------------------------------
+// FFT を使ってスペクトル解析を行うクラス
+//
+// 2021/05/22, Copyright (c) 2021 MIKAMI, Naoki
+//-------------------------------------------------------
+
+#include "FFT_Analyzer.hpp"
+
+namespace Mikami
+{
+ FftAnalyzer::FftAnalyzer(int nFft)
+ : N_FFT_(nFft), fft_(nFft), wBlk_(nFft),
+ xData_(nFft), wData_(nFft), yFft_(nFft/2+1) {}
+
+ void FftAnalyzer::Execute(const Array<float> &xn, Array<float> &absFt)
+ {
+ xData_ = xn; // データを作業領域にコピー
+
+ // 直流分を除去
+ float sum = 0;
+ for (int n=0; n<N_FFT_; n++) sum = sum + xData_[n];
+ float ave = sum/N_FFT_;
+ for (int n=0; n<N_FFT_; n++) xData_[n] = xData_[n] - ave;
+
+ wData_ = wBlk_.Execute(xData_); // 窓掛け
+ fft_.Execute(wData_, yFft_); // FFT の実行
+ for (int n=0; n<=N_FFT_/2; n++) // 絶対値に変換
+ absFt[n] = 100.0f*abs(yFft_[n]);
+ }
+}
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/MyFFT_Analyzer/FFT_Analyzer.hpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,39 @@
+//-------------------------------------------------------
+// FFT を使ってスペクトル解析を行うクラス（ヘッダ）
+//
+// 2021/05/22, Copyright (c) 2021 MIKAMI, Naoki
+//-------------------------------------------------------
+
+#ifndef FFT_ANALYZER_HPP
+#define FFT_ANALYZER_HPP
+
+#include "Array.hpp"
+#include "fftReal.hpp"
+#include "Blackman.hpp"
+
+namespace Mikami
+{
+ class FftAnalyzer
+ {
+ public:
+ // nFft: FFT のデータ点の数
+ FftAnalyzer(int nFft);
+ virtual ~FftAnalyzer() {}
+ void Execute(const Array<float> &xn, Array<float> &absFt);
+
+ private:
+ const int N_FFT_;
+
+ FftReal fft_;
+ BlackmanWindow wBlk_;
+
+ Array<float> xData_; // 解析対象の時系列データ
+ Array<float> wData_; // 窓掛けされたデータ
+ Array<Complex> yFft_; // FFT の結果
+
+ // コピー･コンストラクタおよび代入演算子の禁止のため
+ FftAnalyzer(const FftAnalyzer& );
+ FftAnalyzer& operator=(const FftAnalyzer& );
+ };
+}
+#endif // FFT_ANALYZER_HPP
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MyFFT_Analyzer/UIT_FFT_Real.lib Thu Sep 09 08:52:33 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/UIT_FFT_Real/#0b4975fffc90
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MyFFT_Analyzer/Window.lib Thu Sep 09 08:52:33 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/Window/#823e9a4ab223
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SerialTxRxIntr.lib Thu Sep 09 08:52:33 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/SerialTxRxIntr/#268977533f95
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/XferBase.hpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,49 @@
+//---------------------------------------------------------------------
+// データを PC へ転送するための抽象基底クラス
+//
+// 2021/07/11, Copyright (c) 2021 MIKAMI, Naoki
+//---------------------------------------------------------------------
+
+#include <string>
+#include "SerialRxTxIntr.hpp"
+using namespace Mikami;
+
+#ifndef XFER_ABSTRACT_BASE_HPP
+#define XFER_ABSTRACT_BASE_HPP
+
+class XferBase
+{
+public:
+ // コンストラクタ
+ XferBase(SerialRxTxIntr& rxTx, int size)
+ : SIZE_(size), xn_(size), rxTx_(rxTx) {}
+
+ // データを PC へ転送（0 ～ 10,000 の範囲の値を 2 文字で表すコード化を利用）
+ void ToPC(const float data[])
+ {
+ Convert(data);
+ string str = "";
+ for (int n=0; n<SIZE_; n++)
+ {
+ str += ((xn_[n] >> 7) & 0x7F) + 0x10;
+ str += (xn_[n] & 0x7F) + 0x10;
+ }
+ rxTx_.TxString(str+"\n");
+ }
+
+protected:
+ const int SIZE_; // PC に送るデータの数
+ Array<uint16_t> xn_; // PC に送るデータ
+
+private:
+ SerialRxTxIntr& rxTx_;
+
+ // データを転送する際の形式に変換
+ // data 元のデータ
+ virtual void Convert(const float data[]) = 0;
+
+ // コピー･コンストラクタおよび代入演算子の禁止のため
+ XferBase(const XferBase&);
+ XferBase& operator=(const XferBase&);
+};
+#endif // XFER_ABSTRACT_BASE_HPP
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/XferSpectrum.hpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,36 @@
+//---------------------------------------------------------------------
+// スペクトル解析の結果を PC へ転送するための XferBase の派生クラス
+//
+// 2021/07/11, Copyright (c) 2021 MIKAMI, Naoki
+//---------------------------------------------------------------------
+
+#include "XferBase.hpp"
+
+#ifndef XFER_SPECTRUM_DERIVED_HPP
+#define XFER_SPECTRUM_DERIVED_HPP
+
+class XferSpectrum : public XferBase
+{
+public:
+ // コンストラクタ
+ XferSpectrum(SerialRxTxIntr& rxTx, int size)
+ : XferBase(rxTx, size) {}
+
+private:
+ // スペクトル解析の結果を転送する形式に変換
+ // yAbs FFT の結果の絶対値
+ virtual void Convert(const float yAbs[])
+ {
+ const float MAX = 10000;
+ for (int n=0; n<SIZE_; n++)
+ {
+ float x = yAbs[n];
+ xn_[n] = (x > MAX) ? (uint16_t)MAX : (uint16_t)x;
+ }
+ }
+
+ // コピー･コンストラクタおよび代入演算子の禁止のため
+ XferSpectrum(const XferSpectrum&);
+ XferSpectrum& operator=(const XferSpectrum&);
+};
+#endif // XFER_SPECTRUM_DERIVED_HPP
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/main.cpp Thu Sep 09 08:52:33 2021 +0000
@@ -0,0 +1,110 @@
+//---------------------------------------------------------------------
+// FFT によるスペクトルアナライザ，白色雑音発生器付き (Nucleo-F446RE 用)
+// ● PC 側のプログラム： "CQ_FFT_Analyzer"
+// ● ボーレート
+// ポートの検索時： 9,600 baud
+// ポートの検索後： 460,800 baud
+//
+// ● 入力： A1
+// ● 白色雑音の出力：A2
+//
+// 2021/07/11, Copyright (c) 2021 MIKAMI, Naoki
+//---------------------------------------------------------------------
+
+#include <string>
+#include "Array.hpp"
+#include "DSP_AdcIntr.hpp"
+#include "DSP_Dac.hpp"
+#include "FFT_Analyzer.hpp"
+#include "DoubleBuffer.hpp"
+#include "Coefs_IIR_LP.hpp" // 縦続形 IIR フィルタの係数
+#include "IirCascade.hpp" // 縦続形 IIR フィルタ
+#include "MSeq16.hpp" // M 系列発生器
+#include "XferSpectrum.hpp"
+using namespace Mikami;
+
+#ifndef __STM32F446xx_H
+#error "Use Nucleo-F446RE"
+#endif
+
+const int N_FFT_ = 1024; // FFT の点数
+const int N_FRAME_ = N_FFT_; // １フレーム当たり標本化するデータ数
+const int N_SPC_ = N_FFT_/2 + 1; // 有効なスペクトルの点数
+const int RATIO_ = 10; // オーバーサンプリングの倍率
+const int N_TX_ = 501; // PC に転送するデータ数
+
+DoubleBuffer<float> buf_(N_FRAME_); // AD の結果を保存するバッファ
+DspAdcIntr myAdc_(10.24f*RATIO_, A1); // 標本化周波数： 102.4 kHz
+DspDac myDac;
+IirCascade df1_(ORDER1_, CK1_, G01_); // ダウンサンプリング用 Anti-alias フィルタ
+IirCascade df2_(ORDER2_, CK2_, G02_); // 白色雑音発生用低域通貨フィルタ
+MSeq16 mSeq_; // M 系列信号発生器（N = 16）
+
+// ADC 変換終了割り込みに対する割り込みサービス･ルーチン
+void AdcIsr()
+{
+ static int count = 0;
+
+ float xn = myAdc_.Read(); // AD 変換された値を取得
+
+ float noise = df2_.Execute(mSeq_.Execute());
+ myDac.Write(0.8f*noise); // 白色雑音出力
+
+ float yn = df1_.Execute(xn); // ダウンサンプリング用 Anti-alias フィルタの実行
+
+ if (++count >= RATIO_)
+ {
+ buf_.Store(yn); // ダウンサンプリングされたデータをバッファへ格納
+ count = 0;
+ buf_.IsFullSwitch(); // バッファが満杯であればバッファを切り替える
+ }
+}
+
+int main()
+{
+ SerialRxTxIntr rxTx; // PC との通信用，最初は 9600 baud
+ XferSpectrum tx(rxTx, N_TX_); // PC に転送するためのオブジェクトの生成
+ FftAnalyzer analyzer(N_FFT_); // FFT によるスペクトル解析オブジェクトの生成
+
+ Array<float> sn(N_FFT_, 0.0f); // スペクトル解析の対象となるデータ
+ Array<float> absFt(N_SPC_); // 解析結果：スペクトルの絶対値
+
+ NVIC_SetPriority(ADC_IRQn, 0); // AD変換終了割り込みの優先度が最高
+ NVIC_SetPriority(USART2_IRQn, 1);
+
+ bool ready = false; // スペクトルの計算終了で true
+ bool okGo = false; // "GO" を受信したら true
+
+ myAdc_.SetIntrVec(&AdcIsr); // AD変換終了割り込みの割り当て
+ while (true)
+ {
+ // PC からのコマンドの解析
+ if (rxTx.IsEol()) // 受信バッファのデータが有効になった場合の処理
+ {
+ string str = rxTx.GetBuffer();
+ if (str == "FFTAnalyzer")
+ {
+ rxTx.TxString("ACK\n"); // PC からの "FFTAnalyzer" に対して "ACK" を送信
+ wait_ms(10);
+ rxTx.Baud(460800); // 以降は 460,800 baud
+ }
+ if (str.substr(0, 2) == "GO")
+ okGo = true; // データの転送要求あり
+ }
+
+ if (buf_.IsFull()) // 入力データが満杯の場合，以下の処理を行う
+ {
+ for (int n=0; n<N_FRAME_; n++) sn[n] = buf_.Get(n);
+ analyzer.Execute(sn, absFt); // スペクトル解析の実行
+ ready = true; // スペクトル解析終了
+ }
+
+ // 転送要求がありスペクトル解析が終了している場合にデータを PC へ転送する
+ if (okGo && ready)
+ {
+ tx.ToPC(absFt); // データを PC へ転送
+ ready = false;
+ okGo = false;
+ }
+ }
+}
\ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mbed.bld Thu Sep 09 08:52:33 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/mbed_official/code/mbed/builds/65be27845400 \ No newline at end of file