がれすたさんのDIY日記

電子回路、Python、組み込みシステム開発、自作エフェクターを語るblog

CMSIS DSP:biquad-low-passをやってみる

今回はフィルタのお話です。
今まで作ってたSAIペリフェラルを使ったデジタルエフェクタのプラットフォームがある程度動作するようになったので部屋にこもってオシロスコープとお友達になりながら格闘していました。
んなわけでいくつかエフェクト使った後にフィルタないといろいろできないやん!ってことでCMSIS DSP APIを使って高速フィルタ計算をやらせてみたって記事です。

んなわけで使うのが比較的有名なbiquad filter(双二次IIR)ってやつ
フィルタ係数を変えるとハイパスだったりバンドパスだったりになったり便利な奴です。

ARM CMSIS DSPの図を引用するとこんな感じのフィルタ
タイプで言うと直接型Ⅱの転置構成
f:id:gsmcustomeffects:20170930205453p:plain

んなわけで係数は5つ

というか正確には6つ


\displaystyle
\frac{b_{0} + b_{1}z^{-1}+b_{2}z^{-2}}{a_{0} + a_{1}z^{-1} + a_{2}z^{-2}}

この式ではa0にあたる部分があるのだが基本形はここを1としてa0で割った形で表現する。(ここの数値がクオリティファクターにかかわるがARM CMSIS DSPのBiquadでは1で固定なのでカットオフ点で-3dBになる。
というか定数項は1にして正規化しろって習った気がする
次にこいつを差分方程式で表していく

\displaystyle y[n] = b_{0}  x[n] + b_{1}  x[n-1] + b_{2}  x[n-2] - a_{1}  y[n-1] - a_{2}  y[n-2]

ここでn-1,n-2ってのはディレイ信号でありn-1が一個前の信号n-2が二個前の信号ということになる。

ARMのDocumentではこう解説してる

y[n] = b0 * x[n] + d1;
d1 = b1 * x[n] - a1 * y[n] + d2;
d2 = b2 * x[n] - a2 * y[n];

まずd2の信号から考えると入力にb2をかけたものと出力に-a2をかけたものを足している。
示すとこの部分だ

f:id:gsmcustomeffects:20170930224020p:plain

次にd1の部分を考える。
先ほど求めたd2をまとめて書くとこういう感じ

f:id:gsmcustomeffects:20170930224206p:plain

最後にそいつらをまとめて

f:id:gsmcustomeffects:20170930224434p:plain

そんなわけでこういう式を書けば実現できるわけだがARMさんがCMSIS DSP Libraryと言うのを提供していてそれを使うと高速演算できるので今回はこれを使わせていただく。

CMSIS DSP Lib

てなわけでCMSIS DSPを使っていく
ライブラリの導入はこの記事が詳しいので読んでくれ。

gsmcustomeffects.hatenablog.com

ここにある通り導入が済んだらBiquadフィルタのAPIを使うための準備をしていく
順序はこの通り

  1. インスタンスの生成
  2. 初期化
  3. APIのコール

インスタンスの生成

arm_biquad_cascade_df2T_instance_f32 S;

このように宣言するちなみにこいつはこういう構造を持ってる

  typedef struct
  {
    uint8_t numStages;         /**< number of 2nd order stages in the filter.  Overall order is 2*numStages. */
    float32_t *pState;         /**< points to the array of state coefficients.  The array is of length 2*numStages. */
    float32_t *pCoeffs;        /**< points to the array of coefficients.  The array is of length 5*numStages. */
  } arm_biquad_cascade_df2T_instance_f32;

ステージ数とディレイバッファと係数である。
そのうち後者の二つはポインタ型なので適宜自分で定義してやる必要がある。

float32_t pCoeffs1[5];
float32_t buffer1[2];

定義したら係数を適当なところで代入してやる。
今回はカットオフ1kHzにしたのでこんな感じ。

  pCoeffs1[0] = 0.0010232;//b0
  pCoeffs1[1] = 0.0020464;//b1
  pCoeffs1[2] = 0.0010232;//b2
  pCoeffs1[3] = -1.0*-1.90750;//a1
  pCoeffs1[4] = -1.0*0.911594497;//a2

ここで注意が一つその辺のIIR設計ソフトウエアで設計すると係数が符号なしで出てくるのでブロック線図を見て符号をきちんと読み替える必要がある。
f:id:gsmcustomeffects:20170930225838p:plain

初期化

先ほど宣言した構造体メンバをセットしてくれるAPIがある。

arm_biquad_cascade_df2T_init_f32(&S, 1, pCoeffs1, buffer1);

第二引数の1はフィルタの段数なので2次系なら1、4次にしたかったら2(その数だけバッファ変数bufferの大きさを確保する)

APIコール

float in,out;
arm_biquad_cascade_df2T_f32(&S, &in, &out, 1);

ここで第四引数はブロック数
自分の場合はサンプルバイサンプル処理なので1にしている。
in outに値する部分は処理系に準じて適宜変更してくれていい。

自分ならこんな感じ

f:id:gsmcustomeffects:20170930230408p:plain

実際にやってみる。

ファンクションジェネレーターで1kHzの波形を入れてみる。

f:id:gsmcustomeffects:20170930230543p:plain

カットオフ点なので-3dBぐらいになっている。

信号処理時間

f:id:gsmcustomeffects:20170930231449p:plain

一応1us程度だった。
自分のリアルタイム信号処理環境では48kHzサンプリングなので21us(たぶん19uSecぐらいが限界)までは余裕がある。

キャッシュの有無による処理時間の差異

現在使っているSTM32F767ではIキャッシュとDキャッシュがある。
それに関して考察しているブログがあるので明確な違いはそこを見てほしい。

というわけで先ほどの波形はキャッシュをオンにした状態だがキャッシュをオフにしたときはこのぐらい遅くなる。
f:id:gsmcustomeffects:20170930232605p:plain

二倍近く遅くなっていることがわかる。

まとめ

CMSIS DSP APIでフィルターを動かすことができた。
今後は複数の処理を直列にしたデジタル信号処理をやって行きたい。