２１．メル尺度とMFCCを知ろう - keep learning blog（キープラーニングブログ）

――考えるな。感じるんだ。――

メル尺度とMFCCを知ろう

メル尺度とMFCCを知ろう

一ヶ月近くブログを書かずにいたら、はてなブログから私の登録アドレスにメールが来ました。「あなたが最後に更新してから一ヶ月以上が経ちました。そろそろブログを書いてみませんか？」だそうです。~~大きなお世話です~~確かにもうそんなに更新していなかったんですね。

いちおうこのブログの立ち位置は自分の備忘録用であって、書きたい内容がなければ放置しておいてもいいのですが、本業（研究）の方で論文の草稿を書き終わったので、ようやく時間に余裕ができました。ようするに今やっとブログを書く心の余裕が生まれました。

さて、私の研究内容は機械学習を用いた音声処理です。それなのに、これまで両者を完全に切り離されたものとして記事を書いてきました。特に前回の記事は「犬と猫の画像をCNNで分類する」という内容で、いったい画像が音声と何の関係があるのだと疑問に思われる方もおられると思います。

今回はその間のミッシングリンクを埋める存在、メル尺度（Mel-scale）という考え方をご紹介します。

f:id:yuki0718:20200128023752j:plain

人間の聴覚とメル尺度

音声の話をするとき、人間の生物学的な知覚について考察することは避けられません。人間の可聴範囲（音が鳴っていると感知できる範囲）は20Hzから20kHzと言われていますが、我々の日常生活に溢れている典型的な音、例えば人間の話し声などはだいたい300Hzから700Hzの周波数であり、類稀なる才能を持ったソプラノ歌手であってもせいぜい2.5kHz（2500Hz）が最高音だそうです。

他方、蝙蝠（コウモリ）やイルカなどの生物は、人間には聞き取ることができない超音波（100kHzを超える高周波数の音）を飛ばして、障害物を認識したり仲間と交信したりすることができます。しかし私たちは蝙蝠でもイルカでもありません。人間社会の音声を利用する以上、その分析もまた人間の聴覚に沿って処理されるべきでしょう。

1937年、コンピュータ学者であり心理学者でもあったStanley Smith Stevensらは、臨床実験により人間が二つの周波数の差を等間隔だと感じる心理的なピッチ（基本周波数）の尺度に関する論文「A Scale for the Measurement of the Psychological Magnitude Pitch」を発表しました。これをメル尺度（Mel-scale）と呼びます。

論文によれば、心理音響学（psychoacoustics）的見地から、人間は低周波領域の周波数差には敏感（350Hzと450Hzの違いはすぐ分かる）なのに対し、高周波領域の周波数差には鈍感（10kHzと10.1kHzの違いはほとんど分からない）なのだそうです*1。

つまり、物理量である周波数Hzは、音波を数式で表すのには便利でも、人間の（心理的作用も含む）聴覚に沿った尺度ではないということです。

この周波数差を「心理的に等間隔」な尺度に変換したものがメル尺度（Mel-scale）というわけです。メル尺度は、近似的に以下のような対数関数でスケーリングすることができます。

$\displaystyle Mel \stackrel{\mathrm{def}}{\equiv} \frac{1000}{\log_{e} (1+1000/700)} \log_{e} \left(1+ \frac{Hz}{700} \right)$

この関数は、標準的な音声処理の教科書には必ずといっていいほど登場する式で、ちょうど人間の話し声の最高音である700Hzくらいまでは直線的ですが、そこを過ぎたところから急激に増加が遅くなります。すなわち、MelはHzが大きくなるほどHzの変化に鈍感な尺度といえます。

メル尺度で等量だけ異なる音は、人間の聴覚では同じくらい異なる音として知覚されます。周波数Hzとは違って、Mel=100と110の違いはMel=1000と1010の違いと等価（Equivalent）だといえるわけです。

メル尺度スペクトログラムとフィルタバンク

かなり前に以下の記事で紹介したとおり、現実の連続的な音の波動をマイクなどで電気信号に変えたあと、コンピュータに入力したものが音声データです。音声は時間変化するデータなので、離散フーリエ変換（Discrete Fourier transformation）と呼ばれる数学的な操作を行えば、周波数領域（Frequency domain）に「移す」ことができます。

keep-learning.hatenablog.jp

また、同記事では、音声を数10msec～100msec程度の幅の時間窓に切り取ってそれぞれの窓を離散フーリエ変換する手法、すなわち短時間フーリエ変換（STFT）が音声処理では一般的に用いられているとご説明しました。

STFTの出力結果は、周波数Fの関数であると同時に、窓が切り取った時刻（時間）の関数でもあり、サーモグラフィのような二次元マッピングが可能です。一般的に横軸を時間、縦軸を周波数とし、STFTの振幅Aを $20\log_{10}(A)$ （デシベル単位）に変換して色で表現したグラフを、スペクトログラム（Spectrogram）と呼びます。

f:id:yuki0718:20190728035947p:plain — スペクトログラムの例

実は、この縦軸（周波数軸）をメル尺度で書き直したスペクトログラム（Mel-scale spectrogram）は、機械学習の分野ではポピュラーな入力になっています。スペクトログラムは画像なので、前回の記事でご説明したような畳み込みニューラルネットワーク（CNN）に入力すれば、音声の分類器を作ることができます*2。

Mel-scale spectrogramの導出は、以下のような手順で実行されます。

１．高周波を強調するようなFIRフィルタ（ $b_{0}=1.0$ , $b_{1}=-0.97$ が一般的）をかける。

２．適当な窓長と窓間隔でSTFTをかける。（width=30msec, interval=15msecくらいが一般的）

３．得られたSTFTで各時間窓の絶対振幅（Absolute amplitude）を二乗してエネルギースペクトルにする。

４．各時間窓のエネルギーにメル・フィルタバンクと呼ばれるフィルタをかけてメル尺度にする。

５．メル尺度になったエネルギースペクトルの対数をとってdB単位に変換する。

ポイントになるのは上記４．のところで出てくるメル・フィルタバンク（Mel-fileterbank）です。これは以下のようなメル尺度で等間隔な三角形状のフィルタになります。

f:id:yuki0718:20200128003629p:plain — メル尺度で等間隔な三角フィルタとそのHzスケール

メル尺度で等間隔な繰り返しフィルタをパワースペクトルに掛けるということは、つまりはメル尺度で等間隔な点（Channelと呼ばれます）それぞれの近傍でエネルギーを積分することと等価です。積分されたエネルギーbinは、私たちが全体の音波を知覚したときに「あ、～番目のChannelの音が強く鳴ってるぞ」と感じる心理的インパクトに対応します。

メル・フィルタバンクをPythonで実装

メル・フィルタバンクは、以下に示すような100行に満たないPythonコードで簡単に実装することができます。

import math
import glob
import soundfile as sf
import numpy as np
from scipy import signal as sg
from scipy import fftpack as fp

#■Function for generating Mel-scale filters■
def melFilterBank(Fs, fftsize, Mel_scale, Mel_cutf, Mel_channel, Mel_norm):
    
    #Define Mel-scale parameter m0 based on "1000Mel = 1000Hz"
    m0 = 1000.0 / np.log(1000.0 / Mel_scale + 1.0)
    
    #Resolution of frequency
    df = Fs / fftsize
    
    #Mel-scale filters are periodic triangle-shaped structures
    #Define the lower and higher frequency limit of Mel-scale filers
    Nyq = Fs / 2
    f_low, f_high = Mel_cutf
    if f_low is None:
        f_low = 0
    elif f_low < 0:
        f_low = 0
    if f_high is None:
        f_high = Nyq
    elif f_high > Nyq or f_high <= f_low: 
        f_high = Nyq
    #Convert into Mel-scale
    mel_Nyq = m0 * np.log(Nyq / Mel_scale + 1.0)
    mel_low = m0 * np.log(f_low / Mel_scale + 1.0)
    mel_high = m0 * np.log(f_high / Mel_scale + 1.0)
    #Convert into index-scale
    n_Nyq = round(fftsize / 2)
    n_low = round(f_low / df)
    n_high = round(f_high / df)
    
    #Calculate the Mel-scale interval between triangle-shaped structures
    #Divided by channel+1 because the termination is not the center of triangle but its right edge
    dmel = (mel_high - mel_low) / (Mel_channel + 1)
    
    #List up the center position of each triangle
    mel_center = mel_low + np.arange(1, Mel_channel + 1) * dmel
    
    #Convert the center position into Hz-scale
    f_center = Mel_scale * (np.exp(mel_center / m0) - 1.0)
    
    #Define the center, start, and end position of triangle as index-scale
    n_center = np.round(f_center / df)
    n_start = np.hstack(([n_low], n_center[0 : Mel_channel - 1]))
    n_stop = np.hstack((n_center[1 : Mel_channel], [n_high]))
    
    #Initial condition is a 0 padding matrix
    output = np.zeros((n_Nyq, Mel_channel))
    
    #Repeat every channel
    for c in np.arange(0, Mel_channel):
        
        #Slope of a triangle(growing slope)
        upslope = 1.0 / (n_center[c] - n_start[c])
        
        #Add a linear function passing through (nstart, 0) to output matrix 
        for x in np.arange(n_start[c], n_center[c]):
            #Add to output matrix
            x = int(x)
            output[x, c] = (x - n_start[c]) * upslope
        
        #Slope of a triangle(declining slope)
        dwslope = 1.0 / (n_stop[c] - n_center[c])
        
        #Add a linear function passing through (ncenter, 1) to output matrix 
        for x in np.arange(n_center[c], n_stop[c]):
            #Add to output matrix
            x = int(x)
            output[x, c] = 1.0 - ((x - n_center[c]) * dwslope)
        
        #Normalize area underneath each Mel-filter into 1
        #[ref] T.Ganchev, N.Fakotakis, and G.Kokkinakis, Proc. of SPECOM 1, 191-194 (2005)
        #https://pdfs.semanticscholar.org/f4b9/8dbd75c87a86a8bf0d7e09e3ebbb63d14954.pdf
        if Mel_norm == True:
            output[:, c] = output[:, c] * 2 / (n_stop[c] - n_start[c])
    
    #Return Mel-scale filters as list (row: frequency, column: Mel-channel)
    return output

メル・フィルタバンク（Mel-filterbank）の定義にはいくつか流派というかアレンジがあるので、このコードは一般的な実装と少し異なるかもしれません*3。

具体的には、Configurationパラメータとして、上述したメル尺度の定義式に存在するしきい値Mel_scale（通常700）や、フィルタの端を切るカットオフ周波数Mel_cutf=[fmin, fmax]（通常[0, Nyquist周波数]）を外から指定可能にするためのアレンジを加えています。Mel channelはチャンネル数（通常30～120程度）を表し、fmin～fmaxの領域を等間隔に分割します。

また、音声処理ライブラリ「Slaney’s Auditory Toolbox」に倣って、Mel_norm=Trueに指定すると、各チャンネルのフィルタに囲まれる面積を1に正規化（Normalization）して、エネルギー積分が各チャンネルで等価になるようにしています*4。Falseにしておくと面積ではなく三角フィルタの頂点の高さが1に正規化されます（前出のグラフはこちらで出力）。

なお、この関数からは、行が周波数、列がメル尺度の行列という形でフィルタバンクが出力されます。STFTで各時間窓の絶対振幅を二乗したパワースペクトルを計算し、フィルタバンクとの行列積を計算してから対数をとってdB単位に変換すれば、Mel spectrogramが得られます。

MFCC（Mel-frequency cepstrum coefficient）

Mel-scale spectrogramは人間の心理的な尺度に沿った音声の特徴量表現であり、画像としてCNNに入力すれば音声分類で絶大なパフォーマンスを発揮します。しかし、スペクトログラムのような画像はCNNへの入力には適している反面、情報量が多すぎて冗長（Redundant）であるため、GMMのような古典的なモデルへの入力には向きません。

そこで、ついでと言ってはなんですが、音声の機械学習モデルで古くから入力データとして頻繁に用いられてきた、MFCC（Mel-frequency cepstrum coefficient）と呼ばれる情報の圧縮化手法を紹介しておきます。

MFCCとは、メル・フィルタバンクのチャンネル数（典型的には30～120程度）だけ存在するdB単位のスペクトログラムデータを、DCT（Discrete cosine transformation）でフーリエ逆変換処理し、得られた係数を次数が低いものから有限数個だけ抜粋してしまうという手法です。数学的にはメル尺度スペクトルの包絡関数（Envelope of spectrum）を取得していることになります*5。

MFCCの計算手順は、Mel-scale spectrogramのものと途中まで完全に同じで、最後の「５．対数をとってdB単位に変換する。」の後に「６．DCTを行って次数が小さい方からいくつかの係数を抜粋する。」というプロセスが加わるだけです。

このとき、次数に0次の係数を含むかどうか、何次までとるのか、といった自由度はあるものの、典型的には0次を含む13個程度のDCT係数をMFCCとして抽出します。

MFCCをPythonで実装

MFCCについても、先ほどのメル・フィルタバンク生成関数を利用して、Pythonで簡単に実装することができます。なお、以下のコードからDCT以降の部分をカットすれば、Mel-scale spectrogramが出力されます（MFCCの方が一手間余計にかかっている）。

#■Function for getting MFCC■
def get_MFCC(folder_path, frame_length, frame_shift, Mel_scale, Mel_cutf, Mel_channel, Mel_norm, MFCC_num, MFCC_lifter, Add_deriv):
    
    #Inicialize list
    mfcc = []
    
    #Get .wav files as an object
    files = glob.glob(folder_path + "/*.wav")
    print("Folder:" + folder_path)
    
    #For a progress bar
    nfiles = len(files)
    unit = math.floor(nfiles/20)
    bar = "#" + " " * math.floor(nfiles/unit)
    
    #Repeat every file-name
    for i, file in enumerate(files):
        
        #Display a progress bar
        print("\rProgress:[{0}] {1}/{2} Processing...".format(bar, i+1, nfiles), end="")
        if i % unit == 0:
            bar = "#" * math.ceil(i/unit) + " " * math.floor((nfiles-i)/unit)
            print("\rProgress:[{0}] {1}/{2} Processing...".format(bar, i+1, nfiles), end="")
        
        #Read the .wav file
        data, Fs = sf.read(file)
        #Transform stereo into monoral
        if(isinstance(data[0], list) == True):
            wavdata = 0.5*data[:, 0] + 0.5*data[:, 1]
        else:
            wavdata = data
        
        #Normalize the wave ranging between -1 and 1
        wavdata = (wavdata - np.mean(wavdata))
        wavdata = wavdata / np.amax(np.abs(wavdata))
        
        #Calculate the index of window size and overlap
        FL = round(frame_length * Fs)
        FS = round(frame_shift * Fs)
        OL = FL - FS
        
        #Pass through a pre-emphasis fileter to emphasize the high frequency
        wavdata = sg.lfilter([1.0, -0.97], 1, wavdata)
        
        #Execute STFT
        F, T, dft = sg.stft(wavdata, fs=Fs, window='hamm', nperseg=FL, noverlap=OL)
        Adft = np.abs(dft)[0 : round(FL/2)]**2
        
        #Call my function for generating Mel-scale filters(row: fftsize/2, column: Channel)
        filterbank = melFilterBank(Fs, FL, Mel_scale, Mel_cutf, Mel_channel, Mel_norm)
        
        #Multiply the filters into the STFT amplitude, and get logarithm of it
        meldft = Adft.T @ filterbank
        if np.any(meldft == 0):
            meldft = np.where(meldft == 0, 1e-9, meldft)
        meldft = np.log10(meldft)
        
        #Get the DCT coefficients (DCT: Discrete Cosine Transformation)
        dct = fp.realtransforms.dct(meldft, type=2, norm="ortho", axis=-1)
        
        #Trim the MFCC features from C(0) to C(MFCC_num)
        dct = np.array(dct[:, 0 : MFCC_num])
        
        #Lift up the cepstrum by sine-function depending on channel
        if MFCC_lifter is not None:
            #K.Paliwal, "Decorrelated and liftered filter-bank energies for robust speech recognition." Eurospeech(1999)
            #[ref] https://maxwell.ict.griffith.edu.au/spl/publications/papers/euro99_kkp_fbe.pdf
            lifter = 1 + 0.5 * MFCC_lifter * np.sin(np.pi * np.arange(0, MFCC_num) / MFCC_lifter)
            nframes = dct.shape[0]
            lifter = lifter[np.newaxis, :]
            lifter = np.tile(lifter, (nframes, 1)) #Duplicate in order to Hadamard operation
            dct = dct * lifter
        
        #Compute the 1st and 2nd derivatives
        if Add_deriv == True:
            deriv1 = np.diff(dct, n=1, axis=0)[1:] #Adjust the length into 2nd derivative (1:)
            deriv2 = np.diff(dct, n=2, axis=0)
            dct = dct[2:] #Adjust the length into 2nd derivative (2:)
            dct = np.concatenate([dct, deriv1], axis=1)
            dct = np.concatenate([dct, deriv2], axis=1)
        
        #Add results to lists sequentially
        mfcc.append(dct)
    
    #Finish the progress bar
    bar = "#" * math.ceil(nfiles/unit)
    print("\rProgress:[{0}] {1}/{2} Completed!   ".format(bar, i+1, nfiles), end="")
    print()
    
    #Return the result (as list format not numpy tensor)
    return mfcc

MFCCの計算にもいくつかアレンジがあるので、一般的なMFCCのコードとは少し異なるかもしれません。

具体的には、MFCC_num（典型的には13）個のMFCCを計算するだけでなく、Add_deriv=Trueにしておくと、その時間フレーム一階微分と二階微分も合わせて算出するようにしています。つまり、各時間窓ごとにMFCC_num $\times$ 3個の特徴量が抽出されます。

情報量を圧縮しておいて再び増やすのか、と思われるかもしれませんが、経験上これはパフォーマンス向上にかなり寄与します。これはMFCCのみならずMel-scale spectrogramでも同様です*6。

また、1999年にK.Paliwalらが発表した論文「Decorrelated and liftered filter-bank energies for robust speech recognition」によると、MFCCは高次の係数になるほど減衰する（Attenuate）傾向があるので、高次の係数をリフトアップすると雑音に対してロバストになって音声認識の精度が向上するそうです。

論文中では様々なリフトアップ関数を比較していましたが、私のコードではリフターパラメータのMFCC_lifter（典型的には22）をNone以外に指定すると、sin関数でMFCCをリフトアップするように実装しています。

以上のように、MFCCは古典的な音声処理においてデファクトスタンダード（de facto standard）になっていた（いる）歴史があるため、長い歴史の中で、こういった細かい工夫・アレンジが多数存在しています。

一言に「MFCCを使って機械学習したらこういう結果になりました」と誰かが論文で発表したりネットに書いていたりしても、MFCCを導出する過程の細かい処理が異なると再現性が得られないので、ライブラリを使ってMFCCを計算する際には注意しなければなりません*7。

今日のところは以上です。Pythonでちょっとしたコードを書けば、音声をMel-scale spectrogramで二次元画像に変換したり、MFCCを取得したり（冗長なら全時間フレームで時間平均を取ったり）することで、音声の特徴量を数値化できます。

メル尺度の特徴量は、音声分類のみならず、話者認識、声質変換などで広く用いられている一般的なものです。古典モデルかDNNかを問わず、機械学習には優れた入力データの存在が必須なので、機械学習に興味がある方は（音声自体には興味がなくても）覚えておいて損はないと思います。

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

*1:これは、日常生活に溢れている人間の話し声が300Hz～700Hzという低周波領域であることに起因します。身近にない音を区別する能力が発達した生物など、真っ先に淘汰されるので進化論的にありえないということです。

*2:Acoustic Scene Classificationと呼ばれる音声シーン分類の学術領域では、様々なDNNモデルの活用が提案されていますが、半分以上の論文がMel-scale spectrogramを入力データとしてCNNに入力している印象でした。

*3:各アレンジの詳細ついては2005年にT. Ganchevらのまとめた「Comparative Evaluation of Various MFCC Implementations on the Speaker Verification Task」という論文が大変分かりやすいです。

*4:Pythonの音声処理ライブラリlibrosaでも、メル・フィルタバンクは面積によって正規化するように実装していました。

*5:参考URL: http://wantee.github.io/2015/03/14/feature-extraction-for-asr-mfcc/

*6:もちろん、入力の次元が3倍になるのですから、機械学習モデルに通したときの計算量も3倍、あるいはその累乗で負荷がかかります。お試しで学習させるときはAdd_deriv=Falseにして微分を含めない方が身のためです。

*7:これはMFCCに限らず、あらゆる物理量の計算ライブラリに言えることです。中身を知らないままライブラリを使って結果だけ得るのは、サイエンティストとしてもエンジニアとしても危険な行為です――と殊勝なことを言いつつ、内部実装を完全に理解しないままライブラリを使うことなんて私もザラにあるんですけれども。