ケプストラム分析によるスペクトル包絡推定

ケプストラム分析を用いてスペクトル包絡推定を行いました。参考文献 [1] によると、ケプストラム分析は1964年にNoll らにより提案された方法です。ちなみにケプストラム(cepstrum)という名称は、spectrumのアナグラム(specだけを逆から書くとcepstrum)となっています。

ケプストラム分析

ケプストラム分析について説明します。スペクトル包絡について説明してから、ケプストラム分析の考え方、リフタリングについて説明しようと思います。

スペクトル包絡とは

ソースフィルタモデルで母音を音声合成 で説明しましたように、ソースフィルタモデルではz変換された音声 Y は以下のように表せます。

\displaystyle{
Y(z) = H(z) X(z)
}

ここで、H(z) は声道の伝達関数、X(z) は声帯振動で得られた音源をz変換したものです。

つまり、音声の周波数特性 Y(ω) は以下のように表せるということです。

\displaystyle{
Y(\omega) = H(\omega) X(\omega)
}

H(ω) は声道フィルタの周波数特性、X(ω) は音源の周波数特性です。

声道フィルタの振幅特性 |H(ω)| についてはスペクトル包絡とも呼ばれています。

今回はケプストラム分析を用いて、実際の音声からこのスペクトル包絡を推定します。

ケプストラム分析の考え方

ケプストラム分析の考え方は意外に単純です。「あ」を発音した音声データの対数振幅スペクトルを図1に示します。ちなみに0dBに対応するのはフルスケールの正弦波のパワーです。

図1:「あ」の対数振幅スペクトル
図1:「あ」の対数振幅スペクトル

ソースフィルタモデルから音声の対数振幅スペクトルは以下のように表せます。

\displaystyle{
\begin{align}
\log{(|Y(\omega)|)} &= \log{(|H(\omega) X(\omega)|)}\\
&=  \log{(|H(\omega)|)} +  \log{(|X(\omega)|)}
\end{align}
}

つまり、音声の対数振幅スペクトルは声道フィルタの対数振幅スペクトルと音源の対数振幅スペクトルの加算で表せるということです。

声帯振動波形はパルス列のようになりますので、音源の対数振幅スペクトルも基本周波数の間隔でパルスが立つパルス列となります。つまり、図1が周期的にギザギザになっているのは音源の特性が原因と考えられます。

一方、声道フィルタの対数振幅スペクトルは緩やかな形状を持ちますので、さらにフーリエ変換(DFT)することで声道特性と音源特性が低域と高域に分離できるのではないかというのがケプストラム分析の考え方です。

音声の対数振幅スペクトルをさらにDFTしたものをケプストラムと呼びます。また、ケプストラムの横軸についてはケフレンシー(quefrency はfrequecyのアナグラムです)軸と呼びます。

補足:参考文献[1][2] では、対数振幅スペクトルをDFTではなくIDFTしたものをケプストラムとしています。理由としては、対数振幅スペクトルが偶関数であるため、以下のようにIDFTとDFTが係数を除いて、同じになるからです。

\displaystyle{
\begin{align}
{\rm IDFT}[\log{|Y[k]|}]&= \frac{1}{N} \sum_{k=0}^{N-1} \log{|Y[k]|} e^{j 2\pi kn/N} \\
&=  \frac{1}{N} \sum_{k=0}^{N-1} \log{|Y[N-k]|} e^{-j 2\pi kn/N} \\
&=  \frac{1}{N} \sum_{k=0}^{N-1} \log{|Y[-k]|} e^{-j 2\pi kn/N} \\
&= \frac{1}{N} \sum_{k=0}^{N-1} \log{|Y[k]|} e^{-j 2\pi kn/N}  \\
&= \frac{1}{N} {\rm DFT}[\log{|Y[k]|}]
\end{align}
}

また、スペクトルを逆フーリエ変換したと考えると横軸は時間 [s] とみなせるので、音声分野では逆フーリエ変換したものをケプストラムと呼ぶ説明が多いようです。これ以降の説明ではIDFTしたものをケプストラムとします。

リフタリング

「あ」のケプストラムは図2のようになります。

図2:「あ」のケプストラム
図2:「あ」のケプストラム

音源特性は図2のように基本周期  T_o にピークがあり、声道特性は低ケフレンシーに集中するため、ケプストラムの低ケフレンシー部だけを抽出し、これをDFTすることでスペクトル包絡が推定できます。

ケプストラムの低ケフレンシー部だけを抽出することをリフタリング(liftering は filtering のアナグラムです)と呼びます。

ケプストラムを抽出する範囲については、参考文献やWEBに明確な記載はなかったので、私は基本周期の半分まで抽出しました。

プログラム

ケプストラムと推定したスペクトル包絡を出力するプログラムは以下のcepstrum.py です。

解析するWAVファイルをソースコードと同じディレクトリに入れて、wav_namestartを書き換えて、「python cepstrum.py」と実行すれば、ceps.png(ケプストラムのグラフ) と sp_env.png(スペクトル包絡のグラフ)が出力されます。

import soundfile as sf
import numpy as np
from scipy.fft import rfft, irfft
import scipy.signal as sg
import matplotlib.pyplot as plt

# パラメータ
wav_name = "aiueo.wav"   # 読み込むWAVデータの名前
ceps_name = "ceps.png"   # 出力するPNGデータの名前
spec_name = "sp_env.png" # 出力するPNGデータの名前
window = "hann"          # 窓関数の種類
start = 10000            # DFTする開始点
N = 2**(14)              # DFT点数

# WAVファイルを読み込む
x, fs = sf.read(wav_name)
w = sg.get_window(window, N) # 窓関数の作成
w = w / np.sum(w)            # 窓関数の補正

# 対数振幅スペクトルを求める
eps = np.finfo(np.float64).eps         # マシンイプシロン
spec = rfft(x[start:start+N]*w)        # DFTする
spec_log = np.log10(np.abs(spec)+eps)  # 対数変換
freq  = np.arange(N//2+1) * (fs / N)   # 周波数軸

# ケプストラムを求める
ceps = irfft(spec_log)       # IDFTする
t = np.arange(N)*1000 / fs   # ケフレンシー軸

# ケプストラムのグラフ出力
To_l = int(fs/800) # 基本周波数の推定範囲の上限を800Hzとする
To_h = int(fs/40)  # 基本周波数の推定範囲の下限を40Hzとする
ylim = np.max(ceps[To_l:To_h]) # 基本周期のピーク
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, ceps, c="red")
ax.set_ylim([-0.02,ylim+0.02]) # 基本周期のピーク+0.2まで表示
ax.set_xlim([0,30])            # 30msまで表示
ax.set_xlabel("Quefrency [ms]") 
ax.set_ylabel("log amplitude")
fig.savefig(ceps_name)  # ケプストラムを出力

# リフタリング
To = np.argmax(ceps[To_l:To_h])+To_l  # 基本周期の点数を求める
lifter = To//2                        # 基本周期の半分まで抽出
ceps[lifter:N-lifter+1] = 0           # リフタリング
sp_env = np.real(rfft(ceps))          # DFTして実部だけ取り出す

# スペクトル包絡のグラフ出力
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(freq, 10*spec_log, c="red")
ax.plot(freq, 10*sp_env, c="blue")
ax.set_xlim([0,4000])   # 4kHzまで表示
ax.set_ylim([-50,-10]) 
ax.set_xlabel("Frequency [Hz]") 
ax.set_ylabel("Amplitude [dB]")
fig.savefig(spec_name)  # スペクトル包絡を出力

7~13行目:読み込むWAVファイルや出力ファイルの名前、窓関数の種類、データの解析する位置、FFT点数などのパラメータを決めています。

15~18行目:WAVファイルの読み込みと窓関数の作成を行っています。窓関数の補正も行っています。

20~24行目:窓関数をかけてDFTしています。さらに対数振幅スペクトルを求めています。対数変換する際は、計算エラーが発生しないように微小量(マシンイプシロン)を加算してから対数変換しています。

26~28行目:対数振幅スペクトルをIDFTして、ケプストラムを求めています。

30~41行目:ケプストラムのグラフを出力しています。いい感じにグラフが表示されるようにy軸の上限を基本周期のピーク+0.2にしています。

43~47行目:基本周期の点数を求めて、基本周期の半分までケプストラムを抽出しています。基本周期は1.25ms (800Hz) ~25ms (40Hz) の間の最大値を探すことで求めています。また、リフタリングする際は、ケプストラムはN/2で偶対称となっているので、(N-lifter+1)~(N-1) のケプストラムについても残しておきます。

49~58行目:対数振幅スペクトルとスペクトル包絡のグラフを出力しています。

ケプストラム分析の結果

以下の「あいうえお」を発音した音声をケプストラム分析しました。

あいうえおの音声

変な音声になっているのは、記事末尾に記載しているデータベースの音声から「あ」「い」「う」「え」「お」を抜き出して、つなぎ合わせてるからです。

「あいうえお」をケプストラム分析してスペクトル包絡を推定した結果は図3のようになりました。

図3:「あいうえお」のスペクトル包絡
図3:「あいうえお」のスペクトル包絡

緑の点線は、参考文献 [3] に記載されていた男性33人のフォルマント周波数 F1~F3 の平均です。

緑の点線がある箇所にスペクトル包絡の山がだいたいあるので、上手く推定できているかな? 正直いって、これが上手く推定できていると断定はできませんが、多分上手くできているでしょう。

おわりに

ケプストラム分析を用いてスペクトル包絡の推定を行いました。実は今回ケプストラム分析を紹介したのは、次回の記事でスペクトル包絡を用いて、音声変換がしたかったためです。次回、音声変換について書きたいと思います。

参考文献
[1] 森勢将雅、”音声分析合成”、コロナ社、2018.
[2] 川村新、”音声音響信号処理の基礎と実践”、コロナ社、2021.
[3] 電子情報通信学会、”聴覚と音声”、コロナ社、1973.

使用したデータについて
この記事で信号処理した音声は髙橋弘太研究室の話速バリエーション型音声データベース(SRV-DB)で提供されているATR25文の読み上げ音声を使用させていただきました。

【音声データベースを提供している髙橋弘太研修室のページ】
http://www.it.cei.uec.ac.jp/SRV-DB/

お問い合わせフォーム プライバシーポリシー

© 2021 Setoti All rights reserved.