平面波と球面波とは？〜pythonによる可視化まで〜【音響理論基礎】

今回は音響学で最もよく扱う波である平面波と球面波について概説した後，pythonのanimationで可視化する方法を説明します。

3次元波動方程式からスタートしますので，前回の記事を適宜ご参照ください。

【音響理論基礎】波動方程式の導出

それでは本題です。

平面波と球面波の導出

ここでは平面波と球面波の意味と導出について簡単にお話します。

平面波・球面波はどちらも波動方程式の解となっていますので波動方程式から出発して導出していきます。

準備

今回は時間領域のフーリエ変換を使って導出していきますが，波の伝搬表現の都合上フーリエ変換及び逆フーリエ変換を以下のように定義します。

\begin{align}
f(t) &= \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{-i\omega t} d\omega \notag \\
F(\omega) &= \int_{-\infty}^{\infty} f(t) e^{i\omega t} dt \notag
\end{align}

ここで以前の記事とは符号が違うことに注意して下さい。

※実は符号はどちらでも良いのですが，個人的にこちらの方が慣れているのでこのようにさせて頂きました笑

この時空間方向のフーリエ変換は位相の符号を逆にします。
即ち$x$についてのフーリエ変換，逆フーリエ変換を以下のように定義します。

\begin{align}
f(x) &= \frac{1}{2\pi} \int_{-\infty}^{\infty} F(k_{x}) e^{i k_{x} x} d k_{x} \notag \\
F(k_{x}) &= \int_{-\infty}^{\infty} f(x) e^{-i k_{x} x} dx \notag
\end{align}

今回の記事内では空間方向のフーリエ変換は使いませんが，今後の記事で出す予定なので記載しておきます。

平面波

平面波は名前の通り平面の波が伝搬していくようなイメージです。
3次元波動方程式

$$\Delta p(\boldsymbol{r}, t) – \frac{1}{c^{2}} \frac{\partial^{2} p(\boldsymbol{r}, t)}{\partial t^{2}} = 0$$

について解を考えます。
導出の方法は幾つかありますが，今回は(簡単なので)ヘルムホルツ方程式経由で導出していきます。

定常状態(周波数領域)の音圧を考えてみます。
波動方程式をフーリエ変換するのですが，その前に

\begin{align}
p(\boldsymbol{r}, t) &= \mathscr{F}^{-1}[\mathscr{F}[p(\boldsymbol{r}, t)]] \notag \\
&= \mathscr{F}^{-1}[ P(\boldsymbol{r}, \omega) ] \notag
\end{align}

であることから

\begin{align}
\frac{\partial p(\boldsymbol{r}, t)}{\partial t} &= \frac{1}{2\pi} \frac{\partial}{\partial t} \int_{-\infty}^{\infty} P(\boldsymbol{r}, \omega) e^{-i\omega t} d\omega \notag \\
&= \frac{1}{2\pi} \int_{-\infty}^{\infty} { -i\omega P(\boldsymbol{r}, \omega) e^{-i\omega t} } d\omega \notag
\end{align}

なので両辺フーリエ変換すると

\begin{align}
\mathscr{F} \left[ \frac{\partial p(\boldsymbol{r}, t)}{\partial t} \right] &= -i\omega P(\boldsymbol{r}, \omega) \notag \\
\Rightarrow \mathscr{F} \left[ \frac{\partial^{2} p(\boldsymbol{r}, t)}{\partial t^{2}} \right] &= – \omega^{2} P(\boldsymbol{r}, \omega) \notag \\
\end{align}

とできます。このことから波動方程式の両辺をフーリエ変換すると

$$\mathscr{F} \left[ \Delta p(\boldsymbol{r}, t) – \frac{1}{c^{2}} \frac{\partial^{2} p(\boldsymbol{r}, t)}{\partial t^{2}} \right] = \Delta P(\boldsymbol{r}, \omega) + k^{2} P(\boldsymbol{r}, \omega) = 0$$

のように出来ます。ここで$k=\frac{\omega}{c}$は波数と呼ばれる空間の周波数のようなものです。

また，今後は時間領域であるか周波数領域であるかは大抵文脈で判断できるため，簡単のため上式も以下のように書くことがあります。

$$\Delta p + k^{2} p = 0$$

ちなみにこの式をヘルムホルツ方程式という音響学で非常に重要な式です。

このヘルムホルツ方程式の解を考えてみると

$$p(\boldsymbol{r}, \omega) = \tilde{A}(\omega) e^{i \boldsymbol{k} \cdot \boldsymbol{r}}$$

となることがわかるかと思います。(「2階微分=定数×元の関数」の形を考えたらすぐわかりますね！)

また，ここで$\tilde{A}(\omega)$は$\omega$についての任意関数，$\boldsymbol{k} := [k_{x}, k_{y}, k_{z} ]^{T}$は波数ベクトルです。

波数ベクトルの各成分は各方向への波数を表していて，$c,\omega$から定まる波数$k$に対して$k^{2} = k_{x}^{2} + k_{y}^{2} + k_{z}^{2}$の関係があります。

次に平面波の物理的性質から$\omega$についての任意関数$\tilde{A}(\omega)$を決定していきます。
平面波は単一周波数でのみ意味をもつのでディラックのデルタ関数を用いて

$$\tilde{A}(\omega) = 2 \pi A(\omega) \delta(\omega – \omega_{0})$$

つまり

$$p(\boldsymbol{r}, \omega) = 2 \pi A(\omega) \delta(\omega – \omega_{0}) e^{i \boldsymbol{k} \cdot \boldsymbol{r}}$$

とします。この時間領域での表現は逆フーリエ変換を用いて

\begin{align}
p(\boldsymbol{r}, t) &= \mathscr{F}^{-1}[ P(\boldsymbol{r}, \omega) ] \notag \\
&= \int_{-\infty}^{\infty} A(\omega) \delta(\omega – \omega_{0}) e^{i \boldsymbol{k} \cdot \boldsymbol{r}} e^{-i\omega t} d\omega \notag \\
&= A(\omega_{0}) e^{i \boldsymbol{k} \cdot \boldsymbol{r}} e^{-i\omega_{0} t} \notag
\end{align}

これが各周波数$\omega_{0}$の平面波の表現です。

$A(\omega)$は音の振幅や位相に関わるのですが，今回はあまり気にしないので1にして

$$p(\boldsymbol{r}, t) = e^{i \boldsymbol{k} \cdot \boldsymbol{r}} e^{-i\omega_{0} t}$$

として話を進めます。
平面波や球面波は目で見た方が分かり易いので後程pythonで可視化します。

球面波

球面波も名前の通りある点(今回は原点)から球状に広がっていく綺麗な波です。

平面波の導出ではデカルト座標の波動方程式から出発したのですが，球面波は動径(半径)方向の関数になるので，まずは波動方程式を極座標で表現することを考えます。

ここで問題なのが微分作用素$\Delta$の扱いです。

これを極座標表現するのは愚直にチェーンルールを使って微分していけば簡単にできるのですが，ちょっと面倒なのでここでは飛ばします。

結論をいうとデカルト座標の$\Delta = \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}}$は極座標で以下のようになります。

$$\Delta = \frac{1}{r^{2}} \frac{\partial}{\partial r} \left(r^{2} \frac{\partial}{\partial r}\right) + \frac{1}{r^{2} \sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial}{\partial \theta} \right) + \frac{1}{r^{2} \sin^{2} \theta} \frac{\partial^{2}}{\partial \phi^{2}}$$

今球面波を考えているので音圧は球対称になるはずなので$\theta, \phi$方向の微分は0になります。ここで

$$\frac{1}{r^{2}} \frac{\partial}{\partial r} = \frac{\partial^{2}}{\partial r^{2}} + \frac{2}{r} \frac{\partial}{\partial r}$$

かつ

$$\frac{\partial^{2} (rp)}{\partial r^{2}} = \frac{\partial}{\partial r} \left( p + r\frac{\partial (rp)}{\partial r} \right) = r\frac{\partial^{2} p}{\partial r^{2}} + 2\frac{\partial p}{\partial r}$$

なので

$$\Delta p = \frac{1}{r} \frac{\partial^{2} (rp)}{\partial r^{2}}$$

と変形できます。これより波動方程式の極座標表現は

$$\frac{1}{r} \frac{\partial^{2} (rp)}{\partial r^{2}} – \frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}} = 0$$

両辺$r$を乗じて

$$\frac{\partial^{2} (rp)}{\partial r^{2}} – \frac{1}{c^{2}} \frac{\partial^{2} (rp)}{\partial t^{2}} = 0$$

という形にできます。これはデカルト座標で考えた波動方程式と殆ど同じ形ですね。
平面波と同様に解を考えると以下のようになるのがすぐに分かります。

\begin{align}
rp &= \tilde{A}(\omega) e^{ikr} \notag \\
\Rightarrow p &= \tilde{A}(\omega) \frac{e^{ikr}}{r} \notag
\end{align}

これが球対称な解の形です。球面波を考えると単一周波数で意味を持つので

$$\tilde{A}(\omega) = 2 \pi A(\omega) \delta(\omega – \omega_{0})$$

つまり平面波同様に

$$p(\boldsymbol{r}, \omega) = 2 \pi A(\omega) \delta(\omega – \omega_{0}) \frac{e^{ikr}}{r}$$

となって逆フーリエ変換により時間領域表現

$$p(\boldsymbol{r}, t) = A(\omega_{0}) \frac{e^{ikr}}{r} e^{-i\omega_{0} t}$$

これが各周波数$\omega_{0}$の球面波の表現です。
平面波同様$A(\omega)$を1として

$$p(\boldsymbol{r}, t) = \frac{e^{ikr}}{r} e^{-i\omega_{0} t}$$

を球面波として実装していきます。

pythonによる平面波と球面波の可視化

ではpythonで今まで導出した単純な波を可視化していきます。

平面波・球面波に限らず可視化は物理現象のイメージを捉えるのにとても重要なのでpythonを使う方は是非色んな描画の方法を試してみて下さい。

平面波のコーディング

平面波を描画するpythonのコードは次のようになります。

freq = 1000
omega = 2 * np.pi * freq
rho = 1.293
kappa = 142.0e3
c = np.sqrt(kappa/rho)
k = omega / c
theta = np.pi/2
k_vec = k*np.array([np.cos(theta), np.sin(theta)])
row = 201
col = 201
width = 0.01
y, x = np.mgrid[0:row:1+width, 0:col:1+width] * width
y = y - width*(row/2)
x = x - width*(col/2)
disp = np.array([np.min(x), np.max(x), np.min(y), np.max(y)])
p = np.exp(1j * (x*k_vec[0] + y*k_vec[1])) # plane wave
fig, ax = plt.subplots()
ax.imshow(p.real, cmap=plt.cm.bwr, vmin=-1, vmax=1, interpolation='bicubic', origin='lower', extent=disp)
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_xticks([-1, -0.5, 0, 0.5, 1], minor=False)
ax.set_yticks([-1, -0.5, 0, 0.5, 1], minor=False)

出力される画像は次のようになります。

平面の波面であることがわかりますね。

球面波のコーディング

球面波を描画するpythonのコードは次のようになります。

freq = 1000
omega = 2 * np.pi * freq
rho = 1.293
kappa = 142.0e3
c = np.sqrt(kappa/rho)
k = omega / c
source = np.array([0, -0.8])
row = 201
col = 201
width = 0.01
y, x = np.mgrid[0:row:1+width, 0:col:1+width] * width
y = y - width(row/2) x = x - width(col/2)
disp = np.array([np.min(x), np.max(x), np.min(y), np.max(y)])
r = np.sqrt(np.power(x-source[0],2) + np.power(y-source[1],2))
p = np.exp(1jkr) / (4np.pir) # spherical wave
p[np.isinf(p)] = 0
p[np.abs(p)>1] /= np.abs(p[np.abs(p)>1])
fig, ax = plt.subplots()
ax.imshow(p.real, cmap=plt.cm.bwr, vmin=-1, vmax=1, interpolation='bicubic', origin='lower', extent=disp)
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.set_xticks([-1, -0.5, 0, 0.5, 1], minor=False)
ax.set_yticks([-1, -0.5, 0, 0.5, 1], minor=False)

ここで実行の際に警告が出るかと思いますが，気にしなくていいです(0割りの警告です)。
※もしどうしても気になるのであればsourceの場所を若干ずらすと解決します

また，出力される画像は次のようになります。

こちらも球面状の波面を見ることができます。

animation関数による動画生成

最後にanimationという関数を用いた動画生成のコーディングを少しだけお話しします。

matplotlibに入っているものなので次のようにインポートすれば良いです。

import matplotlib.animation as animation

今回は保存方法や関数の詳細等は省いてコードだけお示しします。
平面波の動画を生成するコードは以下のようになります。

def plot_time_domain(i, fs, P, disp, omega, ax):
    t = i / fs
    p = P * np.exp(-1j * omega * t)
    ax.cla()
    ax.imshow(p.real, cmap=plt.cm.bwr, vmin=-1, vmax=1, interpolation='bicubic', origin='lower', extent=disp)
    ax.set_xlabel("x [m]")
    ax.set_ylabel("y [m]")
    ax.set_xticks([-1, -0.5, 0, 0.5, 1], minor=False)
    ax.set_yticks([-1, -0.5, 0, 0.5, 1], minor=False)
    fig.tight_layout()
    fig.show()

# main
fps = 15
bitrate = 10000
fs = 48000
freq = 1000
nframes = int(fs/freq)
omega = 2 * np.pi * freq
rho = 1.293
kappa = 142.0e3
c = np.sqrt(kappa/rho)
k = omega / c
theta = np.pi/2
k_vec = k*np.array([np.cos(theta), np.sin(theta)])
row = 201
col = 201
width = 0.01
y, x = np.mgrid[0:row:1+width, 0:col:1+width] * width
y = y - width*(row/2)
x = x - width*(col/2)
disp = np.array([np.min(x), np.max(x), np.min(y), np.max(y)])
P = np.exp(1j * (x*k_vec[0] + y*k_vec[1]))
fig, ax = plt.subplots()
ani = animation.FuncAnimation(fig, plot_time_domain, fargs=(fs, P, disp, omega, ax), interval=100, frames=nframes, repeat=False)

少し長いですが，殆ど先程の画像の描画コードと同じです。
変更点は関数を1つ追加したのとanimation関数を用いたことくらいでしょうか。

これによって生成された動画は以下です。

球面波も同じようにして生成出来ますので省略しますが，結果は以下のようになります。

動くとより分かりやすいですね！

まとめ

今回は音響学で重要な平面波と球面波について概説し，導出した後pythonで可視化してみました。

式で見ると難しいように感じるかもしれませんが，慣れればそんなことはないですし可視化すればとっても単純で綺麗な性質をもっています。

参考書としては以下のものを挙げておきます。

音響分野の非常に有名な先生が書いており，基礎から分かりやすく書いてあります。
音響学最初の一冊としてもオススメです。

それではお疲れ様でした！