メル周波数ケプストラム（MFCC）

メル周波数ケプストラム係数（ Mel Frequency Cepstral Coefficient ） †

はじめに †

スペクトラム（spectrum）とは †

• 音声や地震波などの周期性のある信号は、どれだけ複雑な信号であっても、単純な波に分解できる（フーリエの定理）
  • 単純な波...単一の周波数と振幅をもつ正弦波、余弦波
• 上記の定理に従って、ある信号の周波数成分と振幅成分を抽出したものが周波数スペクトラム（スペクトル）です*1。

• スペクトラムの求め方（MATLABによる説明）
  • 以下のような複雑な波を考えます。
    
    • 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/音声の加工 / 純音を作る を参照

• フーリエ変換によって、複合音を構成する各純音の周波数、振幅（パワー）の値を調べます。
  • フーリエ変換の詳細な説明（限定公開）
  • MATLABでフーリエ変換を行うには、以下のコードを実行*2

    %離散フーリエ変換
    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の結果は複素数で与えられます。実数部分と虚数部分の絶対値を振幅スペクトル、絶対値の二乗をパワースペクトル、偏角を位相スペクトルとよびます*3。

• 結果
  
  • それぞれの波の周波数と振幅の値が取得できています。*4
  • フーリエ変換は、複雑な波形の周期成分と振幅成分を取り出す手法であるといえます*5。
  • 音声のサンプリング周波数が変われば、縦軸の値も変わります。
  • 音声を短時間ごとに区切り、各フレームのスペクトラムを時間軸に沿って並べたものがスペクトログラムです。

• 【補記】縦軸をデシベル表記に変えてみます。

  %プロット（デシベル）
  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]);

  
  • 振幅スペクトルは20×log、パワースペクトルは10×log*6がデシベルの値になります。*7
  • ※音圧の単位について Okumura's Blog 波形グラフの縦軸？

逆フーリエ変換とは †

• フーリエ変換の結果から元の波形を再現することを、逆フーリエ変換とよびます。
  • MATLABでdftの値から波を生成するには、上に続けて以下のコードを実行

      newdata = ifft(dft);
      %元の音声データと同じデータ数（442）でプロットする
      plot(time*1000, newdata(1:length(wavdata))); xlabel('時間[ms]'); ylabel('振幅');

• 結果
  
  • 元の波形を再現できています。
• 複雑な波にフーリエ変換を行うことで、時間・空間軸から周波数・空間軸に変換することができ、また逆フーリエ変換を行うことで元の時間・空間軸に再変換できます。

音声波形の特性 †

音声のスペクトラム †

• まず、音声波形のスペクトラムを求めてみます。
  • a.wav を使って以下のコードを実行

      [wav, Fs] = wavread('a.wav');
      t = 0 : 1/Fs : length(wav)/Fs;
      plot(t(1:length(t)-1)*1000, wav); xlabel('時間[ms]'); ylabel('振幅');

• 結果
  

• 母音の定常部（中心部分）40ms（0.04秒）を切り出し、ハニング窓をかけて丸めたもの（フレーム）を分析対象にします。
  • 上に続けて以下のコードを実行

      %中央部分を切り出し
      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);

• 結果
  

• 切り出した音声のスペクトラムを求めます*8。

    %離散フーリエ変換
    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]);

  • フーリエ変換の次数が大きいと解析結果が細かくなります*9。
  • 振幅スペクトルはフーリエ変換の結果（複素数で表現されます*10）の絶対値です。MATLABでは「Adft = abs(dft);」と書いても求められます。その場合、パワースペクトルは「Pdft = abs(dft) .^ 2;」（振幅スペクトルの二乗）と書けます。
• 結果
  
  • 母音/a/に含まれる音声の周波数成分と振幅成分が分かりました。
  • 横軸のスケール合ってる? フォルマント解析 の結果も要確認!

• なお、振幅スペクトルを求める過程で、信号の位相の情報が除去されます。解析をするだけなら問題ないですが、音声の復号を行う場合は位相の情報が必要です。そのため、ここで位相スペクトルも計算しておきます。

    Angdft = atan2(imag(dft), real(dft));

  • MATLABでは「Angdft = angle(dft)」と書いても求められます。

音声の対数スペクトラム †

• 音声信号は、音源特性（基本周波数）と声道特性（共鳴周波数）の畳み込みである、と考えることができます。
  
  • そこで、音声信号の振幅スペクトルをS(ω)、音源信号の振幅スペクトルをG(ω)、声道の伝達特性をH(ω)とすると、以下のように表現できます*11。
    
    • ωは周波数をあらわします。スペクトルは複素数で表現され、絶対値が周波数の大きさ、偏角が波の位相をあらわしています*12
    • τの値は...*13

• 対数の基本的な演算より、両辺の対数をとると、
  
  • 複数のスペクトルの畳み込み状態を、対数スペクトルの足し算で表すことができました*14。

• 対数スペクトルをプロットしてみます。
  • 上に続けて以下のコードを実行

      %対数振幅スペクトル
      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]);

  • 結果
    
    • 負値の意味は?
    • 対数にかけることで、音源に由来する成分と声道に由来する成分が「畳み込まれた」波が、「足し合わされた」波に変換されました。
    • フーリエ変換の結果をよく見ると、大局的な変動（声道の伝達特性）に加えて、規則正しく周波数成分のピークが繰り返され（音源信号）、周期的なパターンを示していることが分かります。すなわち、フーリエ変換の結果も周期性のある信号であるといえます*15。

ケプストラム（cepstrum） †

ケプストラム分析（alanysis） †

• 対数スペクトルをもう一度フーリエ変換にかけて、音源の周波数と声道の周波数を切り分けることを考えます。
  • 振幅スペクトル（音声信号→フーリエ変換→絶対値）に対数をかけて、さらにもう一度フーリエ変換をかけたもの（スペクトラムのスペクトラム）を、ケプストラムとよびます。*16
  • フーリエ変換のフーリエ変換は、フーリエ変換の結果に対して逆フーリエ変換をかけるのと同じ処理になります。
    • ただし、絶対値をとっているので、（実数ケプストラムでは）位相スペクトルの情報は消えます。
  • MATLABで（対数）ケプストラムを生成するには、上に続けて以下のコードを実行します。*17

      %対数振幅スペクトルに逆フーリエ変換を実行して実部のみを取り出す
      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）は周波数成分の逆数（時間の尺度）です。
  • ケフレンシー上のピークは、対数スペクトル上で周期的に出現する周波数成分（倍音）をあらわしています。

• 低ケフレンシー成分を拡大すると、以下のようになります。
  
  • 対数スペクトルのケプストラムを求めたことで、高周期で変動する成分（音源成分）と、畳み込まれていた声道特性の成分が分離できています...*18
  • 基本周期はケフレンシー軸上で5.4[ms]くらいのところにあらわれているので、逆数をとって（1/0.0054）、約185[Hz]が基本周波数（f0）の値です（要説明!）*19。
  • 一番左の値（ケプストラムの0次の係数、C0）は、直流成分をあらわします（フーリエ変換 a_0の抽出）。
  • 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） †

• 上で求めた対数ケプストラムから声道特性の部分（スペクトル包絡）を抽出してみます。
  • ケプストラム係数の低次成分だけを残し、高次成分を0にした上でFFTによってスペクトラムに再変換します。

      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]);

    • 2011/07/20 パワースペクトルの計算式を修正しました。ご指摘ありがとうございました！*20
    • 0次成分+1～120次成分を残しました。
    • このとき、対称形になっている同じ値の部分も残します*21。
  • 結果
    
    • 上がリフタリングなし、下がリフタリング後
    • 基本周波数の成分が除去されて、スペクトル包絡が抽出できています*22。

音声再合成 †

• リフタリングによって音源信号成分を除去した振幅スペクトルから音声波形を生成してみます。
  • リフタリング後の振幅スペクトルと、FFTの分析 の時に求めた位相スペクトルを合わせて、sin波成分とcos波成分に戻し*23、逆フーリエ変換をして実部を取り出したものが音声波形です。

      % 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);

    • 2012/12/18 上記スクリプトの2行目は「sound_fft = AdftSpc .* exp(Angdft*j);」が正しいです（要虚数変換）。修正しました。
  • 結果
    
    • 波形の形は違いますが、声道成分が保存されているため、再合成音声も/a/に聞こえます。

• なお、ここでは音源信号成分を除去しましたが、逆に音源信号のみを残し、声道特性成分を除外することもできます。
  • リフタリングによって残す成分の次数を逆にしてみます。

      %リフタリング
      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);

  • 結果
    

• なお、音声の音源信号成分のみを取り出したいのであれば、ここで紹介したケプストラムを使った方法よりもローパスフィルタを使うほうが簡単です。

• フレームが複数ある連続音声の再合成を行う場合は、各フレームごとに音声を再合成して足しあわせる必要があります。
  • MATLAB Note/音声の加工/フレーム同士の結合 を参照して下さい。
  • リフタリングなどの処理によって振幅スペクトルの形状が変わっている場合、位相スペクトルとの不一致によって再合成音声にノイズが混ざります。これを解決するためにLSEE-MSTFTMアルゴリズムが提案されています。MATLAB Note/音声の加工/LSEE-MSTFTMアルゴリズム を参照して下さい。
    • これらの手法によって 声質を維持したピッチシフト などができます。

メル周波数ケプストラム †

• 対数ケプストラムの低次成分は、音声のスペクトル包絡（声道成分に由来した周波数特性）を表現できていることがわかりました。この対数ケプストラムの低次成分に対して、ヒトの周波数知覚特性を考慮した重み付けをした特徴量を、メル周波数ケプストラム係数（MFCC）と呼びます。

メル尺度（メル対数、メル周波数） †

• 聴覚尺度 / メル尺度 を参照してください。
• 試しに、音声の対数スペクトラム に戻り、対数スペクトラムの横軸（周波数軸）をメル周波数に変換してみます。*24

  %対数振幅スペクトル
  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)])

  • 結果
    
    • メル対数スペクトラム（下段）では、対数スペクトラム（上段）の低周波数成分が引き伸ばされて、高周波数成分が圧縮されていることがわかります。

メル周波数スペクトラム †

• 振幅スペクトルをメル尺度上で等間隔なフィルタバンクにかけて、各帯域のスペクトル成分を取り出します。*25
  • 音声認識では、経験的に20個の帯域のフィルタバンクを使います。
• （メル周波数軸上で等間隔な）各帯域の振幅スペクトルの和をとり、20次元の振幅スペクトルに圧縮します。
• 20次元に圧縮された振幅スペクトルの対数をとり、対数振幅スペクトルにします。*26

メル周波数ケプストラム †

• 上記で求めたメル周波数（スケールで圧縮された対数振幅）スペクトラム（20次元）に対して、離散コサイン変換(DCT) を行い、ケプストラム（20次元）に変換します。*27
  • コサイン変換はフーリエ変換のバリエーションのひとつといえます*28。

メル周波数ケプストラム係数 †

• 続いて、ケプストラムの低次成分（スペクトルの声道成分）12次元を取り出します。*29
  • MFCCの次数（c[1],c[2],...）が、取り出した次元数に対応します。

• 音声認識では、ケプストラムの低次成分を取り出した後、正規化処理（C0 除去，リフタリング処理，ケプストラム平均除去）によってMFCCを求めることが一般的のようです*30
  • ケプストラム平均除去によって、定常な音響系などの特性に対応する成分が除去されます。*31

MFCCの意義と応用 †

メリット †

• ケプストラム
  • 低次成分
    • 個人差の大きいピッチ成分を除去して、音韻の特定にとって重要な声道の音響特性（すなわち、口腔の形）のみを抽出できる
    • フォルマント成分の抽出*32
    • 音声の定常性の検出 ... MFCCの差分（デルタケプストラム）の二乗和 が、閾値以下なら定常（≒母音）、閾値以上なら非定常（≒子音）
  • 高次成分
    • 声道の音響特性を除去して、ピッチ成分を抽出できる
• メル周波数ケプストラム係数（MFCC）
  • ヒトの聴覚上重要な周波数成分が引き伸ばされて、ケプストラム全体における割合が増える
  • メルフィルタバンクを通すことで、メル周波数ケプストラムの特徴量の次元数が減り、計算の負荷が減る

デメリット †

• 雑音のスペクトルが特定の帯域に集中している場合、ケプストラムのすべての係数に影響を及ぼす*33
• ピッチの情報は含まれない*34

動的特徴量（デルタパラメータ）の追加 †

• 母音のような定常的な音韻はメルケプストラムでよく表現できますが、子音や半母音のような動的な音韻を表現するには不十分です。
• そこで、音声認識では、メルケプストラムの時間経過による変化をあらわしたデルタメルケプストラムも使います。
  • 詳細は 動的特徴量（デルタパラメータ） を参照

MFCC解析ツール †

• MFCC解析のツール を参照

検索

更新履歴

最新の20件

edit
限定版