ケプストラム解析

ケプストラムについて

ケプストラム解析は、音声処理やイメージ処理などの分野で種々に適用されている非線形信号処理手法です。

シーケンス x の "複素ケプストラム" は、x のフーリエ変換の複素自然対数を求め、その結果のシーケンスを逆フーリエ変換することにより計算されます。

$x_{}^{ˆ} = \frac{1}{2 π} \int_{- π}^{π} \log [X (e^{j ω})] e^{j ω n} d ω .$

ツールボックスの関数 cceps によりこの演算を行って、入力シーケンスに対する複素ケプストラムを推定し、入力シーケンスと同じ長さの実数シーケンスが返されます。

次に、エコー検出のアプリケーションに cceps を使用します。まず、以下のようにして 100 Hz でサンプリングした 45 Hz の正弦波を作成します。信号の立ち上りから 0.2 秒後に、半分の振幅をもつ信号のエコーを付加します。

t = 0:0.01:1.27;
s1 = sin(2*pi*45*t);
s2 = s1 + 0.5*[zeros(1,20) s1(1:108)];

新しい信号の複素ケプストラムを計算し、プロットします。

c = cceps(s2);
plot(t,c)

Figure contains an axes object. The axes object contains an object of type line.

複素ケプストラムは 0.2 秒の点にエコーを表すピークを示します。

信号 x の "実数ケプストラム" は、場合によっては単にケプストラムと呼ばれるものですが、x のフーリエ変換の大きさの自然対数を求め、そして、その結果のシーケンスの逆フーリエ変換を得ることにより計算されます。

$c_{x} = \frac{1}{2 π} \int_{- π}^{π} \log | X (e^{j ω}) | e^{j ω n} d ω .$

このツールボックスの関数 rceps は、この演算を行ってシーケンスの実数ケプストラムを返します。返されるシーケンスは、入力ベクトルと同じ長さの実数値ベクトルです。

関数 rceps では、入力と同じ実数ケプストラムをもつ固有の最小位相シーケンスが返されます。シーケンスの実数ケプストラムと最小位相再現を共に得るには、[y,ym] = rceps(x) を使用します。ここで、y は実数ケプストラム、ym は x の最小位相再現です。以下の例は、rceps の 1 出力が x と同じ実数ケプストラムをもつ固有の最小位相シーケンスであることを示します。

y = [4 1 5];                 % Non-minimum phase sequence
[xhat,yhat] = rceps(y);
xhat2 = rceps(yhat); 
[xhat' xhat2']

ans = 3×2

    1.6225    1.6225
    0.3400    0.3400
    0.3400    0.3400

逆複素ケプストラム

複素ケプストラムを逆変換するには、関数 icceps を使用します。関数 cceps では、入力のアンラップされた位相がゼロ周波数で連続になるように、データに依存した位相変調が実行されるため、逆変換が複雑になります。位相の変調は整数倍の遅延と等しくなります。出力引数を 2 つ設定すると、cceps では、この遅延の項が返されます。

x = 1:10;
[xhat,delay] = cceps(x)

xhat = 1×10

    2.2428   -0.0420   -0.0210    0.0045    0.0366    0.0788    0.1386    0.2327    0.4114    0.9249

delay = 1

複素ケプストラムを逆変換するには、以下のように元の遅延パラメーターを設定して、icceps を使用します。

icc = icceps(xhat,2)

icc = 1×10

    2.0000    3.0000    4.0000    5.0000    6.0000    7.0000    8.0000    9.0000   10.0000    1.0000

上の例が示すように、複素ケプストラムになんらかの変更を加えると、元の遅延の項が無効となる可能性があります。複素ケプストラムを厳密に逆変換できません。

参考

cceps | icceps | rceps