【python】自己相関関数を用いたピッチ検出【サウンドプログラミング】

こんばんは．1000chです．今週頭から内定先のインターンが始まり，今年はTAとして参加生のサポートをしています．自分が昨年参加したインターンに，このような形で携われるのはなかなか感慨深いものがあります．振り返るとこのブログを始めたきっかけもこのインターンでの出会いだったりして，時が経つのは早いものだなぁと思うばかりです．

さて，今回は音のピッチ検出プログラムを作ってみます．キーワードは自己相関関数．いきますよー．

ピッチ検出アルゴリズム

入力音声のピッチ検出には，様々なアルゴリズムが知られています．くわしくはおぐぐりいただくか，wikipedia先生に聞いてみてください．

今回は手始めということで，簡単そうな「自己相関関数」を利用したピッチ検出を試みます．

相関関数

相関関数r(m)は，ざっくり言うと「信号の周期性」を表現できる値です．信号をmサンプルずらしたとき，元の信号にどれくらい似通っているかを示す値です．相関関数はその導出法により自己相関関数，相互相関関数と分かれています．詳しくはこちらの3.5章を参考です．

ディジタル・サウンド処理入門―音のプログラミングとMATLAB(Octave・Scilab)における実際 (ディジタル信号処理シリーズ)

作者: 青木直史
出版社/メーカー: CQ出版
発売日: 2006/03
メディア: 単行本
購入: 2人クリック: 28回
この商品を含むブログ (10件) を見る

それぞれの定義式はこんな感じ．

違いは総和をとるindexの違いだけですね．自己相関関数の方はmが大きくなる程加算するサンプル数が減るため，小さくなる傾向になります．これはずらすサンプル数が増える程，元の信号から遠ざかる，というイメージに近いのです．

mサンプルずらした際．元信号に似た波形になる時，相関関数は大きくなります．これはすなわち．元信号が周囲mサンプルであることを示します．

つまり，相関関数r(m)を最大化するmを求めれば，信号の周期がわかり，そこから信号の周波数，ピッチを計算できるのです．原理は簡単ですね．

実装してみた

numpy/scipyでの相関関数の実装はこちらを参照．

correlation - numpy reference correlation - scipy reference

さしあたりscipyでの自己相関関数，相互相関関数の利用はこんな感じでできます．

import scipy.signal as sig

# dt:信号
# N:窓サイズ
# max_m:最大ずらすサンプル数

# 自己相関関数
def auto_correlate(dt):
    cor = sig.correlate(dt,dt,mode="full")
    return cor[cor.size/2:]

# 相互相関関数
def correlate(dt):
    return sig.correlate(dt[:N+max_m],dt[:N],mode="valid")

自己相関関数を利用し，信号の周期を求めるプログラムは次のようになります．なお，計算された自己相関関数の最大値を求める処理は少々愚直です．自己相関関数は全くずらさない（m=0）ときが最大のため，そこから下がり始め，再び上がり始めてからの最大値を求めています．

# coding:utf-8
# Last Change:2014_08_27_00:48:01.

import numpy as np
import scipy.signal as sig
import scipy.io.wavfile as scw
import matplotlib.pyplot as plt

# 自己相関関数
# r(m) = sum( x(n)*x(n+m) ) n=0~N-m-1, m=0~N-1
def auto_correlate(dt):
    cor = sig.correlate(dt,dt,mode="full")
    return cor[cor.size/2:]

# 初出peakの検出
# 愚直過ぎる..(´・ω・`)
def calc_first_peak_index(array):
    i = 0
    while array[i] > array[i+1]:
        i += 1
    while array[i] < array[i+1]:
        i += 1
    return i

if __name__ == '__main__':

    src_name = './tone_A4.wav' 

    fs, dt = scw.read(src_name)
    dt = dt / (2 ** 15)
    ts = 1 / fs
   
    # 自己相関係数の計算
    win_size = 2 ** 9 # 512
    N = 2 ** 9
    start = 1000
   
    cor = auto_correlate(dt[start:start+N])
    peak_index = calc_first_peak_index(cor)
    freq = 1 / (peak_index * ts)
    print(freq)

    plt.plot(np.arange(cor.size) * ts, cor)
    plt.show()