RIAAフォノイコライザをIIRデジタルフィルタで実装する ⑤Octaveによる位相検証

かなり間があいてしまったが、この記事は、

RIAAフォノイコライザをIIRデジタルフィルタで実装する　①準備編
RIAAフォノイコライザをIIRデジタルフィルタで実装する　②係数計算編
RIAAフォノイコライザをIIRデジタルフィルタで実装する　③シミュレーション編
RIAAフォノイコライザをIIRデジタルフィルタで実装する　④実装・評価編

の続きで、IIRデジタルフィルタで実装したレコード再生用のRIAAイコライザに関するものです。

IIRデジタルフィルタは、位相特性も含めてアナログフィルタを等価的に再現できる。
とはいえ、扱う信号の周波数がナイキスト周波数（サンプリング周波数の1/2）に近づくにつれ誤差が大きくなる。
前回までのシミュレーション及び実機検証でも96kHz サンプリングにおいて、20kHzでおよそ＋0.5dBのゲイン誤差を確認している。
今回はゲイン誤差に加えて、位相誤差がサンプリング周波数によってどの程度ちがうかを確認するため、フリーで使用できる数値解析用のプログラミング言語Octaveを使って解析した。

①アナログCR RIAAイコライザ回路のLTSpiceでのシミュレーション
データ検証の基準として、アナログ構成のCR型RIAAイコライザをLTSpiceによってシミュレーションする。
シミュレーションした回路と、AC解析結果をそれぞれ図1，図2に示す。


図1．基準とするRIAAイコライザ回路
RIAAフォノイコライザをIIRデジタルフィルタで実装する　①準備編で計算した定数による回路。



図2．AC解析結果（Gain-Phase）
LTSpiceのAC解析シミュレーション。


②OctaveによるIIRデジタルRIAAイコライザのシミュレーション
数値解析用プログラミング言語Octaveを使って、RIAAフォノイコライザをIIRデジタルフィルタで実装する　②係数計算編で算出した係数を使い、RIAAイコライザ回路のシミュレーションを行った。
Octaveは非常に便利な解析用プログラミング言語で、行列計算を基本として、あらゆる数値解析に応用できる。
今回はOctaveのIIRフィルタ解析の機能を使ってシミュレーションを行った。
使用したOctaveは最新版のVer.6.4.0GUI。ダウンロードはこちら。
OctaveでIIRフィルタのシミュレーションを行うには次のようにする。

式1で与えられるIIRフィルタを考える。

  ……式1
ただし、
　、、　、　

このとき、
b = G*conv(b1,b2);
a = conv(a1,a2);
[h,w]=freqz(b,a,4096,'whole',96000);
を行うと、式1で示されるIIRフィルタ出力の角周波数ベクトルｗと応答周波数ベクトルhが得られるので、h,wをプロットすることでゲイン特性と位相特性のグラフが得られる。
ここで、convはベクトルの畳み込み演算（多項式同士の乗算）、freqzはデジタルフィルタの周波数応答を計算する関数、4096はサンプル数、96000はサンプリング周波数。
RIAAフォノイコライザをIIRデジタルフィルタで実装する　②係数計算編で計算したサンプリング周波数96kHzでの双一次変換によるRIAAイコライザの伝達関数は、

……式2
なので、
ｂ1＝［1　0］
ｂ2＝［1　‐0.96777105  0］
ａ１＝［1　0］
ａ２＝［1　-1.866859545　0.86728426］
G＝1
とすれば、このRIAAイコライザの特性がシミュレーションできる。
実際のコードは次のとおり。（コード中ではG=ｂ0＝1）
（注）subplotは複数のグラフを並べて表示する関数（今回はゲインのグラフと位相のグラフ）、semilogxはx軸を対数にする関数。

b0 = 1;
b1 = [1 0];
b2 = [1 -0.96777105  0];
a1 = [1 0];
a2 = [1 -1.866859545 0.86728426];
b = b0*conv(b1,b2);
a = conv(a1,a2);
[h,w]=freqz(b,a,4096,'whole',96000);
subplot(2,1,1);
semilogx(w,20*log10(abs(h)));
grid;
ylabel('Magnitude (dB)');
title('Bilinear Transform RIAA EQ @96kHz');
subplot(2,1,2);
semilogx(w,angle(h)*180/pi);
grid;
ylabel('Phase (Deg)');
xlabel('Normalized Frequency (Hz)');

このコードをOctaveのコマンドウィンドウにコピペすれば、次のグラフが表示される。


図3．双一次変換法による96kHzサンプリングRIAAイコライザ特性

※上記のプログラムリストを、テキストエディタでbilinear96.mと名前をつけてunicode（UTF-8）で保存し、
保存したホルダをOctaveの”現在のディレクトリ”に指定した上で、コマンドウィンドウで
>>bilinear96
としてENTER↓でも実行できる。


同じく192kHzサンプリングでは、

b0 = 1;
b1 = [1 0];
b2 = [1 -0.98375463 0];
a1 = [1 0];
a2 = [1 -1.93124941 0.93135924];
b = b0*conv(b1,b2);
a = conv(a1,a2);
[h,w]=freqz(b,a,4096,'whole',192000);
subplot(2,1,1);
semilogx(w,20*log10(abs(h)));
grid;
ylabel('Magnitude (dB)');
title('Bilinear Transform RIAA EQ @192kHz');
subplot(2,1,2);
semilogx(w,angle(h)*180/pi);
grid;
ylabel('Phase (Deg)');
xlabel('Normalized Frequency (Hz)');


図4．双一次変換法による192kHzサンプリングRIAAイコライザ特性


同様に384kHzサンプリングでは、

b0 = 1;
b1 = [1 0];
b2 = [1 -0.99184419 0];
a1 = [1 0];
a2 = [1 -1.96505172 0.96507966];
b = b0*conv(b1,b2);
a = conv(a1,a2);
[h,w]=freqz(b,a,4096,'whole',384000);
subplot(2,1,1);
semilogx(w,20*log10(abs(h)));
grid;
ylabel('Magnitude (dB)');
title('Bilinear Transform RIAA EQ @384kHz');
subplot(2,1,2);
semilogx(w,angle(h)*180/pi);
grid;
ylabel('Phase (Deg)');
xlabel('Normalized Frequency (Hz)');


図5．双一次変換法による384kHzサンプリングRIAAイコライザ特性


③アナログ回路およびサンプリング周波数別の特性比較
今回のテーマは、RIAAイコライザの理論値（＝アナログ回路のLTSpiceシミュレーション）とIIRフィルタによるデジタルRIAAでサンプリング周波数ごとに、どの程度誤差が出るかを検証することなので、図2～図5を比較する。
100hz、1kHz、10kHz、20kHzでのゲインと位相の値を表1に示す。これは得られたグラフにカーソルを当てて直読した。

表1．各回路の周波数ごとのゲイン、位相
※ゲインは1kHzで-20dBとなるように正規化した。


前回までの記事で、96kHzサンプリングのデジタルRIAAでは、20kHzのゲインにおよそ+0.5dB程度誤差が出ることが確認されているが、今回のシミュレーションでも同じ結果となった。
今回は位相について注目したい。一番下の黄色のマスが、LTSPICEによるアナログRIAA回路でのシミュレーション結果で、これを基準とする。
このときの各回路での位相誤差を図6に示す。


図6．位相特性の比較

図6の一番下の黄色いプロットが基準となるアナログ回路によるRIAAのシミュレーションで、20kHzにおける位相が-85°だ。
サンプリング周波数が384kHz、192kHz、96kHzと減少するにつれ、位相は-76°、-67°、-48°と乖離してゆく。
これが聴感上どの程度の差になって表れるか定量化することは難しいが、IIRによるRIAAイコライザで高域の位相特性まで忠実に再現しようとするなら、サンプリング周波数は高いほど良いということになる。
予想としては「96kHzのサンプリング周波数であれば必要十分」程度に考えていたので、意外な結果と言える。


●2022/02/26追記
RIAAフォノイコライザをIIRデジタルフィルタで実装する　②係数計算編で触れたが、双一次変換によって計算した伝達関数である式２は、
もともと式3として算出されたものに対してｚ^-1項の係数をネグって、z^-2項の係数をz^-1にシフトしたものだ。

　…式３


　…式２（再掲）

これはRIAAフォノイコライザをIIRデジタルフィルタで実装する　②係数計算編で書いた通り、式2を採用することで高域のゲイン誤差が抑えられるということがわかったが、Octaveで解析したところ式3のほうが位相誤差が小さい。
そのため追加評価として式3での解析結果を以下に記しておく。

式3によって計算すると、96k㎐サンプリングのOctaveコードと得られた特性は次のとおり。

b0 = 1;
b1 = [1 0];
b2 = [1 0.03222895 -0.96777105];
a1 = [1 0];
a2 = [1 -1.866859545 0.86728426];
b = b0*conv(b1,b2);
a = conv(a1,a2);
[h,w]=freqz(b,a,4096,'whole',96000);
subplot(2,1,1);
semilogx(w,20*log10(abs(h)));
grid;
ylabel('Magnitude (dB)');
title('Bilinear Transform RIAA EQ @96kHz(2)');
subplot(2,1,2);
semilogx(w,angle(h)*180/pi);
grid;
ylabel('Phase (Deg)');
xlabel('Normalized Frequency (Hz)');


図７．双一次変換法による96kHzサンプリングRIAAイコライザ特性（式3による）


同じく192k㎐サンプリングでは、

b0 = 1;
b1 = [1 0];
b2 = [1 -0.98375463 0];
a1 = [1 0];
a2 = [1 -1.93124941 0.93135924];
b = b0*conv(b1,b2);
a = conv(a1,a2);
[h,w]=freqz(b,a,4096,'whole',192000);
subplot(2,1,1);
semilogx(w,20*log10(abs(h)));
grid;
ylabel('Magnitude (dB)');
title('Bilinear Transform RIAA EQ @192kHz(2)');
subplot(2,1,2);
semilogx(w,angle(h)*180/pi);
grid;
ylabel('Phase (Deg)');
xlabel('Normalized Frequency (Hz)');


図８．双一次変換法による192kHzサンプリングRIAAイコライザ特性（式3による）


同様に384k㎐サンプリングでは、

b0 = 1;
b1 = [1 0];
b2 = [1 0.00815581 -0.99184419];
a1 = [1 0];
a2 = [1 -1.96505172 0.96507966];
b = b0*conv(b1,b2);
a = conv(a1,a2);
[h,w]=freqz(b,a,4096,'whole',384000);
subplot(2,1,1);
semilogx(w,20*log10(abs(h)));
grid;
ylabel('Magnitude (dB)');
title('Bilinear Transform RIAA EQ @384kHz(2)');
subplot(2,1,2);
semilogx(w,angle(h)*180/pi);
grid;
ylabel('Phase (Deg)');
xlabel('Normalized Frequency (Hz)');


図９．双一次変換法による384kHzサンプリングRIAAイコライザ特性（式3による）


以上のシミュレーション結果による各周波数でのゲイン、位相を表2に示す。

表２．各回路の周波数ごとのゲイン、位相（式３による）


以上より、式3による計算方式では、96k㎐サンプリングのゲインにおよぞ-1.4ｄBの誤差が出たが、位相誤差はほぼゼロとなった。
サンプリング192k㎐ではゲイン、位相ともに問題のないレベルで、必要十分といえそうだ。

なので、結論としては、

・リソースの都合で96k㎐サンプリングで行いたい場合は、ゲイン誤差を優先して式２による係数にするか、位相誤差を優先して式３による係数にするか要検討。

・192k㎐サンプリングが許される場合は、式３による係数を採用すればゲイン、位相ともにほぼ誤差のない結果が得られる。


【おまけ】
今回使用した数値解析用プログラミング言語Octaveは、有名なMATLAB（有償）と互換性があり、フリーで使えるツールです。
今回のようなオーディオの解析だけではなくありとあらゆる数値解析やグラフ表示に使うことができるのでたいへん有用です。

・3次元プロットの例


上の3Dプロットのプログラムは次のとおり。

tx = ty = linspace (-8, 8, 41)';
[xx, yy] = meshgrid (tx, ty);
r = sqrt (xx .^ 2 + yy .^ 2) + eps;
tz = sin (r) ./ r;
mesh (tx, ty, tz);
xlabel ("tx");
ylabel ("ty");
zlabel ("tz");
title ("3-D Sombrero plot");