Pythonを使って音声データからスペクトログラムを作成する

概要

Pythonを使って音楽ファイルを読み込み、サウンドスペクトログラムを作る手順を整理しました。大まかな手順は以下です。

  • pydubを使って音楽ファイルを読み込み
  • numpyを使ってフーリエ変換を実施し、スペクトログラムを作成
  • matplotlibとseabornでスペクトログラムを可視化
  • サウンドスペクトログラムとは以下のような横軸時間、縦軸周波数スペクトルをとったヒートマップです。

    スペクトログラムって何?ということについては↓このあたりがイメージつかみやすいかと思います。
    http://blog.media.teu.ac.jp/2015/05/post-0ec1.html
    音声ファイルはスーパーマリオ64のウォーターランドの曲を使わせていただきました。

    https://www.nintendo.co.jp/ds/series/dsi/menu/sound/download.html

    以下のモジュールをインポートしておきます。(ない場合にはpip等でインストールおねがいします)

    In [1]:
    from pydub import AudioSegment
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    import seaborn as sns
    
    %matplotlib inline
    

    なお、pydubを使うにはffmpegが必要になるので、準備しておきましょう。以下記事等を参考にしてみてください。
    Pythonのmatplotlibでgifアニメを作成する

    pydubを使って音楽ファイルを読み込み

    読み込みは以下1行で実行できます。以下例ではm4aを対象にしていますが、mp3やwav等でも同様のコード(m4a部分を変える)で読み込むことができるはずです。

    In [2]:
    sound = AudioSegment.from_file("ウォーターランド.m4a", "m4a")
    

    基本情報は以下で取得できます。
    https://qiita.com/itoru257/items/8af2902d8ce851ae74ea

    In [3]:
    #チャンネル数(1:mono, 2:stereo)
    sound.channels
    
    Out[3]:
    2
    In [4]:
    #サンプルレート(Hz)
    sound.frame_rate
    
    Out[4]:
    44100
    In [5]:
    #再生時間(秒)
    sound.duration_seconds
    
    Out[5]:
    55.72643990929705

    get_array_of_samplesを使うことで、生データをListで受け取ることができます。後々のことを考えて、numpyで読み込み、加えてステレオ音声から片側だけの音を抽出しておきます。

    In [6]:
    samples = np.array(sound.get_array_of_samples())
    sample = samples[::sound.channels]
    

    簡単にプロットするとこんな感じ。

    In [7]:
    plt.plot(sample[::10])
    
    Out[7]:
    [<matplotlib.lines.Line2D at 0x7f29b2aaa518>]

    numpyを使ってフーリエ変換を実施し、スペクトログラムを作成

    まずはスペクトルを計算してみる

    フーリエ変換はnumpyの機能を使って簡単に実施することができます。
    https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html

    In [8]:
    #フーリエ変換自体は以下1行を実行するだけ
    spec = np.fft.fft(sample)
    
    In [9]:
    #スペクトルの値(複素数)が入っている。
    spec
    
    Out[9]:
    array([ 909365.00000000     +0.j        ,
             51629.03247526+667871.8941476j ,
           -461312.87203113-154156.11920285j, ...,
            335244.22709894+308623.48241612j,
           -461312.87203112+154156.11920289j,   51629.03247527-667871.89414759j])

    スペクトル(複素数)の取り扱いを簡単に列挙しておきます。

    In [10]:
    #実数部分を取り出し
    spec.real
    
    Out[10]:
    array([ 909365.        ,   51629.03247526, -461312.87203113, ...,
            335244.22709894, -461312.87203112,   51629.03247527])
    In [11]:
    #虚数部分を取り出し
    spec.imag
    
    Out[11]:
    array([      0.        ,  667871.8941476 , -154156.11920285, ...,
            308623.48241612,  154156.11920289, -667871.89414759])
    In [12]:
    #絶対値を取り出し
    np.abs(spec)
    
    Out[12]:
    array([ 909365.        ,  669864.48180705,  486388.39931611, ...,
            455672.19105605,  486388.39931612,  669864.48180704])
    In [13]:
    #偏角を取り出し
    np.angle(spec)
    
    Out[13]:
    array([ 0.        ,  1.49364597, -2.81909085, ...,  0.74407664,
            2.81909085, -1.49364597])
    In [14]:
    #絶対値と偏角から複素数を復元
    np.exp(1j*np.angle(spec))*np.abs(spec)
    
    Out[14]:
    array([ 909365.00000000     +0.j        ,
             51629.03247526+667871.8941476j ,
           -461312.87203113-154156.11920285j, ...,
            335244.22709894+308623.48241612j,
           -461312.87203112+154156.11920289j,   51629.03247527-667871.89414759j])
    In [15]:
    #逆変換
    np.fft.ifft(spec)
    
    Out[15]:
    array([  1.00000000e+00 -3.23473530e-12j,
             1.00000000e+00 -2.99057462e-12j,
             1.00000000e+00 -1.14700955e-12j, ...,
            -1.21269118e-11 -6.54196579e-12j,
            -9.99257535e-12 -7.04638408e-12j,  -1.04776518e-11 -7.39542712e-12j])
    In [16]:
    #周波数
    np.fft.fftfreq(sample.shape[0], 1.0/sound.frame_rate)
    
    Out[16]:
    array([ 0.        ,  0.0179448 ,  0.03588961, ..., -0.05383441,
           -0.03588961, -0.0179448 ])
    In [17]:
    #スペクトルをプロット表示
    freq = np.fft.fftfreq(sample.shape[0], 1.0/sound.frame_rate)
    plt.plot(freq, np.abs(spec))
    
    Out[17]:
    [<matplotlib.lines.Line2D at 0x7f29b2986c50>]

    負の領域は使わないことや、周波数0の値が2倍されていること、対数メモリのほうが見やすことに留意して以下のように表示します。

    In [18]:
    freq = freq[:int(freq.shape[0]/2)]
    spec = spec[:int(spec.shape[0]/2)]
    spec[0] = spec[0] / 2
    
    #plt.xscale("log")
    plt.yscale("log")
    plt.plot(freq, np.abs(spec))
    
    Out[18]:
    [<matplotlib.lines.Line2D at 0x7f29b28fbef0>]

    どうやら100〜1000Hzのスペクトルが強いことが分かりました。

    スペクトログラムを作成

    上では全体の波形を使って1つのスペクトルを算出しましたが、これを一定の幅を取り、それを刻み幅ずつ動かしながら算出していくことでスペクトログラムが作成できます。

    In [19]:
    #窓幅
    w = 1000
    #刻み
    s = 500
    
    In [20]:
    #スペクトル格納用
    ampList = []
    #偏角格納用
    argList = []
    
    #刻みずつずらしながら窓幅分のデータをフーリエ変換する
    for i in range(int((sample.shape[0]- w) / s)):
        data = sample[i*s:i*s+w]
        spec = np.fft.fft(data)
        spec = spec[:int(spec.shape[0]/2)]
        spec[0] = spec[0] / 2
        ampList.append(np.abs(spec))
        argList.append(np.angle(spec))
    
    #周波数は共通なので1回だけ計算(縦軸表示に使う)  
    freq = np.fft.fftfreq(data.shape[0], 1.0/sound.frame_rate)
    freq = freq[:int(freq.shape[0]/2)]
    
    #時間も共通なので1回だけ計算(横軸表示に使う)
    time = np.arange(0, i+1, 1) * s / sound.frame_rate
    
    #numpyの配列にしておく
    ampList = np.array(ampList)
    argList = np.array(argList)
    

    上記でスペクトログラムのデータを作成することができました。

    matplotlibとseabornでスペクトログラムを可視化

    次に実際に可視化をしていきます。

    In [21]:
    df_amp = pd.DataFrame(data=ampList, index=time, columns=freq)
    
    In [22]:
    #seabornのheatmapを使う
    plt.figure(figsize=(20, 6))
    sns.heatmap(data=np.log(df_amp.iloc[:, :100].T), 
                xticklabels=100, 
                yticklabels=10, 
                cmap=plt.cm.gist_rainbow_r,
               )
    
    Out[22]:
    <matplotlib.axes._subplots.AxesSubplot at 0x7f29b262c438>

    100Hz〜500Hzくらいが高いという、かなりはっきりした傾向がありますね。

    偏角の方も可視化してみます。

    In [23]:
    df_arg = pd.DataFrame(data=argList, index=time, columns=freq)
    
    In [24]:
    plt.figure(figsize=(20, 6))
    sns.heatmap(data=df_arg.iloc[:, :100].T, 
                xticklabels=100, 
                yticklabels=10, 
                cmap=plt.cm.gist_rainbow_r,
               )
    
    Out[24]:
    <matplotlib.axes._subplots.AxesSubplot at 0x7f29b26c8c88>

    こちらはほぼランダムに見えますね。それとも、スペクトルと何か相関があるのでしょうか。

    まとめ

    ちょっと今さから感もありますが、Pythonでサウンドスペクトログラムを作ってみました。引き続き、スペクトルと偏角に関係があるかどうか検討していこうと思います。


    コメントを残す

    メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

    CAPTCHA