はじめに

信号処理の理論ばかり追っていてもどう役に立てれば良いのかわからず、なかなかモチベーションが上がりません。

今回は、信号処理の応用例として、人間の声のピッチ、つまり音程を自己相関関数で推定してみます。

音声データを用意する

なんでも良いので人の声の音声データをWaveファイルで用意します。
Scipyを使うと、Waveファイルを簡単に読み込めます。

import scipy.io.wavfile

fs, data = scipy.io.wavfile.read(wave_file_name)

読み込んだ結果はタプルで返ってきます。fsがサンプリングレート、dataがサンプルデータです。
少々面倒なのですが、モノラルチャンネルのWaveファイルの場合、dataは素直なnumpy.ndarrayによる一次元配列になっているものの、ステレオチャンネル以上のWaveファイルの場合、dataはnumpy.ndarrayの二次元配列になっています。

そこで私は、以下のような処理で、チャンネル数に依存せずに扱えるようにデータ構造を変換して扱っています。

import numpy

if type(data[0]) == numpy.ndarray:
    nch = len(data)
    lpcms = data
else:
    nch = 1
    lpcms = numpy.array([data])

lpcms[0]が１チャンネル目のサンプルデータ、lpcms[1]が２チャンネル目のサンプルデータ、という様な形でデータにアクセスできます。

スライド窓で分割

今回はDFT（離散フーリエ変換）による自己相関関数を使うので、指定した窓長で音声信号を切り出します。

offset = int(fs * 10.0)
window_size = fs // 50
window = lpcms[0][offset:offset+window_size]

サンプリングレートに秒をかけて、切り出し位置のオフセットサンプル数を計算しています。例では10秒の位置を指定しています。

次に、サンプリングレートを50で割って窓長を計算しています。窓長の決め方は重要です。人間が声で出せる音程の下限値は諸説ありますが、おそらく50Hzを下回る事はないはずなので、ここでは最低音程＝最長周期長として、サンプリングレートを50Hzで割った値を周期長として採用します。

ちなみに私が出せる低音程の下限値は約78Hzでした。

DFTによる解析窓長が短いと、以下のようなデメリットがあります。

• 周波数分解能が低い
• 低周波成分を解析できない

窓長を長くすればするほど、上記の問題を解消できますが、窓長を長くする事のデメリットもあります。

• 解析窓内に過渡現象が入り込みやすくなり解析結果がぼやける(特に高周波成分)
• アルゴリズム遅延が増える

特に窓長における解析精度の時間分解能と周波数分解能のトレードオフの関係を、不確定性原理と呼ぶこともあるようです。

また、自己相関関数をDFTで計算する場合、実装がFFT(高速フーリエ変換)であれば高速化アルゴリズムの都合により、窓長のサンプル数を2のべき乗にしないと計算できません。

よって、上記の例では省いていますが、実際には窓長を2のべき乗に切り上げる計算も必要になります。

雑に計算するならば、人間の声の音程を解析したい場合、44.1kHzや48kHzの音声に対して窓長は最低でも1024サンプルを必要とします。
そのため、サンプリングレートが44.1kHzや48kHzと既知であるならば、窓長を1024として構いません。

切り出した音声から音程を推定する

自己相関関数とパワースペクトラムの記事で書いたように、切り出した音声断片についてDFTによる自己相関関数を求めます。

import numpy
import scipy.fftpack

def acf_by_dft(values):
    power_spectrum = numpy.abs(scipy.fftpack.fft(values)) ** 2
    return numpy.real(scipy.fftpack.ifft(power_spectrum))

padded = numpy.r_[window, numpy.zero(len(window))]
acf = acf_by_dft(padded)[:len(window)]

ここで、切り出した窓と同じ長さのゼロ値を連結して自己相関関数を計算しています。ゼロ値を挿入すると、入力信号の後半半分は全てゼロなので、自己相関関数の値の大きさは時刻0をピークに段々と減少していき、最終的にはゼロになります。

このゼロ値連結の小細工によって、音程の周波数推定アルゴリズムが単純になります。

自己相関関数acfの後ろ半分は前半分の対象形で出てくるだけなので、この後の推定処理を簡単にするために捨ててしまいます。

人間の声の音程は、基本周波数と呼ばれる一番低い周波数ピークに現れます。自己相関関数に置き換えると、最初の山のピーク位置が、基本周波数に相当する周期長です。

周期的な信号の自己相関関数におけるピーク（山）は、時刻をずらしていった際に支配的な周期のタイミングで現れます。そして、その手前には半周期ずれのタイミングで相関係数が負のピーク（谷）となる、つまり最小となるタイミングが存在します。

一つ目の山を探索する、という処理は難しいですが、値全体から最小値を見つける、という問題は比較的簡単です。そこで、まず最小値を探索し、それ以降の最大値を探索することで、一つ目の山を見つけ出し、周期を推定します。

min_index = numpy.argmin(acf)
max_index = numpy.argmax(acf[min_index:]) + min_index
basic_cycle_length = max_index
basic_frequency = fs / basic_cycle_length
pitch_reliability = acf[basic_cycle_length] / acf[0]

自己相関関数は、自分自身の信号をずらしながら相関係数を計算した結果です。つまり、時間の関数です。ここでの時間はサンプル数のズレを意味しているので、上記のbasic_cycle_lengthはそのままサンプル数単位の周期長を表しています。

あとは、元のWaveファイルのサンプリングレートから、基本周波数basic_frequencyを求めるだけです。簡単ですね。

最後の自己相関関数の正規化は、基本周期長のタイミングでの自己相関値を相互相関係数化する計算です。相互相関係数に変換することで、値は-1 \sim 1の範囲に必ず収まります。

ノイズ状の信号では、前述の方法で機械的に求めた基本周期における相互相関係数が低い値になると期待できるので、この値に対して閾値処理をすることで、人の声かどうかをある程度弁別可能になります。

応用の際の注意点

本記事で紹介した方法は非常にシンプルで実装も容易ですが、実際に応用する際には注意が必要です。

• 基本周期の精度が整数精度なので正確な基本周期とのずれがある
• 低周波（長周期）の分解能は高いが高周波（短周期）の分解能は低い
• 窓内だけを見た自己相関関数なので切り出し方によっては長周期由来の自己相関値が上がりにくい
• 相互相関係数化した自己相関関数値による人の声判別だけでは不十分

人の声の音声信号だけに解析対象をしぼっても、実に様々な信号が存在するので、実践する際はデータを集めてアルゴリズムを適用し、十分にモデル化できていないデータを見つけて、何が足りていないのかを良く考察する必要があります。

ちなみに、ピッチ（pitch）はもともと長さの間隔を意味するので、ピッチと言った場合には「周波数」ではなく「周期」を連想するのが正しい気がします。