AudioUnitでFuzzFaceを作る。その１：ダイオードシミュレーション

FuzzFaceって知っていますか？

Jimi Hendrixも使っていたりする超有名なエフェクターですね。
これが欲しいのですが、お金がないです。
自作するにも環境がないです。あとパーツを買う元気もないです。

そこでAudioUnitでこれを実装しようという話です。プログラム上で組めればタダですから。
まあニッチな内容なので誰も興味ないでしょうが、自分用のメモとしてここに残しておきます。

まず最初に回路を確認します。

f:id:tl_muk99:20160617171115p:plain こちらから引用

シミュレートするにあたって一番重要になってくるのがトランジスタです。今回はこのトランジスタをどうやってシミュレートするかが肝になります。
トランジスタのシミュレートに関しては電子回路シミュレーション(現代非線形科学シリーズ)を参考にしながら設計していきます。

まず基本的な回路シミュレート方針として、FazzFaceの回路を回路方程式に直してそれを解くことでシミュレーションするという形です。

しかしここで問題があり、トランジスタはそのまま回路方程式に突っ込むことはできません。
そこでトランジスタの等価モデルに直してから回路方程式に突っ込んでいく必要があります。

・バイポーラトランジスタ等価モデル

FazzFaceで使われているものはバイポーラトランジスタになるので、それに合う等価モデルを使う必要があります。そして先ほどの参考文献に挙げられている等価モデルは以下の二つです。

Ebers-Mollモデル
Gummel-Poonモデル

今回はこのEbers-Mollモデルを採用してシミュレートを行います。しかしどうやら厳密なモデルが要求される場合はEbers-Mollモデルを改良したGummel-Poonモデルを使うらしいですが……今回は勘弁してください、まずは実装まで漕ぎ着けましょう。

では実際のEbers-Mollモデルの回路は以下のとおりです。
f:id:tl_muk99:20160623145320p:plain wikipediaより引用
見たとおりトランジスタ一個を二つのダイオードと二つの電流源で表現しています。この二つは回路方程式に組み込めるので一安心です。

・ダイオードシミュレーション(交流-直流変換回路)

まあ、自分のやったことを残すという意味も込めてダイオードのシミュレーションにも少し触れていきます。
ここでシミュレートする内容は交流-直流変換回路です。

すごい簡単な回路になります。ではこれを順繰り計算していきましょう。

まずこの回路から数式を洗い出していきます。

${ \displaystyle { i }_{ D } = { i }_{ C }+{ i }_{ R } \cdots \left( 1 \right) \\ E\left( t \right) ={ V }_{ D }+V \cdots \left( 2 \right) \\ I_{ D }=I_{ S }\left( { e }^{ \frac { { V }_{ D } }{ n{ V }_{ t } } }-1 \right) \cdots \left( 3 \right) \\ V=R{ I }_{ R } \cdots \left( 4 \right) \\ { I }_{ C }=C\frac { dV }{ dt } \cdots \left( 5 \right) }$

洗い出し終わったらどんどんこの数式を解いていきます。
(1),(5)より
${ \displaystyle \frac { dV }{ dt } =\frac { { I }_{ C } }{ C } \\ \frac { dV }{ dt } =\frac { -{ I }_{ R }+{ I }_{ D } }{ C } \cdots \left( 6 \right) }$

(3),(4),(6)より
${ \displaystyle \frac { dV }{ dt } = -\frac { V }{ RC } +\frac { { I }_{ S } }{ C } \left( { e }^{ \frac { E\left( t \right) -V }{ n{ V }_{ t } } }-1 \right) \\ F\left( V,t \right) = -\frac { V }{ RC } +\frac { { I }_{ S } }{ C } \left( { e }^{ \frac { E\left( t \right) -V }{ n{ V }_{ t } } }-1 \right) }$

ここで前進オイラー法を使うことで、
${ \displaystyle V\left( t+h \right) =V\left( t \right) +hF\left( V\left( t \right) ,t \right) }$

今回与える必要のあるパラメータはステップ幅 $h$ と抵抗値 $R$ 、静電容量 $C$ 、更にトランジスタのパラメータとして ${ I }_{ S }$ 、 $n{ V }_{ t }$ を用意する必要があります。今回使うのは以下のとおり。

パラメータ	値
$h$	$1/44100$
$R$	$100kΩ$
$C$	$1nF$
${ I }_{ S }$	$16.8\times { 10 }^{ -9 }$
$n{ V }_{ t }$	$55.5\times { 10 }^{ -3 }$

ここまできたらあとはこの数式をプログラミングするだけです。ではpythonで書いてみましょう。

# -*- coding: utf-8 -*-
import math
import numpy as np
import matplotlib.pyplot as plt

h = 1.0/44100.0 #StepWidth (Reciprical of SamplingRate)
R = 100.0 * (10**3) #Resistance 100kΩ
C = 1.0 * (10**(-6)) #Capacitance 1nF

#Parameter of Diode (Is, nVt)
Is = 16.8 * (10**(-9))
nVt = 55.5 * (10**(-3))

x = np.arange(0, 0.1, h)
f = 220
E = 3 * np.sin(2.0 * np.pi * f * x) #InputWave 3*sin(2πf)
V = np.zeros(len(E)) #OutputWave

for t in range(1, len(E)):
	V[t] = V[t-1] + h*(-V[t-1]/(R+C) + (Is/C)*(math.exp((E[t-1]-V[t-1])/nVt) - 1.0))

plt.xlabel("Time[sec]")
plt.ylabel("Amplitude[V]")
plt.ylim(-3.1, 3.1)
plt.plot(x, E, linewidth=2)
plt.plot(x, V, linewidth=2)
plt.show()