WAVEファイルの読み込みと表示 (2018.4.29, 2020.10.4追記)

2018.5.5, 2018.5.20, 2019.2.15追記, 2020.10.4追記

Summary

PySoundFileやscipyでwavファイルを読み込んでmatplotlibでスペクトログラムを表示するサンプル。scipyは量子化ビット数16ビットまでのファイルしか読み込めない制約があるので、PySoundFileを推奨。スペクトログラムの色味をsoxのようにカッコ良くするにはmatplotlib.pyplot.clim(vmin, vmax)を使って表示される値のレンジを限定する必要がある。 importで指定されたpackageをインストールするには、PyCharmのメニューのFile –> Default Preferences –> Project Interpreterの画面で下部の+を押して検索して追加する。(Anacondaを使わずにPyCharmをインストールする方法の詳細はmacOSでPyCharmを使う(2018.4.30)に移動) “import soundfile as sf”では“PySoundFile”パッケージを使う。pysndfileやSoundFileではないので注意

2020.10.4追記 matplotlib関連のアップデートをしたtest6.pyが最新版です。そちらを参考にしてください。

WAVE readのサンプル

Pythonで音響信号処理というサイトを参考にしてMacでPythonを使ってWAVEファイルを読み込むサンプルコードを試してみた。

  1. test1.py: PySoundFile (import soundfile)というパッケージを使ってwavファイルを読み込む。 もとのサイトのソースでは、サンプル数がfftフレームをフレーム長の半分ずつシフトして処理していて、ちょうど端数サンプルなしに処理できるサンプル数でないと正しく処理できていなかった。最後に処理するフレームがfftサイズ未満の場合には、0を足してfftサイズにしてから処理するように修正した。

downloaded package

  1. test2.py: scipy.io.wavfileでもwavファイルを読み込めるが、量子化ビット数16ビットにまでしか対応していない。matplotlibで直接スペクトログラムの計算もできる。

downloaded package

  1. test3.py: 上記の例のmatplotlibの表示パラメータを ’cmap=plt.cm.gist_heat’に変更。24ビット入力に対応したPySondFileを使うように変更。

downloaded package

  1. test4.py: soxのスペクトログラムっぽく表示してみた。[2018.5.20追記]
    しかし、デフォルトのカラーマップでは色味が明るすぎていまひとつかっこよくない。

downloaded package

ちなみに、soxでもスペクトログラムを出力できる。

$ sox filename -n spectrogram

downloaded package

  1. test5.py: スペクトログラムのカラーマップの表示スケールを調整。[2019.2.15追記]
    matplotlibのspecgramで波形から直接スペクトログラムを計算し表示できる。mode=‘default’ (psd)では 10 log 10になっている。それ以外のmode=’magnitude’などでは20 log 10で表示される。 colormapで色指定するが、デフォルトの表示では色味がかっこよくない。 matplotlib.pyplot.clim(vmin, vmax)を使って表示される値のレンジを限定すればsoxっぽい見かけにできる。（Thanks 小泉くん）

downloaded package

  1. test6.py: plt系を卒業してfig, axをメインで使うように書き直し、さらにfig.add_subplot()からfig.add_axes()に変更。ticker.FuncFormatterを使用。[2020.10.4追記]
    matplotlibの使い方を勉強し、よりsoxに近い表示ができるようになった。この版と参考リンクを理解すれば、かなりmatplotlibを使いこなせるようになるはず。

downloaded package downloaded package

ソースコード

PyCharmでプロジェクトを作成して、必要なパッケージを登録。 （登録方法は後述）

downloaded package

test1.py

sys, math, numpy, scipy.fftpack, matplotlib.pyplot


#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# WAV read test
# https://qiita.com/wrist/items/5759f894303e4364ebfd

import sys
import math
import numpy as np
import scipy.fftpack as fft
import matplotlib.pyplot as plt

import soundfile as sf

if __name__ == '__main__':
    plt.close("all")

    # wavファイル読み込み
    filename = sys.argv[1]
    wav, fs = sf.read(filename)
    print("filename: " + filename)

    # ステレオ2chの場合、LchとRchに分割
    wav_l = wav[:, 0]
    wav_r = wav[:, 1]

    # 入力をモノラル化
    xs = (0.5 * wav_l) + (0.5 * wav_r)

    n_len = len(xs)
    n_fft: int = 128
    n_overlap = 2
    n_shift: int = math.ceil(n_fft / n_overlap)
    #print(n_shift, n_len, n_len // n_fft, n_len % n_fft)

    # 中間バッファ
    zs = np.zeros(n_len)
    Zs = np.zeros(n_fft)

    # 出力バッファ
    ys = np.zeros(n_len)

    # 窓関数
    window = np.hanning(n_fft)
    #print(window)

    # FFT & IFFT
    for start in range(0, n_len - n_shift, n_shift):
        #print(start, min(n_fft, n_len - start))

        # when n_len is multiple of n_fft
        #xs_cut = xs[start: start + n_fft]
        # if n_fft <= (n_len - start):
        #     xs_win = xs_cut * window
        # #print(xs_win)
        #
        #     Xs = fft.fft(xs_win, n_fft)
        #
        #     # some signal processing
        #     Zs = Xs
        #     zs = fft.ifft(Zs, n_fft)
        #
        #     # write output buffer
        #     ys[start: start + n_fft] += np.real(zs)

        # solved issues when n_len is not multiple of n_fft
        if n_fft <= (n_len - start):
            xs_cut = xs[start: start + n_fft]
        else:
            xs_cut = np.append(xs[start: n_len], np.zeros(n_fft-(n_len-start)))

        xs_win = xs_cut * window
        #print(xs_win)

        Xs = fft.fft(xs_win, n_fft)

        # some signal processing
        Zs = Xs
        zs = fft.ifft(Zs, n_fft)

        # write output buffer
        if n_fft <= (n_len - start):
            ys[start: start + n_fft] += np.real(zs)
        else:
            ys[start: n_len] += np.real(zs[:(n_len-start)])

    # 冒頭から10秒分プロット
    fig = plt.figure(1, figsize=(8, 10))
    ax = fig.add_subplot(211)
    ax.plot(xs[:fs*10])
    ax.set_title("input signal (" + filename + ")")
    ax.set_xlabel("time [pt]")
    ax.set_ylabel("amplitude")

    ax = fig.add_subplot(212)
    ax.plot(ys[:fs*10])
    ax.set_title("output signal (" + filename + ")")
    ax.set_xlabel("time [pt]")
    ax.set_ylabel("amplitude")

    plt.show()

test2.py

sys, scipy.io.wavfile, matplotlib.pyplot


#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# WAV read test
# http://ism1000ch.hatenablog.com/entry/2014/07/16/022304

import sys
import scipy.io.wavfile as wio
import matplotlib.pyplot as plt

if __name__ == '__main__':
    plt.close("all")

    # wavファイル読み込み
    filename = sys.argv[1]
    rate, data = wio.read(filename)
    print("filename: " + filename)

    # ステレオ2chの場合、LchとRchに分割
    wav_l = data[:, 0]
    wav_r = data[:, 1]

    # 入力をモノラル化
    xs = (0.5 * wav_l) + (0.5 * wav_r)

    pxx, freq, bins, t = plt.specgram(xs, Fs=rate)
    plt.show()

test3.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# WAV read test
# http://ism1000ch.hatenablog.com/entry/2014/07/16/022304
# https://qiita.com/wrist/items/5759f894303e4364ebfd

import sys
#import scipy.io.wavfile as wio
import matplotlib.pyplot as plt

import soundfile as sf


if __name__ == '__main__':
    plt.close("all")

    # wavファイル読み込み
    filename = sys.argv[1]
    wav, fs = sf.read(filename)
    print("filename: " + filename)

    # ステレオ2chの場合、LchとRchに分割
    wav_l = wav[:, 0]
    wav_r = wav[:, 1]

    # 入力をモノラル化
    xs = (0.5 * wav_l) + (0.5 * wav_r)

    n_fft: int = 128
    n_overlap = 2
    n_shift: int = n_fft // n_overlap

    pxx, freq, bins, t = plt.specgram(xs, NFFT=n_fft, Fs=fs, noverlap=n_shift, cmap=plt.cm.gist_heat)
    plt.show()

test4.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# WAV read test
# http://ism1000ch.hatenablog.com/entry/2014/07/16/022304

import sys
#import scipy.io.wavfile as wio
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

import soundfile as sf


if __name__ == '__main__':
    plt.close("all")

    # wavファイル読み込み
    filename = sys.argv[1]
    wav, fs = sf.read(filename)
    print("filename: " + filename)

    # ステレオ2chの場合、LchとRchに分割
    wav_l = wav[:, 0]
    wav_r = wav[:, 1]

    # 入力をモノラル化
    #xs = (0.5 * wav_l) + (0.5 * wav_r)

    n_fft: int = 128
    n_overlap = 2
    n_shift: int = n_fft // n_overlap

    #c_map = plt.cm.gist_heat
    #c_map = plt.cm.gnuplot
    c_map = plt.cm.CMRmap


# spectrogram 2
    plt.style.use('dark_background')

    plt.subplot(211)
    plt.subplots_adjust(hspace=0)
    #plt.xlabel('time (s)')
    plt.ylabel('freq')

    plt.tick_params(labelbottom='off', labeltop='on')
    plt.locator_params(labelbottom='off', labeltop='on')
    #plt.yscale('log')
    pxx, freq, bins, t = plt.specgram(wav_l, NFFT=n_fft, Fs=fs, noverlap=n_shift, cmap=c_map, scale='dB')

    plt.title('avemaria.wav', y=1.12)

    plt.subplot(212)
    plt.xlabel('time (s)')
    plt.ylabel('freq')
    pxx, freq, bins, t = plt.specgram(wav_r, NFFT=n_fft, Fs=fs, noverlap=n_shift, cmap=c_map, scale='dB')


    plt.subplots_adjust(bottom=0.1, right=0.8, top=0.9)
    cax = plt.axes([0.83, 0.1, 0.03, 0.8])
    plt.colorbar(cax=cax)

    #plt.tight_layout()
    #plt.text(0.1, 1.08, 'avemaria.wav')

    plt.savefig("./hist_test4.png")
    plt.show()

test5.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# WAV read test
# http://ism1000ch.hatenablog.com/entry/2014/07/16/022304

import sys
#import scipy.io.wavfile as wio
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

import soundfile as sf


if __name__ == '__main__':
    plt.close("all")

    # wavファイル読み込み
    filename = sys.argv[1]
    wav, fs = sf.read(filename)
    print("filename: " + filename)

    # ステレオ2chの場合、LchとRchに分割
    wav_l = wav[:, 0]
    wav_r = wav[:, 1]

    # 入力をモノラル化
    #xs = (0.5 * wav_l) + (0.5 * wav_r)

    n_fft: int = 4096
    n_overlap = 2
    n_shift: int = n_fft // n_overlap

    # plt.clim()でspectrogram表示の色味をsoxっぽく
    c_vmin: int = -120
    c_vmax: int = -30

    #c_map = plt.cm.gist_heat
    #c_map = plt.cm.gnuplot
    c_map = plt.cm.CMRmap


# spectrogram 2
    plt.style.use('dark_background')

    plt.subplot(211)
    plt.subplots_adjust(hspace=0)
    #plt.xlabel('time (s)')
    plt.ylabel('freq')

    plt.tick_params(labelbottom='off', labeltop='on')
    plt.locator_params(labelbottom='off', labeltop='on')
    #plt.yscale('log')
    pxx, freq, bins, t = plt.specgram(
        wav_l, NFFT=n_fft, Fs=fs, noverlap=n_shift, cmap=c_map, mode='magnitude', scale='dB')
    plt.clim(c_vmin, c_vmax)

    plt.title('avemaria.wav', y=1.12)

    plt.subplot(212)
    plt.xlabel('time (s)')
    plt.ylabel('freq')
    pxx, freq, bins, t = plt.specgram(
        wav_r, NFFT=n_fft, Fs=fs, noverlap=n_shift, cmap=c_map, mode='magnitude', scale='dB')
    plt.clim(c_vmin, c_vmax)


    plt.subplots_adjust(bottom=0.1, right=0.8, top=0.9)
    cax = plt.axes([0.83, 0.1, 0.03, 0.8])
    plt.colorbar(cax=cax)

    #plt.tight_layout()
    #plt.text(0.1, 1.08, 'avemaria.wav')

    plt.savefig("./hist_test5.png")
    plt.show()

test6.py, .vscode/launch.json

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# WAV read test
# http://ism1000ch.hatenablog.com/entry/2014/07/16/022304
# https://takeshid.hatenadiary.jp/entry/2016/01/10/153503

# matplotlib
# 早く知っておきたかったmatplotlibの基礎知識、あるいは見た目の調整が捗るArtistの話(https://qiita.com/skotaro/items/08dc0b8c5704c94eafb9)
# matplotlibのcolorbarを解剖してわかったこと、あるいはもうcolorbar調整に苦労したくない人に捧げる話(https://qiita.com/skotaro/items/01d66a8c9902a766a2c0)
# matplotlibのめっちゃまとめ(https://qiita.com/nkay/items/d1eb91e33b9d6469ef51)
# [Python] matplotlib: 論文用に図の体裁を整える(https://qiita.com/qsnsr123/items/325d21621cfe9e553c17)
# matplotlibでジャーナルに投稿可能な図を作るためのメモ(https://qiita.com/HajimeKawahara/items/03f29367744d5209be42)
# matplotlibのグラフを高解像度で保存する(https://analytics-note.xyz/programming/matplotlib-save-dpi/)
# matplotlibで複数のグラフを表示するときのサイズを揃える(https://qiita.com/sbajsbf/items/b3ce138de83362bc45b0)
# [Matplotlib] タイトルの設定(https://python.atelierkobato.com/title/)
# matplotlib Contour Demo(https://matplotlib.org/gallery/images_contours_and_fields/contour_demo.html#sphx-glr-gallery-images-contours-and-fields-contour-demo-py) 
# seabornの細かい見た目調整をあきらめない(https://qiita.com/skotaro/items/7fee4dd35c6d42e0ebae)

import sys
import os
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker      # FormatterとLoactorを使うためにTickerモジュールをimport

import soundfile as sf
#import scipy.io.wavfile as wio

# 独自のFormatterの指定（縦軸をkHzにしたい）
@ticker.FuncFormatter
def major_formatter_khz(y, pos):
    return '{:.0f}'.format(y/1000)      # 縦軸はkHzにしたいのでyの値を1000で割って整数部のみ表示

if __name__ == '__main__':
    plt.close("all")

    # wavファイル読み込み
    filename = sys.argv[1]                  # filename = "./avemaria.wav"
    wav, fs = sf.read(filename)
    print("filename: " + filename)

    sup_title = os.path.basename(filename)  # 正味のファイル名を取り出してタイトルに設定したい
    outfig_dpi = 350                        # pdf出力の解像度を指定

    # ステレオ2chの場合、LchとRchに分割
    wav_l = wav[:, 0]
    wav_r = wav[:, 1]

    # 入力をモノラル化
    #xs = (0.5 * wav_l) + (0.5 * wav_r)

    n_fft: int = 4096
    n_overlap = 2
    n_shift: int = n_fft // n_overlap

    # plt.clim()でspectrogram表示の色味をsoxっぽく
    c_vmin: int = -120
    c_vmax: int = -30

    #c_map = plt.cm.gist_heat
    #c_map = plt.cm.gnuplot
    c_map = plt.cm.CMRmap

# spectrogram v2.1
# Object orientedなfig, axを使うお作法を試す

    plt.style.use('dark_background')
    fig = plt.figure()              # Figureオブジェクトを作成

# formatterとlocatorを使う
    #ax1 = fig.add_subplot(2,1,1)                   # 通常はfig.add_subplot()を使うことが多いがadd_axes()を使う
    ax1 = fig.add_axes((0.13, 0.5, 0.7, 0.4))       # figに属するAxesオブジェクトを作成（左チャネル用）
    pxx, freq, bins, t = plt.specgram(
        wav_l, NFFT=n_fft, Fs=fs, noverlap=n_shift, cmap=c_map, mode='magnitude', scale='dB')
    plt.clim(c_vmin, c_vmax)

    ax1.set_ylabel('freq (kHz)')
    ax1.set_ylim(1,fs/2)          # 上グラフの縦軸の0と下グラフの最大メモリの数字が被るので1から表示
    ax1.tick_params(top=True, bottom=False, labelbottom=False, labeltop=True)
    ax1.yaxis.set_major_locator(ticker.MultipleLocator(4000))   # 4000Hzごとに周波数を書く
    ax1.yaxis.set_major_formatter(major_formatter_khz)          # 値を1000で割ってkHzにするformatterを呼ぶ

    #ax2 = fig.add_subplot(2,1,2)                   # 通常はfig.add_subplot()を使うことが多いがadd_axes()を使う
    ax2 = fig.add_axes((0.13, 0.1, 0.7, 0.4))       # figに属するAxesオブジェクトをもう一つ作成（右チャネル用）
    pxx, freq, bins, t = plt.specgram(
        wav_r, NFFT=n_fft, Fs=fs, noverlap=n_shift, cmap=c_map, mode='magnitude', scale='dB')
    plt.clim(c_vmin, c_vmax)

    ax2.margins(y=None, tight=True)
    ax2.set_xlabel('time (s)')
    ax2.set_ylabel('freq (kHz)')
    ax2.tick_params(top=False, bottom=True, labelbottom=True, labeltop=False)
    ax2.yaxis.set_major_locator(ticker.MultipleLocator(4000))   # 4000Hzごとに周波数を書く
    ax2.yaxis.set_major_formatter(major_formatter_khz)          # 値を1000で割ってkHzにするformatterを呼ぶ

    fig.suptitle(sup_title, y=0.99)

    fig.subplots_adjust(bottom=0.1, right=0.82, top=0.9)
    cax = plt.axes([0.85, 0.1, 0.03, 0.8])
    cb = plt.colorbar(cax=cax)
    
# できたFigureを保存してさらに画面表示
    plt.savefig("./hist_test6.pdf", dpi=outfig_dpi)             # 解像度dpiを指定して出力
    plt.show()

{
    // Use IntelliSense to learn about possible attributes.
    // Hover to view descriptions of existing attributes.
    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
    "version": "0.2.0",
    "configurations": [
        {
            "name": "Python: Current File",
            "type": "python",
            "request": "launch",
            "program": "${file}",
            "console": "integratedTerminal"
            "args":[
                "./avemaria.wav"
            ]
        }
    ]
}

PyCharmでのパッケージ登録手順

importで指定されたpackageをインストールするには、PyCharmのメニューのFile –> Default Preferences –> Project Interpreterの画面で下部の+を押して検索して追加する。 “import soundfile as sf”では“PySoundFile”パッケージを使う。pysndfileやSoundFileではないので注意

downloaded package

downloaded package

downloaded package

downloaded package

downloaded package

downloaded package

downloaded package

downloaded package

参考

もともとの参考リンク

matplotlib関係の神リンク

Back to Index