Miyazawa’s Pukiwiki
メル周波数ケプストラム（MFCC）
はすでに存在します。
開始行:
#access
#analog
*メル周波数ケプストラム係数（ Mel Frequency Cepstral Coef...
#contents
**はじめに [#i52f1c12]
***スペクトラム（spectrum）とは [#u28deb5a]
-音声や地震波などの周期性のある信号は、どれだけ複雑な信号...
--単純な波...単一の周波数と振幅をもつ正弦波、余弦波
-上記の定理に従って、ある信号の周波数成分と振幅成分を抽出...
-スペクトラムの求め方（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、周波...
coswav_1 = 0.9 * cos(2 * pi * 200 * time); %振幅0.9、周波...
sinwav_2 = 1.8 * sin(2 * pi * 260 * time); %振幅1.8、周波...
coswav_2 = 1.4 * cos(2 * pi * 320 * time); %振幅1.4、周波...
wavdata = 1.4 + [sinwav_1 + coswav_1 + sinwav_2 + coswav_...
plot(time*1000, wavdata); xlabel('時間[ms]'); ylabel('振...
}}
---正弦波、余弦波の詳細は [[MATLAB Note/音声の加工 / 純音...
--フーリエ変換によって、複合音を構成する各純音の周波数、...
---'''[[フーリエ変換の詳細な説明（限定公開）:http://showe...
---MATLABでフーリエ変換を行うには、以下のコードを実行((参...
#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までをf...
%プロット（表示するのはサンプリング周波数の半分まで）
subplot(2,1,1); plot(fscale(1:fftsize/2), Adft(1:fftsize/...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル'); xlim([0 1...
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft(1:fftsize/...
xlabel('周波数[Hz]'); ylabel('パワースペクトル'); xlim([0...
}}
---fftの結果は複素数で与えられます。実数部分と虚数部分の...
--結果
#ref(spectrum_2.png);
---それぞれの波の周波数と振幅の値が取得できています。((左...
---フーリエ変換は、''複雑な波形の周期成分と振幅成分を取り...
---&color(red){音声のサンプリング周波数が変われば、縦軸の...
---[[音声を短時間ごとに区切り:http://shower.human.waseda....
--【補記】縦軸をデシベル表記に変えてみます。
#geshi(matlab){{
%プロット（デシベル）
subplot(2,1,1); plot(fscale(1:fftsize/2), 20*log10(Adft(1...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル[dB]'); xlim(...
subplot(2,1,2); plot(fscale(1:fftsize/2), 10*log10(Pdft(1...
xlabel('周波数[Hz]'); ylabel('パワースペクトル[dB]'); xli...
}}
#ref(spectrum_2dB.png,,80%);
---振幅スペクトルは&color(red){20};×log、パワースペクトル...
---※音圧の単位について [[Okumura's Blog 波形グラフの縦軸...
***逆フーリエ変換とは [#l34aa4d5]
-フーリエ変換の結果から元の波形を再現することを、逆フーリ...
--MATLABでdftの値から波を生成するには、上に続けて以下のコ...
#geshi(matlab){{
newdata = ifft(dft);
%元の音声データと同じデータ数（442）でプロットする
plot(time*1000, newdata(1:length(wavdata))); xlabel('時...
}}
--結果
#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]'); y...
}}
--結果
#ref(cepstrum_1.jpg);
-母音の定常部（中心部分）40ms（0.04秒）を切り出し、[[ハニ...
--上に続けて以下のコードを実行
#geshi(matlab){{
%中央部分を切り出し
center = fix(length(wav) / 2); %中心のサンプル番号
cuttime = 0.04; %切り出す長さ[s]
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(cut...
time = t(center-fix(cuttime/2*Fs) : center+fix(cuttime/...
subplot(2,1,1); plot(time*1000, wavdata); ylabel('振幅');
%ハニング窓をかける
han_window = 0.5 - 0.5 * cos(2 * pi * [0 : 1/length(wav...
wavdata = han_window(1:length(wavdata))' .* wavdata;
subplot(2,1,2); plot(time*1000, wavdata); xlabel('時間[...
% 切り出したフレームの再生
wavplay(wavdata, Fs);
}}
--結果
#ref(cepstrum_2.jpg);
--切り出した音声のスペクトラムを求めます((スペクトラムを...
#geshi(matlab){{
%離散フーリエ変換
fftsize = 2048; %フーリエ変換の次数, 周波数ポイントの数
dft = fft(wavdata, fftsize); % 実部がcos波成分、虚部がs...
%パワースペクトル（実部と虚部をそれぞれ二乗して合計）
Pdft = (real(dft) .^ 2) + (imag(dft) .^ 2);
%振幅スペクトル（パワースペクトルの平方根）
Adft = sqrt(Pdft);
%周波数スケール（サンプリング定理より、OHzからサンプリ...
fscale = linspace(0, Fs, fftsize); %0～Fsまでをfftsi...
%プロット（表示するのはサンプリング周波数の半分まで）
subplot(2,1,1); plot(fscale(1:fftsize/2), Adft(1:fftsiz...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル'); xlim([0...
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft(1:fftsiz...
xlabel('周波数[Hz]'); ylabel('パワースペクトル'); xlim(...
}}
---フーリエ変換の次数が大きいと解析結果が細かくなります((...
---振幅スペクトルはフーリエ変換の結果（複素数で表現されま...
--結果
#ref(cepstrum_3.png);
---母音/a/に含まれる音声の周波数成分と振幅成分が分かりま...
---&color(red){横軸のスケール合ってる? [[フォルマント解析...
--なお、振幅スペクトルを求める過程で、信号の位相の情報が...
#geshi(matlab){{
Angdft = atan2(imag(dft), real(dft));
}}
---MATLABでは「Angdft = angle(dft)」と書いても求められま...
***音声の対数スペクトラム [#d33d7f39]
-音声信号は、音源特性（基本周波数）と声道特性（共鳴周波数...
#ref(cepstrum_7.jpg);
--そこで、音声信号の振幅スペクトルをS(ω)、音源信号の振幅...
#mimetex(|S(\omega)| = |G(\omega)| |H(\omega)| = \int_{\i...
---ωは周波数をあらわします。スペクトルは複素数で表現され...
---τの値は...((畳み込みとは：[[初心者用 畳み込み（たたみ...
--[[対数の基本的な演算:http://ja.wikipedia.org/wiki/%E5%A...
#mimetex(log|S(x)| = log|G(x)| + log|H(x)|);
---複数のスペクトルの畳み込み状態を、対数スペクトルの足し...
-対数スペクトルをプロットしてみます。
--上に続けて以下のコードを実行
#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:ff...
xlabel('周波数[Hz]'); ylabel('対数振幅スペクトル'); xli...
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft_log(1:ff...
xlabel('周波数[Hz]'); ylabel('対数パワースペクトル'); x...
}}
--結果
#ref(cepstrum_4.png);
---&color(red){負値の意味は?};
---対数にかけることで、音源に由来する成分と声道に由来する...
---フーリエ変換の結果をよく見ると、大局的な変動（声道の伝...
**ケプストラム（cepstrum） [#qee1e455]
***ケプストラム分析（alanysis） [#x3efddeb]
-対数スペクトルをもう一度フーリエ変換にかけて、音源の周波...
--振幅スペクトル（音声信号→フーリエ変換→絶対値）に対数を...
--フーリエ変換のフーリエ変換は、フーリエ変換の結果に対し...
---ただし、絶対値をとっているので、（実数ケプストラムでは...
--MATLABで（対数）ケプストラムを生成するには、上に続けて...
#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]...
--結果
#ref(cepstrum_5.png);
---横軸の[[ケフレンシー（quefrency）:http://ja.wikipedia....
---ケフレンシー上のピークは、対数スペクトル上で周期的に出...
--低ケフレンシー成分を拡大すると、以下のようになります。
#ref(cepstrum_6.jpg);
---対数スペクトルのケプストラムを求めたことで、高周期で変...
---基本周期はケフレンシー軸上で5.4[ms]くらいのところにあ...
---一番左の値（ケプストラムの0次の係数、C0）は、直流成分...
---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にした...
#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:fft...
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:fft...
xlabel('周波数[Hz]'); xlim([0,5000]);
}}
---&color(red){2011/07/20 パワースペクトルの計算式を修正...
---0次成分+1～120次成分を残しました。
---このとき、対称形になっている同じ値の部分も残します((参...
--結果
#ref(cepstrum_9_refer.png);
---上がリフタリングなし、下がリフタリング後
---基本周波数の成分が除去されて、スペクトル包絡が抽出でき...
//---スペクトルの絶対値（振幅スペクトル）から作った実数ケ...
***音声再合成 [#f411ed6a]
-リフタリングによって音源信号成分を除去した振幅スペクトル...
--リフタリング後の振幅スペクトルと、[[FFTの分析:http://sh...
#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_...
--結果
#ref(cepstrum_8.png);
---波形の形は違いますが、声道成分が保存されているため、再...
-なお、ここでは音源信号成分を除去しましたが、逆に音源信号...
--リフタリングによって残す成分の次数を逆にしてみます。
#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:f...
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:f...
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://spee...
--リフタリングなどの処理によって振幅スペクトルの形状が変...
---これらの手法によって [[声質を維持したピッチシフト:http...
**メル周波数ケプストラム [#o8763891]
-対数ケプストラムの低次成分は、音声のスペクトル包絡（声道...
***メル尺度（メル対数、メル周波数） [#i7714922]
-[[聴覚尺度 / メル尺度:http://shower.human.waseda.ac.jp/~...
-試しに、[[音声の対数スペクトラム:http://shower.human.was...
#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:ffts...
xlabel('周波数[mel]'); ylabel('対数振幅スペクトル'); xlim...
subplot(2,1,2); plot(mellog(fscale(1:fftsize/2)), Pdft_lo...
%横軸をメルスケールに変換してプロット
subplot(2,1,2); plot(mellog(fscale(1:fftsize/2)), Adft_lo...
xlabel('周波数[mel]'); ylabel('メル対数振幅スペクトル'); ...
}}
--結果
#ref(cepstrum_10.png);
---メル対数スペクトラム（下段）では、対数スペクトラム（上...
***メル周波数スペクトラム [#x7e138e5]
-振幅スペクトルを[[メル尺度上で等間隔なフィルタバンク:htt...
--音声認識では、経験的に20個の帯域のフィルタバンクを使い...
-（メル周波数軸上で等間隔な）各帯域の振幅スペクトルの和を...
-20次元に圧縮された振幅スペクトルの対数をとり、対数振幅ス...
***メル周波数ケプストラム [#zf07f33d]
-上記で求めたメル周波数（スケールで圧縮された対数振幅）ス...
--コサイン変換はフーリエ変換のバリエーションのひとつとい...
***メル周波数ケプストラム係数 [#ne220052]
-続いて、ケプストラムの低次成分（スペクトルの声道成分）12...
--MFCCの次数（c[1],c[2],...）が、取り出した次元数に対応し...
-音声認識では、ケプストラムの低次成分を取り出した後、正規...
--'''ケプストラム平均除去によって、定常な音響系などの特性...
**MFCCの意義と応用 [#l5922372]
***メリット [#g7eac6ba]
-ケプストラム
--低次成分
---個人差の大きいピッチ成分を除去して、音韻の特定にとって...
---フォルマント成分の抽出((Jonathan Darch, Ben Milner, Sa...
---音声の定常性の検出 ... &color(red){MFCCの差分（デルタ...
--高次成分
---声道の音響特性を除去して、ピッチ成分を抽出できる
-メル周波数ケプストラム係数（MFCC）
--ヒトの聴覚上重要な周波数成分が引き伸ばされて、ケプスト...
--メルフィルタバンクを通すことで、メル周波数ケプストラム...
***デメリット [#se2d3797]
-雑音のスペクトルが特定の帯域に集中している場合、ケプスト...
-ピッチの情報は含まれない((メル周波数スペクトラムを使って...
***動的特徴量（デルタパラメータ）の追加 [#qe7a16fb]
-母音のような定常的な音韻はメルケプストラムでよく表現でき...
-そこで、音声認識では、メルケプストラムの時間経過による変...
--詳細は &pgid(,動的特徴量（デルタパラメータ）); を参照
***MFCC解析ツール [#h459b3a5]
-&pgid(,MFCC解析のツール); を参照
終了行:
#access
#analog
*メル周波数ケプストラム係数（ Mel Frequency Cepstral Coef...
#contents
**はじめに [#i52f1c12]
***スペクトラム（spectrum）とは [#u28deb5a]
-音声や地震波などの周期性のある信号は、どれだけ複雑な信号...
--単純な波...単一の周波数と振幅をもつ正弦波、余弦波
-上記の定理に従って、ある信号の周波数成分と振幅成分を抽出...
-スペクトラムの求め方（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、周波...
coswav_1 = 0.9 * cos(2 * pi * 200 * time); %振幅0.9、周波...
sinwav_2 = 1.8 * sin(2 * pi * 260 * time); %振幅1.8、周波...
coswav_2 = 1.4 * cos(2 * pi * 320 * time); %振幅1.4、周波...
wavdata = 1.4 + [sinwav_1 + coswav_1 + sinwav_2 + coswav_...
plot(time*1000, wavdata); xlabel('時間[ms]'); ylabel('振...
}}
---正弦波、余弦波の詳細は [[MATLAB Note/音声の加工 / 純音...
--フーリエ変換によって、複合音を構成する各純音の周波数、...
---'''[[フーリエ変換の詳細な説明（限定公開）:http://showe...
---MATLABでフーリエ変換を行うには、以下のコードを実行((参...
#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までをf...
%プロット（表示するのはサンプリング周波数の半分まで）
subplot(2,1,1); plot(fscale(1:fftsize/2), Adft(1:fftsize/...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル'); xlim([0 1...
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft(1:fftsize/...
xlabel('周波数[Hz]'); ylabel('パワースペクトル'); xlim([0...
}}
---fftの結果は複素数で与えられます。実数部分と虚数部分の...
--結果
#ref(spectrum_2.png);
---それぞれの波の周波数と振幅の値が取得できています。((左...
---フーリエ変換は、''複雑な波形の周期成分と振幅成分を取り...
---&color(red){音声のサンプリング周波数が変われば、縦軸の...
---[[音声を短時間ごとに区切り:http://shower.human.waseda....
--【補記】縦軸をデシベル表記に変えてみます。
#geshi(matlab){{
%プロット（デシベル）
subplot(2,1,1); plot(fscale(1:fftsize/2), 20*log10(Adft(1...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル[dB]'); xlim(...
subplot(2,1,2); plot(fscale(1:fftsize/2), 10*log10(Pdft(1...
xlabel('周波数[Hz]'); ylabel('パワースペクトル[dB]'); xli...
}}
#ref(spectrum_2dB.png,,80%);
---振幅スペクトルは&color(red){20};×log、パワースペクトル...
---※音圧の単位について [[Okumura's Blog 波形グラフの縦軸...
***逆フーリエ変換とは [#l34aa4d5]
-フーリエ変換の結果から元の波形を再現することを、逆フーリ...
--MATLABでdftの値から波を生成するには、上に続けて以下のコ...
#geshi(matlab){{
newdata = ifft(dft);
%元の音声データと同じデータ数（442）でプロットする
plot(time*1000, newdata(1:length(wavdata))); xlabel('時...
}}
--結果
#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]'); y...
}}
--結果
#ref(cepstrum_1.jpg);
-母音の定常部（中心部分）40ms（0.04秒）を切り出し、[[ハニ...
--上に続けて以下のコードを実行
#geshi(matlab){{
%中央部分を切り出し
center = fix(length(wav) / 2); %中心のサンプル番号
cuttime = 0.04; %切り出す長さ[s]
wavdata = wav(center-fix(cuttime/2*Fs) : center+fix(cut...
time = t(center-fix(cuttime/2*Fs) : center+fix(cuttime/...
subplot(2,1,1); plot(time*1000, wavdata); ylabel('振幅');
%ハニング窓をかける
han_window = 0.5 - 0.5 * cos(2 * pi * [0 : 1/length(wav...
wavdata = han_window(1:length(wavdata))' .* wavdata;
subplot(2,1,2); plot(time*1000, wavdata); xlabel('時間[...
% 切り出したフレームの再生
wavplay(wavdata, Fs);
}}
--結果
#ref(cepstrum_2.jpg);
--切り出した音声のスペクトラムを求めます((スペクトラムを...
#geshi(matlab){{
%離散フーリエ変換
fftsize = 2048; %フーリエ変換の次数, 周波数ポイントの数
dft = fft(wavdata, fftsize); % 実部がcos波成分、虚部がs...
%パワースペクトル（実部と虚部をそれぞれ二乗して合計）
Pdft = (real(dft) .^ 2) + (imag(dft) .^ 2);
%振幅スペクトル（パワースペクトルの平方根）
Adft = sqrt(Pdft);
%周波数スケール（サンプリング定理より、OHzからサンプリ...
fscale = linspace(0, Fs, fftsize); %0～Fsまでをfftsi...
%プロット（表示するのはサンプリング周波数の半分まで）
subplot(2,1,1); plot(fscale(1:fftsize/2), Adft(1:fftsiz...
xlabel('周波数[Hz]'); ylabel('振幅スペクトル'); xlim([0...
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft(1:fftsiz...
xlabel('周波数[Hz]'); ylabel('パワースペクトル'); xlim(...
}}
---フーリエ変換の次数が大きいと解析結果が細かくなります((...
---振幅スペクトルはフーリエ変換の結果（複素数で表現されま...
--結果
#ref(cepstrum_3.png);
---母音/a/に含まれる音声の周波数成分と振幅成分が分かりま...
---&color(red){横軸のスケール合ってる? [[フォルマント解析...
--なお、振幅スペクトルを求める過程で、信号の位相の情報が...
#geshi(matlab){{
Angdft = atan2(imag(dft), real(dft));
}}
---MATLABでは「Angdft = angle(dft)」と書いても求められま...
***音声の対数スペクトラム [#d33d7f39]
-音声信号は、音源特性（基本周波数）と声道特性（共鳴周波数...
#ref(cepstrum_7.jpg);
--そこで、音声信号の振幅スペクトルをS(ω)、音源信号の振幅...
#mimetex(|S(\omega)| = |G(\omega)| |H(\omega)| = \int_{\i...
---ωは周波数をあらわします。スペクトルは複素数で表現され...
---τの値は...((畳み込みとは：[[初心者用 畳み込み（たたみ...
--[[対数の基本的な演算:http://ja.wikipedia.org/wiki/%E5%A...
#mimetex(log|S(x)| = log|G(x)| + log|H(x)|);
---複数のスペクトルの畳み込み状態を、対数スペクトルの足し...
-対数スペクトルをプロットしてみます。
--上に続けて以下のコードを実行
#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:ff...
xlabel('周波数[Hz]'); ylabel('対数振幅スペクトル'); xli...
subplot(2,1,2); plot(fscale(1:fftsize/2), Pdft_log(1:ff...
xlabel('周波数[Hz]'); ylabel('対数パワースペクトル'); x...
}}
--結果
#ref(cepstrum_4.png);
---&color(red){負値の意味は?};
---対数にかけることで、音源に由来する成分と声道に由来する...
---フーリエ変換の結果をよく見ると、大局的な変動（声道の伝...
**ケプストラム（cepstrum） [#qee1e455]
***ケプストラム分析（alanysis） [#x3efddeb]
-対数スペクトルをもう一度フーリエ変換にかけて、音源の周波...
--振幅スペクトル（音声信号→フーリエ変換→絶対値）に対数を...
--フーリエ変換のフーリエ変換は、フーリエ変換の結果に対し...
---ただし、絶対値をとっているので、（実数ケプストラムでは...
--MATLABで（対数）ケプストラムを生成するには、上に続けて...
#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]...
--結果
#ref(cepstrum_5.png);
---横軸の[[ケフレンシー（quefrency）:http://ja.wikipedia....
---ケフレンシー上のピークは、対数スペクトル上で周期的に出...
--低ケフレンシー成分を拡大すると、以下のようになります。
#ref(cepstrum_6.jpg);
---対数スペクトルのケプストラムを求めたことで、高周期で変...
---基本周期はケフレンシー軸上で5.4[ms]くらいのところにあ...
---一番左の値（ケプストラムの0次の係数、C0）は、直流成分...
---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にした...
#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:fft...
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:fft...
xlabel('周波数[Hz]'); xlim([0,5000]);
}}
---&color(red){2011/07/20 パワースペクトルの計算式を修正...
---0次成分+1～120次成分を残しました。
---このとき、対称形になっている同じ値の部分も残します((参...
--結果
#ref(cepstrum_9_refer.png);
---上がリフタリングなし、下がリフタリング後
---基本周波数の成分が除去されて、スペクトル包絡が抽出でき...
//---スペクトルの絶対値（振幅スペクトル）から作った実数ケ...
***音声再合成 [#f411ed6a]
-リフタリングによって音源信号成分を除去した振幅スペクトル...
--リフタリング後の振幅スペクトルと、[[FFTの分析:http://sh...
#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_...
--結果
#ref(cepstrum_8.png);
---波形の形は違いますが、声道成分が保存されているため、再...
-なお、ここでは音源信号成分を除去しましたが、逆に音源信号...
--リフタリングによって残す成分の次数を逆にしてみます。
#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:f...
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:f...
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://spee...
--リフタリングなどの処理によって振幅スペクトルの形状が変...
---これらの手法によって [[声質を維持したピッチシフト:http...
**メル周波数ケプストラム [#o8763891]
-対数ケプストラムの低次成分は、音声のスペクトル包絡（声道...
***メル尺度（メル対数、メル周波数） [#i7714922]
-[[聴覚尺度 / メル尺度:http://shower.human.waseda.ac.jp/~...
-試しに、[[音声の対数スペクトラム:http://shower.human.was...
#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:ffts...
xlabel('周波数[mel]'); ylabel('対数振幅スペクトル'); xlim...
subplot(2,1,2); plot(mellog(fscale(1:fftsize/2)), Pdft_lo...
%横軸をメルスケールに変換してプロット
subplot(2,1,2); plot(mellog(fscale(1:fftsize/2)), Adft_lo...
xlabel('周波数[mel]'); ylabel('メル対数振幅スペクトル'); ...
}}
--結果
#ref(cepstrum_10.png);
---メル対数スペクトラム（下段）では、対数スペクトラム（上...
***メル周波数スペクトラム [#x7e138e5]
-振幅スペクトルを[[メル尺度上で等間隔なフィルタバンク:htt...
--音声認識では、経験的に20個の帯域のフィルタバンクを使い...
-（メル周波数軸上で等間隔な）各帯域の振幅スペクトルの和を...
-20次元に圧縮された振幅スペクトルの対数をとり、対数振幅ス...
***メル周波数ケプストラム [#zf07f33d]
-上記で求めたメル周波数（スケールで圧縮された対数振幅）ス...
--コサイン変換はフーリエ変換のバリエーションのひとつとい...
***メル周波数ケプストラム係数 [#ne220052]
-続いて、ケプストラムの低次成分（スペクトルの声道成分）12...
--MFCCの次数（c[1],c[2],...）が、取り出した次元数に対応し...
-音声認識では、ケプストラムの低次成分を取り出した後、正規...
--'''ケプストラム平均除去によって、定常な音響系などの特性...
**MFCCの意義と応用 [#l5922372]
***メリット [#g7eac6ba]
-ケプストラム
--低次成分
---個人差の大きいピッチ成分を除去して、音韻の特定にとって...
---フォルマント成分の抽出((Jonathan Darch, Ben Milner, Sa...
---音声の定常性の検出 ... &color(red){MFCCの差分（デルタ...
--高次成分
---声道の音響特性を除去して、ピッチ成分を抽出できる
-メル周波数ケプストラム係数（MFCC）
--ヒトの聴覚上重要な周波数成分が引き伸ばされて、ケプスト...
--メルフィルタバンクを通すことで、メル周波数ケプストラム...
***デメリット [#se2d3797]
-雑音のスペクトルが特定の帯域に集中している場合、ケプスト...
-ピッチの情報は含まれない((メル周波数スペクトラムを使って...
***動的特徴量（デルタパラメータ）の追加 [#qe7a16fb]
-母音のような定常的な音韻はメルケプストラムでよく表現でき...
-そこで、音声認識では、メルケプストラムの時間経過による変...
--詳細は &pgid(,動的特徴量（デルタパラメータ）); を参照
***MFCC解析ツール [#h459b3a5]
-&pgid(,MFCC解析のツール); を参照
ページ名:
既存のページ名で編集する
