ケプストラム分析を用いてスペクトル包絡推定を行いました。参考文献 [1] によると、ケプストラム分析は1964年にNoll らにより提案された方法です。ちなみにケプストラム（cepstrum）という名称は、spectrumのアナグラム（specだけを逆から書くとcepstrum）となっています。
ケプストラム分析
ケプストラム分析について説明します。スペクトル包絡について説明してから、ケプストラム分析の考え方、リフタリングについて説明しようと思います。
スペクトル包絡とは
ソースフィルタモデルで母音を音声合成 で説明しましたように、ソースフィルタモデルは声帯の振動を音源、声道による影響をフィルタとして考えることができます。
数式で説明すると、音声 Y は以下のように表せるということです。
$$
Y(\omega) = H(\omega) X(\omega)
$$
ここで、Xは周波数領域における声帯振動、Hは声道フィルタの周波数特性です。
声道フィルタの振幅特性 \(|H(\omega)|\) についてはスペクトル包絡と呼ばれています。今回はケプストラム分析を用いて、実際の音声からこのスペクトル包絡を推定します。
ケプストラム分析の考え方
ケプストラム分析の考え方は意外に単純です。「あ」を発音した音声データの対数振幅スペクトルを図に示します。ちなみに0dBに対応するのはフルスケールの正弦波のパワーです。
ソースフィルタモデルから音声の対数振幅スペクトルは以下のように表せます。
$$
\begin{align}
\log{(|Y(\omega)|)} &= \log{(|H(\omega) X(\omega)|)}\\[5pt]
&= \log{(|H(\omega)|)} + \log{(|X(\omega)|)}
\end{align}
$$
つまり、音声の対数振幅スペクトルは声道フィルタの対数振幅スペクトルと音源の対数振幅スペクトルの加算で表せるということです。
声帯振動波形はパルス列のようになりますので、音源の対数振幅スペクトルも基本周波数の間隔でパルスが立つパルス列となります。つまり、図が周期的にギザギザになっているのは音源の特性が原因と考えられます。
一方、声道フィルタの対数振幅スペクトルは緩やかな形状を持ちますので、さらに対数振幅スペクトルをフーリエ変換（DFT）することで声道特性と音源特性が低域と高域に分離できるのではないかというのがケプストラム分析の考え方です。
音声の対数振幅スペクトルをさらにDFTしたものをケプストラムと呼びます。また、ケプストラムの横軸についてはケフレンシー（quefrency はfrequecyのアナグラムです）軸と呼びます。
補足：参考文献[1][2] では、対数振幅スペクトルをDFTではなくIDFTしたものをケプストラムとしています。理由としては、対数振幅スペクトルが偶関数であるため、以下のようにIDFTとDFTが係数を除いて、同じになるからです。
\newcommand{\sb}[1]{_{#1}}
\begin{align}
{\rm IDFT} [\log{|Y[k]|}] &= \frac{1}{N} \sum\sb{k=0}^{N-1} \log{|Y[k]|} e^{j 2\pi kn/N} \\
&= \frac{1}{N} \sum\sb{k=0}^{N-1} \log{|Y[N-k]|} e^{-j 2\pi kn/N} \\
&= \frac{1}{N} \sum\sb{k=0}^{N-1} \log{|Y[-k]|} e^{-j 2\pi kn/N} \\
&= \frac{1}{N} \sum\sb{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したものをケプストラムとします。
リフタリング
「あ」の対数振幅スペクトルをDFTしたケプストラムは以下の図のようになります。
音源特性は図のように基本周期にピークがあり、声道特性は低ケフレンシーに集中するため、ケプストラムの低ケフレンシー部だけを抽出し、これをDFTすることでスペクトル包絡が推定できます。
ケプストラムの低ケフレンシー部だけを抽出することをリフタリング（liftering は filtering のアナグラムです）と呼びます。
ケプストラムを抽出する範囲については、参考文献やWEBに明確な記載はなかったので、私は基本周期の半分まで抽出しました。
プログラム
ケプストラムと推定したスペクトル包絡を出力するプログラムは以下です。
解析するWAVファイルをソースコードと同じディレクトリに入れて、wav_nameとstartを書き換えて、「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) # スペクトル包絡を出力
20～24行目：対数変換する際は、計算エラーが発生しないように微小量（マシンイプシロン）を加算してから対数変換しています。
30～41行目：ケプストラムのグラフを出力しています。いい感じにグラフが表示されるようにy軸の上限を基本周期のピーク+0.2にしています。
43～47行目：基本周期の点数を求めて、基本周期の半分までケプストラムを抽出しています。基本周期は1.25ms (800Hz) ～25ms (40Hz) の間の最大値を探すことで求めています。また、リフタリングする際、ケプストラムはN/2で偶対称となっているので、(N-lifter+1)～(N-1) のケプストラムについても残しておきます。
ケプストラム分析の結果
「あいうえお」をケプストラム分析してスペクトル包絡を推定した結果を以下の図に示します。
緑の点線は、参考文献 [3] に記載されていた男性33人のフォルマント周波数 F1～F3 の平均です。
緑の点線がある箇所にスペクトル包絡の山がだいたいあるので、上手く推定できているかな？ 正直いって、これが上手く推定できていると断定はできませんが、多分上手くできているでしょう。
おわりに
ケプストラム分析を用いてスペクトル包絡の推定を行いました。実装した感想としては以下の点で扱いが難しいと感じました。
- 求めたスペクトル包絡がギザギザ（特に「あ」や「う」）なためピークの位置でフォルマント周波数が求められない
- リフタリングする範囲を決めるのに困る
実装は比較的簡単でしたので、簡易的に求めるのには良いかもしれません。
■参考文献
[1] 森勢将雅、”音声分析合成”、コロナ社、2018．
[2] 川村新、”音声音響信号処理の基礎と実践”、コロナ社、2021.
[3] 電子情報通信学会、”聴覚と音声”、コロナ社、1973．