EXCELでスピーカーシミュレーション（７） - シミュレーションのトビラ

７．EXCELでシミュレーションする

７．EXCELでシミュレーションする

　実際にEXCELで具体的な数値を入力して出力音圧の周波数特性をシミュレーションしてみます。ちなみに、各定数の具体的な数値は、ドライバーユニットのカタログに載っていたり、寸法を測って求めることができるものもありますが、計算や実験をしないとわからないものもあります。そのあたりは別途説明したいと思いますので、ここでは機械的に入力していきます。
　音響インピーダンスには3種類あると書きましたが、ここでは、「機械インピーダンス」を用いて計算していきます。

７－１．ドライバーユニットの音圧

　まず、ドライバーユニットの特性を求めます。
　最初に、準備段階として、周波数を縦方向に20Hzから20kHzまで列記します(かなり面倒です)。定数とかを入力する欄を空けておきたいので、H列当たりに入力します。セルH1は見出しとして「f」とでもしておき、数値はセルH2から入力します。周波数の間隔は、JIS Z 8601の標準数列を使うのが良いと思います。測定機もこの数列の間隔で測定しているのが多いので、実測値と比較するときに便利です。今回はR80で並べてみます。

　次に、隣のI列に $j \omega$ を入力します。 $j$ は虚数単位です。EXCELは複素数を扱えます。複素数の実部を0にして虚部を $\omega$ とすれば $j \omega$ となります。
　具体的には複素数を定義する関数のCOMPLEX()を用いて、セルI2に
=COMPLEX(0,2*PI()*H2,"j")

と入力します。最初の引数は実部ですので0、次の引数は虚部ですので $\omega=2 \pi f$ 、最後の引数は虚数単位の選択で" $j$ "を入力します。そして、20kHzまでオートフィルします。ちなみに、20Hzでは、125.663706143592jとなっていればOKです。
　次は、力 $F$ を入力します。力 $F$ は、
$F=\dfrac{BLE}{Z+Z_0}=\dfrac{BLE}{R_c+j \omega L_c}$
でしたので、まず定数を入力しておきます。数値はとりあえず下図のようにしておきましょう。とある8オームのドライバーユニットの実測値や換算値です。アンプの電圧は8オームユニットに対して1Wの入力になるときの電圧です。

　これらを使って、J列に $F$ を入力します。 $j \omega$ を含むので、複素数の計算方法で入力する必要があります。実数のときとは違い特殊な関数で入力します。

演算	実数	複素数
和	A1+B1	IMSUM(A1,B1)
差	A1-B1	IMSUB(A1,B1)
乗	A1*B1	IMPRODUCT(A1,B1)
除	A1/B1	IMDIV(A1,B1)

　ということで、セルJ2には、
=IMDIV($B$2*$B$3*$B$4,IMSUM($B$5,IMPRODUCT(I2,$B$6)))
と入力してオートフィルします。

定数のところは$で絶対参照にしておかないと、オートフィルのときに参照が狂うのでご注意。20Hzでは6.60631418561728-0.0598911188399136jとなっていればOKです。
　続いて各インピーダンスを入力していきます。まず、定数を入力します。こんなかんじです。

　質量 $m$ は振動板だけでなく、エッジとダンパーの半分くらいとかボビンとかボイスコイルとか振動するものすべての重さを足します。でもそれだとぶっ壊さないと測れないので、実際には実験値から計算で求めます。スチフネス $s$ と機械抵抗 $r$ も実験から求められます。
　放射質量 $m_a$ と放射抵抗 $ｒ_a$ は本当は周波数の関数なんですが、ka<<1の周波数では近似式として $m_a=\frac{8}{3}\rho_0(\frac{S_d}{\pi})^{\frac{3}{2}}$ 、 $r_a=\rho_0c\frac{k^2S_d^2}{2\pi}$ を用いることができます。ただし、感覚的にはka<0.5くらいの周波数までならばこの近似式でだいたい合うのですが、それ以上の周波数だと $m_a$ は急激に低下し、 $r_a$ は急激に上昇するので、合わなくなります。
　空気密度 $\rho_0$ は1.293、振動板の直径はエッジの半分くらいまでを足して212mmくらいなので、振動板面積を $\frac{\pi d^4}{4}$ で計算しています。
　周波数に依存するインピーダンスは、 $\frac{BL^2}{Z+Z_0}$ と $m$ と $m_a$ と $s$ です。
　 $\frac{BL^2}{Z+Z_0}$ は、セルK2に、
=IMDIV( ($B$2*$B$3)^2,IMSUM( $B$5,IMPRODUCT(I2,$B$6) ) )

と入力し、 $m$ と $m_a$ は足し合わて、セルL2に、
=IMPRODUCT(I2,$B$7+$B$10)

と入力し、 $s$ は、セルM2に、
=IMDIV($B$8,I2)

と入力して、それぞれオートフィルします。
　次に、粒子速度 $V$ を計算します。計算式は、
$V =\dfrac{F}{Z_{total}}$
でしたので、セルN2に、
=IMDIV( J2,IMSUM(K2,L2,M2,$B$9+$B$11) )

と入力します。 $r$ と $r_a$ を絶対参照で加えるのを忘れないようにしましょう。
　最後に音圧を求めます。計算式は、
$p= j \omega \rho_0 \dfrac{S_dV}{2 \pi x} e^{-jkx}$
でした。
　 $x$ は測定距離で、周波数特性は1mで測定しますから1にします。
　 $k$ は $\frac{\omega}{c}$ ですので、音速340m/sを使って、セルO2に、
=2*PI()*H2/340

と入力してオートフィルします。
　そして、音圧の式として、セルP2に、
=IMPRODUCT( I2,1.293,IMDIV( IMPRODUCT( $B$12,N2),2*PI()*1 ),IMEXP(COMPLEX(0,-O2*1,"j") ) )
と入力してオートフィルします。

　 $e^a$ の $a$ の部分は、実数の場合はEXP()関数を使いますが、複素数の場合はIMEXP()となります。
　これで音圧の周波数特性が求まりました。しかし、これで求まった音圧 $p$ は、複素数で表現されています。これでは分からないので、絶対値をとって実数にします。絶対値はIMABS()関数で計算できますので、セルQ2に、
=IMABS(P2)
と入力してオートフィルします。これで音圧の値が求まりました。単位はPaです。一般的にはデシベルの方が慣れていると思いますので、20uPaを0dBとして計算します。セルR2に、
=20*LOG10(Q2/(20*10^-6))
と入力してオートフィルします。
　そして、H列の周波数とR列の音圧をグラフにとると、このようになります。

　また、P列の複素数の音圧の偏角をとると、位相遅れを計算できます。偏角はIMARGUMENT()関数で求めることができますが、結果はラジアンで出てきます。度にしたければ、PI()/180で割ます。セルS2に、
=IMARGUMENT(P2)/(PI()/180)
と入力してオートフィルします。DEGREES()関数を使った方がはやいかも。
　そして、H列の周波数とS列の位相とをグラフにとって、第二軸に加えると、このようになります。

　-180度で+180度に折り返しているのは、位相が一周して更に遅れていることを表します。これは距離 $x$ を1mに設定して、その時点の音圧を求めているので、距離による位相遅れが周波数に比例して大きくなることに起因します。ドライバーユニットの振動系の純粋な位相遅れを見たければ、 $e^{-jkx}$ の部分を外して計算すればOKです。 $e^{-jkx}$ は絶対値には影響を与えず、距離による位相遅れにのみ影響を与える項なので、このようになるはずです。

　位相遅れはいろんな情報を持っているのですが、例えば位相変化を角周波数で微分すると群遅延という値が得られます。群遅延はザックリ言うと、周波数のずれに対して音の発生が遅れる度合いがどのくらいかを示したもので、要は群遅延が大きい周波数の音は遅れて聞こえるってことです。定義としては、
$T_{gd}=-\dfrac{d\phi}{d\omega}$ となります。単位は(sec)で、 $\phi$ は位相角です。
　群遅延も求めておきましょう。セルT2は空白にして、セルT3に
=-(IMARGUMENT(P3)-IMARGUMENT(P2))/(2*PI()*(H3-H2))
と入力し、オートフィルします。グラフにとると、このような感じです。

　折り返しのところで計算がおかしくなってしまうので、距離遅延の $e^{-jkx}$ を抜いて計算すると、

こんな感じになって見やすくなります。
　ともあれ、一応の目的を達成できましたので、定数の値をいろいろいじって遊んでみてください。グラフの曲線が変わると思います。

７－２．密閉型エンクロージャの音圧

　次は密閉型エンクロージャに入れたときの音圧をシミュレーションしてみます。
　といっても、実際にはドライバーユニットの等価回路にキャビネット内の空気のスチフネス $s_c$ と機械抵抗 $r_c$ が直列に入るだけです。
　ドライバーユニットの計算シートに続けて入力していきます。まず、定数を下図のように入力します。

　キャビネットのサイズは61.4cm×32.7×33.1cmで。細かいですが、実際作ったのがこれくらいだったので。キャビネット内空気の機械抵抗は適当です。
　次に、キャビネットのスチフネスを入力します。空洞内のスチフネスは、
$s_{Am}=\dfrac{\gamma P_0}{V_c}Sd^2$
です。比熱比 $\gamma$ は1.4で、大気圧 $P_0$ は101,300Paです。そして、セルU2に、
=IMDIV( $B$14,I2 )
と入力してオートフィルします。
　あとは粒子速度 $V$ を求めます。セルV2に、
=IMDIV( J2,IMSUM( K2,L2,M2,U2,$B$9+$B$11+$B$15 ) )
と入力します。 $r_c$ を絶対参照で加えるのを忘れないようにしましょう。
　そして、音圧 $p$ は、セルW2に、
=IMPRODUCT( I2,1.293,IMDIV( IMPRODUCT( $B$12,V2 ),2*PI()*1),IMEXP( COMPLEX(0,-O2*1,"j" ) ) )
と入力し、X列とY列には絶対値とデシベル値を、Z列とAA列には位相と群遅延を入力します。
　グラフはこんな感じです。

　キャビネット内の空気で振動板の動きが制動されるので、ほんの少し低域が出なくなります。その代わり群遅延は若干改善します。

７－３．位相反転型エンクロージャの音圧

　次は位相反転型エンクロージャに入れたときの音圧をシミュレーションしてみます。追加する定数は下図の通りです。

　ポートの直径は15.6cmで長さは10cmとしておきます。ポートの開口面積は $\frac{\pi D_p^2}{4}$ で計算します。
　ポート内の空気の質量は、ポートの体積に空気の密度をかければいいのですが、ポートのような筒形状の場合、開口部周囲の空気を若干巻き込んで動くらしく、その分を長さで補正すると良いようです。長さの追加分は、直径×0.73だそうですので、 $\rho_0 S_p (L_p+0.73*D_p)$ という感じです。
　ポート内の空気の機械抵抗ですが、これは筒の中を動く流体の粘性抵抗となるので、一応近似式があり、 $8 \pi \mu_0 L_p$ です。 $\mu_0$ は空気の粘度で、1.822×10 $^{-5}$ Pa・Sです。
　ポートの放射質量は、振動板の放射質量の $S_d$ を $S_p$ に変えるだけです。放射抵抗は適当です。
　最後に、振動板面積 $S_d$ とポート開口面積 $S_p$ の比の二乗を計算しておきます。ポートの各インピーダンスは、この係数を乗じることになります。この点について少し説明します。
　位相反転型はキャビネット内の空気を介して、ドライバーユニットの振動系とポート内の振動系が結合した構造をしています。このように、異なる質量を空気バネで結合する場合、等価回路としてはトランスが介在する形となります。

　というのは、ポートの開口面積は振動板よりも小さいので、振動板によってキャビネット内の空気に与えた力がポートに加わるときには、小さい面積に大きな圧力がかかるため、ポートに加わる力は振動板とポート開口部の面積比の分だけ大きくなります。そうなると、機械回路としては、変圧比＝振動板面積/ポート開口面積のトランスが入った状態と等価に考えることができます。そして、ドライバーユニット側と一次側、ポート側を二次側と考えたときに、機械回路を一次側換算すると、トランスの計算のように、ポートの各インピーダンスは変圧比の二乗が掛かった状態になるのです。これが、振動板面積 $S_d$ とポート開口面積 $S_p$ の比の二乗を乗じることになる理由です。
　

　ちなみに、二乗を乗じるのは等価回路の計算に機械インピーダンスを用いているからです。比音響インピーダンスを用いる場合、比音響インピーダンスはもともと断面積で除しているので、ポートの各インピーダンスは振動板面積 $S_d$ とポート開口面積 $S_p$ の比 $\frac{S_d}{S_p}$ を乗じます。また、音響インピーダンスを用いる場合、音響インピーダンスは機械インピーダンスに対して断面積の二乗を除しているので、ポートの各インピーダンスは何も乗じずそのまま計算できます。このように、3種類の音響インピーダンスのどれを用いるかによって等価回路の考え方が変わるので要注意です。
　さて、EXCELに戻ります。位相反転型エンクロージャでは、振動板の振動速度 $V_1$ と、ポートの振動速度 $V_2$ を求める必要がありました。式は下記です。
$V_1 = \dfrac{z_2+z_3}{(z_1+z_2)(z_2+z_3)-z_2^2}F$
$V_2 = \dfrac{z_2}{(z_1+z_2)(z_2+z_3)-z_2^2}F$
　EXCELの入力作業を分かりやすくするために、式に合わせて各インピーダンスを $z_n$ の形にしましょう。まず、セルAC1を $z_1$ として、セルAC2に、
=IMSUM( K2,L2,M2,$B$9+$B$11 )
と入力し、セルAD1を $z_2$ として、セルAD2に、
=IMSUM(U2,$B$15)
と入力し、セルAE1を $z_3$ として、セルAE2に、
=IMPRODUCT($B$23,IMSUM(AB2,$B$20+$B$22))
と入力します。 $z_3$ だけは、 $(\frac{S_d}{S_p})^2$ をIMPRODUCTで掛ける必要があることに注意です。
　そして、列AF、AGをそれぞれ、V1、V2とします。セルAF2には、
=IMDIV(IMPRODUCT(IMSUM(AD2,AE2),J2),IMSUB(IMPRODUCT(IMSUM(AC2,AD2),IMSUM(AD2,AE2)),IMPOWER(AD2,2)))
と入力し、セルAG2には、
=IMDIV(IMPRODUCT(AD2,J2),IMSUB(IMPRODUCT(IMSUM(AC2,AD2),IMSUM(AD2,AE2)),IMPOWER(AD2,2)))
と入力します。実数の指数は^記号で入力できましたが、複素数の指数は、IMPOWER()という関数を用います。
　音圧は、トータルだけ分かればいいのですが、せっかくなので振動板から出ている音圧 $p_1$ と、ポートから出ている音圧 $p_2$ も個別に出してみます。
　AH、AI、AJの列はそれぞれ $p_1$ の音圧、デシベル値、位相とします。セルAH2には、
=IMPRODUCT(I2,1.293,IMDIV(IMPRODUCT($B$12,AF2),2*PI()*1),IMEXP(COMPLEX(0,-O2*1,"j")))
と入力し、セルAI2には、
=20*LOG10(IMABS(AI2)/(20*10^-6))
と入力し、セルAJ2には、
=IMARGUMENT(AI2)*(180/PI())
と入力します。
　あとは同じ要領で、列AK、AL、AMにはポートの音圧、デシベル値、位相を入力し、列AN、AO、APには、トータルの音圧、デシベル値、位相を入力します。群遅延は列AQにトータルのものだけ求めておきましょう。
　こんな感じになっていればOKです。

　各グラフはこのようになります。

　ところで、EXCELで位相角のアンラップをする関数ってないんですね。見にくい・・・IF分使って無理やりできんこともないでしょうけど、面倒。

７－４．ドロンコーン型エンクロージャの音圧

　最後にドロンコーン型エンクロージャに入れたときの音圧をシミュレーションしてみます。追加する定数は下図の通りです。

　ドロンコーン型は位相反転型に対してスチフネスが追加になるだけなので、計算式はほぼコピーで使えます。
　ドロンコーンの振動板の直径はドライバーユニットの振動板と同じでもいいのですが、少し小さめに20cmくらいにして、質量も25gくらいにします。スチフネスは2000くらいにしましょうか。機械抵抗は適当で。放射質量はポート開口面積をドロンコーンの振動板面積に入れ替える。変圧比も忘れずに。
　列ASには、ドロンコーンの質量と放射質量を合わせて、セルAR2に、
=IMPRODUCT(I2,$B$26+$B$29)
と入力し、列ATにはドロンコーンのスチフネスを入れます。セルAS2に、
=IMDIV($B$27,I2)
と入力します。等価回路では $z_3$ が変わるだけなので、列ATに $z_3'$ とでもして、セルAT2に
=IMPRODUCT($B$31,IMSUM(AR2,$B$28+$B$30,AS2))
と入力します。変圧比を掛けるのを忘れずに。
　あとは位相反転型と同じ要領で、列AU、AVに振動速度V1'、V2'を入力します。続いて列AW、AX、AYにはドライバーユニットの振動板による音圧、デシベル値、位相を入力し、列AZ、BA、BBにはドロンコーンの振動板による音圧、デシベル値、位相を入力し、列BF、BD、BE、BFにはトータルの音圧、デシベル値、位相、群遅延を入力します。
　こんな感じになっていればOKです。