【AdC2015】物理数学：ベッセル関数とFM音源
アドベントカレンダーとは以下略
ロボット技術研究会でもアドベントカレンダーをやっています。そちらの方も合わせてご利用くださいませ。titech-ssr.blog.jp
今日のお題
FM音源を数学的に解析します。具体的には、物理数学でよく出てくる特殊関数の一つ、(第一種)ベッセル関数を使って、最も単純な構成のアルゴリズムで生成された波形が倍音を持つことを示します。
なぜFM音源なのかといえば…
- 個人的に興味があるから
shitaro-happy-physics.hatenablog.jp
shitaro-happy-physics.hatenablog.jp
- この前行われたロボット技術研究会後期研究報告会でFM音源が紹介されたから
です。個人的に、そしてサークル的にホットな話題だから取り上げてみました。
FM音源と倍音
FM音源とは
FM音源とは、周波数変調(Frequency Modulation)を用いた*1音作りのことです。簡単に言うと、元の波形を別の波形で揺さぶって変化させて出力する音源です。主にYAMAHAがFM音源モジュールを作成しています。
YAMAHAでは、波形を生成する回路のことをオペレータと呼びます。オペレータの内、変調されるもの(元の波形)をキャリア、変調するもの(揺さぶる波形)をモジュレータと呼びます。
FM音源ではいくつかのモジュレータを組み合わせて、キャリアを通して波形を出力することで音を生成します。いくつかのモジュレータとキャリアの組み合わせのことをアルゴリズムと呼びます。
以下の議論では、簡単のため、使用するアルゴリズムは正弦波を出力するモジュレータ１つと正弦波を出力するキャリアをつなげた単純なものとします。
FM音源の数式
今回のアルゴリズムによって生成される波形は次のように表せます。
\begin{align*}
s(t) = a_c \sin(\omega_c t + a_m \sin(\omega_m t))
\end{align*}
振幅\( a_c\)と\( a_m \)はそれぞれキャリア振幅、モジュレータ振幅と呼ばれるパラメータです。角振動数\(\omega_c \)と\( \omega_m\)はそれぞれキャリア振動数、モジュレータ振動数と呼ばれるパラメータです。
特に、モジュレータをいじるだけで倍音が作れるので、初期のパソコンや携帯電話、ゲーム機などのスペック的に余裕がないものに多く使われていたそうです。今回のテーマは、本当に倍音が現れるのか、どう現れるのかを物理数学の知識で見ていきたいと思います。
ベッセル関数(Bessel function)
実はFM音源の数式はベッセル関数で展開することができます。ここではそのベッセル関数の導入をします。
ベッセル微分方程式
ベッセル微分方程式は次の微分方程式で表されます。
\begin{align*}
x^2 \frac{\mathrm{d}^2 y}{\mathrm{d} x^2}
+
x \frac{\mathrm{d} y}{\mathrm{d} x}
+
(x^2 - n^2)y
=
0
\end{align*}
例えば円柱座標系におけるラプラス方程式を変数分離したときの、動径方向に関する微分方程式がこのベッセル微分方程式になります。ベッセル微分方程式は、太鼓の膜振動や円柱型導波管の電磁波の伝播などでよく見かける有名な微分方程式なので一度は解いたことがあると思います。
この微分方程式の解の一つは第一種ベッセル関数と呼ばれる特殊関数で与えられます。
\begin{align*}
J_n (x) = \left( \frac{x}{2} \right)^n \sum_{k=0}^{\infty} \frac{(-1)^k}{k!\Gamma(n+k+1)}\left( \frac{x}{2} \right)^{2k}
\end{align*}
第一種ベッセル関数は振動しながら減衰していく振る舞いをします。
天下り的ですが、この第一種ベッセル関数の母関数表示は次で与えられます。
\begin{align*}
\exp \left( \frac{x(z-z^{-1})}{2} \right) = \sum_{n=-\infty}^{\infty} J_n (x) z^n
\end{align*}
第一種ベッセル関数の母関数表示で\(z=e^{iy}\)とおくと次の式が得られます。
\begin{align*}
e^{ix \sin(y)} = \sum_{n=-\infty}^{\infty} J_n (x) e^{inx}
\end{align*}
この式の両辺の実部と虚部を比較すると次の関係式が得られます。
\begin{align*}
\cos(x \sin(y)) &= \sum_{n=-\infty}^{n=\infty} J_n (x) \cos (ny) \\
\sin(x \sin(y)) &= \sum_{n=-\infty}^{n=\infty} J_n (x) \sin (ny)
\end{align*}
注目すべきは、\(\cos(\sin)\)と\(\sin(\sin)\)が第一種ベッセル関数と三角関数の積の級数で表されているところです。このように、初等関数のいくつかはベッセル関数の級数で表されることが知られています。
ベッセル関数でFM音源の倍音に迫る
三角関数の加法定理を使って最初のsinを展開すると\(\cos(\sin)\)と\(\sin(\sin)\)の項が現れます。
\begin{align*}
s(t) = a_c (\sin(\omega_c t)\cos(a_m \sin(\omega_m t)) + \cos(\omega_c t)\sin(a_m \sin(\omega_m t)) )
\end{align*}
ところで\(\cos(\sin)\)と\(\sin(\sin)\)は第一種ベッセル関数と三角関数の積の級数で表せるのでした。
\begin{align*}
\cos(a_m \sin(\omega_m t)) &= \sum_{n=-\infty}^{n=\infty} J_n (a_m) \cos (n \omega_m t) \\
\sin(a_m \sin(\omega_m t)) &= \sum_{n=-\infty}^{n=\infty} J_n (a_m) \sin (n \omega_m t)
\end{align*}
よって、波形は次のように級数でまとめることができます。
\begin{align*}
&s(t)\\
=&
a_c
\sum_{n=-\infty}^{n=\infty} J_n (a_m) (\sin(\omega_c t)\cos(n \omega_m t) + \cos(\omega_c t)\sin(n \omega_m t)) \\
=&
a_c
\sum_{n=-\infty}^{n=\infty} J_n (a_m) \sin ( (\omega_c + n \omega_m)t )
\end{align*}
つまり、FM音源の波形は「振幅は第一種ベッセル関数」、「振動は正弦波」である波の重ねあわせで与えられることが分かります。しかも、振動は\(\sin ( (\omega_c + n \omega_m)t )\)であることから、基音(振動数\(\omega_c\))を中心に\(\omega_c \pm \omega_m,\ \omega_c \pm 2\omega_m, \ldots\)のピークが現れることが分かります。よって、FM音源の波形には倍音が現れることが分かりました。
モジュレータでキャリア*2を揺さぶるだけで複数の倍音が現れるので、制御パラメータが少なくてもいろんな音を表現できるのですね！
└（՞ةڼ◔）」＜Besselを感じるゥ〜！
