- 追加された行はこの色です。
- 削除された行はこの色です。
#freeze
#access
#analog
*メル周波数ケプストラム係数（ Mel Frequency Cepstral Coefficient ） [#d96a8b16]
#contents
**はじめに [#i52f1c12]
***スペクトラム（spectrum）とは [#u28deb5a]
-音声や地震波などの周期性のある信号は、どれだけ複雑な信号であっても、単純な波に分解できる（フーリエの定理）
--単純な波...単一の周波数と振幅をもつ正弦波、余弦波
-上記の定理に従って、ある信号の周波数成分と振幅成分を抽出したものが[[周波数スペクトラム（スペクトル）:http://ja.wikipedia.org/wiki/%E5%91%A8%E6%B3%A2%E6%95%B0%E3%82%B9%E3%83%9A%E3%82%AF%E3%83%88%E3%83%AB]]です((日本語では「スペクトル」が一般的ですが、 cepstrum が spectrum のアナグラムであることを意識して、ここではスペクトラムとしました（[["spectrum"と「スペクトラム」と「スペクトル」 - 教えて！goo:http://oshiete1.goo.ne.jp/qa1287977.html]]）))。
-スペクトラムの求め方（MATLABによる説明）
--以下のような複雑な波を考えます。
#ref(spectrum_1.jpg);
---MATLABで上の波を生成するには、以下のコードを実行
#geshi(matlab){{
Fs = 8820; % サンプリング周波数 8820Hz
time = 0 : 1 / Fs : 0.05;
sinwav_1 = 1.2 * sin(2 * pi * 130 * time); %振幅1.2、周波数130[Hz]の正弦波
coswav_1 = 0.9 * cos(2 * pi * 200 * time); %振幅0.9、周波数200[Hz]の余弦波
sinwav_2 = 1.8 * sin(2 * pi * 260 * time); %振幅1.8、周波数260[Hz]の正弦波
coswav_2 = 1.4 * cos(2 * pi * 320 * time); %振幅1.4、周波数320[Hz]の余弦波
wavdata = 1.4 + [sinwav_1 + coswav_1 + sinwav_2 + coswav_2];
plot(time*1000, wavdata); xlabel('時間[ms]'); ylabel('振幅');
}}
---正弦波、余弦波の詳細は [[MATLAB Note/音声の加工 / 純音を作る:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/index.php?MATLAB%20Note/%E9%9F%B3%E5%A3%B0%E3%81%AE%E5%8A%A0%E5%B7%A5#o225da40]] を参照
--フーリエ変換によって、複合音を構成する各純音の周波数、振幅（パワー）の値を調べます。
---'''[[フーリエ変換の詳細な説明（限定公開）:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki/index.php?%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B%EF%BC%88FFT%EF%BC%89]]'''
---MATLABでフーリエ変換を行うには、以下のコードを実行((参考文献：[[MATLAB fft:http://dl.cybernet.co.jp/matlab/support/manual/r13/toolbox/matlab/ref/?/matlab/support/manual/r13/toolbox/matlab/ref/fft.shtml]] [[複素共役:http://ja.wikipedia.org/wiki/%E8%A4%87%E7%B4%A0%E5%85%B1%E5%BD%B9]] [[フーリエ変換の実際:http://www009.upp.so-net.ne.jp/hachinami/note004/index.htm]]))
#geshi(matlab){{
%離散フーリエ変換
fftsize = 2048; %フーリエ変換の次数, 周波数ポイントの数
dft = fft(wavdata, fftsize);
%振幅スペクトル（dftの絶対値、conj:複素共役）
Adft = sqrt(dft .* conj(dft)); %Adft = abs(dft); としても同じ
%パワースペクトル（dftの絶対値の二乗）
Pdft = abs(dft) .^ 2;
%周波数スケール（サンプリング定理より、OHzからサンプリング周波数まで）
fscale = linspace(0, Fs, fftsize); %0Hz～8820Hzまでをfftsize個に分割
%プロット（表示するのはサンプリング周波数の半分まで）
subplot(2,1,1); plot(fscale(1:fftsize/2), Adft(1:fftsize/2));
xlabel('周波数[Hz]'); ylabel('振幅スペクトル'); xlim([0 1000]);
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft(1:fftsize/2));
xlabel('周波数[Hz]'); ylabel('パワースペクトル'); xlim([0 1000]);
}}
---fftの結果は複素数で与えられます。実数部分と虚数部分の絶対値を振幅スペクトル、絶対値の二乗をパワースペクトル、偏角を位相スペクトルとよびます((参考文献：[[MATLABプログラミング:http://arx.ee.utsunomiya-u.ac.jp/azuma/courses/04Elab2/main.pdf]]、[[フーリエ変換の分かりやすい説明:http://bluedb.org/archives/73]]、''[[小野測器 - デジタル信号処理の基礎 - ８ ”「いろいろなパワースペクトル」”:http://www.onosokki.co.jp/HP-WK/eMM_back/03_05_29.htm]]''))。
--結果
#ref(spectrum_2.png);
---それぞれの波の周波数と振幅の値が取得できています。((左端の大きな値はスペクトルの直流成分で、音声のエネルギー（音圧）と相関の高い値です。))
---フーリエ変換は、''複雑な波形の周期成分と振幅成分を取り出す手法である''といえます((より詳細な解説： [[デジタル信号処理の例：Matlab によるFFT:http://www-akaz.ist.osaka-u.ac.jp/~pak/lecture/HumanInfEng/FFTexample.pdf]] [[特別研究I: MATLAB演習（2005 年度版）:http://winnie.kuis.kyoto-u.ac.jp/~kitahara/memos/matlab-2005.pdf]]))。
---&color(red){音声のサンプリング周波数が変われば、縦軸の値も変わります。};
---[[音声を短時間ごとに区切り:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#o5c30647]]、各フレームのスペクトラムを時間軸に沿って並べたものが[[スペクトログラム:http://ja.wikipedia.org/wiki/%E3%82%B9%E3%83%9A%E3%82%AF%E3%83%88%E3%83%AD%E3%82%B0%E3%83%A9%E3%83%A0]]です。
--【補記】縦軸をデシベル表記に変えてみます。
#geshi(matlab){{
%プロット（デシベル）
subplot(2,1,1); plot(fscale(1:fftsize/2), 20*log10(Adft(1:fftsize/2)));
xlabel('周波数[Hz]'); ylabel('振幅スペクトル[dB]'); xlim([0 1000]);
subplot(2,1,2); plot(fscale(1:fftsize/2), 10*log10(Pdft(1:fftsize/2)));
xlabel('周波数[Hz]'); ylabel('パワースペクトル[dB]'); xlim([0 1000]);
}}
#ref(spectrum_2dB.png,,80%);
---振幅スペクトルは&color(red){20};×log、パワースペクトルは&color(red){10};×log((&color(red){2010/12/18 修正しました。'''[[ご指摘:http://shower.human.waseda.ac.jp/~m-kouki/pass/20101215_mfcc_refer.txt]]'''ありがとうございました!};))がデシベルの値になります。((参考： ''[[周波数特性 (ディジタル信号処理):http://ufcpp.net/study/dsp/frequency.html]]''))
---※音圧の単位について [[Okumura's Blog 波形グラフの縦軸？:http://oku.edu.mie-u.ac.jp/~okumura/blog/node/2426]]
***逆フーリエ変換とは [#l34aa4d5]
-フーリエ変換の結果から元の波形を再現することを、逆フーリエ変換とよびます。
--MATLABでdftの値から波を生成するには、上に続けて以下のコードを実行
#geshi(matlab){{
newdata = ifft(dft);
%元の音声データと同じデータ数（442）でプロットする
plot(time*1000, newdata(1:length(wavdata))); xlabel('時間[ms]'); ylabel('振幅');
}}
--結果
#ref(spectrum_3.jpg);
---元の波形を再現できています。
--複雑な波にフーリエ変換を行うことで、時間・空間軸から周波数・空間軸に変換することができ、また逆フーリエ変換を行うことで元の時間・空間軸に再変換できます。
**音声波形の特性 [#yc248b2b]
***音声のスペクトラム [#n36c0d49]
-まず、音声波形のスペクトラムを求めてみます。
--&ref(a.wav); を使って以下のコードを実行
#geshi(matlab){{
[wav, Fs] = wavread('a.wav');
t = 0 : 1/Fs : length(wav)/Fs;
plot(t(1:length(t)-1)*1000, wav); xlabel('時間[ms]'); ylabel('振幅');
}}
--結果
#ref(cepstrum_1.jpg);
-母音の定常部（中心部分）40ms（0.04秒）を切り出し、[[ハニング窓:http://ja.wikipedia.org/wiki/%E7%AA%93%E9%96%A2%E6%95%B0#.E3.83.8F.E3.83.B3.E7.AA.93]]をかけて丸めたもの（[[フレーム:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#gcaedeef]]）を分析対象にします。
--上に続けて以下のコードを実行
#geshi(matlab){{
%中央部分を切り出し
center = fix(length(wav) / 2); %中心のサンプル番号
cuttime = 0.04; %切り出す長さ[s]
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(cuttime/2*Fs));
time = t(center-fix(cuttime/2*Fs) : center+fix(cuttime/2*Fs));
subplot(2,1,1); plot(time*1000, wavdata); ylabel('振幅');
%ハニング窓をかける
han_window = 0.5 - 0.5 * cos(2 * pi * [0 : 1/length(wavdata) : 1]);
wavdata = han_window(1:length(wavdata))' .* wavdata;
subplot(2,1,2); plot(time*1000, wavdata); xlabel('時間[ms]'); ylabel('振幅');
% 切り出したフレームの再生
wavplay(wavdata, Fs);
}}
--結果
#ref(cepstrum_2.jpg);
--切り出した音声のスペクトラムを求めます((スペクトラムを求める前に、高域強調の処理（プリエンファシス）を行うことが一般的ですが、ここでは飛ばしました。プリエンファシスのやり方は[[音声波形のフレーム分割および前処理 / 高域強調:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#q962d911]]を参照してください。))。
#geshi(matlab){{
%離散フーリエ変換
fftsize = 2048; %フーリエ変換の次数, 周波数ポイントの数
dft = fft(wavdata, fftsize); % 実部がcos波成分、虚部がsin波成分
%パワースペクトル（実部と虚部をそれぞれ二乗して合計）
Pdft = (real(dft) .^ 2) + (imag(dft) .^ 2);
%振幅スペクトル（パワースペクトルの平方根）
Adft = sqrt(Pdft);
%周波数スケール（サンプリング定理より、OHzからサンプリング周波数まで）
fscale = linspace(0, Fs, fftsize); %0～Fsまでをfftsize個に分割
%プロット（表示するのはサンプリング周波数の半分まで）
subplot(2,1,1); plot(fscale(1:fftsize/2), Adft(1:fftsize/2));
xlabel('周波数[Hz]'); ylabel('振幅スペクトル'); xlim([0,5000]);
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft(1:fftsize/2));
xlabel('周波数[Hz]'); ylabel('パワースペクトル'); xlim([0,5000]);
}}
---フーリエ変換の次数が大きいと解析結果が細かくなります(([[fftsize = 1024のときと8192のときの比較:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/index.php?plugin=attach&pcmd=open&file=fft_point1.png&refer=%E3%83%A1%E3%83%AB%E5%91%A8%E6%B3%A2%E6%95%B0%E3%82%B1%E3%83%97%E3%82%B9%E3%83%88%E3%83%A9%E3%83%A0%EF%BC%88MFCC%EF%BC%89]]。また、フレーム窓長が違っても解析結果の細かさに差が出ます（[[例:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/index.php?plugin=attach&pcmd=open&file=fft_point2.png&refer=%E3%83%A1%E3%83%AB%E5%91%A8%E6%B3%A2%E6%95%B0%E3%82%B1%E3%83%97%E3%82%B9%E3%83%88%E3%83%A9%E3%83%A0%EF%BC%88MFCC%EF%BC%89]]）))。
---振幅スペクトルはフーリエ変換の結果（複素数で表現されます((フーリエ変換結果の実部・虚部の意味については [[フーリエ変換の実際（信号処理 ランダム・ウォーク）:http://www009.upp.so-net.ne.jp/hachinami/note004/index.htm]] や [[4.1 フーリエ級数(Fourier Series)/変換(Fourier Transform):http://www2.ttcn.ne.jp/~ymizu/hkiso2.html]]、[[第4 回MATLAB 講習会-信号の表現-:http://www.fl.ctrl.titech.ac.jp/~old/home/html/course/SS/01SS/matlab/matlab4.pdf]] が参考になります。))）の絶対値です。MATLABでは「Adft = abs(dft);」と書いても求められます。その場合、パワースペクトルは「Pdft = abs(dft) .^ 2;」（振幅スペクトルの二乗）と書けます。
--結果
#ref(cepstrum_3.png);
---母音/a/に含まれる音声の周波数成分と振幅成分が分かりました。
---&color(red){横軸のスケール合ってる? [[フォルマント解析:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/index.php?MATLAB%20Note/%E9%9F%B3%E5%A3%B0%E3%81%AE%E5%88%86%E6%9E%90#y6a12906]] の結果も要確認!};
--なお、振幅スペクトルを求める過程で、信号の位相の情報が除去されます。解析をするだけなら問題ないですが、音声の復号を行う場合は位相の情報が必要です。そのため、ここで位相スペクトルも計算しておきます。
#geshi(matlab){{
Angdft = atan2(imag(dft), real(dft));
}}
---MATLABでは「Angdft = angle(dft)」と書いても求められます。
***音声の対数スペクトラム [#d33d7f39]
-音声信号は、音源特性（基本周波数）と声道特性（共鳴周波数）の畳み込みである、と考えることができます。
#ref(cepstrum_7.jpg);
--そこで、音声信号の振幅スペクトルをS(ω)、音源信号の振幅スペクトルをG(ω)、声道の伝達特性をH(ω)とすると、以下のように表現できます((参考文献：[[インパルス応答の測定:http://www.kyushu-id.ac.jp/~samejima/aip/impulse_response.pdf]]))。
#mimetex(|S(\omega)| = |G(\omega)| |H(\omega)| = \int_{\infty}^{-\infty} |G(\omega)||H(\omega-\tau)|d\tau );
---ωは周波数をあらわします。スペクトルは複素数で表現され、絶対値が周波数の大きさ、偏角が波の位相をあらわしています((参考文献：[[スペクトル:http://www.asp.c.dendai.ac.jp/courses/basic/kihon_spectrum01.pdf]]))
---τの値は...((畳み込みとは：[[初心者用 畳み込み（たたみこみ）解説:http://www.ice.tohtech.ac.jp/~nakagawa/laplacetrans/convolution1.htm]]))
--[[対数の基本的な演算:http://ja.wikipedia.org/wiki/%E5%AF%BE%E6%95%B0#.E5.9F.BA.E6.9C.AC.E7.9A.84.E3.81.AA.E6.BC.94.E7.AE.97]]より、両辺の対数をとると、
#mimetex(log|S(x)| = log|G(x)| + log|H(x)|);
---複数のスペクトルの畳み込み状態を、対数スペクトルの足し算で表すことができました((参考文献：[[ケプストラム分析:http://vision.kuee.kyoto-u.ac.jp/lecture/dsp/dictionary/cepstrum.htm?keepThis=true&TB_iframe=true&height=400&width=500]]))。
-対数スペクトルをプロットしてみます。
--上に続けて以下のコードを実行
#geshi(matlab){{
%対数振幅スペクトル
Adft_log = log10(abs(dft));
%対数パワースペクトル
Pdft_log = log10(abs(dft) .^ 2);
%プロット
subplot(2,1,1); plot(fscale(1:fftsize/2), Adft_log(1:fftsize/2));
xlabel('周波数[Hz]'); ylabel('対数振幅スペクトル'); xlim([0,5000]);
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft_log(1:fftsize/2));
xlabel('周波数[Hz]'); ylabel('対数パワースペクトル'); xlim([0,5000]);
}}
--結果
#ref(cepstrum_4.png);
---&color(red){負値の意味は?};
---対数にかけることで、音源に由来する成分と声道に由来する成分が「畳み込まれた」波が、「足し合わされた」波に変換されました。
---フーリエ変換の結果をよく見ると、大局的な変動（声道の伝達特性）に加えて、規則正しく周波数成分のピークが繰り返され（音源信号）、周期的なパターンを示していることが分かります。すなわち、フーリエ変換の結果も周期性のある信号であるといえます((参考文献：B. P. Bogert, M. J. R. Healy, & J. W. Tukey: The quefrency alanysis of time series for echoes: cepstrum, pseudo-autocovariance, cross-cepstrum, and saphe cracking, Proceedings of the Symposium on Time Series Analysis, pp.209-243, 1963. レイ・D・ケント, リチャーズ・リード, 音声の音響分析, 海文堂, 1996．))。
**ケプストラム（cepstrum） [#qee1e455]
***ケプストラム分析（alanysis） [#x3efddeb]
-対数スペクトルをもう一度フーリエ変換にかけて、音源の周波数と声道の周波数を切り分けることを考えます。
--振幅スペクトル（音声信号→フーリエ変換→絶対値）に対数をかけて、さらにもう一度フーリエ変換をかけたもの（スペクトラムのスペクトラム）を、ケプストラムとよびます。((文献によっては、「パワースペクトルに対数をかけてフーリエ変換する」と説明されているものもあります。 ))
--フーリエ変換のフーリエ変換は、フーリエ変換の結果に対して逆フーリエ変換をかけるのと同じ処理になります。
---ただし、絶対値をとっているので、（実数ケプストラムでは）位相スペクトルの情報は消えます。
--MATLABで（対数）ケプストラムを生成するには、上に続けて以下のコードを実行します。((&color(red){2012/07/10 ケフレンシー軸の計算方法を修正しました。};))
#geshi(matlab){{
%対数振幅スペクトルに逆フーリエ変換を実行して実部のみを取り出す
cps = real(ifft( Adft_log ));
%ケフレンシー軸、0[s]～切り出した音声の長さまでをサンプル数で分割
quefrency = linspace(0, cuttime, cuttime*Fs);
%ケプストラムをプロット
plot(quefrency(1:fftsize/2)*1000, cps(1:fftsize/2));
xlabel('ケフレンシー[ms]'); ylabel('対数振幅ケプストラム');
}}
// quefrency = time - min(time); %ケフレンシー軸、0[s]～切り出した音声の長さ[s]まで
--結果
#ref(cepstrum_5.png);
---横軸の[[ケフレンシー（quefrency）:http://ja.wikipedia.org/wiki/%E3%82%B1%E3%83%97%E3%82%B9%E3%83%88%E3%83%A9%E3%83%A0#quefrency]]は周波数成分の逆数（時間の尺度）です。
---ケフレンシー上のピークは、対数スペクトル上で周期的に出現する周波数成分（倍音）をあらわしています。
--低ケフレンシー成分を拡大すると、以下のようになります。
#ref(cepstrum_6.jpg);
---対数スペクトルのケプストラムを求めたことで、高周期で変動する成分（音源成分）と、畳み込まれていた声道特性の成分が分離できています...((参考文献：[[ケプストラム:http://hil.t.u-tokyo.ac.jp/~sagayama/applied-acoustics/2009/B3-Cepstrum.pdf]]))
---基本周期はケフレンシー軸上で5.4[ms]くらいのところにあらわれているので、逆数をとって（1/0.0054）、約185[Hz]が基本周波数（f0）の値です（&color(red){要説明!};）((この方法でf0を求める手法を、ケプストラムピッチ抽出法といいます。ケプストラムピッチ抽出法は高精度であるが、雑音に弱いといわれていて、改良法が提案されています。詳細は MATLAB マルチメディア信号処理 下，池原 雅章, 真田 幸俊, 島村 徹也，培風館 p.40～ を参照。))。
---一番左の値（ケプストラムの0次の係数、C0）は、直流成分をあらわします（'''[[フーリエ変換 a_0の抽出:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki/index.php?%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B%EF%BC%88FFT%EF%BC%89#hf56c43d]]'''）。
---0次の係数を覗いた残りの値は、左右対称の関係にあります。
>> cps(1:6)
ans =
-1.4016
0.3868
0.0788
0.0445
0.0272
0.0952
>> cps( length(cps)-4 : length(cps) )
ans =
0.0952
0.0272
0.0445
0.0788
0.3868
***リフタリング（liftering） [#aa1863e8]
-上で求めた対数ケプストラムから声道特性の部分（スペクトル包絡）を抽出してみます。
--ケプストラム係数の低次成分だけを残し、高次成分を0にした上でFFTによってスペクトラムに再変換します。
#geshi(matlab){{
cps_lif = cps;
cps_lif(120:length(cps)-118) = 0; %リフタリング
%リフタリング前のパワースペクトルをプロット
dftSpc = fft(cps, fftsize); %リフタリング前のケプストラムの離散フーリエ変換（対数振幅スペクトルに戻す）
AdftSpc = abs(10 .^ dftSpc); %対数振幅スペクトルを振幅スペクトルに戻す
PdftSpc = AdftSpc .^ 2; %振幅スペクトルをパワースペクトルに変換
subplot(2,1,1); plot(fscale(1:fftsize/2), PdftSpc(1:fftsize/2));
xlabel('周波数[Hz]'); xlim([0,5000]);
%リフタリング後のパワースペクトルをプロット
dftSpc = fft(cps_lif, fftsize);
AdftSpc = abs(10 .^ dftSpc);
PdftSpc = AdftSpc .^ 2;
subplot(2,1,2); plot(fscale(1:fftsize/2), PdftSpc(1:fftsize/2));
xlabel('周波数[Hz]'); xlim([0,5000]);
}}
---&color(red){2011/07/20 パワースペクトルの計算式を修正しました。ご指摘ありがとうございました！};((&color(red){'''[[ご指摘:http://shower.human.waseda.ac.jp/~m-kouki/pass/20110601_mfcc_refer.zip]]'''ありがとうございました!};))
---0次成分+1～120次成分を残しました。
---このとき、対称形になっている同じ値の部分も残します((参考文献：[[MATLABによる音声処理 9 母音の音色の分析（法政大伊藤先生）:http://www.slp.k.hosei.ac.jp/~itou/lecture/2009/SpeechProcessing/06_text.pdf]]))。
--結果
#ref(cepstrum_9_refer.png);
---上がリフタリングなし、下がリフタリング後
---基本周波数の成分が除去されて、スペクトル包絡が抽出できています((したがって、ケプストラムの低次成分を取り出す処理は、[[対数スペクトルに対するLPC分析:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#y6a12906]] と似た処理であるといえます。音声認識では、LPC対数スペクトルに対してケプストラムを求めて、より高精度に声道成分の抽出を行う場合もあります（IT Text 音声認識システム, オーム社, 2006.）))。
//---スペクトルの絶対値（振幅スペクトル）から作った実数ケプストラムは位相の情報が削除されているので、再変換しても元のスペクトルは再現されません。元通りに変換するには、[[複素ケプストラム:http://www.mathworks.com/access/helpdesk_ja_JP/help/toolbox/signal/cceps.html]]を使う必要があります。
***音声再合成 [#f411ed6a]
-リフタリングによって音源信号成分を除去した振幅スペクトルから音声波形を生成してみます。
--リフタリング後の振幅スペクトルと、[[FFTの分析:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/66.html#n36c0d49]] の時に求めた位相スペクトルを合わせて、sin波成分とcos波成分に戻し(([[信号処理工学II 第３回：音声信号処理その２～処理編～:http://apple.ee.uec.ac.jp/~spe2/handout/spe2_handout3.pdf]]（電気通信大長井先生）を参考にしました))、逆フーリエ変換をして実部を取り出したものが音声波形です。
#geshi(matlab){{
% sin波成分とcos波成分に戻す
sound_fft = AdftSpc .* exp(Angdft*j);
% 逆フーリエ変換をして実部を取り出す
sound = real(ifft(sound_fft, fftsize));
% プロット
subplot(2,1,1); plot(wavdata); title('元の波形');
subplot(2,1,2); plot(sound(1:fftsize/2)); title('再合成した波形');
% 再合成した波形の再生
wavplay(sound, Fs);
}}
---&color(red){2012/12/18 上記スクリプトの2行目は「sound_fft = AdftSpc .* exp(Angdft''*j'');」が正しいです（要虚数変換）。修正しました。};
--結果
#ref(cepstrum_8.png);
---波形の形は違いますが、声道成分が保存されているため、再合成音声も/a/に聞こえます。
-なお、ここでは音源信号成分を除去しましたが、逆に音源信号のみを残し、声道特性成分を除外することもできます。
--リフタリングによって残す成分の次数を逆にしてみます。
#geshi(matlab){{
%リフタリング
cps_lif_P = cps;
cps_lif_P(2:120) = 0;
cps_lif_P(length(cps)-118:length(cps)-1) = 0;
%リフタリング前のパワースペクトルをプロット
dftSpc_P = fft(cps, fftsize); %リフタリング前のケプストラムの離散フーリエ変換（対数振幅スペクトルに戻す）
AdftSpc_P = abs(10 .^ dftSpc_P); %対数振幅スペクトルを振幅スペクトルに戻す
PdftSpc_P = AdftSpc_P .^ 2; %振幅スペクトルをパワースペクトルに変換
subplot(2,2,1); plot(fscale(1:fftsize/2), PdftSpc_P(1:fftsize/2));
xlabel('周波数[Hz]'); xlim([0,5000]); title('リフタリング前のパワースペクトル');
%リフタリング後のパワースペクトルをプロット
dftSpc_P = fft(cps_lif_P, fftsize);
AdftSpc_P = abs(10 .^ dftSpc_P);
PdftSpc_P = AdftSpc_P .^ 2;
subplot(2,2,3); plot(fscale(1:fftsize/2), PdftSpc_P(1:fftsize/2));
xlabel('周波数[Hz]'); xlim([0,5000]); title('リフタリング後のパワースペクトル');
% sin波成分とcos波成分に戻す
sound_fft_P = AdftSpc_P .* exp(Angdft*j);
% 逆フーリエ変換をして実部を取り出す
sound_P = real(ifft(sound_fft_P, fftsize));
% プロット
subplot(2,2,2); plot(wavdata); title('元の波形');
subplot(2,2,4); plot(sound_P(1:fftsize/2)); title('再合成した波形');
% 再合成した波形の再生
wavplay(sound_P, Fs);
}}
--結果
#ref(cepstrum_11.png,,80%);
--なお、音声の音源信号成分のみを取り出したいのであれば、ここで紹介したケプストラムを使った方法よりもローパスフィルタを使うほうが簡単です。
-フレームが複数ある連続音声の再合成を行う場合は、各フレームごとに音声を再合成して足しあわせる必要があります。
--[[MATLAB Note/音声の加工/フレーム同士の結合:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#ae15abec]] を参照して下さい。
--リフタリングなどの処理によって振幅スペクトルの形状が変わっている場合、位相スペクトルとの不一致によって再合成音声にノイズが混ざります。これを解決するためにLSEE-MSTFTMアルゴリズムが提案されています。[[MATLAB Note/音声の加工/LSEE-MSTFTMアルゴリズム:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#b7c530cc]] を参照して下さい。
---これらの手法によって [[声質を維持したピッチシフト:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/75.html#y517f110]] などができます。
--[[MATLAB Note/音声の加工/フレーム同士の結合:http://speechresearch.fiw-web.net/75.html#ae15abec]] を参照して下さい。
--リフタリングなどの処理によって振幅スペクトルの形状が変わっている場合、位相スペクトルとの不一致によって再合成音声にノイズが混ざります。これを解決するためにLSEE-MSTFTMアルゴリズムが提案されています。[[MATLAB Note/音声の加工/LSEE-MSTFTMアルゴリズム:http://speechresearch.fiw-web.net/75.html#b7c530cc]] を参照して下さい。
---これらの手法によって [[声質を維持したピッチシフト:http://speechresearch.fiw-web.net/75.html#y517f110]] などができます。
**メル周波数ケプストラム [#o8763891]
-対数ケプストラムの低次成分は、音声のスペクトル包絡（声道成分に由来した周波数特性）を表現できていることがわかりました。この対数ケプストラムの低次成分に対して、ヒトの周波数知覚特性を考慮した重み付けをした特徴量を、メル周波数ケプストラム係数（MFCC）と呼びます。
***メル尺度（メル対数、メル周波数） [#i7714922]
-[[聴覚尺度 / メル尺度:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/index.php?%E8%81%B4%E8%A6%9A%E5%B0%BA%E5%BA%A6#lbf3c463]] を参照してください。
-試しに、[[音声の対数スペクトラム:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/66.html#d33d7f39]] に戻り、対数スペクトラムの横軸（周波数軸）をメル周波数に変換してみます。((参考文献：[[単旋律楽曲における機械学習を用いた楽器音の音源同定:http://hawaii.naist.jp/~mizuki-i/data/papers/sigchallenge07.pdf]]))
#geshi(matlab){{
%対数振幅スペクトル
Adft_log = log10(abs(dft));
%対数パワースペクトル
%Pdft_log = log10(abs(dft) .^ 2);
%プロット
subplot(2,1,1); plot(fscale(1:fftsize/2), Adft_log(1:fftsize/2));
xlabel('周波数[mel]'); ylabel('対数振幅スペクトル'); xlim([0,5000]);
subplot(2,1,2); plot(mellog(fscale(1:fftsize/2)), Pdft_log(1:fftsize/2));
%横軸をメルスケールに変換してプロット
subplot(2,1,2); plot(mellog(fscale(1:fftsize/2)), Adft_log(1:fftsize/2));
xlabel('周波数[mel]'); ylabel('メル対数振幅スペクトル'); xlim([0,mellog(5000)])
}}
--結果
#ref(cepstrum_10.png);
---メル対数スペクトラム（下段）では、対数スペクトラム（上段）の低周波数成分が引き伸ばされて、高周波数成分が圧縮されていることがわかります。
***メル周波数スペクトラム [#x7e138e5]
-振幅スペクトルを[[メル尺度上で等間隔なフィルタバンク:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/index.php?MATLAB%20Note/%E9%9F%B3%E5%A3%B0%E3%81%AE%E5%88%86%E6%9E%90#x126ec36]]にかけて、[[各帯域のスペクトル成分を取り出します:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/index.php?MATLAB%20Note/%E9%9F%B3%E5%A3%B0%E3%81%AE%E5%88%86%E6%9E%90#qeaeca52]]。((詳細は [[MATLAB Note/音声の分析/フィルタバンク解析/メルフィルタバンク:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#x126ec36]] 以下を参照して下さい。))
--音声認識では、経験的に20個の帯域のフィルタバンクを使います。
-（メル周波数軸上で等間隔な）各帯域の振幅スペクトルの和をとり、20次元の振幅スペクトルに圧縮します。
-20次元に圧縮された振幅スペクトルの対数をとり、対数振幅スペクトルにします。((詳細は [[MATLAB Note/音声の分析/フィルタバンク解析/メルフィルタバンク解析（振幅スペクトルを分割するやり方）:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#v962bcc1]] を参照して下さい。))
***メル周波数ケプストラム [#zf07f33d]
-上記で求めたメル周波数（スケールで圧縮された対数振幅）スペクトラム（20次元）に対して、[[離散コサイン変換(DCT):http://wapedia.mobi/ja/%E9%9B%A2%E6%95%A3%E3%82%B3%E3%82%B5%E3%82%A4%E3%83%B3%E5%A4%89%E6%8F%9B]] を行い、ケプストラム（20次元）に変換します。((詳細は [[MATLAB Note/音声の分析/フィルタバンク解析/メルフィルタバンク解析（振幅スペクトルを分割するやり方）:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#v962bcc1]] を参照して下さい。))
--コサイン変換はフーリエ変換のバリエーションのひとつといえます((参考：[[離散コサイン変換（Discrete Cosine Transform ; DCT）:http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/advanced/dct/]]))。
***メル周波数ケプストラム係数 [#ne220052]
-続いて、ケプストラムの低次成分（スペクトルの声道成分）12次元を取り出します。((詳細は [[MATLAB Note/音声の分析/フィルタバンク解析/メルフィルタバンク解析（振幅スペクトルを分割するやり方）:http://shower.human.waseda.ac.jp/~m-kouki/pukiwiki_public/73.html#v962bcc1]] を参照して下さい。))
--MFCCの次数（c[1],c[2],...）が、取り出した次元数に対応します。
-音声認識では、ケプストラムの低次成分を取り出した後、正規化処理（C0 除去，リフタリング処理，ケプストラム平均除去）によってMFCCを求めることが一般的のようです((参考文献：岩野公司, 小島要, 古井貞煕: マルチバンド音声認識のためのLDAに基づく帯域重み推定手法, 電子情報通信学会技術研究報告. SP, 音声, 106(78), pp.13-18, 2006.、[[西村義隆, 篠崎隆宏, 岩野公司, 古井貞煕, 周波数帯域ごとの重みつき尤度を用いた雑音に頑健な音声認識, 信学技報, 2003.:http://www.furui.cs.titech.ac.jp/publication/2003/sp2003-116.pdf]]))
--'''ケプストラム平均除去によって、定常な音響系などの特性に対応する成分が除去されます。'''((IT Text 音声認識システム, オーム社, p.14, 2006. より引用))
**MFCCの意義と応用 [#l5922372]
***メリット [#g7eac6ba]
-ケプストラム
--低次成分
---個人差の大きいピッチ成分を除去して、音韻の特定にとって重要な声道の音響特性（すなわち、口腔の形）のみを抽出できる
---フォルマント成分の抽出((Jonathan Darch, Ben Milner, Saeed Vaseghi, "MAP prediction of formant frequencies and voicing class from MFCC vectors in noise," Speech Communication, Vol.48, Issue.11, pp.1556-1572, 2006. によれば、直接フォルマント周波数に対応するMFCCの値はないが、GMMを使ったモデルによって高い相関を得られた))
---音声の定常性の検出 ... &color(red){MFCCの差分（デルタケプストラム）の二乗和}; が、閾値以下なら定常（≒母音）、閾値以上なら非定常（≒子音）
--高次成分
---声道の音響特性を除去して、ピッチ成分を抽出できる
-メル周波数ケプストラム係数（MFCC）
--ヒトの聴覚上重要な周波数成分が引き伸ばされて、ケプストラム全体における割合が増える
--メルフィルタバンクを通すことで、メル周波数ケプストラムの特徴量の次元数が減り、計算の負荷が減る
***デメリット [#se2d3797]
-雑音のスペクトルが特定の帯域に集中している場合、ケプストラムのすべての係数に影響を及ぼす((岩野公司, 小島要, 古井貞煕: マルチバンド音声認識のためのLDAに基づく帯域重み推定手法, 電子情報通信学会技術研究報告. SP, 音声, 106(78), pp.13-18, 2006. より引用))
-ピッチの情報は含まれない((メル周波数スペクトラムを使って、ピッチ情報を残したまま音声認識に利用する研究もある。谷口 勝則，村上 仁一，池原 悟，"モーラ情報を用いたフィルタバンクによる孤立単語認識，" 電子情報通信学会技術研究報告. NLC, 言語理解とコミュニケーション，102(527)，pp.63-68，2002． を参照。))
***動的特徴量（デルタパラメータ）の追加 [#qe7a16fb]
-母音のような定常的な音韻はメルケプストラムでよく表現できますが、子音や半母音のような動的な音韻を表現するには不十分です。
-そこで、音声認識では、メルケプストラムの時間経過による変化をあらわしたデルタメルケプストラムも使います。
--詳細は &pgid(,動的特徴量（デルタパラメータ）); を参照
***MFCC解析ツール [#h459b3a5]
-&pgid(,MFCC解析のツール); を参照
