電気回路のシミュレーション（Julia Advent Calendar 2022）

Julia Advent Calendar 2022の18日目の記事です。毎年Juliaで音に関係したことを書いていて、これまでにこのブログで書いたアドカレ記事は以下の通りです。

これまでは音の分析についてばかり書いていました。今年は音を扱うのは変わりませんが、少し毛色を変えて音の加工をする電気回路のシミュレーションについて書いてみます。

はじめに
ACMEのチュートリアル
この回路に楽音を通してみる

はじめに

さて、いまからさかのぼること約25年、Propellerhead社からReBirth RB-338が登場したときにはとても驚きました。RB-338は、RolandのTR-808とTB-303という2種類の音源（リズムマシンとシンセサイザ）を再現したソフトウェア・シンセサイザです。すでにYamaha VL-1やVP1のような物理モデル音源もあるにはありましたが、高い演算パワーが必要なのでハードウェアで実現するのが精一杯という印象でした（そしてとても高価でした）。それがパソコンの速度向上のおかげで、ソフトウェアでアナログシンセを実現できる時代になったのでした。

いま考えればアナログシンセの電子回路をそのままシミュレーションしているわけではないとは思うものの、当時Verilog HDLなどを触れる環境にあったので、RB-338は電子回路をリアルタイムでシミュレーションしているのではないかと思って胸を熱くしていました。ゆくゆくは自分もそういったアナログシンセを作ってみたいと夢想したりもしましたが、現在はさらにコンピュータが高速化したので、リアルタイムで電子回路をシミュレーションできるツールも登場しています。そういったものの例としては、LiveSPICEという、SPICEで作ったモデルがVST プラグインとしてDAWから呼び出せるオープンソースのツールがあります。演奏音のリアルタイム加工に使っても大きな違和感は生じないくらいの低遅延で動作します。

そんなわけで今回は、Juliaで使える回路シミュレータを紹介します。リアルタイム処理は難しいですが、入出力の関係を見たりするのには十分かと思います。

楽器用のエフェクタ（音の加工をする装置）には、山彦やホールのような響きをつけたり、シュワシュワ・ワオワオ言うような効果をつけたり、ギュイーンと迫力のある音にしたりと、さまざまな種類のものがあります。今日は、回路シミュレータACME.jl（Analog Circuit Modeling and Emulation for Julia）のチュートリアルに沿って、入力音にやわらかな歪み（ひずみ）を加えつつ高音域を削るような電気回路を見ていきます。

ACMEのチュートリアル

ACMEでは電子回路を以下のコードのように記述します。まず

voltagesourceは電源（単位：ボルト）
resistorは抵抗器（単位：オーム）
capacitorはコンデンサ（単位：ファラド）
diodeはダイオード（isで飽和電流を指定）

という具合に、それぞれパーツを準備します。次に、それらのパーツがどのようにつながっているかを⟷記号で示してあげます。抵抗やコンデンサは電流の方向は関係ないので[1]と[2]で両端を示しますが、電源やダイオードは方向が重要なので[+]と[-]で指定します（ダイオードの場合は[+]がアノード、[-]がカソードです）。

using ACME

circ = @circuit begin
    j_in = voltagesource()
    r1 = resistor(1e3)
    c1 = capacitor(47e-9)
    d1 = diode(is=1e-15)
    d2 = diode(is=1.8e-15)
    j_out = voltageprobe()
    j_in[+] ⟷ r1[1]
    j_in[-] ⟷ gnd
    r1[2] ⟷ c1[1] ⟷ d1[+] ⟷ d2[-] ⟷ j_out[+]
    gnd ⟷ c1[2] ⟷ d1[-] ⟷ d2[+] ⟷ j_out[-]
end

これを回路図にしてみると以下のようになるでしょうか。上のコードとの対応がわかりやすいように[1][2][+][-]を書き入れてあります。

この回路を音に対して実行できるモデルにします。音楽CDと同じく1秒間に44100回の処理ができるように、標本化周期を1/44100にしました。

model = DiscreteModel(circ, 1/44100)

1000 Hzの正弦波信号を作って、このモデルに突っ込んでみましょう。また、その結果をグラフ化して確認してみます。

x = sin.(2π * 1000 / 44100*(0:44099)');
y = run!(model, x);

using Plots
pyplot();
t = [0:199] ./ 44100 .* 1000;
plot(t, x[1:200], xlabel="Time (ms)", ylabel="Amplitude",
  label="Input (1000 Hz puretone)");
plot!(t, y[1:200], label="Output")

水色が入力の1000 Hz正弦波、橙色が出力波形です。ダイオードによるクリッピングで上下がつぶされて平坦になっていることがわかります。また、波形が時間方向に少しずれて、ピークを中心に左右対称ではなくなっています。

これを周波数分析してみましょう。ここで使っているplot_spectrumは以前作った自作関数です（関連記事：パワースペクトル計算の2つの方法 (Julia Advent Calendar 2019) - 丸井綜研）。

plot_spectrum(vec(y), 44100)
savefig("spectrum_1kHz_sinewave.png")

1000 Hzのところが一番強く出ていますが、その整数倍である2000, 3000, 4000, ... Hzなどにもわずかな信号が入っており、信号が少々歪んでしまうのが分かります。特に奇数倍の3000, 5000, 7000, ... Hzが比較的強く出ているのが、同じくらいの飽和電流のダイオードを対称に並べていることの効果です。（全く同じダイオードを対称配置にすると奇数次の成分だけが出るようになります）

この回路に楽音を通してみる

録音した楽器音を回路に通すのも簡単です。ギターでも歌声でも良いですが、sound.wavというファイルにステレオで楽音が入っているとします。左右チャンネル別々に回路モデルを適用してファイルに保存するというのが次のコードです。

using FileIO: load, save
import LibSndFile
using SampledSignals

snd = load("sound.wav")
model = DiscreteModel(circ, 1/snd.samplerate);
outputL = run!(model, map.(Float64, snd.data[:,1]'));
outputR = run!(model, map.(Float64, snd.data[:,2]'));
snd_out = SampleBuf(map.(PCM16Sample, [outputL' outputR']), snd.samplerate)
save("output.wav", snd_out)