倒頻譜

倒頻譜的範例

倒頻譜（cepstrum），顧名思義，就是將頻譜（spectrum）的英文前四個字母反過來寫。倒頻譜是為了某些時候，為了計算方便，將原來信號的頻譜先轉成類似分貝的單位，再作逆傅里叶变换，把它視為一種新的訊號做處理。倒頻譜有複數倒頻譜，及實數倒頻譜。

倒頻譜被定義在1963的論文（Bogert等）。定義如下：

• 字義：倒頻譜（信號）是信號頻譜取對數的傅立葉變換後的新頻譜（信號），有時候會稱頻譜的倒頻譜。
• 數學上：信號的倒頻譜 = IFT ( log ( | FT (信号) | ) + j2πm )（m為實數）
• 演算法：信号 -> 傅立叶变换 -> 取绝对值 -> 取对数 -> 相位展开 -> 逆傅立叶变换 -> 倒频谱

複數倒頻譜擁有頻譜大小跟相位的資訊，實數倒頻譜只有頻譜大小的資訊，各有各的不同應用。

複數倒頻譜與實數倒頻譜

複數倒頻譜

$\hat{x}\left[n\right]=\underset{-\frac{1}{2}}{\overset{\frac{1}{2}}{\int }}\hat{X}\left(F\right)e^{j2\pi F}dF$
其中$\hat{X}\left[F\right]=\log |X(F)|+j\arg [X(F)]$
可能遭遇的問題
1. $\log 0=-\mathrm{\infty }$
2. $\arg [X[n]]$有無限多的解
當輸入是實數時,因為$\log |X(F)|$偶對稱，$\arg [X(F)]$奇對稱,所以複數倒頻譜的值為實數

實數倒頻譜

$C\left[n\right]=\underset{-\frac{1}{2}}{\overset{\frac{1}{2}}{\int }}\log |X(F)|e^{j2\pi Fn}dF$
可能遭遇的問題
1. $\log 0=-\mathrm{\infty }$

應用

• 倒頻譜可以被視為在不同頻帶上變化速率的資訊，倒頻譜一開始被發明在地震或炸彈產生的地震回音，現今也被使用在分析雷達訊號，以及訊號處理等問題。
• 自相關倒頻譜(autocepstrum)被定義為倒頻譜的自相關性，自相關倒頻譜有時在分析處理回傳訊號時比倒頻譜還準確。
• 倒頻譜在處理人聲訊號以及音樂訊號有非常好的效果，例如梅爾頻率倒頻譜(Mel-Frequency Cepstrum)，用來做聲音的辨認，偵測音高等。近年來梅耳倒頻譜也被應用在音樂資訊的回覆。
• 倒頻譜在聲學中可以將聲帶震動的影響去除。
• 倒頻譜用在處理多路徑問題時(如聲波的迴音、電磁波的折、反射等)，如果將其他路徑干擾視為雜訊，為了消除雜訊，利用倒頻譜，不需測量每條多路徑的延遲時間，可以利用傳送多次信號，觀察其他路徑在倒頻譜上的效果，並且加以濾除。
• 語音大致上是由音高、聲帶脈衝、聲門波形所組成，我們可以利用倒頻譜將這三種元素在倒頻域上分開，以利於做語音訊號的分析。
• 倒頻譜的微分適用於影像處理上的圖形辨認(pattern recognition)。
• 倒頻譜與同型聲音理論(homomorphic sound theory)有關。

倒頻譜觀念

頻譜圖上的獨立變數是頻率，而倒頻譜圖上的獨立變數為倒頻率(quefrency)，倒頻率是一種時間的度量單位。舉個例子，聲音訊號取樣速率等於44100赫茲，在倒頻譜上有個很大的值在倒頻率等於100，代表實際上在44100/100=441赫茲有很大的值，這值出現在倒頻譜上因為頻譜上週期性出現，而頻譜上出現的週期與倒頻譜很大的值出現的位置有關。

倒濾波器

濾波器(filter)常使用在頻譜上，用來保存或刪除我們所要或不要的資訊，經過上面的許多討論，不難猜到，倒濾波器(lifter)就是在倒頻譜上所使用的濾波器。低通的倒濾波器跟低通濾波器有點類似，它可以藉由在倒頻譜上乘以一個window係數，使倒頻譜上的高倒頻率被壓抑，如此依來，當信號轉回時域空間時會變成一個較平滑的信號。

計算倒頻譜的方法

直接計算IDTFT(反離散時間傅立葉變換)

$\hat{x}\left[n\right]=\underset{-\frac{1}{2}}{\overset{\frac{1}{2}}{\int }}\hat{X}\left(F\right)e^{j2\pi F}dF$
問題: $\hat{X}\left(F\right)$ 可能會無限大, 且對於arg(x[n])有無限多個解

利用Z轉換的零點與極點

先對信號做Z轉換, 並整理一下係數, 讓他變成下面的形式
$X\left(Z\right)=\frac{A{Z}^r\prod _{k=1}^{m_i}(1-a_k{Z}^{-1})\prod _{k=1}^{m_0}(1-b_{k}Z)}{\prod _{k=1}^{P_i}(1-c_k{Z}^{-1})\prod _{k=1}^{P_0}(1-d_{k}Z)}$
其中$\left|a_k\right|,\left|b_k\right|,\left|c_k\right|,\left|d_k\right|\le 1$

分子:
第一項A是係數
第二項$Z^r$是延遲
第三項是位於單位圓內的零點
第四項是位於單位圓外的零點

分母:
第一項是位於單位圓內的極點
第二項是位於單位圓外的極點

對$X\left(Z\right)$取log變成$\hat{X}\left(Z\right)$
$\hat{X}\left(Z\right)=logX\left(Z\right)=\log A+r\log Z+\sum _{k=1}^{m_i}\log (1-a_k{Z}^{-1})+\sum _{k=1}^{m_0}\log (1-b_{k}Z)-\sum _{k=1}^{P_i}\log (1-c_k{Z}^{-1})-\sum _{k=1}^{P_0}\log (1-d_{k}Z)$
假設r=0, 因為這只是延遲, 並不會破壞波形
根據Z轉換所得到的系數, 我們可以利用泰勒展開得到Z的反轉換
$\hat{x}\left[n\right]=\{\begin{array}{cc}\log A& \text{if }n=0\\ -\sum _{k=1}^{m_i}\frac{{a_k}^n}{n}+\sum _{k=1}^{P_i}\frac{{c_k}^n}{n}& \text{if }n>0\\ \sum _{k=1}^{m_0}\frac{{b_k}^{-n}}{n}-\sum _{k=1}^{P_0}\frac{{d_k}^{-n}}{n}& \text{if }n<0\end{array}$

注意事項
1.$\hat{x}\left[n\right]$總是IIR(無限脈衝響應)
2.對於FIR(有限脈衝響應)的情況, $c_k=0,d_k=0$

利用Z轉換與微分

$Z\cdot {\hat{X}}^{\prime }\left(Z\right)=Z\cdot \frac{X^{\prime }\left(Z\right)}{X\left(Z\right)}$
$Z{X}^{\prime }\left(Z\right)=Z{\hat{X}}^{\prime }\left(Z\right)\cdot X\left(Z\right)$
對其做Z的反轉換
$nx[n]=\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}k\hat{x}\left[k\right]x[n-k]$
故
$x[n]=\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}\frac{k}{n}\hat{x}\left[k\right]x[n-k]forn\ne 0$

分別對於x[n]的四種不同的狀況做延伸
1.對於x[n]是因果(causal)和最小相位(minimum phase) i.e. $x[n]=\hat{x}\left[n\right]=0,n<0$
對於$x[n]=\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}\frac{k}{n}\hat{x}\left[k\right]x[n-k]forn\ne 0$
可得出
$x[n]=\sum _{k=0}^{\mathrm{\infty }}\frac{k}{n}\hat{x}\left[k\right]x[n-k]forn>0$
故
$x[n]=\hat{x}\left[n\right]x[0]+\sum _{k=0}^{n-1}\frac{k}{n}\hat{x}\left[k\right]x[n-k]$
2.對於x[n]是最小相位(minimum phase)
$\hat{x}\left[n\right]=\{\begin{array}{cc}0& \text{if }n<0\\ \frac{x[n]}{x[0]}-\sum _{k=0}^{n-1}\frac{k}{n}\hat{x}\left[k\right]\frac{x[n-k]}{x[0]}& \text{if }n>0\\ \log A& \text{if }n=0\end{array}$
3.對於x[n]是反因果(anti-causal)且最大相位(maximum phase) i.e. $x[n]=\hat{x}\left[n\right]=0,n>0$
$\begin{array}{cc}x[n]& =\sum _{k=n}^0\frac{k}{n}\hat{x}\left[k\right]x[n-k]forn<0\\ & =\hat{x}\left[n\right]x[0]+\sum _{k=n+1}^0\frac{k}{n}\hat{x}\left[k\right]x[n-k]\\ \end{array}$
4.對於x[n]是最大相位(maximum phase)
$\hat{x}\left[n\right]=\{\begin{array}{cc}0& \text{if }n>0\\ \frac{x[n]}{x[0]}-\sum _{k=n+1}^0\frac{k}{n}\hat{x}\left[k\right]\frac{x[n-k]}{x[0]}& \text{if }n<0\\ \log A& \text{if }n=0\end{array}$

特性

1. 複數倒頻譜至少以$\frac{1}{n}$的速度衰退
$|\hat{x}\left[n\right]|=c|\frac{{\alpha }^n}{n}|-\mathrm{\infty }<n<\mathrm{\infty }$
其中 $\alpha =max(a_k,b_k,c_k,d_k)$
2. 如果X(Z)沒有在單位圓以外的零點和極點, 則
$\hat{x}\left[n\right]=0foralln<0$
因為$b_k,d_k=0$
3. 如果X(Z)沒有在單位圓以內的零點和極點, 則
$\hat{x}\left[n\right]=0foralln>0$
因為$a_k,c_k=0$
4. 如果x[n]是有限長度, 則$\hat{x}\left[n\right]$是無限長度

同態解卷積的應用(Application of Homomorphic Deconvolution)

同態解卷積有非常多應用面，尤其是在聲學工程和語音分析方面的實用性

(1) 回聲的均衡化

$y[n]=x[n]+\alpha x[n-N_p]$

其中$y[n]$是接收到的訊號，$x[n]$是原始訊號，$N_p$是延遲的樣本數，$\alpha$是衰減係數

令$p[n]$是脈衝響應，描述原始訊號與回聲訊號之間的關係

$p[n]=\delta [n]+\alpha \delta [n-N_p]$，其中$\delta$是單位脈衝函數

$y[n]=x[n]+\alpha x[n-N_p]=x[n]*p[n]$

系統函數 $P(Z)=1+\alpha Z^{-N_p}$

透過對系統函數進行對數轉換，簡化回聲成分的分析和處理

$\hat{P}(Z)=log(1+\alpha Z^{-N_p})=\sum _{k=1}^{\mathrm{\infty }}(-1{)}^{k+1}\frac{{\alpha }^k}{k}{Z}^{-k{N}_p}$

將$\hat{P}(Z)$轉換到時域

$\hat{p}[n]=\sum _{k=1}^{\mathrm{\infty }}(-1{)}^{k+1}\frac{{\alpha }^k}{k}\delta (n-k\cdot N_P)$

(2) 聲學工程

$y[n]=x[n]*h[n]$，其中$y[n]$是合成音樂，$x[n]$是原始音樂，$h[n]$是脈衝響應(例如建築物空間的影響)

(3) 語音分析

透過在complex cepstrum domain中進行濾波，分離這些成分，使得對語音訊號的理解和處理更為精確。

$s[n]=g[n]*v[n]*p[n]$，其中$s[n]$是語音波，$g[n]$是全局波形，$v[n]$是聲道脈衝，$p[n]$是音高，*是卷積

(4) 地震信號分析

(5) 任意波傳播的多路徑分析

梅爾頻率倒頻譜

梅爾頻率倒頻譜是倒頻譜的一種應用，梅爾頻率倒頻譜常應用在聲音訊號處理，對於聲音訊號處理比倒頻譜更接近人耳對聲音的分析特性，而梅爾頻率倒頻譜與倒頻譜的差別在於:

1. 梅爾頻率倒頻譜的頻帶分析是根據人耳聽覺特性所設計，人耳對於頻率的分辨能力，是由頻率的"比值"決定，也就是說，人耳對200赫茲和300赫茲之間的差別與2000赫茲和3000赫茲之間的差別是相同的。
2. 梅爾頻率倒頻譜是針對訊號的能量取對數，而倒頻譜是針對訊號原始在頻譜上的值取對數。
3. 梅爾頻率倒頻譜是使用離散餘弦轉換，倒頻譜是用離散傅立葉變換。
4. 梅爾頻率倒頻譜係數足夠描述語音的特徵。

梅爾頻率倒頻譜係數(MFCCs)的推導步驟：

1. 將信號做傅立葉變換
2. 頻譜上的值取絕對值再平方成為能量，在乘上頻譜上對應的梅爾頻率倒頻譜三角重疊窗(window)的係數。
3. 對每個梅爾頻率取對數。
4. 作離散餘弦轉換。
5. 求得梅爾頻率倒頻譜係數。

梅爾頻率倒頻譜應用

• 梅爾頻率倒頻譜係數常利用在辨認語音技術上，例如辨認電話中說話的人的身份。
• 利用每種樂風、或樂器在梅爾頻域上有不同特性來分析音樂的種類與類型，並且可以加以分類。

雜訊敏感性

梅爾頻率倒頻譜係數很容易被外來的雜訊所破壞，因此有些研究結果指出，在求梅爾頻率倒頻譜係數時，在作離散餘弦轉換前，提升適當的能量(大約2或3倍)，以減少雜訊在低能量成份的影響。

梅爾頻率倒頻譜優點

相較於原始的倒頻譜

• 有絕對值平方

卷積

倒頻譜領域上的一項重要的特性為二信號卷積之產生，其產生之程序為二倒頻譜值(cepstra)之相加：

$x_1*x_2\to x{\text{'}}_1+x{\text{'}}_2$

微分倒頻譜(differential cepstrum)

定義

${\hat{x}}_d(n)=Z^{-1}\frac{X^{\prime }(Z)}{X(Z)}$ 或 ${\hat{x}}_d[n]=\underset{-\frac{1}{2}}{\overset{\frac{1}{2}}{\int }}\frac{X^{\prime }(F)}{X(F)}{e}^{i2\pi F}dF$

$(\frac{d}{dZ}{\hat{X}}_d(Z)=\frac{d}{dZ}logX(Z)=\frac{X^{\prime }(Z)}{X(Z)})$
If $x(n)=x_1(n)*x_2(n)$
$X(Z)=X_1(Z)X_2(Z)$
$X^{\prime }(Z)={X_1}^{\prime }(Z)X_2(Z)+X_1(Z){X_2}^{\prime }(Z)$
$\frac{X^{\prime }(Z)}{X(Z)}=\frac{{X_1}^{\prime }(Z)}{X_1(Z)})+\frac{{X_2}^{\prime }(Z)}{X_2(Z)})$
$\therefore {\hat{x}}_d(n)={\hat{x}}_{1d}(n)+{\hat{x}}_{2d}(n)$
優點: (a)沒有模糊的相位 (b)可以處理延遲問題

特性

(1)微分倒頻譜在shift和scaling時，結果不改變。
ex: $y[n]=AX[n-r]$
$\Rightarrow {\hat{y}}_d(n)=\{\begin{array}{c}{\hat{x}}_d(n),n\ne 1\\ -r+{\hat{x}}_d(1),n=1\end{array}$
(proof):
$Y(z)=A{z}^{-r}X(z)$
$Y(z)=A{z}^{-r}{X}^{\prime }(z)-rA{z}^{-r-1}X(z)$
$\frac{Y^{\prime }(z)}{Y(z)}=\frac{X^{\prime }(z)}{X(z)}-r{z}^{-1}$
(2)複數倒頻譜$\hat{C}[n]$ 與 微分倒頻譜 ${\hat{x}}_d[n]$和原訊號x[n]有關
$\hat{C}(n)=\frac{-{\hat{x}}_d(n+1)}{n},n\ne 0$ diff cepstrum
$-(n-1)x(n-1)=\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}{\hat{x}}_d(n)x(n-k)$ recursive formula
$\Rightarrow$複數頻譜做得到的事情, 微分倒頻譜也做得到
(3)如果x[n]是最小相位(minimum phase),則${\hat{x}}_d[n]=0$,當$n\le 0$
minimum phase 意思為 no poles 或 zeros 在單位圓外
(4)如果x[n]是最大相位(maximum phase),則${\hat{x}}_d[n]=0$,當$n\ge 2$
maximum phase 意思為 no poles 或 zeros 在單位圓內
(5)如果x(n)為有限區間,則${\hat{x}}_d[n]$為無限區間

• 複數倒頻譜的衰減率反比於n
• 微分倒頻譜的衰減率下降

$\therefore {\hat{x}}_d(n+1)=n\hat{c}(n)\propto n\frac{1}{n}=1$

範例

• $x[0]=1,x[1]=0.5$ ,otherwise 0 , Find its cepstrum.

$x[n]\stackrel{Ztransform}{⟶}X(Z)\stackrel{log}{⟶}\hat{X}(Z)\stackrel{Z^{-1}}{⟶}\hat{x}[n]$

step 1. Z transform: $X(Z)=1+0.5Z^{-1},pole=-0.5$
step 2. log: $\hat{X}(Z)=\sum _{k=1}^{m_i}log(1-(-0.5Z^{-1}))$
step 3. reverse Z transform: $\hat{x}[n]=\sum _{n=0}^N-\frac{-0.5^n}{n},n>0$

• $\hat{x}[0]=1$ ,otherwise 0 , Find its inverse cepstrum.

$\hat{x}[n]\stackrel{Ztransform}{⟶}\hat{X}(Z)\stackrel{exp}{⟶}X(Z)\stackrel{Z^{-1}}{⟶}x[n]$

step 1. Z transform: $\hat{X}[n]=Z^{-1}$
step 2. exp: $e(\frac{1}{z})=\sum _{n=0}^{\mathrm{\infty }}\frac{\frac{1}{z^n}}{n!}$
step 3. reverse Z transform: $x[n]=\{\begin{array}{c}\frac{1}{n!},n\ge 0\\ 0,otherwise\end{array}$

• Suppose that an IIR filter is $H(Z)=\frac{2z^3-4z^2-z+2}{2z^2-2z+1}$

$x[n]\stackrel{Ztransform}{⟶}X(Z)\stackrel{log}{⟶}\hat{X}(Z)\stackrel{Z^{-1}}{⟶}\hat{x}[n]$

step 1. Z transform: $H(Z)=\frac{(-2)(z)(z-\frac{\sqrt{2}}{2}{z}^{-1})(z+\frac{\sqrt{2}}{2}{z}^{-1})(1-\frac{1}{2}z)}{(1-\frac{1+j}{2}{z}^{-1})(1-\frac{1-j}{2}{z}^{-1})}$
step 2. log: $\hat{H}(Z)=log(-2)+3log(z)+log(1\pm \frac{\sqrt{2}}{2}{z}^{-1})+log(1-\frac{1}{2}z)-log(1-\frac{1\pm j}{2}{z}^{-1})$
step 3. reverse Z transform: $\hat{h}[n]=\{\begin{array}{c}log(-2),n=0\\ -\frac{{(\frac{\sqrt{2}}{2})}^n+{(\frac{-\sqrt{2}}{2})}^n}{n}+\frac{{(\frac{1+j}{2})}^n+{(\frac{1-j}{2})}^n}{n},n>0\\ \frac{{(\frac{1}{2})}^{-n}}{n},n<0\end{array}$

${\displaystyle }$

參考文獻

1. B. P. Bogert, M. J. R. Healy, and J. W. Tukey: "The quefrency analysis of time series for echoes: cepstrum, pseudo-autocovariance, cross-cepstrum, and saphe cracking". Proceedings of the Symposium on Time Series Analysis (M. Rosenblatt, Ed) Chapter 15, 209-243. New York: Wiley, 1963.
2. D. G. Childers, D. P. Skinner, R. C. Kemerait, "The Cepstrum: A Guide to ProcessingTemplate:Wayback," Proceedings of the IEEE, Vol. 65, No. 10, October 1977, pp. 1428-1443.
3. Jian-Jiun Ding, Advanced Digital Signal Processing class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2008
4. Jian-Jiun Ding, Advanced Digital Signal Processing class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2024