FFTNetを読んで，メルケプストラム係数（Mel Cepstral Coefficients）が何かわからなかったので調べてみると，メル周波数ケプストラム係数（MFCC）ばかりが検索上位に来てしまい困った．なので，今後のためにもメルケプストラム係数のメモと，Pythonで出力する方法をまとめる．
また，メルケプストラムと表現している文献が多いので，以下ではメルケプストラムと省略する．

参考

メルケプストラムとは

統計的声質変換を行うための知識と手法では以下のように説明されている．

スペクトル包絡は音色をよく表現できる音響特徴ですが、次元が大きく扱いづらいという問題があります。スペクトル包絡は周波数方向に周期性があるので、周波数方向に逆フーリエ変換することで、表現能力をあまり損なわずに低次元な特徴量を得ることができます。これがケプストラムです。（ケプストラム【cepstrum】はスペクトル【spectrum】のアナグラムです。）

人間の聴覚を考慮することで、ケプストラムよりも表現能力を保持させた、メルケプストラムという音響特徴量があります。統計的声質変換にはこの特徴量がよく用いられます。サンプリング周波数によりますが、スペクトル包絡の次元数は512や1024と大きい一方、メルケプストラムの次元数は60ほどに小さくしても元と同程度の品質で復元できます。（別の概念でメル周波数ケプストラム係数MFCCというものがありますが、これはメルケプストラムとは異なるので注意。）

つまり

スペクトル包絡だとデータが大きいので，表現力をある程度保ちつつ情報量を圧縮したのがケプストラム
人間の聴覚特性を考慮したケプストラムがメルケプストラム

Pythonでの実装

メルケプストラムはSPTK（Signal　Processing　Toolkit）を使うと簡単に求められる．しかし，SPTKはLinux環境でないと動作しないため，pysptkを利用する．

pysptkでメルケプストラムを出力

pysptk.mcepで計算できるようだが，イマイチ使い方がわからなかったので，今後調べる．他に，pysptk.sp2mc(powerspec, order, alpha)という関数で出力が可能．pysptkのリポジトリを参考にそれぞれの入力を示す．

powerspec

スペクトル包絡（spectral envelope）を入力すれば良い．スペクトル包絡はWORLDを使えば算出可能．
こちらもPythonラッパーがあるのでこれを利用する．

order

メルケプストラムの次数．機械学習等のモデルパラメータに合わせて設定を行う

alpha

all-pass constantと呼ばれるもの．SPTK Reference manualを参照すると，サンプリング周波数に対して

16kHz：0.42
10kHz：0.35
8kHz：0.31

に設定すると良いとのこと

出力の確認

試しにLJ Speechを使って，メルケプストラムを出力してみる．

alphaの設定

SPTKのマニュアルには16kHzまでしか書いていなかったが，LJ Speechのサンプリング周波数は22.5kHzである．そこで，既知の数値から適切な値を予測した．サンプリング周波数に対して，alphaの値が対数っぽく変化してるなと考え，少し乱暴に計算してみた所$$alpha = \frac{1}{2.9}\log fs$$でそれっぽい値が算出できた．ただし，$fs[kHz]$とする．この式でalphaを求めると，0.46となったので，とりあえずこれを使ってメルケプストラムを求める

サンプルコード

サンプルコードを以下に示す．メルケプストラムの次数は25とした．

mcep.py

import numpy as np
import librosa
import pyworld as pw
import pysptk as sptk
import matplotlib.pyplot as plt

wav, fs = librosa.load('LJ001-0001.wav')
wav = wav.astype(np.float)

#スペクトル包絡算出にWORLDを使用
_f0, t = pw.dio(wav, fs)
f0 = pw.stonemask(wav, _f0, t, fs)
sp = pw.cheaptrick(wav, f0, t, fs)
ap = pw.d4c(wav, f0, t, fs)

#メルケプストラムに変換前のスペクトル包絡を可視化
plt.figure()
plt.imshow(np.log10(sp).T, origin='lower')#対数を取ることで見やすくなる
plt.colorbar()
plt.title('spectral envelope')

#メルケプストラムの算出及び可視化
mcep = sptk.sp2mc(sp, order=25, alpha=0.46)
plt.figure()
#スペクトル包絡に合わせてアスペクト比の調節
plt.imshow(mcep.T, aspect=mcep.shape[0]/mcep.shape[1]/3.7, origin='lower')
plt.colorbar()
plt.title('Mel-cepstrum coefficients')

#メルケプストラムからスペクトル包絡に変換及び可視化
sp_from_mcep = sptk.mc2sp(mcep, alpha = 0.46, fftlen = 1024)
plt.figure()
plt.imshow(np.log10(sp_from_mcep).T, origin='lower')
plt.colorbar()
plt.title('spectral envelope from mcep')

#WORLDでそれぞれ音声波形を合成し，保存
synthesized = pw.synthesize(f0, sp, ap, fs)
librosa.output.write_wav('./world.wav', synthesized, fs)

syn_mcep = pw.synthesize(f0, sp_from_mcep, ap, fs)
librosa.output.write_wav('./world_mcep.wav', syn_mcep, fs)

出力

WORLDから算出したスペクトル包絡