FFT によるスペクトル解析器
Dependencies: Array_Matrix mbed SerialTxRxIntr UIT_FFT_Real DSP_ADDA
Revision 0:91cc5a03f0ca, committed 2021-01-08
- Comitter:
- MikamiUitOpen
- Date:
- Fri Jan 08 02:27:46 2021 +0000
- Commit message:
- 1
Changed in this revision
diff -r 000000000000 -r 91cc5a03f0ca Array_Matrix.lib --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Array_Matrix.lib Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/Array_Matrix/#d3aa1ddb57e1
diff -r 000000000000 -r 91cc5a03f0ca DSP_ADDA.lib --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DSP_ADDA.lib Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/DSP_ADDA/#a1dcee67c67e
diff -r 000000000000 -r 91cc5a03f0ca DoubleBuffer.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DoubleBuffer.hpp Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,58 @@ +//-------------------------------------------------------- +// ダブル・バッファの template クラス +// 内部のバッファは通常の配列を使用 +// +// 2019/11/22, Copyright (c) 2019 MIKAMI, Naoki +//-------------------------------------------------------- + +#ifndef DOUBLE_BUFFER_2DARRAY_HPP +#define DOUBLE_BUFFER_2DARRAY_HPP + +template<class T, int N> class DoubleBuffer +{ +public: + // コンストラクタ + explicit DoubleBuffer(T initialValue) + : ping_(0), pong_(1), index_(0), full_(false) + { + for (int k=0; k<2; k++) + for (int n=0; n<N; n++) buf_[k][n] = initialValue; + } + + // データを格納 + void Store(T data) { buf_[ping_][index_++] = data; } + + // 出力バッファからデータの取り出し + T Get(int n) const { return buf_[pong_][n]; } + + // バッファが満杯でバッファを切り替える + void IfFullSwitch() + { + 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: + T buf_[2][N]; // 標本化したデータのバッファ + int ping_, pong_; // バッファ切替用 + int index_; // 入力データのカウンタ + bool full_; // 満杯の場合 true + + // コピー・コンストラクタおよび代入演算子の禁止のため + DoubleBuffer(const DoubleBuffer&); + DoubleBuffer& operator=(const DoubleBuffer&); +}; +#endif // DOUBLE_BUFFER_2DARRAY_HPP +
diff -r 000000000000 -r 91cc5a03f0ca IIR_Filter/Biquad.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/IIR_Filter/Biquad.hpp Fri Jan 08 02:27:46 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
diff -r 000000000000 -r 91cc5a03f0ca IIR_Filter/Coefs_IIR_LP.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/IIR_Filter/Coefs_IIR_LP.hpp Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,48 @@ +//----------------------------------------------------- +// 縦続形 IIR フィルタの次数と係数の定義 +// +// 2020/11/08, Copyright (c) 2020 MIKAMI, Naoki +//----------------------------------------------------- + +#include "Biquad.hpp" + +// FFT アナライザで使うフィルタ +// 標本化周波数が 100 kHz の場合,5 kHz 以上で +// 少なくとも 60 dB 減衰させる LPF + +// 低域通過フィルタ +// 連立チェビシェフ特性 +// 次数 : 10 次 +// 標本化周波数:100.00 kHz +// 遮断周波数 : 4.80 kHz +// 通過域のリップル: 0.50 dB +// 阻止域の減衰量 :60.00 dB +const int ORDER1_ = 10; +const Biquad::Coefs CK1_[] = { + { 1.824022E+00f, -8.381143E-01f, -9.384721E-01f, 1.0f}, // 1段目 + { 1.856865E+00f, -9.028289E-01f, -1.790351E+00f, 1.0f}, // 2段目 + { 1.884825E+00f, -9.577163E-01f, -1.875764E+00f, 1.0f}, // 3段目 + { 1.898856E+00f, -9.845865E-01f, -1.896090E+00f, 1.0f}, // 4段目 + { 1.905948E+00f, -9.962457E-01f, -1.901699E+00f, 1.0f}}; // 5段目 +const float G01_ = 1.221772E-03f; // 利得定数 + + + +// 白色雑音発生器で使うフィルタ +// 標本化周波数が 100 kHz の場合, +// 100/(4π) ≒7.96 kHz 以上で +// 少なくとも 60 dB 減衰させる LPF + +// 低域通過フィルタ +// 連立チェビシェフ特性 +// 次数 : 6 次 +// 標本化周波数:100.00 kHz +// 遮断周波数 : 5.00 kHz +// 通過域のリップル: 0.50 dB +// 阻止域の減衰量 :60.00 dB +const int ORDER2_ = 6; +const Biquad::Coefs CK2_[] = { + { 1.789361E+00f, -8.093387E-01f, -6.031744E-01f, 1.0f}, // 1段目 + { 1.825484E+00f, -8.899602E-01f, -1.690793E+00f, 1.0f}, // 2段目 + { 1.869667E+00f, -9.677499E-01f, -1.803334E+00f, 1.0f}}; // 3段目 +const float G02_ = 1.404132E-03f; // 利得定数 \ No newline at end of file
diff -r 000000000000 -r 91cc5a03f0ca IIR_Filter/IirCascade.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/IIR_Filter/IirCascade.hpp Fri Jan 08 02:27:46 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
diff -r 000000000000 -r 91cc5a03f0ca MSeq16.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MSeq16.hpp Fri Jan 08 02:27:46 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
diff -r 000000000000 -r 91cc5a03f0ca MyFFT_Analyzer/Blackman.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MyFFT_Analyzer/Blackman.hpp Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,50 @@ +//------------------------------------------------------------------- +// Blackman 窓による窓掛け +// ゼロ詰め(zero-padding)の機能を持つ +// +// 2020/11/07, Copyright (c) 2020 MIKAMI, Naoki +//------------------------------------------------------------------- + +#ifndef HAMMING_WINDOW_HPP +#define HAMMING_WINDOW_HPP + +#include "mbed.h" +#include "Array.hpp" + +namespace Mikami +{ + class BlackmanWindow + { + public: + // コンストラクタ + BlackmanWindow(uint16_t nData, uint16_t nFft) + : N_(nData), NFFT_(nFft), w_(nData) + { + const float PI2L = 6.283185f/(float)nData; + const float PI4L = 2*PI2L; + for (int k=0; k<nData; k++) + w_[k] = 0.42f - 0.5f*cosf(k*PI2L) + 0.08f*cosf(k*PI4L); + } + + // デストラクタ + virtual ~BlackmanWindow() {} + + // 窓掛けを実行 + void Execute(const Array<float> &x, Array<float> &y) + { + for (int n=0; n<N_; n++) y[n] = x[n]*w_[n]; + for (int n=N_; n<NFFT_; n++) y[n] = 0; + } + + private: + const int N_; + const int NFFT_; + + Array<float> w_; + + // コピー・コンストラクタおよび代入演算子の禁止のため + BlackmanWindow(const BlackmanWindow& ); + BlackmanWindow& operator=(const BlackmanWindow& ); + }; +} +#endif // HAMMING_WINDOW_HPP \ No newline at end of file
diff -r 000000000000 -r 91cc5a03f0ca MyFFT_Analyzer/FFT_Analyzer.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MyFFT_Analyzer/FFT_Analyzer.cpp Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,32 @@ +//------------------------------------------------------- +// FFT を使ってスペクトル解析を行うクラス +// +// 2020/11/07, Copyright (c) 2020 MIKAMI, Naoki +//------------------------------------------------------- + +#include "FFT_Analyzer.hpp" + +namespace Mikami +{ + FftAnalyzer::FftAnalyzer(int nFft) + : N_FFT_(nFft), fft_(nFft), wHm_(nFft, 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; + + wHm_.Execute(xData_, wData_); // 窓掛け + + 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
diff -r 000000000000 -r 91cc5a03f0ca MyFFT_Analyzer/FFT_Analyzer.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MyFFT_Analyzer/FFT_Analyzer.hpp Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,39 @@ +//------------------------------------------------------- +// FFT を使ってスペクトル解析を行うクラス(ヘッダ) +// +// 2020/11/07, Copyright (c) 2020 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 wHm_; + + 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
diff -r 000000000000 -r 91cc5a03f0ca SerialTxRxIntr.lib --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SerialTxRxIntr.lib Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/SerialTxRxIntr/#deeef404ff49
diff -r 000000000000 -r 91cc5a03f0ca UIT_FFT_Real.lib --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UIT_FFT_Real.lib Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/MikamiUitOpen/code/UIT_FFT_Real/#0b4975fffc90
diff -r 000000000000 -r 91cc5a03f0ca Xfer.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Xfer.hpp Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,58 @@ +//--------------------------------------------------------------------- +// スペクトル解析の結果を PC へ転送するためのクラス +// +// 2020/02/05, Copyright (c) 2020 MIKAMI, Naoki +//--------------------------------------------------------------------- + +#include <string> +#include "Array.hpp" +#include "SerialRxTxIntr.hpp" +using namespace Mikami; + +#ifndef XFER_CONVERT_TOPC_HPP +#define XFER_CONVERT_TOPC_HPP + +class Xfer +{ +public: + // コンストラクタ + Xfer(SerialRxTxIntr& rxTx, int size) + : SIZE_(size), xn_(size), rxTx_(rxTx) {} + + // スペクトル解析の結果を転送する形式に変換 + // yAbs FFT の結果の絶対値 + // magn 倍率 + void Convert(const float yAbs[], float magn = 1.0f) + { + const float MAX = 10000; + for (int n=0; n<SIZE_; n++) + { + float x = magn*yAbs[n]; + xn_[n] = (x > MAX) ? MAX : (uint16_t)x; + } + } + + // データを PC へ転送(0 ~ 10,000 の範囲の値を 2 文字で表すコード化を利用) + void ToPC() + { + string str = ""; + for (int n=0; n<SIZE_; n++) + { + div_t a = div(xn_[n], 100); + str += a.quot + 0x10; + str += a.rem + 0x10; + } + rxTx_.TxString(str+"\n"); + rxTx_.TxString("EOT\n"); + } + +private: + const int SIZE_; // PC に送るデータ数 + Array<uint16_t> xn_; // PC に送るデータ + SerialRxTxIntr& rxTx_; + + // コピー・コンストラクタおよび代入演算子の禁止のため + Xfer(const Xfer&); + Xfer& operator=(const Xfer&); +}; +#endif // XFER_CONVERT_TOPC_HPP \ No newline at end of file
diff -r 000000000000 -r 91cc5a03f0ca main.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main.cpp Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,113 @@ +//--------------------------------------------------------------------- +// FFT によるスペクトルアナライザ,白色雑音発生器付き (Nucleo-F446RE 用) +// +// ● ST-Link Firmware の V2.J37.M26 で動作確認 +// +// ● ST-Link Firmware のアップグレードには stsw-link007_V2-37-26.zip +// に含まれている "ST-LinkUpgrade.exe" を使う +// +// ● PC 側のプログラム: "F446_FFT_Analyzer" +// ● ボーレート: 460800 baud +// ● 受信データの文字列の終了マーク: "\r" +// +// ● 入力: A1 +// ● 白色雑音の出力:A2 +// +// 2021/01/08, Copyright (c) 2021 MIKAMI, Naoki +//--------------------------------------------------------------------- + +#include "mbed.h" +#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 "Xfer.hpp" +using namespace Mikami; + +const int N_FFT_ = 1024; // FFT の点数 +const int N_FRAME_ = N_FFT_; // 1フレーム当たり標本化するデータ数 +const int N_SPC_ = N_FFT_/2 + 1; // 有効なスペクトルの点数 +const int RATIO_ = 10; // オーバーサンプリングの倍率 + +DoubleBuffer<float, N_FRAME_> buf_(0); // AD の結果を保存するバッファ +DspAdcIntr myAdc_(10*RATIO_, A1); // 標本化周波数: 100 kHz +DspDac myDac; +IirCascade df1_(ORDER1_, CK1_, G01_); // ダウンサンプリング用 Anti-alias フィルタ +IirCascade df2_(ORDER2_, CK2_, G02_); // 白色雑音発生用 +MSeq16 mSeq_; + +// ADC 変換終了割り込みに対する割り込みサービス・ルーチン +void AdcIsr() +{ + static int count = 0; + + float xn = myAdc_.Read(); + + 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_.IfFullSwitch(); // バッファが満杯であればバッファを切り替える + } +} + +int main() +{ + SerialRxTxIntr rxTx(32, 460800); // PC との通信用 + Xfer tx(rxTx, N_SPC_); // 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 + float magn = 1; // 転送する際の倍率 + + myAdc_.SetIntrVec(&AdcIsr); // AD変換終了割り込みの割り当て + while (true) + { + // PC からのコマンドの解析 + if (rxTx.IsEol()) // 受信バッファのデータが有効になった場合の処理 + { + string str = rxTx.GetBuffer(); + if (str == "FFTAnalyzer") + rxTx.TxString("ACK\n"); // PC からの "FFTAnalyzer" に対して "ACK" を送信 + if (str.substr(0, 2) == "GO") + { + okGo = true; // データの転送要求あり + float x = (float)atoi(str.erase(0, 2).c_str()); + magn = powf(10, x/20.0f); + } + } + + if (buf_.IsFull()) // 入力データが満杯の場合,以下の処理を行う + { + for (int n=0; n<N_FRAME_; n++) sn[n] = buf_.Get(n); + analyzer.Execute(sn, absFt); // スペクトル解析の実行 + tx.Convert(absFt, magn); // スペクトル解析の結果を転送する形式に変換 + ready = true; // スペクトル解析終了 + } + + // 転送要求がありスペクトル解析が終了している場合にデータを PC へ転送する + if (okGo && ready) + { + tx.ToPC(); // データを PC へ転送 + ready = false; + okGo = false; + } + } +} \ No newline at end of file
diff -r 000000000000 -r 91cc5a03f0ca mbed.bld --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mbed.bld Fri Jan 08 02:27:46 2021 +0000 @@ -0,0 +1,1 @@ +https://os.mbed.com/users/mbed_official/code/mbed/builds/65be27845400 \ No newline at end of file