その中で「自己相関関数」というものが紹介されており，非常に面白そうだったので実装してみることにします．
1.自己相関関数とは
「あー」というような定常な音の場合，同じような音の波形が繰り返されています．
もっとも，実際の音声では完全にイコールにはならないと思いますが，概ねこのような性質を持っています．こういった関数を周期Tの周期関数と呼びます．例えばサイン関数は周期2πの周期関数です．
この周期Tを求める手法として上記の本では「自己相関関数」を用いています．関数fに対する自己相関関数R_ffを，今回次のように定義します．
R_{ff}(\tau) = \frac{(f(t),f(t+\tau))}{||f(t)||\,\,||f(t+\tau)||}\\
ただし，実装するときにどうするかをはっきりしたいので，関数の内積とノルムは次のように定義します．
(f(t),f(t+\tau)) = \int_{t_0}^{t_1} f(t)f(t+\tau)dt \\
||f(t)|| = \sqrt{\int_{t_0}^{t_1}f(t)^2 dt }\\
||f(t+\tau)|| =\sqrt{ \int_{t_0}^{t_1}f(t+\tau)^2 dt} \\
関数f(t)とf(t+τ)の相関係数みたいなものなので，関数のノルムで割るかどうかや，関数の平均値を減じるかどうかは人それぞれだと思います．
さてこの自己相関関数R_ffを実際に計算していくとこのようになります．
τを変化していき，τ=nT ( n = 0,1,2,... )となる点で似たパターンとなるため，R_ffが極大値を取ることがわかります．ブログ用．自己相関関数の値の変化． pic.twitter.com/Qmul8V82lH— もりとにー (@TonyMooori) 2016年3月7日
2.抽出アルゴリズム
ソースコードを貼る前に，実際にプログラムを組んで「あー」という音を読み込ませて，その音声に対する自己相関関数を求めるとこんな感じになります．
恐らくフラクタル図形のように「あー」という音の中にも似た部分があるようです．それはそれで面白いのですが……．
まぁ自分の目で見ればどれが目的の周期なのか一目瞭然なのですが，できれば自動的に判定したいところです．
とりあえずFFTで高周波成分を除去するハイパスフィルタを実装するとこんな感じになります．緑色の部分がハイパスフィルタで高周波成分を除去したもの．移動平均でも良いと思います．
その周期で音声の波形を10個ほど取り出すとこんな感じになった．
この繰り返されている波形データだけで音声認識とかできないか考え中なところです．
3.プログラム
コメントは多めにしたつもりではありますがうまくいかない場合はコメントを参考にしつつ頑張って下さい．
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #coding:utf-8 | |
| import numpy as np | |
| import pandas as pd | |
| import matplotlib.pyplot as plt | |
| import wave | |
| from scipy import signal | |
| """ | |
| 参考 | |
| [1] 波形を見る - 人工知能に関する断創録 http://aidiary.hatenablog.com/entry/20110519/1305808715 | |
| [2] 雨宮好文 監修, 佐藤幸男 著『信号処理入門』 オーム社 | |
| """ | |
| def read_wave(path): | |
| """ waveファイルを読み込む関数 """ | |
| wf = wave.open(path , "r" ) | |
| fs = wf.getframerate() # サンプリング周波数 | |
| x = wf.readframes(wf.getnframes()) | |
| x = np.frombuffer(x, dtype= "int16") / 32768.0 # -1 - +1に正規化 | |
| wf.close() | |
| # fsはサンプリング周波数 | |
| return (x,fs) | |
| def HPF(x,n_dim): | |
| """ ハイパスフィルタ，n_dimは打ち切る周波数 """ | |
| c = np.fft.rfft(x) # 高速フーリエ変換 | |
| c[n_dim:] = 0.0 # 高周波成分除去 | |
| return np.fft.irfft(c) | |
| def cross_correlation(v1,v2,low,high,tau): | |
| """ ベクトルv1,v2の相互相関関数 """ | |
| # v1の要素番号lowからhightまで | |
| # v2の要素番号low+tauからhigh+tauまでを取り出す | |
| x1 = v1[low:high] | |
| x2 = v1[(low+tau):(high+tau)] | |
| # ノルムを1にする | |
| x1 = x1 / np.linalg.norm(x1) | |
| x2 = x2 / np.linalg.norm(x2) | |
| # x1,x2の内積が相関係数となる | |
| return x1.T.dot(x2) | |
| if __name__=="__main__": | |
| # データ読み取り | |
| print("reading voice.wav") | |
| x,fs = read_wave("./voice.wav") | |
| # 自己相関関数を計算 | |
| print("calculating R_ff") | |
| R = np.array([ cross_correlation(x,x,0,fs,tau) for tau in range(fs) ]) | |
| # 周期を計算してく | |
| print("calculating T") | |
| # まずハイパスフィルタで高周波成分を除去 | |
| R_HPF = HPF(R,len(R) / 20) | |
| # 自己相関関数はτが小さいほうが相関が高いため，τが小さい部分を取り出す | |
| R_HPF[:fs/5] | |
| # 極大値の要素番号を取得 | |
| local_max, = signal.argrelmax(R_HPF) | |
| # 極大値の中で0.2よりも大きい部分の要素番号を取得 | |
| threshold = 0.2 # しきい値 | |
| local_max = local_max[ R_HPF[ local_max ] > threshold ] | |
| # この時 local_max は(うまく行けば) | |
| # array([ 5, 646, 1291, 1938, 2582, 3228, 3872, 4520]) | |
| # のような等差数列になっているので各項の差を取る | |
| diff = np.diff(local_max) | |
| # diffは(うまく行けば)こうなっている | |
| # array([641, 645, 647, 644, 646, 644, 648]) | |
| # 外れ値があるかもしれないので，平均値ではなく中央値を使うと良い | |
| T = (int)(np.median(diff)) | |
| # 周期Tで取り出して表示 | |
| for i in range(10): | |
| plt.plot(x[(i*T):(i*T+T)]) | |
| plt.title("T = %d" % T ) | |
| plt.show() | |
