PythonでWave音源の周波数と音量(デシベル)を取得する方法

前回に引き続き音声関連のプログラム紹介となります。
今回は、Pythonで、Wave音源の周波数と音量(デシベル)を取得する方法になります。
Audio機器なんかの試験で、再生されている音の周波数や音量を見たいことが筆者はあって、簡易的な試験で機材を安く少なくするために、こういったプログラムが欲しいなーってことがあるので試しにやってみました。
※ものによっては直接パソコンに機器からの音声を直接マイクラインに入力すると、過電流となりパソコンが壊れることがありますので、入力電圧を確認してから試してください。私は一度やってしまい煙が出ました。。。

サンプルコード
Wave音源の取得
周波数取得処理
音量(デシベル)の出力

サンプルコード

コードの全貌です。

import sys
import wave 
import numpy as np
import math

FLAME_SIZE  = 8192
    
def read_wavefile(file_path):
    wf        = wave.open(file_path, 'rb')
    frames   = wf.getnframes()      # フレーム数を取得
    Channels = wf.getnchannels()
    sampWidth = wf.getsampwidth()
    framerate = wf.getframerate()
    print( str(frames) + " Frames")
    print( str(Channels) + " channels")
    print( str(sampWidth*8) + " bit" )
    print( str(framerate) + " hz" )
    
    print("\n")
    wf.setpos(int(frames/2))   # 全フレームの真ん中に位置を移動
    buf = wf.readframes(FLAME_SIZE)  # セットした位置からFLAME_SIZEフレームを取得
    wf.close()
    return buf

def get_freq(data):
    data = np.frombuffer(data, dtype= "int16") / 32768.0
    spectrum = np.fft.fft(data)
    maxvalue = 0
    maxidx = 0
    tmpvalue = 0
    flist = np.fft.fftfreq(FLAME_SIZE, d=1.0/44100)             # 周波数リスト

    for j in range( int( len(spectrum) / 2 ) ):
        tmpvalue = spectrum[j]
        if tmpvalue > maxvalue:
            maxvalue = tmpvalue
            maxidx = j

    print( str(int(flist[maxidx])) + "Hz" )

def get_db(data):
    squaressum = 0
    # 累積二乗和を算出
    for i in range(FLAME_SIZE):
        squaressum += data[i] * data[i]
        
    # 平均平方根を算出
    meansquare = squaressum / (FLAME_SIZE/2)

    # 二乗平均平方を取得  (入出力信号レベル)
    rms = math.sqrt(meansquare)
    decibel = 20 * math.log10(rms)
    print( str( int(decibel) ) + "db" )

if __name__ == '__main__':
    args = sys.argv
    
    if len(args) != 1:
        #保存実行
        wave_data = read_wavefile( args[1] )
        get_freq(wave_data)
        get_db(wave_data)

今日はまず実行結果を見せます。
今回の測定に使用した音源データは以下の3つです。
↓のサイトからダウンロードさせていただきました。
http://www.op316.com/tubes/tips/wav2.htm

1khz-0db-30sec.wav
400hz-0db-30sec.wav
5khz-6db-20sec.wav

前回の↓で録音したデータを解析使用と思ったのですが、環境がPCのマイクしかなくノイズが入ってしまうので、ダウンロードさせた頂きました。

では実行結果です。それぞれ、周波数についてはそれっぽい値が取れています。
※すいません。デシベルについては、合っているのか不明です。。。
WindowsMediaPlayerで音源を再生して、iphoneのアプリで計測した感じではそれっぽい値になりました。

続いてはプログラムの解説ですが、解説は次ページにて記載します。

Wave音源の取得

音源の取得は以下の関数で行っています。
前回の書き込みの時と同様に引数で、開くファイルのパスを渡しています。
処理の流れとしてはファイルの制御と同じように、以下のような流れで処理を行います。

①ファイルのオープン
② データの取得
③ ファイルのクローズ

def read_wavefile(file_path):
    wf        = wave.open(file_path, 'rb')                   <- ファイルをオープン
    frames   = wf.getnframes()                               <- トータルのフレーム数を取得 
    Channels = wf.getnchannels()                             <- チャンネル数を取得
    sampWidth = wf.getsampwidth()                            <- データの幅取得
    framerate = wf.getframerate()                            <- フレームレートを取得
    print( str(frames) + " Frames")
    print( str(Channels) + " channels")
    print( str(sampWidth*8) + " bit" )
    print( str(framerate) + " hz" )
    print("\n")

    wf.setpos(int(frames/2))                                <- データを取得する位置を指定
    buf = wf.readframes(FLAME_SIZE)                         <- 指定サイズデータを取得
    wf.close()                                              <- ファイルをクローズ
    return buf

周波数取得処理

以下の関数で周波数を取得しています。
申し訳ありませんが、私はこういった計算処理がちょっと苦手で理解できていないので、ざっくり説明です。。。
あまり詳しくないので、コードで簡単に注釈しておきます。
詳細が知りたい方は、別の専門的なサイトを見てみてください。。。

def get_freq(data):
    data = np.frombuffer(data, dtype= "int16") / 32768.0   <- データが正規化
    spectrum = np.fft.fft(data)                            <- データのスペクトル?を取得
    maxvalue = 0
    maxidx = 0
    tmpvalue = 0
    flist = np.fft.fftfreq(FLAME_SIZE, d=1.0/44100)        <- 周波数のリストを生成

    for j in range( int( len(spectrum) / 2 ) ):
        tmpvalue = spectrum[j]
        if tmpvalue > maxvalue:                     <- 最も大きい値となっている部分を抽出
            maxvalue = tmpvalue          
            maxidx = j

    print( str(int(flist[maxidx])) + "Hz" )         <- 周波数を出力

私は、フーリエ変換についてよく理解できていませんが、各周波数帯を並べて、どこにデータが一番あるのかを見てそこがこの音声の周波数だ！って言っています。
うーん。原理が理解できない。。。

続いて音量の出力です。

音量(デシベル)の出力

以下で、音量を取得しています。
流れ的には以下の通りです。

①累積2乗和を算出
②平均平方根を算出
③2乗平均平方を取得 <-これが入力信号レベル

最後に、 ③で取った値のlog10を取って20をかければ、デシベルになるらしいです。

def get_db(data):
    squaressum = 0
    # 累積二乗和を算出
    for i in range(FLAME_SIZE):
        squaressum += data[i] * data[i]
        
    # 平均平方根を算出
    meansquare = squaressum / (FLAME_SIZE/2)

    # 二乗平均平方を取得  (入出力信号レベル)
    rms = math.sqrt(meansquare)
    decibel = 20 * math.log10(rms)
    print( str( int(decibel) ) + "db" )

私には、なぜこれでデシベルが出せるのか全く理解できませんが、これでだせます。

これで、音声の周波数と音量を取得できました。
うまく組み込めば長時間の無人試験なんかで、音がとまってないか・異音がなってないかといったことに活用できます。
私は、Pythonではありませんが、C# .net で同様の環境を構築し、無人の音声確認試験環境を構築し試験の効率化をしていたことがあります。

みなさんも一度やってみてください。

ちなみに２つ以下と組み合わせれば、リアルタイムで波形を出力しながら、周波数と音量を確認するなんてことも可能ですなので、やってみてください。

2021/12/21 追記
質問がありましたので、ステレオの音声をLRに分けて周波数を取得するサンプルプログラムを追記します。

デジタル音声では、LとRのデータは連続していますので、sampleWitdh単位でデータを分割すれば、LとRのチャンネルにデータを分けることができます。
今回の場合では、16Bit音声ですので、２Byte単位でデータを分離していきます。
そこまでできれば分離したデータをそれぞれ、単純にFFTすれば周波数が取得できます。

import sys
import wave 
import numpy as np
import math

FLAME_NUM  = int(8192)
DATA_NUM    = int(FLAME_NUM*4)  # Stereo 16Bit *  FLAME_NUM Frames
    
def read_wavefile(file_path):
    wf        = wave.open(file_path, 'rb')
    frames   = wf.getnframes()      # フレーム数を取得
    Channels = wf.getnchannels()
    sampWidth = wf.getsampwidth()
    framerate = wf.getframerate()
    print( str(frames) + " Frames")
    print( str(Channels) + " channels")
    print( str(sampWidth*8) + " bit" )
    print( str(framerate) + " hz" )
    
    print("\n")
    wf.setpos(int(frames/2))   # 全フレームの真ん中に位置を移動
    buf = wf.readframes(FLAME_NUM)  # セットした位置からFLAME_SIZEフレームを取得
    wf.close()
    
    lbuf = bytearray()
    rbuf = bytearray()
    #データをLChとRChに分離
    for i in range(0 ,DATA_NUM):
        if i%4 == 0 or i%4 == 1:
            lbuf.append(buf[i])
        if i%4 == 2 or i%4 == 3:
            rbuf.append(buf[i])
            
    return lbuf,rbuf

def get_freq(data):
    ndata = np.frombuffer(data, dtype= "int16") / len(data)
    spectrum = np.fft.fft(ndata)
    maxvalue = 0
    maxidx = 0
    tmpvalue = 0
    flist = np.fft.fftfreq(FLAME_NUM, d=1.0/44100)             # 周波数リスト
    
    for j in range( int( len(spectrum)/2 ) ):
        tmpvalue = spectrum[j]
        if tmpvalue > maxvalue:
            maxvalue = tmpvalue
            maxidx = j

    print( str(int(flist[maxidx])) + "Hz" )


if __name__ == '__main__':
    args = sys.argv
    
    if len(args) != 1:
        #保存実行
        
        wave_ldata,wave_rdata = read_wavefile( args[1] )
        get_freq(wave_ldata)
        get_freq(wave_rdata)