dsPIC33FJ64GP802 (10) --- DSP 命令と DSP ライブラリ (8)

 「楽音」の「音名」でいうと「A」(880 Hz)と、「A#」(932.33 Hz) を中心周波数とする 1/12 オクターブ・バンド・フィルタふたつを IIRTransposed() 関数で実現し、内部で発生させたリニア周波数スイープのサイン波をフィルタリングした結果を下に示します。
 フィルタ係数は 12 kHz サンプリングで中央周波数 880 Hz と 932.33 Hz になる値を求めており、実際にはサンプリング周波数 48 kHz で動作させているので、中心周波数はそれぞれ 4 倍の、3.52 kHz と 3.729 Hz になっています。

 青色の線が 880 Hz (実際は 3.52 kHz)、赤色の線が 932.33 Hz (実際は 3.729 Hz) のフィルタの出力です。
 2 チャネル・ステレオ DAC で両者を同時に出力していますが、WaveSpectra では L/R どちらかを切り替えて 1 チャネルずつの観測しかできません。
 そのため、1 チャネルごとに観測して結果を .WSO ファイルに残し、後で 2 チャネル分の結果を読み込んでグラフ表示しています。
 ピーク周波数付近を、リニア周波数スケールで拡大した結果を下に示します。
 隣接バンドとの境界の周波数でゲインは -3 dB 程度、隣接バンドの中心周波数付近でゲインは約 -20 dB 程度になっています。

 プログラムで使っている IIR フィルタ係数の配列の内容を下に示します。 配列は X-データ領域に置くものとし、アライメントは自然な値の 2 バイトとしています。

// 1/12 octave band filter coeffs for IIRTransposed()
// center freq. = 880 Hz (A) @ fs = 12 kHz
// exported from Filter component of PSoC Creator
fractional _XDATA(2) IIRT_coef_880[] = {
  
// biquad section 1
  Q15( 0.00329756736755371/2), //  num10/2
  Q15( 0.0/2                ), //  num11/2
  Q15( 1.76804947853088/2),    // -den11/2
  Q15(-0.00329756736755371/2), //  num12/2
  Q15(-0.973729372024536/2),   // -den12/2

// biquad section 2
  Q15(  0.0366625785827637/2), //  num20/2
  Q15( -0.0733251571655273/2), //  mum21/2
  Q15(  1.7899694442749/2),    // -den21/2
  Q15(  0.0366625785827637/2), //  num22/2
  Q15( -0.987082719802856/2),  // -den22/2

// biquad section 3
  Q15(  0.0190212726593018/2), //  num30/2
  Q15(  0.0380425453186035/2), //  num31/2
  Q15(  1.76904010772705/2),   // -den31/2
  Q15(  0.0190212726593018/2), //  num32/2
  Q15( -0.98647665977478/2),   // -den32/2
}; // fractional IIRT_coef_880[]]
// 1/12 octave band filter coeffs for IIRTransposed()
// center freq. = 932.328 Hz (A#) @ fs = 12 kHz
// exported from Filter component of PSoC Creator
fractional _XDATA(2) IIRT_coef_932[] = {

// biquad section 1
  Q15( 0.00348567962646484/2), //  num10/2
  Q15( 0.0/2),                 //  num11/2
  Q15( 1.74199986457825/2),    // -den11/2
  Q15(-0.00348567962646484/2), //  num12/2
  Q15(-0.972188949584961/2),   // -den12/2

// biquad section 2
  Q15( 0.0359559059143066/2), //  num20/2
  Q15(-0.0719120502471924/2), //  num21/2
  Q15( 1.76560187339783/2),   // -den21/2
  Q15( 0.0359559059143066/2), //  num22/2
  Q15(-0.986324310302734/2),  // -den22/2

// biquad section 3
  Q15( 0.021756649017334/2),  //  num30/2
  Q15( 0.043513298034668/2),  //  num31/2
  Q15( 1.74227499961853/2),   // -den31/2
  Q15( 0.021756649017334/2),  //  num32/2
  Q15(-0.985674142837524/2),  // -den32/2
}; // fractional IIRT_coef_932[]]

 このフィルタ係数は、PSoC Creator の「Filter」コンポーネント付属の設計プログラムにより求めたものを以下の点を除いて「そのまま」利用しているものです。

  • 1/2 してから Q15 に量子化
  • 順序の変更および分母の符号反転
    (分子_0、分子_1、分子_2、 分母_1、分母_2
    →分子_0、分子_1、-分母_1、分子_2、-分母_2)

 伝達関数のペアリング/オーダリング/スケーリングは Filter コンポーネントが行っており、リニア周波数スイープする単一周波数の (ほぼ) フルスケール正弦波入力に対するフィルタリングで、この係数を利用してのクリッピング等のトラブルはありませんでした。
 PSoC Creator で回路図上に Filter コンポーネントを配置し、ダブル・クリックして現れるフィルタ設計画面で下のように入力します。

880 Hz のフィルタに対しては

  • Center = 0.88
  • Bandwidth = 0.05083786

 932.328 Hz のフィルタに対しては、

  • Center = 0.932328
  • Bandwidth = 0.0538608

と入力します。
 必要な値を入力したあとは、「Final Coefficients」タブをクリックすると、フィルタ係数の計算結果が表示されます。 880 Hz の場合の結果を下に示します。

 これをコピー・アンド・ペーストし、前述のような変更を加えて使用します。
 IIRCanonic() 用のフィルタ係数については、「自前」でスケーリングした結果を下に示します。

// 1/12 octave band filter coeffs for IIRCanonic()
// center freq. = 880 Hz @ fs = 12 kHz
fractional _XDATA(2) IIRC_coef_880[] = {

#define IIRC_prescale  (0.011000)
#define IIRC_postshift (0)

// biquad section 1
  Q15(-0.97371/2),   // -den12/2
  Q15( 1.76786/2),   // -den11/2
  Q15(-0.0076923/2), //  num12/2
  Q15( 0.0000000/2), //  num11/2  
  Q15( 0.0076923/2), //  num10/2  

// biquad section 2
  Q15(-0.98707/2),   // -den22/2
  Q15( 1.78979/2),   // -den21/2
  Q15(-0.023810/2),  //  num22/2
  Q15( 0.000000/2),  //  num21/2  
  Q15( 0.023810/2),  //  num20/2  

// biquad section 3
  Q15(-0.98646/2),   // -den32/2
  Q15( 1.76886/2),   // -den31/2
  Q15(-1.14000/2),   //  num32/2
  Q15( 0.00000/2),   //  num31/2  
  Q15( 1.14000/2),   //  num30/2  
  
// scaling factors
  Q15(0.5 * 0.011), // initialGain
  0                 // finalShift
}; // fractional IIRC_coef_880[]]

(2018 年 10 月 11 日追記: スケーリング・ファクタの扱いを変更しました。)
 IIRCanonic() の「initialGain」は、フィルタ係数の配列とは別の IIRCanonicStruct に設定しなければならないのが、ちょと不便なので、「initialGain」と「finalShift」の値をフィルタ係数配列の末尾に追加しておいて、一体として扱い、IIRCanonicStruct の初期化の際に参照して構造体を設定するようにしました。
 「自前」のスケーリングは、GNU Octave の「freqz()」関数で各段の周波数応答をプロットしてピーク・ゲインを「目視」しながら、ピーク・ゲインを 0 dB 以下に凹ませるように「手動」でスケーリング係数を調整して行いました。