Notebook

Day01::Pythonではじめたい音声認識

音声認識

ことのいきさつ

Twitterで面白そうなタグ #100DaysOfCode というものを見つけました.三日坊主克服のために分野を問わず,進捗をブログにまとめます.

音声処理

音声認識 (機械学習プロフェッショナルシリーズ)

音声認識 (機械学習プロフェッショナルシリーズ)

時系列分野に対する知識が疎いと感じたため読みます.


人体と音

音の高さ(音程)

人間が認識できる音の高さは一般に20Hzから20000Hzと言われており,人間の耳が知覚する音の高さと,実際の周波数との関係は線形ではなく,高音域になるほど音程の変化に対する知覚は小さくなると言われています.

人間の知覚する音の高さを測る尺度としては,メル尺度(mel scale)がよく使われています.これは,1000Hzの音を基準とし,その音に対し$n$倍の高さに知覚される周波数を$n\times 1000$melと表記することで表現されます.例えば,以下のファントの式が周波数とmelの変換によく用いられます. $$ m=\frac{1000}{\log_{10}2}\log_{10}\left(\frac{f}{1000}+1\right) $$

音の大きさ(音圧)

音の物理的強度は音圧と呼ばれ,単位はPa.melとおなじように図るための尺度が存在して,音圧が$p$であるときの音圧レベル(Sopund Pressure Level:SPL)$L_p$は $$ L_p=10\times\log_{10}\frac{p^{2}}{p_0^{2}}=20\times\log_{10}\frac{p}{p_0} $$ と表され,単位はdB SPLです.ただし$p_0=20\times10^{-6}$(Pa)で,人間の耳が聞こえる最小の音圧に対応しています.

音の大きさの感じ取り方も周波数に依存し,人間が同じ音の大きさであると感じる音の音圧レベルと周波数間の関係をプロットした曲線は等ラウドネス曲線と呼ばれます.

Wikipediaより

人間の感じる音の大きさの単位はホン(phon)であり,これは1000Hzの音の音圧レベルをホンと定めています.ホンが等しい音どうしは人間には同じ音の大きさとして聞こえます.例えば1000Hzの20dBの音は80Hzの50Hzの音の大きさと等しく,どちらも20phonの音といえます.

他にも人間の耳には興味深い効果があり,その例としては,

  • 大きな音の鳴った前後の音が聞こなくなる経時マスキング効果
  • 大きな音の鳴っている時,その周波数前後の音が聞こえなくなる同時マスキング効果
  • 対面の相手の音声認識において無意識に口唇の動きで発音を読み取るマガーグ効果
  • 騒音の中でも遠くの音を集中して聞こうとすることで聞こえるようになるカクテルパーティ効果

などがあります.


人間の発声の過程は,肺からの息で押し出された声帯の振動から声門波を発生させ,それが声道を通り口唇から放射されるという過程で,声帯の振動周波数は基本周波数と呼ばれます.有声音の場合は声帯は振動するが,無声音の場合には振動しないことが特徴.

声道は,複数の太さの音響管をつないだ構造と近似することができて,その構造故に複数の共鳴周波数があり,その周波数の音が特に強く発声されます.この周波数はフォルマント周波数と呼ばれている.また,フォルマント周波数が複数存在する場合は周波数が低い方から順に,第一フォルマント(F1),第二フォルマント(F2)...などと呼ばれています.また,舌と歯をつけて発生するような子音などは摩擦音として高周波音域の音として観察されます.

音声の持つ特徴のうち特に重要なものは音韻音素と呼ばれます.音韻はある言語における音声認識の最小単位の音の集合であり,音素は音韻同様音声を構成する集合であるが,こちらは音の物理的に特徴により分類したものである.

世界中の言語の音素を表現するために国際音声字母(IPA:International Phonetic Alphabet)が定められており,下表は日本語の子音のIPAにおける音素分類です.

調声器官 口唇有声 口唇無声 歯有声 歯無声 口蓋有声 口蓋無声 声帯無声
摩擦音 f z s sh h
破擦音 dz tz ch
破裂音 b p d t g k
半母音 w r y
鼻音 m n

対して日本語の母音(aiueo)は声道や舌の形を変えることでフォルマント周波数を変えることで生成されます.

また日本語の場合には上表の他に撥音(ん/ɴなど),拗音(ゃゅょ/fya,tjaなど),促音(っ/ːなど)の音素が存在します.

音声の分析

AD変換

音声自体は空気の振動であり,マイクロホンにより電圧の連続的な変化として扱われるため,計算機上で処理を行うためには連続値を離散値へと変換させる必要があります(AD変換).

AD変換の際はまず,アナログ信号を一定間隔で取り出し離散的な時間間隔のデータとします.この操作は標本化(Sampling)と呼ばれています.

さて,一般的なCD等の音源を見ると44.1kHzで標本化されていることが多いが,これは人間の可聴域の上界が20kHz程度であることとサンプリング定理が成り立つことからそう定められています.


サンプリング定理

あるアナログ信号が0HzからW[Hz]までの音で構成されているとしたとき,2W[Hz]以上の周波数で標本化を行えば元の信号を完全に復元できる


標本化したデータに対して,次に音の大きさを離散化させます.この処理は量子化(Quantization)と呼ばれており,多くは後に2進データとして符号化するために2の累乗個に分割されます.こう言った場合には,分割数は量子化ビット数などと呼ばれます.すなわち,8bitの量子化ビット数は標本一つあたり8bitの情報が与えられ,音を256の細かさに分類して表現することができるということを意味します.一般には16bit(6.5万分割),32bit(42億分割)の量子化ビット数が用いられます.

Pythonでの実装例

音声は効果音ラボ落ち着いた女性の「ご静聴ありがとうございました」を利用します.

以下作業ディレクトリ直下に音声を置いたと仮定してまずは依存パッケージとかを入れます.

brew install ffmpeg --with-libvorbis --with-ffplay --with-theora
pip install pydub numpy matplotlib

以下Jupyter Notebook上で行います.

wavデータで読み込みたいので変換をします.

import pydub
sound = pydub.AudioSegment.from_mp3("info-lady1-goseichouarigatou1.mp3")
sound.export("voice.wav", format="wav")

wavデータの読み込みには標準パッケージのwaveを利用します.そして,得られたバイナリデータをNumPyのfrom_buffer関数でNumPy配列にします.

import wave
import numpy as np

def get_wave_profile(file):
    wf = wave.open(file, "rb")
    # モノラルかステレオか
    nc = wf.getnchannels() 
    # サンプルの大きさ(量子化ビット数)
    ss = wf.getsampwidth()
    # サンプリングレート
    sr = wf.getframerate()
    # フレーム数
    fl = wf.getnframes()
    # 生データ(バイナリ)
    data = wf.readframes(fl)
    # NumPy配列に変換
    data = np.frombuffer(data, dtype="int%d" % (8 * ss))
    # 正規化
    data = data / data.max()
    print(file)
    print("\tチャネル数\t\t%s" % nc)
    print("\t量子化バイト数\t\t%s Byte" % ss)
    print("\tサンプリングレート\t%s Hz" % sr)
    print("\tフレーム数\t\t%s" % fl)
    print("\t秒数\t\t\t%.3f 秒" % (fl / sr))
    wf.close()
    params = {
        "channels": nc,
        "width": ss,
        "rate": sr,
        "length": fl
    }
    return data, params
  
data, params = get_wave_profile("voice.wav")
#voice.wav
#  チャネル数             1
#  量子化バイト数           2 Byte
#  サンプリングレート     48000 Hz
#  フレーム数             98816
#  秒数                  2.059 秒

波形を表示してみましょう

import matplotlib.pyplot as plt
plt.figure(figsize=(11, 5))
plt.plot(range(len(data)), data / data.max())
plt.show()

音声編集ソフトとかでよく見る形ですね.

ついでに,関数から音を生成する関数を作ります.

import struct
'''
func: 関数
fs:   サンプリング周波数
t:    秒数
name: 出力ファイル名
'''
def create_wave(func, fs, t, name):
    point = np.arange(0,fs*t)
    wave = func(point/fs)

    wave = [int(x * 32767.0) for x in wave]
    binwave = struct.pack("h" * len(wave), *wave)

    w = wave.Wave_write(name)
    # 1:モノラル, 2:16bitの量子化, fs:サンプリング周波数
    # バイトでのサンプルサイズ, 圧縮形式(現状PythonではNoneだけがサポートされている)
    p = (1, 2, fs, len(binwave), 'NONE', 'not compressed')
    w.setparams(p)
    w.writeframes(binwave)
    w.close()

検証のために,単純なSin波と少し複雑な合成Sin波を生成させます.

create_wave(lambda x: np.sin(2 * np.pi * 440 * x), 44100, 0.3, "sin.wav")
create_wave(
    lambda x: 0.2 * np.sin(2 * np.pi * 680 * x) + \
              0.3 * np.sin(2 * np.pi * 440 * x) + \
              0.5 * np.sin(2 * np.pi * 1880 * x),
    44100,
    0.3, 
    "sins.wav"
)

sin.wav(少し長めに出力しています)

sins.wav(少し長めに出力しています)

以下,この音声を用いて音声解析を行います.

高音域強調

高音域になるほど音圧は減衰するため,それを補うため高音域を強調します.一般には1オクターブあたり6dBの強調が行われます.ハイパスフィルタを用いて高音域強調を行う場合を考えると,一次有限インパルス応答をハイパスフィルタとして用いる場合,その式は $$ H(z)=1-\alpha z^{-1}~~(z=e^{j\omega},\omega=2\pi f/f_0) $$ となる($f_0$はサンプリング周波数)ので,その対数ゲインを考えるれば,単位をdBにするため両辺を10倍することで $$ 10\log{|H|^{2}}=10\log(1+\alpha^{2}+2\alpha\cos\omega) $$ となります.6dB/1octの高音強調を達成させるために,$\alpha=0.97$程度の値がハイパスフィルタによく用いられます.

音声フレーム

音声のような時系列信号を解析するためには,一般に一定の時間間隔で特徴量を出力させます.その際の間隔ごとに分析される対象は音声フレームと呼ばれ,間隔自体はフレーム間隔やフレーム周期などと呼ばれます.

フレーム周期は短すぎても計算量が増加し,長すぎても時系列的な特徴が失われてしまうことから,適用な長さに定める必要があります.

人間の一般的な音素の長さがおよそ30ms以上であるというデータをもとに,10ms程度のフレーム間隔として取られることが多いです.

つぎにその間隔で信号をとる場合,単純な矩形窓で取る方法が愚直に考えられますが,それだと両端の不連続な波形が高周波に対してのノイズを生む危険性があります.

img

そのため,フレームの中心から端にかけて減衰するハミング窓やハニング窓を音声データに対して適用し,高音域のノイズを抑えています.

ハミング窓

$$ W_H(n)=0.54-0.46\cos\left(\frac{2n\pi}{N-1}\right)(n=0,1,2,...,N-1) $$

ハニング窓

$$ W_N(n)=0.5-0.5\cos\left(\frac{2n\pi}{N-1}\right)(n=0,1,2,...,N-1) $$ いずれも,0未満,N以上の時には0を取ります.この二つの窓の違いによる分析性能の変化はさほどないものの一般的にはハミング窓の方が良く用いられます.

窓自体の長さはフレーム長と呼ばれます.フレーム長は一般に20ms-80ms程度の値を取ることが多いです.たとえばフレーム周期10ms,フレーム長50msの場合には0-50ms,10-60ms,20-70msのデータが特徴量として利用されます.

音声の特徴量抽出

ほとんどの音声認識アルゴリズムでは,実際の音声から,音声認識にとって重要な特徴量を抽出する処理が行われます.

短時間フーリエ解析

フーリエ解析のなかで,この音声信号処理のように窓関数をずらしながらフーリエ解析を行うことを短時間フーリエ解析といいます.

信号$x_N[n]$が周期$N$の信号であるとすると,以下の式で表される離散フーリエ変換(DFT)$X_N[k]$が存在します. $$ X_N[k]=\sum_{n=0}^{N-1}x_N[n]e^{-j2\pi nk/N} (0\leq k<N) $$ その逆変換ももちろん存在して $$ x_N[n]=\frac{1}{N}\sum_{k=0}^{N-1}X_N[k]e^{j2\pi nk/N}(0\leq n<N) $$ となる.実際の処理では,高速化させた高速フーリエ変換(FFT)を用いて信号解析を行います.

FFTではビンの数$N$は2のべき乗である必要があり,そのため窓幅とサンプリング周波数を調節して$N$が条件を満たすような値にする必要があります(例えば窓幅32ms,サンプリング周波数16kHz→$N=512$など).高速フーリエ変換を用いることで,計算量を$O(N^{2})$から$O(N\log_2N)$にまで抑えることが可能になります.

$X_N[k]$は複素スペクトルと呼ばれ,そのゲインを用いて,音の大きさを表すパワースペクトルを $$ \begin{align} S_k&=\frac{1}{N}|X_N[k]|^{2}\\ &=\frac{1}{N}\left|\sum_{n=0}^{N-1}x_N[n]e^{-j2\pi nk/N}\right| \end{align} $$ と定義することができます.フォルマント解析などに用いられるスペクトログラムは軸を時間と周波数とし,色をその周波数,時間におけるパワースペクトルの大きさとして表現したものです.人間が音声を認識する際,音声のゲインの情報を専ら使い,位相の情報はそこまで使われないことが知られており,この事実からスペクトログラムに,音声認識で必要な情報が多く含まれているといえます.

Pythonでの実装

NumPyには,高速フーリエ変換を行うための関数が実装されているため,それを利用します.

まずは,単純な離散信号をフーリエ変換する関数を定義します.(この関数は短時間ではなく渡された信号全体に関してフーリエ変換を行います.)

def fft(X, params):
    # scipyでも以下のフーリエ変換の処理は可能です
    X = np.fft.fft(X)
    # 結果Xに対応する周波数を求めます
    freqs = np.fft.fftfreq(len(X), d=1 / params["rate"])
    gain = np.abs(X)
    # 位相はあまり使いませんが一応返り値として返します.
    phase = np.angle(X)
    return freqs, gain, phase

生成した2つのSin波に対して,先ずはfft関数を用いてフーリエ変換してみます.

sin_data = get_wave_profile("sin.wav")
f, g, p = fft(*sin_data)
n = len(f)
# FTは共役な複素数が帰るという特性上,正負で鏡像となる(偶対称)
# FFTの場合,前半部分が正となるのでその部分だけを利用する
f, g, p = f[:n//2], g[:n//2], p[:n//2]
fig, axs = plt.subplots(2, 1, figsize=(12, 7), sharex=True)
axs[0].plot(f, g)
# np.unwrap で位相が -pi ~ piの範囲に限定されなくなる.
axs[1].plot(f, np.unwrap(p))
axs[0].set_ylabel("Amp.")
axs[1].set_ylabel("Phase")
axs[1].set_xlabel("Frequency")
plt.show()

これを2つの波に対して適用してみましょう.

sin.wav

sins.wav

しっかりと,入力周波数部分で振幅が大きな値を示しています.

では,次に短時間フーリエ解析です.NumPyには波形を$N$分割した際に各々のハニング/ハミング関数が幾つになるかを返すnp.hanning(N),np.hamming(N)が用意されています.これを利用して,ハミング窓を利用した短時間フーリエ解析を行います.

# res: フレーム長
# intv : フレーム周期
def stft_hamming(X, params, res=512, intv=32):
    f = None
    gs = [] # (フレーム数, フレーム長)の行列
    for i in range(len(X) // intv):
        w = X[i * intv:i * intv + res]
        # 端は切り捨てる
        if len(w) != res:
            continue
        # ハミング関数.要素ごとにかければ端点の連続性を仮定できる.
        hamm = np.hamming(len(w))
        f, g, _ = fft(w * hamm, params)
        gs.append(g[:res//2])
    return f[:res//2], np.array(gs)

FFT部分は前に定義したfft関数の使い回しで,for文を使った単純な短時間フーリエ解析です.

これを用いて,sin波と音声波形に対して短時間フーリエ解析を行いましょう.

sin_data = get_wave_profile("sin.wav")
f, gs = fft_hamming(*sin_data)

# 音声に対してSTFTを行う場合はこっち
# v_data = get_wave_profile("voice.wav")
# f, gs = stft_hamming(*v_data)

plt.figure(figsize=(7, 3))
plt.plot(f, 20 * np.log10(gs[10] / gs.max()))
plt.xlabel("Frequency [Hz]")
plt.ylabel("Power [dB]")
plt.show()

こうすることで,ある一つの窓の音圧が表示されます.ただし,ここでのdBは最大音圧に対する相対値であって絶対的な圧力から定義されるdB SPLではありません.

ついでに,アニメーションとして出力します.

def stft_anim_hamming(X, params, res=2048, intv=128):
    ims = []
    fig = plt.figure(figsize=(7, 3))
    f, gs = stft_hamming(X, params, res, intv)
    for g in gs:
        im = plt.plot(f[:res//2], 20 * np.log10(g[:res//2] / gs.max()), c="b")
        plt.xlim((0, 10000))
        plt.ylim((-100, 10))
        plt.xlabel("Freq [Hz]")
        plt.ylabel("Log Power [dB]")
        ims.append(im)
    ani = anim.ArtistAnimation(fig, ims, interval=1000*intv/params["rate"])
    ani.save("output_hamm.mp4", writer="ffmpeg")
    
# Jupyter Notebook上でMP4動画を再生できます
from IPython.display import Video
v_data = get_wave_profile("voice.wav")
stft_anim_hamming(*v_data)
Video("output_hamm.mp4", width=400, embed=True)

メルフィルタバンク

音声認識の観点からは,単純なビン幅による分割ほど細かな分割は要らず,隣接したスペクトルを纏めて,バリアンスを削減しても,音声認識上は問題は少ないといえます.

分散を小さくさせたいわけですが,高周波ほど人間の耳の分解能が低くなるという特性を反映させるために,周波数をメル尺度に変換して,その上で等間隔となるようにバンドパスフィルタを適用させます.すなわち周波数尺度に直すと高周波ほど幅が広い三角形状のバンドパスフィルタを入力に対して適用させます.

各々のメルフィルタに適用した出力を,用いることで次元数を削減することが可能となります.この出力の対数を次のケプストラム特徴量抽出で用います.

Pythonでの実装

メルフィルタバンクを可視化させる

# f_min, f_max: 周波数スケールでの最小値と最大値
# n: 分割数
def mel_filter_show(f_min, f_max, n=20):
    # 周波数 -> メル変換
    mel = lambda x: 1000 * (np.log10(x / 1000 + 1))/ np.log10(2)
    # メル -> 周波数変換
    invmel = lambda y:1000 * (np.power(10, (y * np.log10(2)/ 1000)) - 1)
    # メルスケールでの最大値と最小値
    (m_min, m_max) = mel(f_min), mel(f_max)
    # メルスケールでの等分割を周波数スケールに直す
    mp = invmel(np.linspace(m_min, m_max, n + 2))
    t = np.arange(f_min, f_max)
    plt.xlabel("Freq [Hz]")
    for i in range(1, n + 1):
        mf1 = lambda x: 1 / (mp[i] - mp[i - 1]) * (x - mp[i - 1])
        mf2 = lambda x: 1 / (mp[i + 1] - mp[i]) * (mp[i + 1] - x)
        # バンドパスフィルタの定義
        mf = lambda x: np.maximum(np.minimum(mf1(x), mf2(x)), 1e-80)
        # 周波数スケールでの可視化
        plt.plot(t, mf(t, mp))
    plt.show()
    
mel_filter_show(100, 10000, 10)

これで,メルフィルタバンクがどの周波数帯の音に対して強く働くかが可視化されます.

メルフィルタバンクの実装

上記のコードを流用しメルフィルタバンクを適用させる関数を作成します.次に適用する関数のことを考慮し,対数パワースペクトルを出力させます.

# 引数
# f: 周波数の格納されたベクトル
# spec: (フレーム数 Nf, フレーム長 fl)のゲインの行列
# f_min, f_max: 周波数スケールでの最小値と最大値
# n: 分割数
# 返り値
# mp: メルフィルタバンクの周波数(shape=(n,))
# mspec: 各々のフレームのメルフィルタバンクにおける対数パワースペクトルの大きさ(shape=(Nf , n))
def mel_filter(f, spec, f_min, f_max, n=20):
    mel = lambda x: 1000 * (np.log10(x / 1000 + 1))/ np.log10(2)
    invmel = lambda y:1000 * (np.power(10, (y * np.log10(2)/ 1000)) - 1)
    (m_min, m_max) = mel(f_min), mel(f_max)
    mp = invmel(np.linspace(m_min, m_max, n + 2))
    t = np.arange(f_min, f_max)
    mspec = np.zeros((spec.shape[0], n))
    for i in range(1, n + 1):
        mf1 = lambda x: 1 / (mp[i] - mp[i - 1]) * (x - mp[i - 1])
        mf2 = lambda x: 1 / (mp[i + 1] - mp[i]) * (mp[i + 1] - x)
        mf = lambda x: np.maximum(np.minimum(mf1(x), mf2(x)), 1e-80)
        # i番目のフィルタが周波数に対して作用する度合い
        coef = mf(f)
        # @ はドット積(spec.dot(coef))と同義
        mspec[:, i - 1] = np.log10(spec @ coef)
    return mp, mspec
v_data = get_wave_profile("voice.wav")
f, gs = fft_hamming(*v_data)
mp, ms = mel_filter(f, gs, 0, 25000, 20) # 20 or 80

plt.figure(figsize=(7, 3))
plt.plot(f, np.log10(gs[180]))
plt.plot(mp[1:-1], ms[180], marker='o')
plt.show()

これで,あるフレームにおけるメルフィルタバンクに変換した値が出力されます.青線で記した単なる対数パワースペクトルに比べると,バリアンスが大幅に減少していることがわかります.

$N​$=20の場合

$N$=80の場合

ケプストラム特徴量

これまで生の波形をフーリエ変換し,そのゲインをメル尺度でメルフィルタバンクを用いて次元を減らし,その対数をとりました.音声を微細構造部分と包絡成分(ノイズを取り除いたような滑らかな部分)に分けると,微細な方は声門の波形を,包絡成分は声道のインパルス応答を表していることが知られています.音声認識では主に後者の方が必要となることが多いので,これをメルフィルタバンクの出力から取り出します.

声門波に対するインパルス応答$h$を考える際,声門波によるフィルタリング$v$との畳み込みにより,観測信号$y$が記述できます.すなわち,畳み込み和 $$ y[n]=\sum_{m=-\infty}^{\infty}v[m]h[n-m]=v[n]*h[n] $$ で観測信号が記述されます.畳み込み和のフーリエ変換は積に直すことができて $$ Y[k]=V[k]H[k] $$ 今回のメルフィルタバンクの出力であるパワースペクトルは対数を取っているので, $$ \log Y[k]=\log V[k]+\log H[k] $$ となり,観測信号が包絡部分$H$と声門のフィルタ部分$V$に分けることができます.

人間の聴覚特性は,パワースペクトルの対数と対応していると言われています.したがってここで,パワースペクトル$S[k]$の対数を考えると, $$ \log S[k]=\log|Y[k]|^{2}=2\log|V[k]|+2\log|H[k]| $$ となります.では,この対数パワースペクトルに逆離散コサイン変換を適用します.パワースペクトルの対称性を利用すると, $$ \begin{align} c[n]&=\frac{1}{N}\sum_{k=0}^{N-1}\log S_k\exp\left(j\frac{2\pi kn}{N}\right)\\ &=\frac{1}{N}\sum_{k=0}^{N-1}\log S_k\cos\left(\frac{2\pi kn}{N}\right)\\ &=\frac{1}{N}\sum_{k=0}^{N-1}\log V_k\cos\left(\frac{2\pi kn}{N}\right)+\frac{1}{N}\sum_{k=0}^{N-1}\log H_k\cos\left(\frac{2\pi kn}{N}\right) \end{align} $$ と,変換できます.この$c$は,ケプストラム(cepstrum:spectrumのアナグラム)と呼ばれます.ケプストラムの単位はケフレンシー(quefrency:frequencyのアナグラム)です.

ケプストラム-ケフレンシー図を書くと包絡成分$H_k$が低いケフレンシー部分に,声門波成分$V_k$が高ケフレンシー部分に現れます.そのため,低ケフレンシー部分を取り出すローパスリフタリング(liftering:filteringのアナグラム)をして,低ケフレンシー部分を取り出します.一般にはケフレンシーの低い方から10-15個の値が用いられます.

特にメルフィルタバンクの対数スペクトラルパワーのケプストラムはメルケプストラムや,メル周波ケプストラム係数(MFCC:Mel Frequency Cepstral Coefficient)と呼ばれます.

from scipy.fftpack.realtransforms import dct
def get_cepstrum(mp, mspec):
    ceps = dct(ms, type=2, norm="ortho", axis=1)
    return ceps

v_data = get_wave_profile("voice.wav")
f, gs = stft_hamming(*v_data, res=1024, intv=64)
mel_data = mel_filter(f, gs, 0, 25000, 80)
ceps = get_cepstrum(*mel_data)
cep = ceps[180]
plt.plot(range(len(cep)), cep)
plt.ylim((y.min() - 5, y.max() + 5))
# 低い方から12個までを特徴量として利用します
thres = 12
plt.vlines(thres - 1, -100, 100, color="r")
print(cep[:thres])
# [45.15669001  1.60017792 -1.15822766  1.2903854  -2.15728379  1.93866199
# -1.27352706 -1.64276638 -0.81609222 -0.20153028 -0.99046816  0.63004835]

これで,ようやくメルケプストラムが抽出されます.お疲れ様でした...

差分特徴量

音声の動的な特徴を捉えることも音声認識には有効な手段であることが知られています.そのためには,2通りの方法が考えられます.

一方は,隠れマルコフモデル(HMM)や超短期記憶(LSTM)などといった時系列のデータの特徴を学習するモデルを使う手段で,もう一方は,フレーム間の特徴量の推移を新たな特徴量とする手段です.

具体的には隣接フレーム間で,ケプストラムの差分を取りそれを新たな特徴量とさせます.これはデルタケプストラムや,差分ケプストラムなどと呼ばれています. $$ \Delta c_n(t)=\frac{\sum_{k=-K}^{K}kc_n(t+k)}{\sum_{k=-K}^{K}k^{2}} $$ 一般には$K=2\sim5$の値を用います.また差分の差分を新たに特徴量として加える場合があります.


まとめ

音声解析では,生の信号で扱われることは少なく,様々な分析を経てなるべく特徴を保持したまま次元削減して利用しています.具体的には

  • AD変換: マイクでの録音の際に離散化されます.16kHz前後が多い.
  • 高域強調:パワーの弱い高音域を強調します.
  • STFT: 信号に対し,窓関数をずらしてフーリエ解析を行います.
  • メルフィルタバンク: 次元を削減してバリアンスを抑えます.
  • ケプストラム変換: 信号の包絡部分をケプストラム領域に逆離散コサイン変換をして抽出します.
  • 差分特徴量: 時系列変化の特徴を追加で考慮します.

という過程を経て利用されます.

誤植等があればご報告お願いいたします.

参考資料

余談

100日こんな調子は絶対無理なので2日目以降はボチボチやります.