１５．自分の音痴度を測ってみよう

――音楽の美は、その一瞬の短さにおいて生命に似ている。――

三島 由紀夫

自分の音痴度を測ってみよう

秋口に差し掛かり、けっこう外の陽気も北欧らしくなってきました。日本は秋の風物詩の台風が例年にない甚大な被害をもたらしているようですね。欧州で暮らしてみて、地震や台風といった災害がない国がとてもうらやましいと思いました。

さて、前にラストステージとかなんとか言っておきながら、「Pythonで音声処理」の記事をまた書きます。別に求められてないでしょうけど、この前受けた（地獄の）集中講義で覚えたアルゴリズムの中に面白いものがあったので、一つピックアップしたいと思います。

ちなみに、この写真はこの前遊びに行ったデンマーク人のお友達（御年88歳のおばあさん）のおうちです。妻が毎週通っている裁縫教室でお友達になったので、私もおうち訪問に誘ってもらいました。

ピッチ推定（Pitch estimation）アルゴリズム

以前、第７回の記事では、ギターの音には固有周波数（Natural Frequency）と呼ばれる量が存在し、その整数倍の音を作り出して重ね合わせれば人工的にギター風の音を生成できることをご説明しました。

keep-learning.hatenablog.jp

今回やりたいことは、まさにこの逆です。すなわち、固有周波数を出発点として、ギター音のような音声データを生成した前回の記事に対し、今回は音声データを出発点として、そこに含まれる音の固有周波数がどの程度なのか推定（Estimate）することを考えます。

音大生などであれば親しみがあると思いますが、人間の声などの固有周波数は音楽用語でピッチ（Pitch）と呼びます。このピッチ推定（Pitch estimation）というタスクは、極めてベーシックな命題であるにも関わらずたいへん奥が深く、あらゆる条件で完璧にピッチを抽出できるようなアルゴリズムは未だに存在しません*1。

Harmonic model（高調波モデル）

ピッチ推定の方法は、たいへん古くから研究されていることもあり、音波の自己相関関数を用いる方法や離散フーリエ変換によるピリオドグラム法など、様々な方法が考案されています。ここでは、直感的に分かりやすく、実装もさほど大変ではないHarmonic modelについて紹介したいと思います。

今、サンプリング周波数で録音された長さNの音声データが手元にあったとします。ピッチ（固有周波数）はようするに周期の逆数ですから、理想的な周期性を持つ音声データは、以下のようなL次の三角関数からなる高調波（Harmonic）の重ね合わせで記述できると推測されます。

しかし、現実の音声データが本当に三角関数だけで記述できることは稀であり、実際はピッチに関係ないノイズなどの誤差（Error）が加わります。すなわち、

のように記述されるということです。ここで、表現を簡便にするために、無次元角周波数を導入して、

この形を仮定したピッチ推定手法のことを、Harmonic model（高調波モデル）と呼びます。シグマ和の範囲を規定するパラメータLは、理想とする周期波の複雑さを表すパラメータであり、音声データの性質に応じて3～10くらいの値が用いられます。

Harmonic modelの行列表記

以上で数学的には十分な定義が与えられました。しかし、プログラムで実装することを考えると、こうした数式は行列積の形で記述しておくと便利*2です。そこで、一見遠回りするようにも見えますが、以下のような行数N、列数2Lの行列を定義します。

さらに、三角関数の係数を規定する長さ2Lのベクトルも定義します。

これらを用いることで、Harmonic modelの定義式は以下のようにきわめてシンプルな行列方程式の形で記述することができます。

なお、太字になったとは、各要素が音声データの各時刻（サンプリング点）に対応する長さNのベクトルです。

最小二乗誤差

モデルを仮定したので、次は三角関数の係数と肝心のピッチを表す無次元角周波数をどうやって推定するのか考えなければなりません。

そこで、直感的ではありますが、実際の音声データと理想的な三角関数との誤差であるの二乗和を最小にするのはどうでしょうか。音声データに含まれるピッチを推定するうえで、理想との誤差が一番小さい解を探すことには正当性がありそうです（最尤推定）。すなわち、

を最小化するようにとを決めます。これがHarmonic modelで定式化された命題です。

一般に、このような一次線形結合のノルム最小化問題は、数学の世界では線形回帰問題と呼ばれており、以下の正規方程式（Normal equation）と完全に等価であることが知られています。

この式で登場するは累乗ではなく、行列の転置記号です。をこのように選べば、誤差の二乗和が最小になります。この方程式の解を元の二乗誤差に代入すればの成分が消去でき、のみに依存する式になります。

ここで、ノルムの二乗はベクトルに左からそのベクトルの転置を掛け算したものですから、整理すると以下の表式になります。

第一項はに依存しません。また、第二項にはという負の係数が付いているため、「二乗誤差を最小化するようなを見つける」という元の命題は、「を外した第二項を最大化するようなを見つける」という命題に置き換えることができます。これを数学的に表現すると、

と記述されます。は、「右に書かれた関数を最大化するような引数（Argument）を以上以下の範囲内で求めよ」という命題を表すための数学記号で、何か特別な計算をしているわけではありません。

Harmonic summation近似

ここまで順調に定式化できたので、プログラムに実装することは容易です。やるべきことは、上で定義したを行列の形で構築し、ピッチ推定を行いたい音声データと一緒に上述した式、

に当てはめて、行列計算ライブラリなどで計算させます。これを様々なに対して実行し、最も大きい値を与えるをピッチの推定値として出力すれば終わりです。

しかしながら、言うは易く行うは難し。実はこのHarmonic modelの計算、非常に計算負荷が高いことが知られています。特に逆行列を計算するところの負荷が重く、10秒程度の音声データに実行した場合でも数時間程度の計算時間がかかります。コンピュータサイエンスの世界では、このように原理的・数学的に可能であっても、コンピュータの処理量の関係で現実的（Practical）に不可能になるケースが多々あります。

どうにかして、この問題を私のヘボPCでも数分程度の計算時間で終了するくらいにスケールダウンできないでしょうか。実は、音声データの長さNが周期（ピッチの逆数）に比べて十分に大きい場合には、三角関数の直交性から以下のように近似できます。

この近似式を使うことで、計算負荷の高い逆行列の部分が単位行列になって消えてしまうので、係数を取り払った目的関数（最大値を求める対象の関数）は驚くほど簡単になります。

このNが大きい場合のみ近似的に成立するモデルをHarmonic summation modelと呼び、その計算負荷の軽さからピッチ推定のベースラインモデルとして広く用いられています。

Pythonで実装

例のごとく、いきなりですがソースコード載せます。おもしろいのはプログラムの実践部分なので、興味のない方はソースコードを読み流して次の節に進んでください。まずはピッチ推定のための行列計算を行うサブ関数から紹介します。

### Imports ###
import sys
import math
import time
import soundfile as sf
import IPython.display as ipd
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal as sg

### For estimating the pitch based on harmonic summation model ###
def get_fundamental_freq(wavdata, Fs, f0_min, f0_max, resolution, L_max):
    
    #Parameter check
    N = len(wavdata)
    if f0_min < Fs / N:
        print('The f_min should be more than the frequency resolution(Fs / N).')
        sys.exit()
    if f0_max > Fs:
        print('The f_max should be less than sampling rate Fs.')
        sys.exit()
    
    #Define the range of frequency
    f_list = np.arange(f0_min, f0_max, resolution)
    
    #Generate the z-vector and harmonic summation list
    z = np.arange(0, N)
    Harmonic_sum = []
    Z_list = []
    
    #Repeat for f one by one
    for i in range(len(f_list)):
        
        #Construct the zc and zs matrix
        zc = []
        zs = []
        for L in range(1, L_max + 1):
            zc.append(np.cos(2*np.pi*f_list[i]*L*z/Fs))
            zs.append(np.sin(2*np.pi*f_list[i]*L*z/Fs))
        
        #Stack the zc and zs to construct Z-matrix
        zc = np.array(zc).T
        zs = np.array(zs).T
        Z = np.concatenate([zc, zs], axis=1) #row=N,column=2L
        
        #Compute the harmonic summation
        Zx = Z.T @ wavdata
        Harmonic_sum.append(Zx.T @ Zx)
        Z_list.append(Z)
    
    #Search for argument max (maximize the harmonic summation value)
    Harmonic_sum = np.array(Harmonic_sum)
    arg_i = np.argmax(Harmonic_sum)
    
    #Fundamental frequency
    arg_f = f_list[arg_i]
    #Maximum value
    Harmonic_max = Harmonic_sum[arg_i]
    #Z(0) N x 2L matrix
    Z0 = Z_list[arg_i]
    
    return arg_f, Harmonic_max, Z0

returnで返されるarg_fというのが、引数wavdataという音声データから推定されたピッチ（固有周波数）です。当然ながら、例えば歌を唄う場合などにはメロディに合わせてピッチが変動するので、音声を時間セグメント（Time segment）ごとに区切ってこの関数に代入してやる必要があります。それを行っているメイン関数はこちらです。

if __name__ == "__main__":
    
    #Set up
    segTime = 0.5             #Segment time (seconds)
    overlap_ratio = 0.5       #Overlap ratio of segments (0 < overlap < 1)
    f0_min, f0_max = 110, 350 #Range of fundamental frequency (Hz)
    f_delta = 1               #Resolution of frequency (Hz)
    L_max = 3                 #Degree of Harmonic summation (integer)
    
    #Read the sound data
    [wavdata, Fs] = read_wav_file('./data/record1.wav')
    #[wavdata, Fs] = generate_sinwave(f0=440) #for test use
    
    #Crop the audio
    wavdata = wavdata[round(0.5*Fs) : len(wavdata)-round(0.5*Fs)]
    wavLen = len(wavdata)
    ipd.display(ipd.Audio(data=wavdata, rate=Fs))
    
    #Calculate the parameter
    segLen = round(segTime * Fs)
    N_seg = math.floor((wavLen - overlap_ratio * segLen) / ((1-overlap_ratio) * segLen))
    stride = math.floor((1-overlap_ratio) * segLen)
    time_vector = (segLen/2 + stride * np.arange(0, N_seg)) / Fs
    
    #Do the analysis
    f0_vector = np.zeros((N_seg)) #F0 is defined as a vector
    for i in range(N_seg):
        #Measure time
        start = time.time()
        
        #Trim the wave-data into each segment
        crop_wav = wavdata[round(i*stride) : round(i*stride + segLen)] #Crop the wave-data for each segment
        [arg_f, Harmonic_max, Z0] = get_fundamental_freq(crop_wav, Fs, f0_min, f0_max, f_delta, L_max) #Call my function
        f0_vector[i] = arg_f #Divide the sampling rate (Fs) by the tau to get fundamental frequency
        
        #Report time
        print('Time={:.2f}sec: F0 estimates={:.1f}Hz, Process_time={:.1f}sec'.format((i*stride+segLen/2)/Fs, f0_vector[i], time.time() - start))
        
        #Preserve the second segment for plotting
        if i == round(N_seg)/2:
            X = crop_wav
            Z = Z0
    
    #Construct the Harmonic model S from Z
    S1 = Z.T @ X
    S2 = np.linalg.inv(Z.T @ Z)
    S = Z @ S2 @ S1
    
    #Plot the results
    plot_wavs(X, S, Fs, xmin=0, xmax=segTime/10)
    plot_STFT(wavdata, Fs, time_vector, f0_vector, f0_min, f0_max, win_size=segLen, win_overlap=segLen*0.75)
    print('Averaged F0 estimates: {:.1f}Hz'.format(np.average(f0_vector)))

plot_wavs等は単なるグラフ表示用の自己定義関数なので、この記事では割愛しています。メイン関数では、SetupのところでSegmentの時間長さや周波数の値域・分解能等をパラメータ指定します。動作が始まると各セグメントごとに行列計算を行い、推定されたピッチの推移がスペクトログラムに表示されます。10秒くらいの音声データを0.5秒のセグメントで分割しても数分で終了します。

アカペラの音痴指数を計算

まずは手始めに私のアカペラで「ラ」の音を収録してみました。スマホのマイクで録ったので音質*3はあまり良くないです。（音量注意！お聞き苦しい声が5秒間続きます。）

この音声データを上述したプログラムで各セグメントに分割してピッチ推定した結果がこちらです。

アカペラで出したラの音のSpectrogramとピッチ推移

スペクトログラム（Spectrogram）に重ねてプロットされている青い点が、私の出した「ラ」の推定ピッチです。出力結果によると私のピッチは219Hzだそうです。音楽に詳しい方ならご存知のとおり、一般に女声ソプラノでラの音は440Hzと定義されています。私は男なので1オクターブ下がって220Hzが「ラ」の正解ピッチです。パーフェクトピッチ（絶対音感）からは程遠いものの、救いようがないほど音痴ではないようです。

では続いて、今度はアカペラで歌を唄ってみました。曲は軽快なJ-POPでもQUEEN*4でも何でも良かったのですが、やはり正解の音程が分かる曲でないと音痴度が測れないので、日本人なら誰でも知っている国歌「君が代」を30秒くらい独唱しました。

うわあ改めて聞くと恥ずかしい。あまり歌は聞かないでください。さて、これをプログラムに通した結果はこちらです。

アカペラで歌った君が代のSpectrogramとピッチ推移

0.5秒という短いセグメントで分割しているので、おおむね私の声の音程変化を追うことができています。（ときどき倍音を捉えてしまっていたり、息継ぎのところで変な周波数を拾ったりしていますが、明らかな間違いは一見して分かります。）この結果を、ネットから拾った君が代の楽譜から求めた正しいピッチと一緒にExcelに入力したところ、以下のようになりました。

私の君が代 vs 正しい君が代

うーん、最初の方はちゃんと正しい音程で歌えていたものの、後半へ行くにつれて若干ピッチが上振れしてますね。ちなみに、すべての推定ピッチと正しいピッチとの差を二乗して平均した分散*5は、平方根を取ると5.2Hzでした。すなわち、私のアカペラは平均的に5.2Hzくらい音を外しているということです。

これが良いのか悪いのか分かりませんが、ドの音とレの音はだいたい16Hzくらい離れているので、5.2Hzはピアノの白い鍵同士の間隔の1/3くらいです。まあ「聞くに堪えないほどではない」といったところでしょうか。カラオケ採点機なら80点くらいを付けそうです。

以上、数式が多かったのでちょっと長くなりました。今回はピッチ推定のアルゴリズムの一例として、Harmonic model（高調波モデル）とその近似Harmonic summation modelをご紹介し、具体的にプログラムを実装して自分の声を分析してみました。このように、音声は気軽に扱えるデータなので、理論を勉強したら簡単に実データで遊べるのがおもしろいです。

稚文をお読みいただきありがとうございました。

*1:どのモデルも一長一短というところです。

*2:例えば、Pythonには「numpy」という行列演算に特化した強力なライブラリがあります。

*3:そもそも私の声質が良くないという説もあります。

*4:QUEENならWe are the championsとDon't stop me nowが好きです。フレディ・マーキュリーの声は疲れているときに聞くと何だか元気が出ます。

*5:まさに音痴指数とも呼ぶべき数値です。これが大きいほど音痴で、小さいほど音程がしっかり合っている歌い手ということになります。