[vDSP][信号処理]オーディオ・音声分析への道9 ケプストラム分析01

オーディオ・音声分析への道9 ケプストラム分析01

今回は、ケプストラム分析によって、音信号のフォルマント成分と周波数成分の分離法について書きます。
今回も前回迄同様、vDSPを使用してプログラムを書きますので、iOS用のプログラムでもそのまま使用可能なものを想定します。
今回は、Miyazawa’s Pukiwiki 公開版を参考にしました。
使用するオーディオデータは、このサイトに載っている a.wav を使用します。

ケプストラム

ケプストラムとは、FFT分析によって得られたスペクトラルを更にFFTする事によって得られる"スペクトラル"の事で、1963年にB. P. Bogertらによって定義されました。

ケプストラムの名前の由来は、SpectrumのSpec部分をひっくり返して、Cepsとし、Cepstrumとなっております。スペクトラム領域とケプストラム領域では、様々な信号処理法の名前が変わります。たとば、ケプストラムをフィルタリングする場合は、FilteringのFil部分を逆から読み、Liftering : リフタリングと呼びます。
また、スペクトラムの横軸はFrequency（周波数）ですが、ケプストラムの横軸はQuefrencyとなります。

詳しくは、wikipediaを参照

フォルマント

フォルマントまたはホルマント（英: formant）とは、言葉を発している人の音声のスペクトルを観察すると分かる、時間的に移動している複数のピークのこと。

wikipedia

フォルマントとして得られたピークは、主に音声の母音を表す指標になります。つまり、類似したフォルマント成分を持つ音声信号は、同様の母音を持つことになります。

vDSPによるプログラム

手順

1. オーディオデータからFFTサイズ分のデータを取得
2. 窓関数をかける
3. FFTを行う
4. 対数スペクトラルをとる
5. iFFTを行う
6. リフタリングを行い、低Quefrencyと高Quefrencyを分離する
7. 低Quefrencyを再度FFTする事でフォルマント成分を取得
8. 高Quefrencyを再度FFTする事で周波数成分を取得

手順1~2

前のこの記事で紹介したライブラリを使用します。

初めに、必要なメモリ領域を確保します。

main.c

//for FFT analysis
    float *buffer = malloc(sizeof(float)*fft.fftsize);
    float *magnitude = malloc(sizeof(float)*fft.ffthalfsize);
    float *mag_log = malloc(sizeof(float)*fft.ffthalfsize);
    float *cepstrum = malloc(sizeof(float)*fft.fftsize);

前のこの記事では、オーディオデータからFFTサイズ分のデータを取得する手法は以下の通りにしました。

main.c

for(i=0;i<frame_num*OVERLAP -1 ;i++){
        //copy data
        for(j=0;j<FFTSIZE;j++)
            buffer[j] = data[i*FFTSIZE/OVERLAP + j];

// ここにFFT処理

}

FFTの処理は、上の最初の for loop 内に行います。

手順3 FFT

FFTは以下の通り

main.c

FFT(buffer,&fft);

この強度(amplitude)を計算します。

main.c

// sqrt(real*real + img*img)
vDSP_zvabs(&fft.splitComplex, 1, magnitude, 1, fft.ffthalfsize);

5フレーム目のフレームをプロットすると以下のようになります。

手順4 対数スペクトラルをとる

強度の対数をとります。対数を取る為には以下の関数を使用します。

main.c

// Vector convert to decibels, power, or amplitude.
extern void vDSP_vdbcon(
    const float *__vDSP_A,
    vDSP_Stride  __vDSP_IA,
    const float *__vDSP_B,
    float       *__vDSP_C,
    vDSP_Stride  __vDSP_IC,
    vDSP_Length  __vDSP_N,
    unsigned int __vDSP_F)
        __OSX_AVAILABLE_STARTING(__MAC_10_4, __IPHONE_4_0);

以下の様に使います。
最後の引数は、パワースペクトラルの場合は1を、amplitudeスペクトラルの場合は1をいれます。

main.c

float scaling_db = 1.0;
vDSP_vdbcon(magnitude,1,&scaling_db,mag_log,1,fft.ffthalfsize,0);

結果をプロットすると以下の通りです。(5フレーム目)

上の対数スペクトラルには、ただの強度を取ったスペクトラルと比べ、より小さいエネルギーを持ったFFTビンも確認できます。
ここから、細かい波と大まかな波両方が明確に現れています。フォルマント成分は、この大まかな波の成分、スペクトラルの輪郭部分を抽出する事によって大まかに得られます。
つまり、このスペクトラルを音信号として見たとき、低周波数成分と高周波数成分を分離する事ができれば、フォルマント成分と周波数成分を分離する事ができたと言えます。これにはフィルタリング（cepstrumではリフタリング）法にて達成できます。FFTによるフィルタリングはこの記事で述べました。

手順5 iFFTを行う

リフタリング第一歩として対数スペクトラルを再度FFTします。
この時、iFFTにて実行するやり方を下に示します。

対数スペクトラルをiFFTするには、以下の様にします。

main.c

    //0 padding
    vDSP_vclr(fft.splitComplex.realp,1,fft.ffthalfsize);
    vDSP_vclr(fft.splitComplex.imagp,1,fft.ffthalfsize);
        //insert mag_log to complex
        for(j=0;j<fft.ffthalfsize;j++){
            fft.splitComplex.realp[j] = mag_log[j];
            fft.splitComplex.imagp[j] = 0;
        }

splitComplex を再使用する時は、念のため0値で初期化します。
対数スペクトラルの値を Complex型 の実数部に代入し、虚数部には0を入れます。

そして、iFFTをし、結果を Complex型 から float型 に変換します。

main.c

iFFT(buffer,&fft);
vDSP_ztoc(&fft.splitComplex, 1, (COMPLEX *)cepstrum, 2, fft.ffthalfsize);

iFFTの結果(FFTサイズの半分の画像)

Cepstrumの0番目の値について

ここで、普通に全ての値をプロットしても上手くいきません。
というのも一番最初の値の扱いに注意する必要があるからです。
一番初めの値、0番目のビンには、他と比べて非常に大きな値が入っています。初めの10個の数値を書き出すと以下の様になります。

cepstrum[0] = -195920.968750
cepstrum[1] = 16780.914062
cepstrum[2] = 4551.083496
cepstrum[3] = 2393.166748
cepstrum[4] = 1619.474609
cepstrum[5] = 4564.550293
cepstrum[6] = 2171.340576
cepstrum[7] = 3664.994141
cepstrum[8] = 3249.794189
cepstrum[9] = 1745.072510

見てみると、0番目と1番目のビンでは100倍ほど値が異なります。(正負を無視して)
上の画像では、0番目の値を無視してプロットしたものを載せています。もし、0番目の値を無視しない場合、画像は以下のようになってしまいます。

さて、なぜでしょう。

さっそく調べてみます。
FFTでは前半の半FFTサイズと後半の半FFTサイズでは対称形になっています。Cepstrumの結果をFFTサイズ分プロットすると、以下のようになります。(0番目をのぞく)

では、後半の最後の10個の値を書き出すと以下の通りになります。

cepstrum[2038] = 340.088654
cepstrum[2039] = 1745.073120
cepstrum[2040] = 3249.794189
cepstrum[2041] = 3664.994385
cepstrum[2042] = 2171.340576
cepstrum[2043] = 4564.550293
cepstrum[2044] = 1619.474731
cepstrum[2045] = 2393.166504
cepstrum[2046] = 4551.083496
cepstrum[2047] = 16780.916016

おや？

最初の5つと最後の5つを比べます。

cepstrum[0] = -195920.968750 //非対称！！！
cepstrum[1] = 16780.914062 <- 2047と対称
cepstrum[2] = 4551.083496　<- 2046と対称
cepstrum[3] = 2393.166748　<- 2045と対称
cepstrum[4] = 1619.474609　<- 2044と対称

cepstrum[2043] = 4564.550293　<- 5と対称
cepstrum[2044] = 1619.474731　<- 4と対称
cepstrum[2045] = 2393.166504　<- 3と対称
cepstrum[2046] = 4551.083496　<- 2と対称
cepstrum[2047] = 16780.916016　<- 1と対称

こうすると、0番目のビンだけが対称型になっていないことが分かります。では、0番目は何を表しているのでしょうか。

これは、対数スペクトラル全体の強度を表しています。対数スペクトラルを音信号としてとらえるならば、対数スペクトラルには巨大なDCオフセットがのっていると言えます。

DCオフセットとは？

対数スペクタルをプロットした画像を見ると、スペクトラルは一つの大きな山になっています。
そこで、対数スペクトルを正規化して確認してみます。つまり、0を中心とした信号に変換します。

方法は、対数スペクトラルの値の平均値をとり、その値を対数スペクトラルから除算します。

main.c

float average = 0;
for(j=0;j<fft.ffthalfsize;j++){
    average +=mag_log[j];
}

average = average/fft.ffthalfsize;    

for(j=0;j<fft.ffthalfsize;j++){
    mag_log[j] -= average;
}

そして、再度Cepstrumを求め、結果を数値で出します。

cepstrum[0] = 4.803051
cepstrum[1] = 16876.628906
cepstrum[2] = 4455.369629
cepstrum[3] = 2488.880371
cepstrum[4] = 1523.761230

0番目の数値が小さくなっている事がわかります。

プロットする際に、0番目を除去しなくても問題なくできます。(FFTサイズの半分)

赤い丸で囲まれたピークはf0(基本周波数)を示します。
基本周波数の求め方は前回の記事、自己相関関数を使うやり方が良いそうです。

ただし、0番目のビンは移行使わない事がほとんどです。
わざわざ正規化する必要も無く、そのまま無視して行きます。

まとめ

今回は手順5まで、Cepstrumを求める所までを書きました。
長くなってきたので、続きは次回にしたいと思います。

Cepstrumについては、Miyazawa’s Pukiwiki 公開版を参考にさせていただきました。