#はじめに
本記事は、[厚木の民 Advent Calendar 2019][1] の22日目の記事です。
[1]:https://adventar.org/calendars/4560

こんにちは、ホリエといいます。
本日は冬至だそうですね。ゆずサワーを飲んでいます。おいしい

さてみなさん、
ふと**「音」を「見て」みたくなる**ことってありませんか？？ありますよね！
そんなときのために、今回はPythonを使った「音の可視化」方法を紹介していきます。

##だれ向けの記事？

1. [厚木の民][1]
2. 音とプログラミングに興味がある人
3. 周波数ビンごとの情報（位相とか抑圧ゲインとか）をPythonで可視化したい人

3 でソースが見たいだけという方は後半のほうまでスクロールしてください。

ちなみに私はPythonへの造詣が深くない人間なので、至らない点は
コメント欄で優しくご指導ご指摘してしていただけると嬉しいです。
余談ですが造詣は「ぞうけい」と読むのをここで初めて知りました。とても恥ずかしかったです

#目次

• 音を「見る」？？★
• スペクトログラム★
• Pythonでスペクトログラムを描く★
• STFTとimshowを利用した自作スペクトログラム表示関数

サクッと読むときは★印の部分だけ読んでください！

#音を「見る」？？
無色透明な音ですが、その姿を見る方法はいくつかあります。
なじみ深いところだと、

• 音楽プレイヤーなどに入っているボリュームメーター
• （音）オシロスコープ
• クントの実験
  などがあります。クントの実験は科学館などで
  見たことがある人も多いのではないでしょうか？

今回は、音の高さ成分（周波数成分）の時間変化をひと目で確認できる
スペクトログラムというグラフについて紹介しようと思います！

#スペクトログラム
説明はいいからどんなもんか見せろよ！うおお！
という方も多いと思うので、はじめにグラフを見せますね。

https://ja.m.wikipedia.org/wiki/%E3%83%95%E3%82%A1%E3%82%A4%E3%83%AB:Spectrogram_of_violin.png

謎のしましまが見えますね…？
このグラフの縦軸、横軸、色はそれぞれ、

• 縦軸 … 周波数（音の高さ）。上方向にいくほど高い音になる
• 横軸 … 時間
• 色 …… エネルギーの大きさ。上の図だとオレンジ色の部分がエネルギーが大きい部分

を表しています。声のスペクトログラムは「声紋」なんて言われたりもしています。
ちなみに、つよいオタクエンジニアが音声のスペクトログラムを見ると、
音を聞かずに誰（女性声優）の声かがわかるらしいです。こわいですね

たとえば時報の音（ポン↓・ポン↓・ポーン↑）
https://ultra.zone/file/jihou-sine-2f.mp3
を表示してみると、

こんな感じのスペクトログラムになります。
ポン・ポン と比べて、最後のポーンは音が高いので、
スペクトログラムの線が上にあることがわかりますね。
文字に起こすと語彙力ゼロの説明みたい

Pythonでスペクトログラムを描く

さてさて、このスペクトログラムですが、[フリーソフト][2]を使って簡単に描画することができます。
[2]:https://ja.wikipedia.org/wiki/%E6%B3%A2%E5%BD%A2%E7%B7%A8%E9%9B%86%E3%82%BD%E3%83%95%E3%83%88%E3%82%A6%E3%82%A7%E3%82%A2
が、あえて、これをpythonでやっていこうと思います。

import wave
import matplotlib.pyplot as plt
import numpy as np

def wavread(filepath):
    """
    .wav形式のファイルからデータを取得する

    Parameters
    ----------
    filepath : str
        ファイルのパス

    Returns
    -------
    fs : float
        サンプリング周波数
    data : tuple
        音声データ
    """

    wf = wave.open(filepath, "rb")      # ファイルを開く
    fs = wf.getframerate()              # サンプリング周波数の取得
    # getnframesで得たデータ点数分readし，np.frombufferでバイナリから数値化
    wavdata = np.frombuffer(wf.readframes(wf.getnframes()), dtype="int16")
    wf.close()

    # データのプロパティ表示
    print("\n-----input data infomation-----")
    print(f"filepath       : {filepath}")
    print(f"sampling rate  : {fs}[Hz]\n")
    return fs, wavdata


def main():
    fs, wavdata = wavread("wavファイルのパス")
    nfft = 1024  # 高速フーリエ変換点数

    plt.figure()
    plt.specgram(wavdata, NFFT=nfft, Fs=fs)
    plt.show()

if __name__ == '__main__':
    main()

出力スペクトログラムがこちら

このプログラムでは、

1. .wav形式の音声ファイルを音声入出力ライブラリwaveを使って数値データとして読み込み
2. グラフ描画ライブラリmatplotlib内のspecgramという関数に渡し、スペクトログラムの描画
   をしています。Pythonだとこんなに短く簡単に書けてしまうんですね。ライブラリ is つよい

というわけで、ふと手持ちの音楽ファイルや自分の録音した声のスペクトログラムを見てみたくなったら
上記のコードを実行してみてもらえればと思います。

ちなみにAudacityやAuditionなどのソフトを入れたほうがたぶん10倍くらい楽です。Auditionは有料だけど無料体験もできるよ！

#ここからちょっと専門的な話
さすがにライブラリの関数使って終わってしまうのはアレなので
さっきの関数の中身まで掘り下げていこうと思います。

#短時間フーリエ変換（STFT: Short Term Furier Transform）
ものものしい名前ですね。。。

スペクトログラムですが、これを描くためには内部でこんな処理をしています。

1. 時間信号を一定の幅で取り出す
2. 高速フーリエ変換によって周波数スペクトルを得る
3. 1~2を切り出し箇所を進めながら繰り返し適用し、得られたスペクトル列を2次元配列化して3次元カラーマッピングする

#サンプリング周波数？
私たちがスマホで音楽を聴くときも、電話をするときも、聞こえてくる音は一度デジタルなデータとして扱われているわけですね。
サンプリング周波数は、連続した信号（音波：空気の疎密波 や、音波をマイクで電気化した信号など）を、1秒間にどのくらいの間隔でデジタルデータとして取り込むか、を表しています。

この図だと緑色の線が連続信号、青い点々がサンプリングしたあとの信号を表しています。
サンプリング周波数を高くすればするほど、より細かい時間周期の波（=高い周波数の音）をデータとして取り込むことができるようになります。

サンプリングには、**サンプリング周波数の半分の周波数未満の信号は、
標本化した後、もとの状態に完全に戻すことができる。**という性質があります。
（気になる人は標本化定理で調べてみてください）

人間が音として知覚できる一番高い音の周波数は約20000 Hzと言われているので、40000 Hzよりも大きいサンプリング周波数で音信号をサンプリングしておけば、人が聞こえる帯域の音は標本化しても完全に再現できる、というわけですね！（※量子化は今回考えないことにする）こういう理屈があって、CD音楽などのサンプリング周波数は44100 Hzに設定されています。

#フーリエ変換？FFT？

Wikipediaのこのアニメーションが直感的でウルトラわかりやすいのですが、
スペクトログラムで用いているフーリエ変換の役割は、
「切り出してきた音声区間を無限に続く周期信号とみなして、その周期信号の周波数成分を取得する」って感じです。

フーリエ変換によって、周波数スペクトル（解析区間にどんな周波数が含まれているか？という情報）を得ることができます。
コンピュータの世界では、FFTという昔のえらい人が考えた計算方法で、すばやく周波数成分を得ることができるようになっています。すごい

#FFT点数？
FFTを計算する区間（FFT点数）が長くなれば長くなるほど、より細かく周波数成分を解析することができます。
しかし、解析区間が長いと、その区間内での周波数成分の細かい変動なんかがわからなくなってしまいます。（時間分解能が下がる）
時間-周波数解析では、周波数方向の解析の細かさと時間方向のそれにトレードオフの関係があります。

広い区間で周波数解析を行えば、より細かい周波数情報を得ることができる一例です。
下段の周波数成分のグラフを見ると、解析する区間が短くなっていく（オレンジ→緑→赤）につれて
スペクトル（青い点々）が疎になっていることがわかると思います。

#実装

def stft(audio_in, nfft=1024):
    """
    STFTの実装

    Parameters
    ----------
    audio_in : np.array
        処理対象の音声データ(モノラル)
    nfft : int
        FFT点数

    Returns
    -------
    output : np.ndarray
        出力配列 今回はスペクトログラム
    """

    # 周波数解析窓 今回はハミング窓を使用
    window = np.hamming(nfft)


    # 入力音声データ長がFFT点数の倍数でなければ、不足分を0埋め
    if len(audio_in) // nfft != 0:
        # 二次元配列の行方向、終端に不足分の0を追加
        audio_in = np.pad(audio_in, [0, int(nfft - len(audio_in) % nfft)], 'constant')


    # 出力配列の初期化(切り出しループ回数 * スペクトル数)
    output = np.zeros(len(audio_in) // nfft * (nfft // 2 + 1))


    # 信号切り出しループ （信号長/FFT点数長）回繰り返す
    for start in range(len(audio_in)//nfft): 
            

        # FFT点数ぶん入力信号を切り出し
        frame = audio_in[start * nfft:start * nfft + nfft]


        # 周波数解析窓をかけ、FFTして周波数スペクトルを得る
        frame_w = frame * window
        spec = np.fft.rfft(frame_w, nfft)


        # 振幅スペクトル（周波数スペクトルの絶対値）を計算
        amp_spec = abs(spec)
        
        
        # 出力配列に代入
        output[start * (nfft // 2 + 1):start * (nfft // 2 + 1) + (nfft // 2 + 1)] = amp_spec


    # 1次元配列を2次元配列化
    output = np.reshape(output, (len(output)//(nfft // 2 + 1), (nfft // 2 + 1)))
    return output


def imshow_heatmap(data):
    # 2次元配列をヒートマップ表示
    aspect = data.shape[1] / data.shape[0] # 縦横比が正方形になるように調整
    plt.figure()
    plt.imshow(data, interpolation='nearest', cmap='jet' ,aspect=aspect)
    plt.colorbar(fraction=0.030)
    plt.gca().invert_yaxis() # y軸を反転
    plt.show()

fs, wav = wavread(r"C:\hogehoge.wav")
specgram = stft(wav, nfft=1024)
imshow_heatmap(20*np.log10(specgram.T))

ここまで説明したものを実装したものがこちらです。パワースペクトルをdB表示したものになっています。
2次元配列のカラーマッピングには、seaboanのheatmapも使えます。見た目が非常におしゃれになるので良い...すき...なのですが地獄のように描画が重いです。
imshowは描画領域のアスペクト比を指定する必要があるのが若干めんどくさいですが、体感1000倍くらい高速なので、目的に応じて使い分けてみてください。

#あとがき
初めて記事書きました。泣くほど大変でしたがいい経験になった気がしています。
機会があったらはフーリエ変換の位相のところとか周波数解析窓のところとかをもっと詳しく書きたいかもです。機会があったら。