« 2022年1月 | トップページ | 2022年3月 »

2022年2月

2022年2月26日 (土)

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に示す。

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


Ltspice_eq_gpneg
図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フィルタを考える。

Formura_iir  ……式1
ただし、
B1 、B2、 A1、 A2

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

96kiir_riaa……式2
なので、
b1=[1 0]
b2=[1 ‐0.96777105  0]
a1=[1 0]
a2=[1 -1.866859545 0.86728426]
G=1
とすれば、このRIAAイコライザの特性がシミュレーションできる。
実際のコードは次のとおり。(コード中ではG=b0=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のコマンドウィンドウにコピペすれば、次のグラフが表示される。

Bl96k
図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)');

Bl192k
図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)');

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


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

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

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

Samplingvsphase
図6.位相特性の比較

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


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

Formula_org …式3


96kiir_riaa …式2(再掲)

これは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)');

Riaa96k2
図7.双一次変換法による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)');

Riaa192k2
図8.双一次変換法による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)');

Riaa384k2
図9.双一次変換法による384kHzサンプリングRIAAイコライザ特性(式3による)


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

表2.各回路の周波数ごとのゲイン、位相(式3による)
Fullgp

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

なので、結論としては、

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

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


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

・3次元プロットの例
Octave3dgraph

上の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");

| | コメント (0)

2022年2月16日 (水)

PIC16F1705とBME280による温湿度気圧計

Pic_bme280
写真1.ブレッドボードに組んだ温湿度気圧計
PIC16F1705で温湿度気圧センサBME280モジュールを駆動し、LCD(AQM0802)に表示する。


今回はPICマイコンの開発環境MPLAB X IDEでMCC(Microchip Code Configurator)を使って、温湿度気圧センサBME280をI2Cで駆動して、測定した温湿度気圧をI2CのLCDパネルに表示する温湿度気圧計の試作をした。
PICマイコンのレジスタ設定等をGUIで行えるMCCを使う練習をするのが今回の主な目的だ。

これまで、PICの開発はMCCを使わずにやっていたが、最近のPICマイコンは機能が増え、レジスタ設定が非常に複雑になってきたため、手作業でレジスタ設定を行うことが困難になり始めていた。
Microchip側としてもそのあたりを考慮して、GUIで簡単にレジスタ設定ができるMCCをリリースした。
MCCを使うと、レジスタ設定は非常に簡単に行えるようになるが、その反面レジスタ分類ごとにファイルが生成されるため、慣れないとどこで何を設定しているのかがわかりにくく、全体として見通しが悪くなる。そのためいままで敬遠してきたが、今回使ってみて、慣れてくるととても便利だということがわかった。

【温湿度気圧計仕様】
今回試作した温湿度気圧計の仕様と開発環境は次のとおり。

PICマイコン      PIC16F1705-I/P
温湿度気圧センサ     BME280(I2C)
LCD           AQM0802(I2C)
開発環境         MPLAB X IDE v5.35
             XC8 v2.31
             MCC v4.0.2 ※
測定間隔         270.6秒
電源電圧         2.6~4.2V
消費電流         0.2~0.25mA

※MCCはMPLAB X IDEにプラグインとしてインストールする。


【回路】
16f1705_bme280sch_20220216104701

図1.回路図


図1に今回製作した回路図を示す。
温湿度気圧センサBME280は、3.3VLDOを搭載したモジュールを使用した。(写真2)

Bme280module
写真2.BME280モジュール
amazonやaliexpressなどで購入可能。

表示用LCDは秋月のAQM0802を使用した。これは安価でコンパクトなので、実験によく使用する。
ただ、電源電圧が変わると表示のコントラストが変わってしまう。そのため、今回はPICマイコンのA/Dコンバータ(ADC)で電池電圧をモニタし、電圧が3.5V以下に下がるとLCDコントラストを濃くするように調整している。
ADCには電池電圧を100kΩx2で分圧して入力しているが、ここで電流を消費しないように、測定時のみMOSFET AO3406で抵抗分圧をONしている。
主要部品の動作電圧範囲は、
PIC1705  2.3~5.5V
BME280   1.71~3.6V(LDOなしのBME280単体)
AQM0802  2.7~5.5V
となっており、BME280モジュールは3.3VLDO搭載なので、リチウムポリマ電池(最大4.2V)をそのまま使用した。
BME280とAQM0802はI2Cインタフェイスなので、SCLとSDAにプルアップ抵抗が必要だが、今回はBME280モジュールとAQM0802内部にそれぞれ10kΩが搭載されており、合成でそれぞれ5kΩのプルアップ抵抗になっている。

【ソフトウェア】
今回のソースファイルをプロジェクトホルダごとアーカイブしてこの記事の末尾にリンクしておく。
PICマイコンは、測定、表示が済むとSLEEPに入り、WDT(ウォッチドッグタイマ)によってSLEEPから目覚めて測定表示、再びSLEEP、を繰り返す。測定時は目印として一瞬LEDを点灯する。
WDTの最大時間は、MCCでポストスケーラーを最大値の8388608に設定したとき270.6秒となる。つまり、4分30秒に一度測定、表示更新をする。
WDTでSLEEPから目覚めたときは、SLEEPの直後の命令から復帰する。リセットではないので注意が必要だ。
LCDのコントラストが電池電圧によって変わってしまい、電圧が3Vを切るとほとんど見えなくなってしまうため、ADCで電池電圧をモニタし、3.5Vを堺にコントラスト設定を変えている。これはmain.cで電圧低下フラグBattLowを立ててi2c_lcd.cの113行目からの記述で行っている。

【動作】
SLEEP状態での消費電流はおよそ0.2mAで、測定時には瞬間的に1mA程度(LED点灯含む)となる。平均電流はおおむね0.2~0.25mA程度。
400mAhのリチウムポリマ電池での実動作では72日程度の動作を確認した。
SLEEP状態での消費電流は、BME280モジュールに搭載されたLDOのIqが大きいのではないかと思っていたが、実際に確認してみるとIqでの消費はきわめて小さく、LCD(AQM0802)の消費電流が支配的だった。

【ソース】

ダウンロード - 20220214bmeaqm_1705adc_000ok.x.zip


そういうわけでMCCはとても便利なので、どんどん使いましょう(^-^)
MCCの使い方は、トランジスタ技術2021年9月号別冊付録「PIC開発マニュアル2021」を参考にしました。

| | コメント (0)

« 2022年1月 | トップページ | 2022年3月 »