« RIAAフォノイコライザをIIRデジタルフィルタで実装する ①準備編 | トップページ | RIAAフォノイコライザをIIRデジタルフィルタで実装する ③シミュレーション編 »

2020年10月25日 (日)

RIAAフォノイコライザをIIRデジタルフィルタで実装する ②係数計算編(2022/2/26修正追記)

前回はRIAAフォノイコライザの概念と、アナログで実装する場合の復習をした。
今回はいよいよRIAAフォノイコライザをデジタルフィルタに実装するための計算に入る。

 

1.デジタルフィルタとはなにか

アナログ信号が時間的に連続な信号であるのに対して、デジタル信号は時間的に離散している信号だ。つまりアナログでサイン波は時間軸に対してなめらかな正弦波だが、これをデジタルで表現すると、サンプリング周期ごとに値が1つずつ出てくる。たとえばCDに1khzのサイン波が入っているとすれば、CDのサンプリング周波数は44.1kHz、分解能は16bitなので、1/44100秒ごとに-32768(1000000000000000)~32767(0111111111111111)の範囲の値がひとつ出てきて、これを並べてやると飛び飛びの点だけど、1kHzのサイン波が点描されるというわけだ。つまりデジタルで表現されるデータは、サンプリング周期ごとに出てくる数を並べた数列になっている。

毎朝ベランダの温度を記録しているとする。日によっては急に寒くなったり暑くなったりすることがあるので、たとえば半年間、毎朝の温度をプロットしてグラフを描くと、かなり細かくギザギザしてしまう。これを、全体のおおきな変化傾向だけ見たいのでギザギザを減らしたいとしたらどうしたらいいだろうか。たとえば当日を含む過去10日間を足して1/10すれば、常にその日を起点とする過去10日間の平均が得られて、一日の影響度は1/10に抑えられるので、グラフの細かいギザギザが減る。
これを移動平均フィルタといい、デジタルフィルタの考え方の基本だ。(細かなギザギザがなくなったということは、高い周波数成分がカットされたことになるので、これはローパスフィルタだと考えることができる。)
この例を上のサンプリング時間ごとの値という考え方に当てはめれば、一日に1個のデータが出るので、サンプリング周期は一日ということになる。気温のデータ分解能は、-50~50度の範囲を8ビット(0~255)で表わすとすれば、100/256=0.39つまり、0.39度の分解能で表現される。

つまり、CDから16ビットの音楽データが1/44100秒ごとに1個出てくるのも、8ビットのベランダの温度データが一日に1個出てくるのも、現象としては同じなのでとくにむずかしいことはない。

上の移動平均フィルタの例では、最新のデータを起点とする過去10個の入力データを10で割って平均したものを出力とした。言い換えると、過去10個のデータそれぞれを1/10倍して足し合わせているともいえる。過去のデータは10個に限らず何百個でも何千個でもいいし、それぞれにかける係数も1/10に限らず必要に応じて決めてやると、さまざまな特性のフィルタが構成できる。このような構成のフィルタをFIRフィルタという。FIRフィルタの特徴は直線位相特性を持ち、群遅延が全周波数で一定となること、つまりあらゆる信号に対して常にサンプリング時間の遅れでデータが出力されるということだ。

FIRフィルタでは常に過去の入力データのみに係数を掛けて足し合わせているが、そこに過去の出力データにも係数を掛けて足し合わせる、という処理を加えたらどうなるだろうか。これをIIRフィルタといい、FIRにくらべて少ないデータ数で構成でき、アナログで構成したフィルタ回路と等価的に置き換えることができる。もちろん位相特性も再現される。

以上がものすごく大ざっぱなデジタルフィルタの説明だ。くわしくは専門書を参照してほしい。

 

さて今回のミッションは、アナログレーコード再生に使うRIAAイコライザ回路のデジタル化だ。RIAAイコライザも一種のフィルタなので、デジタルフィルタで再現することができる。

前回の記事で書いたが、レコードは逆RIAA特性で録音され、再生時はその逆の特性で再生する。つまり、

Trec・・・式16

Tplay・・・式17

式16で示される逆RIAAの特性で録音されたレコードを式17で示されるRIAA特性で再生されることで、特性が1(フラット)になるということだ。特性というのは振幅特性と位相特性の両方が含まれているので、当然のことながら位相も逆特性になっている。

(注)録音時の式16は、おそらく実際には高域の制限がかかっているので正確ではないが、便宜上このように考える

位相特性も含めて式17と等価なデジタルフィルタを再現するには、IIRフィルタで構成する。

 

2.IIRフィルタの構造

IIRフィルタのブロック図を図8に示す。

Biquad_block

図8.IIR(BiQuad)フィルタブロック図
IIRフィルタの中でも2つ過去までのデータを使う2次IIRフィルタを特にBiQuad回路と呼び、よく使われる

 

図8に示すのは2次IIRフィルタで、現在の入力データをx(n)、現在の出力データをy(n)とすると、2つ過去のデータx(n-2)、y(n-2)まで使用されている。原理的には n-3,n-4…というさらに過去のデータを使用する高次のIIRフィルタも構成可能だが、現実的には図示の2次IIRを必要に応じて多段構成にすることが多い。

図中z-1はz変換の遅延演算子で、z-1を乗じるごとにひとつ過去のデータとなることを示している。x(n)はz-1を一回経るとx(n-1)(ひとつ過去のデータ)、もう一度z-1を経るとx(n-2)となっている。出力データy(n)も同様で、y(n-1)はひとつ過去の出力、y(n-2)はふたつ過去の出力だ。それぞれのデータにそれぞれの係数-a1,-a2,b0,b1,b2が掛けられ、加算されて最新の出力データとして出力される。すなわち出力データy(n)は、

Yn・・・式18

この式はb0~b2、a1,a2が決まればこのままマイコンのプログラムやFPGAに実装できる。その場合、ひとつ過去、ふたつ過去のデータはバッファに格納しておいて、常に最新の入力データx(n)とともに演算して出力データy(n)を演算で求める。この形式の数式を差分方程式という。

式18から伝達関数を求めるには、Z変換を使う。Z変換は次のようになっている。

X_ztransfer

X2_ztransfer_20201024182201

x(n)のz変換をX(z)としたとき、2つ過去のx(n-2)のz変換はX(z)・z-2となる。とりあえずここでは、Z変換はそういう規則の書き換えだと思っておいてほしい。詳しくは専門書を参照のこと。

それでは式18の伝達関数Y(z)/X(z)を求めよう。

Yzxz

上記のように式18の両辺をZ変換して整理し、Y(z)/X(z)を求める。

Yzxz_0・・・式19

これが式18の伝達関数だ。ただし、b0を1に正規化することが多いので、実用的には

Yzxz_1・・・式20

としておくほうが便利かもしれない。
この場合の差分方程式は、

Xysabun_1・・・式21

こうすることで、最新の入力データx(n)をそのまま使えるので、伝達関数やコードが見やすくなる。
この場合、式20右辺の先頭に掛かっているb0は定数なので、省いてしまっても特性に影響はない。

伝達関数は、具体的には

Riaa_iir_ex0_20201025134601 ・・・式22

のように表わされる。これを見れば係数がすべてわかるので、直ちにハードウェアに実装できる。
この例では分子のz-2の項の係数は0なので書いていない。
慣れてくると、「ああこれは96kHzサンプリングのRIAAですね。」と、一瞥しただけでわかるようになる。(嘘)

 

3.フィルタ係数の計算

以上に書いたように、デジタルフィルタの設計とは、各係数を計算して求めることに他ならない。インターネットで探せばRIAAの係数はいろいろなところで拾ってこられるので、計算ができなくても作ることはできる。ただせっかくなので、今回は係数を求める計算もやってみる。

いきなり計算に入る前に、式17のRIAA伝達関数を見てほしい。この式では分母がsの2次式、分子がsの1次式になっている。伝達関数の見方として、分母=0を満たすsを「極」、伝達関数=0を満たすsを「ゼロ点」という。この式を見ると分母が2次式なので極が2つ、分子を評価すると一次式なのでゼロ点が1つだということがわかる。この極とゼロ点の個数はデジタルフィルタに置き換えても継承されるので、これから設計しようとしているIIRによるRIAA回路も分母が2次式、分子が1次式になる。(zの場合はマイナスなので、2次ならz-2だ)
つまりこれから求めようとしているRIAA回路は式22のような形になるはずだ。これを頭に入れておいてほしい。

注)この記事を書いている2020/10/24時点で、理解できていないところがあって、それも含めて以下の文章を書きます。いずれ理解が深まって解決したら、その時に追記する予定です。

 

3-1.双一次S-Z変換による係数計算

 

何度も出てきているようにアナログRIAAイコライザ回路の伝達関数をラプラス演算子で表わすと式17になる。

Tplay・・・式17

これをZ変換のz演算子で表現することが、このミッションの目的だ。
つまりRIAAの伝達関数を、ラプラス変換のsドメインからZ変換のzドメインに置き換えることだ。

Z変換は両側ラプラス変換を離散化したもので、Tをサンプリング周期とすると、Z変換におけるzとラプラス変換sは次のように対応している。

Z_vs_e_1

そしてこの置き換えから、次の式をラプラス変換領域の伝達関数に代入すればZ変換領域に変換できるとされている。

S_z_henkankousiki・・・式23

これは双1次s-Z変換法という、一般的によく使われる方法だという。だがまてよ、

Z_vs_e_2・・・式24

であるならば、

Z_vs_e_3・・・式25

となるはずだ。どうも納得がいかない、とあちこち調べてみたらこれは、e^xのマクローリン展開による近似式だそうだ。つまり

Maclaurin・・・式26

e^xのマクローリン展開は式26だ。そしてこの最初の2項1+xを使うと

Maclaurin_z_s・・・式27

となり、これをsのついて解くと式23になる。
それにしてもe^xのマクローリン展開をたった2項の1+xで打ち切るって、乱暴すぎないだろうか。誤差とかだいじょうぶなのかな。
ちなみに式25に出てくるs=logzは複素数の対数であり、これは多価関数なのでそのまま扱うのが困難なのだそうだ。

それでは教科書どおり式23で話を進める。

式17のRIAA伝達関数に式23を代入して、Z領域の伝達関数に変換してみる。ここではサンプリング周期Tに対し、サンプリング周波数Fsを、Fs=1/Tとして扱う。サンプリング周波数Fsは具体的には44100、48000,96000などのなじみのある数値だ。

Z_fs_s

Tplays_z Riaa_s_z_20201025130501 ・・・式28

かなり恐ろしげな式になってしまったが、恐れることはない。Fs,T1,T2,T3はすべて定数なので、電卓かEXCELで計算できる。
ただここで注目してほしいのは、分子にz-2の項ができてしまっている点だ。この章の冒頭で書いたとおり、いま求めているRIAA特性は極が2つゼロ点が1つの伝達関数になるので、分母が2次式、分子が1次式になるはずだが、上の計算では分子が2次式になってしまっている。これを書いている2020/10/25時点でわからない点だ。
ともあれ、式28で数値計算するとどうなるのか見てみよう。

T1~T3はRIAAの規格、サンプリング周波数は96kHzにしてみる。

Riaa96000input

この定数でEXCELで計算すると、

Zenkeisan0 ・・・式29

確認のためEXCELでシミュレーションしたところ、次のようになった。
理論値(dB)がRIAAの規格から求めた理論値で、1kHzを0に正規化したもの、正規化(dB)がシミュレーション結果、エラー(dB)が理論値との差分だ。


表3-a.伝達関数 式29のシミュレーション結果
20220226

10kHz で-0.397dB、20㎐で-1.479dBと、高域でゲイン誤差が若干大きめに出る。
ところで先人たちの計算をみると、次のようになっていることが多い。

Z2sakujo…式22

つまり、z-1の項をネグってz-2の項をz-1にシフトしている。これはすでに前の記事で式22として示したのと同じものだ。
確認のためEXCELでシミュレーションしたところ、次のようになった。

 

表3-b.伝達関数 式22のシミュレーション結果
Simulation0_20201025141901


20kHzが最大誤差で0.537dB となっているが、実用許容範囲だと思う。
シミュレーション結果を見ると、10k㎐以上の高域で、上で計算した表3-aよりも誤差が改善している。
また、この結果は下のインパルス不変変換法による係数でのシミュレーション結果にも近い。
●2022/2/26追記  Octaveによるシミュレーションを行ったところ、式22の係数ではゲイン誤差が抑えられるが、位相誤差が大きくなり、式29の係数では逆にゲイン誤差が大きくなるが位相誤差が抑えられているということがわかった。
RIAAフォノイコライザをIIRデジタルフィルタで実装する ⑤Octaveによる位相検証参照。

3-2.インパルス不変変換法による係数計算

もう一つ、一般的な変換方法として、インパルス不変変換法がある。

次のような1次ローパスフィルタ(積分回路)があるとする。

Sekibun_dentatsu・・・式30

これをインパルス不変変換法によってS-Z変換するには

Impulse_fuhen_henkan・・・式31

という置き換えを行う。ただしゲイン項(定数項)は省略してある。

今回の課題であるRIAA回路は次のように分解できる。

Riaa_bunkai・・・式32

つまりこれは積分回路2つとその逆数のハイブースト回路の積算であると見れば、式31を使って

Impulse_fuhen_sekisan・・・式33

と表現できる。

サンプリング周波数Fs=96kHz、T1=3180us、T2=318us、T3=75usとして計算すると、

Impulsefuhen_keisuu・・・式34

となる。

EXCELによるシミュレーション結果を表4に示す。

表4.インパルス不変変換、96kHzサンプリング
Impulse_fuhen_error

計算方法はシンプルだが、先ほどの双一次変換法よりも若干誤差が小さくなった。

 

☆次回は計算した係数やネットで入手した係数のシミュレーションなどについて書く予定です。

 

【参考文献】

はじめて学ぶディジタル・フィルタと高速フーリエ変換」 三上直樹 著

実践ディジタル・フィルタ設計入門」 岩田利王 著

 

【参考サイト】

デジタルフィルタードットコム-DIGITALFILTER.COM
理論や実装方法の記事、それに便利なソフトを配布してくれているありがたいサイトです。

Wyne Stegall HomePage
オーディオのページに、インパルス不変変換法によるRIAAイコライザのC言語での実装や、そのほかにもオーディオ技術関係の記事が多数あり、勉強になります。

Familie Beis HomePage
デジタル信号処理の詳しい説明があり、とても勉強になります。無料配布しているデジタルフィルタの設計解析ソフトが便利です。


| |

« RIAAフォノイコライザをIIRデジタルフィルタで実装する ①準備編 | トップページ | RIAAフォノイコライザをIIRデジタルフィルタで実装する ③シミュレーション編 »

コメント

コメントを書く



(ウェブ上には掲載しません)


コメントは記事投稿者が公開するまで表示されません。



« RIAAフォノイコライザをIIRデジタルフィルタで実装する ①準備編 | トップページ | RIAAフォノイコライザをIIRデジタルフィルタで実装する ③シミュレーション編 »