● なんと、アナログ方式のミュージックシンセサイザーを SPICE でシミュレートしてみよう、という屈折した
試みです。 その過程で SPICE の使い方のヒントもいろいろ得られるだろうという目論見はあります。
まずは VCO, VCF, EG (ADSR) + VCA を作りました。
ピッチ (音程) とキー押し制御データはテキストファイル wm_pitch.txt, wm_key.txt から読み込み
ます (V1, V3)。 V2 は VCO 出力鋸波・矩形波のデューティーファクターを指定します。 B2 は
ノイズ源です。 V4 は VCF にビブラートをかけると共に、VCF のカットオフ周波数を VCO 出力
周波数 (基本波) の 2**(7/12) 倍にしています。 VCF のパラメーター Rq がフィルターの Q を
指定します。 エンビロープの ADSR も、それぞれパラメーターで指定します。 scl は VCA の振幅係数
です。 結果はファイル output.wav に書き込まれます。 (ADSR := Attack, Decay, Sustain,
Release)
synth.zip 中の wm_test.asc がトップレベルの回路です。
最初、wm_test.asc 中の VCA パラメーター scl (出力振幅指定) を 3 にしたところ、どうも信号が
クリップしているような音がします。 説明書によると
"The .wav analog to digital converter has a full scale range of -1 to
+1
Volt or Amp(ere)."
ですから、単なる設定ミスでした。 (この説明書の表現も、実に屈折しています。どこにアナログ信号があるのかな。)
● もう少しシミュレーション速度が上げられないものかと思って、まずは VCF に使っている 1pole OPamp を
lib/sub/ にある最も簡単そうな opamp.sub に変えてみました。 ところが予想に反して、かえって遅くなっ
てしまいました。 この理由は、1pole OPamp はコードモデルの OTA を使っているからだと思われます。
B 電圧源の Laplace 式で VCF を表現しようと試したのですが、カットオフ周波数の電圧制御は無理の
ようです。 他にディレーラインもしくは S/H (サンプル・アンド・ホールド回路) を使い、z 変換を利用
したフィルターも可能です。 しかし出力波形が不連続になるのが気になりますし、それほど速くもあり
ません。 結局、この用途には OP アンプを
使用したステート変数フィルターが最適のようです。 High-pass, Band-pass, Low-pass 出力が一挙に
得られるという利点もあります。 他の方法と言うと、いっそのこと LC フィルターを使ってしまう方法は
残っていますが。
z 変換を利用した IIR フィルターについては
http://groups.yahoo.com/group/LTspice/files/%20Tut/ 中の
Z-Transform%20Applications/
で紹介されています。 この例では S/H を自前で作っていますが、現在は A デバイスの S/H が
使えます。 synth.zip 中の z_filter/Draft1.asc をご覧下さい。
● 前記の例は４次のバターワースフィルターです。 ちょっと脇道にそれるようですが、s-z 変換を
利用して Q も自由に設定できる ２次の IIR LPF を構成してみました。 係数の具体的な計算方法は
www.staff.newcastle.ac.uk/oliver.hinton/
.. /Chapter5.pdf
を参照しました。 一応の知識があるものとして基礎部分をかなり飛ばして書いてありますから、
いきなり読んでもわからないと思います。 基礎部分は書籍等にゆずります。 双一次 (bilinear)
近似を使った直接構成法とも言われる方法です。
計算には数式処理プログラム MAXIMA を使いました。 このような紹介があります。
http://phe.phyas.aichi-edu.ac.jp/~cyamauch/maxima/
http://www.interq.or.jp/mars/cherry/maxima/
ダウンロード先はこちらです。 http://maxima.sourceforge.net/download.shtml
普通のプログラミング言語とは少し雰囲気が違います。 "=" は代入操作を表すのではなく、
あくまでも式の両辺が等しいという意味に使います。 つまり方程式の宣言になります。
代入操作は ":" で表します。
まずは例として示されている fc/fs = 1/10 の２次バターワース LPF の計算をたどってみます。
以下の説明中の青文字の部分は MAXIMA に与えたコマンドです。
fc/fs = 1/10 (fc := カットオフ周波数、fs := サンプリング周波数)
H(s) = 1 / (s^2 + sqrt(2)*s + 1)、これの周波数を変換します。
Ωc = 2π*fc/fs = 2π/10 = 0.6283
ωac = K * tan(Ωc/2) = 1 * tan(0.6283/2) = 0.3249 rad/s
H(s) = 1/( (1/0.3249)^2*s^2 + (sqrt(2)/0.3249)*s + 1 )
fpprec : 6 $
s : (1-x)/(1+x) $ この部分が bilinear
z-transformation の働きをします。 x を z^-1 の意味に使います。
r1 : 1/( (1/0.3249)^2*s^2 + (1.414/0.3249)*s + 1 ) $
pfeformat : true $
fortran(combine(expand(r1)));
ほんとうはもう少し通分してほしいのですが、しかたないので手で通分
します。 次のように求まりました。
1/(
5.12117*x**2 / (x**2+2*x+1)
-18.9466*x / (x**2+2*x+1)
+13.8254 / (x**2+2*x+1)
+(x**2+2*x+1) / (x**2+2*x+1)
)
これを x の次数別に整理して繁分数を普通の分数に直します。 分母の１次の係数の符号に注意。
分子： x**2 + 2*x + 1
分母： (5.12117 + 1)*x**2 - (18.9466 - 2)*x + (13.8254 + 1)
分母の 0 次の係数を 1.0 にするため、(13.8254 + 1) で除算します。
fortran(expand((x^2 + 2*x +
1)/(13.8254+1)));
fortran(expand(( (5.12117+1)*z^2 - (18.9466-2)*z + (13.8254+1)
)/(13.8254+1)));
これでようやく係数が求まりました。
分子： 0.0675*x**2 + 0.1349*x + 0.0675
分母： 0.4129*x**2 - 1.14308*x + 1.0
x を z^-1 に戻します。 丸め誤差等を無視すれば ../Chapter5.pdf に示されているの
と同じ値です。 MAXIMA の出力は次数の大きいほうから書いてあることに注意して下さい。
synth.zip 中の z_filter/Draft2.asc をご覧
下さい。 fs = 10kHz にしています。 fs/2 で出力が 0 になります。 fs が低いので、せい
ぜい 2*fc くらいまでが使える範囲です。
(注) 後日、rat() が繁分数の簡約計算に使えることに気づきました。 keepfloat : true; と設定しておくと
いいでしょう。 また、私の目的には fortran() を多用するより display2d : false; のほうが
向いています。
● 次に、欲張って fc と Q を可変にしました。
H(s) = 1 / (s^2 + (1/Q)*s + 1) として周波数変換をします。
Ωc = 2πfc/fs
ωac = K * tan(2πfc/fs/2) = 1 * tan(2πfc/fs/2) = tan(πfc/fs) rad/s
wac = tan(pi*fc/fs) と置きます。
H(s) = 1 / ( (s/wac)^2 + (1/Q)*(s/wac) + 1 )
分母分子に wac^2 を乗じます。
H(s) = wac^2 / (s^2 + (wac/Q)*s + wac^2)
s : (1-x)/(1+x) $
r1 : wac^2 / (s^2 + (wac/Q)*s +
wac^2) $
pfeformat : true $
combine(expand(r1));
2
wac
---------------------------------
2
x - 2 x + 1 wac - wac x 2
------------ + ----------- + wac
2 Q x + Q
x + 2 x + 1
分母に関して通分します。
(x^2 - 2*x + 1)*Q / (Q*(x^2+2*x+1))
(wac - wac*x)*(x+1) / (Q*(x^2+2*x+1))
wac^2*(Q*(x^2+2*x+1)) / (Q*(x^2+2*x+1))
次数ごとにまとめます。
fortran(combine(expand(
(x^2 - 2*x + 1)*Q
+(wac - wac*x)*(x+1)
+wac^2*(Q*(x^2+2*x+1)) )));
Q*wac**2*x**2 - wac*x**2 + Q*x**2
+2*Q*wac**2*x - 2*Q*x
+Q*wac**2 + wac + Q
fortran(factor(Q*wac**2*x**2 -
wac*x**2 +Q*x**2));
fortran(factor(+2*Q*wac**2*x -
2*Q*x));
fortran(factor(+Q*wac**2 + wac +
Q));
a2 : (Q*wac**2-wac+Q)*x**2 $
a1 : +2*Q*(wac-1)*(wac+1)*x $
a0 : +Q*wac**2+wac+Q $
(これが H(z) の分母になります。)
分母分子に Q*(x^2+2*x+1) を乗じます。
分子： wac**2*Q*(x^2+2*x+1)
fortran(expand(wac**2*Q*(x^2+2*x+1)));
b2 : Q*wac**2*x**2 $
b1 : +2*Q*wac**2*x $
b0 : +Q*wac**2 $
分母の 0 次の係数を 1.0 にします。
fortran(a1/a0) $ 2*Q*(wac-1)*(wac+1)*x/(Q*wac**2+wac+Q)
fortran(a2/a0) $ (Q*wac**2-wac+Q)*x**2/(Q*wac**2+wac+Q)
fortran(b0/a0) $ Q*wac**2/(Q*wac**2+wac+Q)
fortran(b1/a0) $ 2*Q*wac**2*x/(Q*wac**2+wac+Q)
fortran(b2/a0) $ Q*wac**2*x**2/(Q*wac**2+wac+Q)
x を z^-1 に戻します。 synth.zip 中の
z_filter/Draft3.asc をご覧下さい。 fs = 44.1kHz にしています。 10kHz くらいまでは充分使えそうです。
(注) 後日、rat() が繁分数の簡約計算に使えることに気づきました。 keepfloat : true; と設定しておくと
いいでしょう。 また、私の目的には fortran() を多用するより display2d : false; のほうが
向いています。
ちなみに一次 IIR LPF の係数の計算式をメモしておきます。 TS は遅延線または S/H による遅延時間です。
* .param fs=44.1k fc=1k
.param TS=1/fs wac=tan(pi*fc/fs)
.param B0=wac/(wac+1)
.param B1=wac/(wac+1)
.param A1=(wac-1)/(wac+1)
● Draft3 の回路ですと、パラメーターで fc と Q を設定できますが、電圧制御にはなって
いませんでした。 z_filter/Draft4.asc を見てください。 カットオフ周波数 fc や Q、そし
てフィルターの係数 A1, A2, B0 〜 B2 には全てノード電圧を割り当てました。 V2, V3 は、
本来なら外部から与えられる制御電圧を作っています。 B5 〜 B12 は回路図としては
描かれていませんが、SPICE で言うところの B 電圧源そのものです。
● 脱線はこのくらいにして、アナログ回路に話を戻します。 VCF を多重帰還回路で表現しま
した。 V2 が f0 を、V3 が Q を指定しています。 周波数 * Q^2 が大きくなると誤差がやや
増えますが、GBW 10MHz の OP アンプでも f0 = 10kHz, Q=10 くらいまでなら何とか使え
ます。 これは BPF です。 synth.zip 中の mf-bpf.asc をご覧下さい。
次は多重帰還 LPF です。 http://focus.ti.com/lit/an/slod006b/slod006b.pdf の
16.3.2.2 Multiple Feedback Topology の項を参考にして、数式処理プログラム MAXIMA も
使いました。 主な関係式は：
A(s) = -(R2/R1)/(1 + wc*C1*(R2 + R3 + R2*R3/R1)s +
wc^2*C1*C2*R2*R3*s^2)
A0 = -R2/R1
a1 = wc*C1*(R2 + R3 + R2*R3/R1)
b1 = wc^2*C1*C2*R2*R3
Q = sqrt(b1)/a1
2*pi*f0 = 1/sqrt(C1*C2*R2*R3)
途中経過は省略して SPICE 語に翻訳すると：
.param R1=10k f0=1k Q=4
.param A0=-1/Q R2=-A0*R1 k=Q/(Q+1) R3=k*R2
.param C1=Q/(2*pi*f0*(k*Q+Q+k)*R1)
.param C2=C1*(k*Q+Q+k)**2/k
つまり R1, f0, Q を指定すれば R2, R3, C1, C2 が求まります。 (A0=-1/Q のほかに C2/C1
比を最小にする、という
条件をつけました。)
グラフの (V(o)/V(n))*(V(f)/V(in)) は回路のループ利得を表し、この例では
十分に余裕があります。
Q を大きく設定すると R2/R1 比が小さくなって U1 の負荷が重くなりすぎる
ことがあります。 R3/R2 比はそれほど過激な値にはなりません。 synth.zip
中の mf-lpf.asc をご覧下さい。
多重帰還 LPF は Q を高く設定すると C2/C1 比がとても大きくなって、あまりよくありません。
もどる
 |