第2章　音に触れる（お試し版）

本章では，最も単純かつ重要な波形である「正弦波」を計算機で生成することから始めます。そして，アナログ／ディジタル変換 (A/D変換）について学びます。

なお，本章において出力される音は，右のQRコードのページでまとめて聞くことができます。

2.1 音を数式で表現する

今後の作業に備えて，波形を描く関数 plot_wave をここで定義しておきます。（関数 plot_wave は，第3章以降では徐々に機能を増やしていきます。ここでは，最も基本的な部分だけを示します。）

In [2]:
def plot_wave(time, amplitude, xtitle = 'Time (s)', \
              ytitle = 'Amplitude (arb.)', hold = False, \
              color = 'blue', marker = ',', legend = '', linestyle = '-'):
    ''' 時間軸データと波形データを受け取り，波形を描く関数の定義
    引数 time:      時刻の離散データ
         amplitude: 上記の時刻データに対応する瞬時振幅値
         xtitle:    横軸のラベル（暗黙値は 'Time (s)')
         ytitle:    縦軸のラベル（暗黙値は 'Amplitude (arb.)'）
         hold:      False（暗黙値）ならば描画し，True ならばデータを保持する。
         color:     グラフの色（暗黙値はブルー）
         marker:    マーカの種類（暗黙値は pixel）
         legend:    凡例の文字列（暗黙値は「なし」）
         linestyle: 線の種類（暗黙値は実線）
    
    '''

    if (marker != ','):  # マークが pixelでない場合には，
        linestyle = ''   # データのプロットのみを行い線では結ばない。

    plt.plot(time, amplitude, color = color, marker = marker, \
             linestyle = linestyle, label = legend) # この１行で描画される。

    if (legend != ''): # 凡例のデータが渡された場合には，凡例を書く関数を呼ぶ。
        plt.legend()

    if (hold == False): # hold == False のときは，plt.show() を呼んで描画する。
        plt.xlabel(xtitle)
        plt.ylabel(ytitle)
        plt.show()

2.1.1 時間の関数としての音の数式表現

高校の物理では，音が伝搬する様子を観察したでしょう（横軸は位置 $x$ のグラフでした）。一方，以下では，時間の関数としての音の表現を考えていきます（横軸は時刻 $t$ のグラフです）。音が伝搬しているときに，ある位置においては「時間とともに圧力が変化」しています。もし，その位置に耳があれば，時間とともに鼓膜が押される／引かれるが繰り返されるので，音として聞こえるのです。その位置にマイクロホンを置いたと考えて，100 Hz の音を 0.1 s 間だけ録音したのがアウトプット2.1の波形です（インプット2.1のプログラムにより描きました）。

In [3]:
# インプット2.1
f = 100.0                       # 周波数が100 Hz の音を考える。
t = np.arange(0, 0.1, 0.1e-3)   # 横軸の範囲を 0～0.1 s とする。
                                # アナログ波形に見えるように 0.1 ms ごとに波形を描く。
p = np.cos(2.0 * np.pi * f * t) # この式が，本節の肝である。
plt.title("Output 2.1")
plot_wave(t, p)

アウトプット2.1に示したのは，代表的な音である余弦波です（このような単一周波数のみを含む音は，澄んだ音色に聞こえるので，音響学の分野では「純音」と呼ばれます）。この純音の数式表現を考えましょう。

高校生のとき，三角関数を定義するのに「単位円上で点を回したこと」を覚えていますか？$\cos (\alpha)$ は「点 (1, 0) から出発して反時計回りに回る点の $x$ 軸への射影」として与えられました。このとき，$\alpha$ はその点と原点を結ぶ線分が， $x$ 軸となす角です．

周波数が $f$ [Hz] の余弦波は，時刻 $t$ [s] の関数として $\cos(2 \pi f t)$ と表されます。周波数が $f$ [Hz] であることは，1 [s] に $f$ 個の波が含まれることを意味します。余弦関数は1周期が2$\pi$ [rad] ですから，点 (1, 0) から出発して $t$ [s] 後には $\alpha＝ 2 \pi f t$ [rad] だけ回るので，このような式で表されます。

【重要】 周波数が $f$ [Hz] の余弦波の数式表現は $\cos(2 \pi f t)$

ところで，周波数が$f$ [Hz] の波は，$1/f$ [s] ごとに繰り返します（アウトプット2.1では，100 Hz の波は 0.01 s ごとに繰り返しています）。この繰り返しに要する時間は周期 $T_0$と呼ばれます。明らかに周期 $T_0 = 1/f$ [s] です。なお，数式の中に出現する $2 \pi f$ は角周波数 $\omega$と呼ばれ，単位は [rad/s] です。

2.1.2 振幅と位相を含んだ数式表現

音の性質を記述するためには，３つの変量：周波数 $f$, 振幅 $A$, 位相 $\theta$ が重要です。これらを含む，より一般的な音の数式表現は， $$ f(t) = A \cos(2 \pi f t + \theta) $$ です。左辺の $f$ は時間の関数 (function) の f で，右辺の $f$ は周波数（frequency）の f であることに注意してください。

振幅 $A$ は，波形の最大値・最小値を決めることは明らかです。以下では，位相 $\theta$ について考えます。ここでは，$f$ = 10 [Hz]，$A$ = 2 として， $\theta$ だけが異なる以下の３つの波を描画してみます（インプット2.2とアウトプット2.2）。 $$ f_0(t) = A \cos(2 \pi f t) ,\ \ \ f_1(t) = A \cos(2 \pi f t + \pi/2), \ \ \ f_2(t) = A \cos(2 \pi f t - \pi/2), $$

In [ ]:
# インプット2.2
f = 10.0              # 周波数は10 Hz の音を考える。
A = 2.0               # 振幅は 2 とする。
theta = np.pi / 2.0   # 位相は± π/2 ずれた場合をとりあげる。

t = np.arange(0, 0.2, 1e-3) # 横軸の範囲を 0.2 s 間として，
                            # アナログ波形に見えるように 1 ms ごとに波形を描く。

f0 = A * np.cos(2.0 * np.pi * f * t)
f1 = A * np.cos(2.0 * np.pi * f * t + theta)
f2 = A * np.cos(2.0 * np.pi * f * t - theta)

plt.title("Output 2.2")
# hold = True では，図を描かずデータのみ送る。
plot_wave(t, f0, hold = True,  color = 'b', linestyle = '-',  legend = '$f_0(t)$')
plot_wave(t, f1, hold = True,  color = 'r', linestyle = '--', legend = '$f_1(t)$')
plot_wave(t, f2, hold = False, color = 'g', linestyle = ':',  legend = '$f_2(t)$')
# hold = False で，これまでに送られたデータを描画する。

アウトプット2.2から明らかなように，位相項がつくと波形が左右にずれます。では位相項が $\theta$ である波 $f_1(t)$ は，時間軸上で $f_0(t)$ からどれだけずれるかを考えましょう。そのために，数式を以下のように変形します。

$\begin{eqnarray*} f_1(t) = A \cos (2\pi f t + \theta) = A \cos \left( 2 \pi f(t + \theta / 2\pi f) \right) = A \cos \left (2 \pi f(t + \theta / \omega) \right) \\ \end{eqnarray*}$

この式のグラフは，$ f_0(t) = A \cos(2 \pi f t) $ に比べて「左に $\theta / \omega $ だけずれる」ことを高校数学で勉強したことを覚えていますよね。 このように時間軸上で左（＝より先の時刻）にずれることを位相が進むといいます。 逆に，位相が $-\theta$ のときは，グラフは時間軸上で右（＝より後の時刻）にずれ，位相が遅れるといいます。

ところで，高校数学で勉強したとおり $f_2(t) = A \cos(2 \pi f t - \pi/2) = A \sin(2 \pi f t) $であり，そのことは上図からも確かめられます。すると，もう正弦波／余弦波といった区別は不要で，ある周波数の波は「余弦波に位相項がついたもの」で統一的に表現できることに気付きます（このことは，後で「位相スペクトル」を学ぶときに重要になります）。

【重要】 正弦波は余弦波よりも $\pi/2$ だけ位相が遅れた波（逆に，余弦波は正弦波よりも $\pi/2$ だけ位相が進んだ波）

[確認課題]（確認課題は，書籍本体には含まれません。サポートページのipynbノートブックに収録しています。）

（Colabメニューバーにある [+コード] ボタンを押して新しいセルを作り，計算すること）

(1) 空気中の音速 $c$ [m/s] は温度 $T$ [°C]により変化し，$c= 331+ 0.6 T$ で近似的に与えられる。$T=20$ のとき，$c$ を求めなさい。

(2) 周波数 $f = 1000$ [Hz] のとき，波長 $\lambda$ [m] を求めなさい。ただし，音速 $c = 340$ [m/s] とする。（重要：1000 Hz は，音響学の種々の場面で基準として用いられる周波数です。そこで，1000 Hz の波の波長を瞬時に導けるようにしておくことは，回折などの音響事象を考察する際に有益です。）

(3) 周波数 $f = 1000$ [Hz] の正弦波について，時間 $t$ の関数として数式表現を示しなさい。sin を使った場合と，cos を用いた場合の両方を示すこと。（紙と鉛筆で示してくれれば十分です。もし可能であれば，[+テキスト] ボタンで新しいセルを作り，\$ と \$ で囲んだ数式を書いて，Shift+Enterで表示してみてください。）

2.2 正弦波を生成して聞いてみる

音のディジタルデータを生成するのは簡単です。アナログ波形について，一定の時間間隔（標本化周期と呼ぶ）ごとに値を調べて配列に納めるだけです。なお，標本化周期 $T_{\rm s}$ [s] は，標本化周波数 $f_s$ [Hz] （1秒間に値を調べる回数）の逆数として与えられます。早速，以下の課題に取り組んでみましょう。

[課題] 周波数: 1000 Hz，振幅: 1, 継続時間: 1 s の余弦波を生成して聞いてみなさい。また波形を描きなさい。なお，標本化周波数: 44.1 kHz とする。

インプット2.3により音を聞くことができます（再生ボタン▶を押すときには，音量に注意してください）。これは，時刻を表す配列 t を 1/$f_s$ ごとに「粗くとる」ことを除けば，インプット2.2までと全く同様です。

In [7]:
# インプット2.3
f  = 1000.0     # 1000 Hz の音を生成する。
A  = 1.0        # 振幅は 1 とする。
fs = 44100.0    # 標本化周波数は 44.1 kHz とする。

t = np.arange(0, 1, 1/fs)  # 横軸の範囲 1 s 間を，1/fs [s] ごとに区切る時刻を配列 t に納める。
sampledWave = A * np.cos(2.0 * np.pi * f * t)  # 標本値を配列 sampledWave に納める。

audio = Audio(sampledWave, rate = fs)  # 引数 rate に標本化周波数を指定します。
audio                                  # オーディオオブジェクト audio を再生します。
Out[7]:

次に，インプット2.4により波形を確認します（アウトプット2.4では，最初の100個のデータだけプロットします）。

In [6]:
#インプット2.4（標本化された系列の描画）
plt.title("Output 2.4")
plot_wave(t[0:100], sampledWave[0:100], marker = 'o', color = 'r')

【備考】これまでにA/D変換を勉強したことのある人は，この課題を与えられたとき，「あれ？量子化精度は？」という疑問が浮かんだと思います。実は，量子化精度を問題にするのは，「最終的に音をWAVファイルなどに記録するときだけ」と考えても不都合はないのです。すなわち，計算機上で音の処理をする場合は，もっぱら「実数値（浮動小数点演算）」で計算を行います（無限に量子化精度を向上させた，と考えても結構です）。ちなみに 上記の Audio 関数には，振幅を気にせずにデータを渡しても大丈夫です（ちょどよい振幅に調整した後で，D/A変換されて再生されます）。

2.3 A/D変換について確認する

前節では，量子化精度は気にしなくてよいことを述べました。しかし，やはり一度は A/D 変換についてしっかり勉強しておくことは重要でしょう。

A/D (analog to digital) 変換は，アナログ信号に対して標本化と量子化の2つのステップを経ることで，ディジタル信号を得る変換です。私達がコンピュータで音を扱うために，必須の変換です。マイクロホンで音圧→電圧の変換を行った上で，A/D変換を行います。それゆえ，以下のグラフの縦軸は「電圧」に関するものと考えてください。

ここでは，50 Hz で振幅 10000 の cos波をA/D変換してみます。まずアナログ（とみなせる）波形を 0.1 s 間分だけ描きます（インプット2.5とアウトプット2.5）。

In [9]:
#インプット2.5
f = 50.0                    # 周波数は 50 Hz の音を考える。
A = 10000.0                 # 振幅は 100000 とすう。

t = np.arange(0, 0.1, 0.1e-3) # 横軸の範囲を 0～0.1 s とする。
                              # アナログ波形に見えるように 0.1 ms ごとに波形を描く。
ft = A * np.cos(2.0 * np.pi * f * t)

plt.title("Output 2.5")
plot_wave(t, ft)

2.3.1 標本化

まず標本化のステップです。標本化とは，標本化周期 $T_{\rm s}$ ごとにアナログ波形の値を調べることです。このとき，標本化周波数 $f_{\rm s}$ は $f_{\rm s} = 1 / T_{\rm s} $ で与えられます。代表的な $f_{\rm s}$ として，CD の 44.1 kHz は知っておきましょう。

標本化定理に基づけば，周波数$f$の信号を正しく標本化するためには，標本化周波数 $f_{\rm s}$ は $f_{\rm s} > 2f$ とする必要があります。そこで，ここでは 信号周波数（50 Hz）の2倍よりも十分に高い$f_{\rm s} = 500$ [Hz] で標本化してみます。プログラムと出力を，インプット2.6とアウトプット2.6に示します。

In [10]:
# インプット2.6
fs = 500.0  # 標本化周波数は 500 Hz とする。

sampledT = np.arange(0, 0.1, 1/fs)                   # 横軸の範囲 0.1 s 間を，1 / fs [s] ごとに区切る時刻を配列 sampledT に納める。
sampledWave = A * np.cos(2.0 * np.pi * f * sampledT) # 標本化周期ごとの値を調べて配列 sampledWave に納める。

plt.title("Output 2.6")
plot_wave(t, ft, hold = True)                               # アナログ波形（のつもり）の描画
plot_wave(sampledT, sampledWave, marker = 'o', color = 'r') # 標本化された系列の描画

「お試し版」はここまでです。

書籍本体では，この後に以下の内容が掲載されています。

2.3.2 不適切な標本化

2.3.3 量子化

2.3.4 よろしくない量子化

2.3.5 量子化誤差を聞いてみる

　最後の例題として「4 bit量子化した音を聞いてみる」を行います。波形と音色の関係について，ちょっと不思議に感ずるかもしれません。

さらにサポートページの第2章には，音圧レベルに関する例題が用意されていますのでご参照ください。