倒频谱

<image>

倒频谱（cepstrum），顾名思义，就是将频谱（spectrum）的英文前四个字母反过来写。倒频谱是为了某些时候，为了计算方便，将原来信号的频谱先转成类似分贝的单位，再作逆傅里叶变换，把它视为一种新的讯号做处理。倒频谱有复数倒频谱，及实数倒频谱。

倒频谱被定义在1963的论文（Bogert等）。定义如下：

- 字义：倒频谱（信号）是信号频谱取对数的傅立叶变换后的新频谱（信号），有时候会称频谱的倒频谱。

- 算法：信号 -> 傅立叶变换 -> 取绝对值 -> 取对数 -> 相位展开 -> 逆傅立叶变换 -> 倒频谱

复数倒频谱拥有频谱大小跟相位的资讯，实数倒频谱只有频谱大小的资讯，各有各的不同应用。

## 复数倒频谱与实数倒频谱

### 复数倒频谱

$\widehat {x}\left[ n\right] =\int _{-\frac {1} {2}}^{\frac {1} {2}}\widehat {X}\left( F\right) e^{j {2 \pi }F}dF$ 其中$\widehat {X}\left[ F\right] = \log |X(F)| + j \arg[X(F)]$ 可能遭遇的问题1. $\log 0 = - \infty$ 2. $\arg[X[n]]$有无限多的解当输入是实数时,因为$\log |X(F)|$偶对称，$\arg[X(F)]$奇对称,所以复数倒频谱的值为实数

### 实数倒频谱

$C \left[ n\right] =\int _{-\frac {1} {2}}^{\frac {1} {2}} \log |X(F)| e^{j {2 \pi }Fn}dF$ 可能遭遇的问题1. $\log 0 = - \infty$

## 应用

- 倒频谱可以被视为在不同频带上变化速率的资讯，倒频谱一开始被发明在地震或炸弹产生的地震回音，现今也被使用在分析雷达讯号，以及讯号处理等问题。
- 自相关倒频谱(autocepstrum)被定义为倒频谱的自相关性，自相关倒频谱有时在分析处理回传讯号时比倒频谱还准确。
- 倒频谱在处理人声讯号以及音乐讯号有非常好的效果，例如梅尔频率倒频谱(Mel-Frequency Cepstrum)，用来做声音的辨认，侦测音高等。近年来梅耳倒频谱也被应用在音乐资讯的回复。
- 倒频谱在声学中可以将声带震动的影响去除。
- 倒频谱用在处理多路径问题时(如声波的回音、电磁波的折、反射等)，如果将其他路径干扰视为噪声，为了消除噪声，利用倒频谱，不需测量每条多路径的延迟时间，可以利用传送多次信号，观察其他路径在倒频谱上的效果，并且加以滤除。
- 语音大致上是由音高、声带脉冲、声门波形所组成，我们可以利用倒频谱将这三种元素在倒频域上分开，以利于做语音讯号的分析。
- 倒频谱的微分适用于影像处理上的图形辨认(pattern recognition)。
- 倒频谱与同型声音理论(homomorphic sound theory)有关。

## 倒频谱观念

频谱图上的独立变数是频率，而倒频谱图上的独立变数为倒频率(quefrency)，倒频率是一种时间的度量单位。举个例子，声音讯号取样速率等于44100赫兹，在倒频谱上有个很大的值在倒频率等于100，代表实际上在44100/100=441赫兹有很大的值，这值出现在倒频谱上因为频谱上周期性出现，而频谱上出现的周期与倒频谱很大的值出现的位置有关。

## 倒滤波器

滤波器(filter)常使用在频谱上，用来保存或删除我们所要或不要的资讯，经过上面的许多讨论，不难猜到，倒滤波器(lifter)就是在倒频谱上所使用的滤波器。低通的倒滤波器跟低通滤波器有点类似，它可以借由在倒频谱上乘以一个window系数，使倒频谱上的高倒频率被压抑，如此依来，当信号转回时域空间时会变成一个较平滑的信号。

## 计算倒频谱的方法

### 直接计算IDTFT(反离散时间傅立叶变换)

$\widehat {x}\left[ n\right] =\int _{-\frac {1} {2}}^{\frac {1} {2}}\widehat {X}\left( F\right) e^{j {2 \pi }F}dF$ 问题: $\widehat {X}\left( F\right)$ 可能会无限大, 且对于arg(xn)有无限多个解

### 利用Z转换的零点与极点

先对信号做Z转换, 并整理一下系数, 让他变成下面的形式$X\left( Z \right)=\cfrac{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| \leq 1$

分子:

第一项A是系数第二项$Z^r$是延迟第三项是位于单位圆内的零点第四项是位于单位圆外的零点

分母:

第一项是位于单位圆内的极点第二项是位于单位圆外的极点

对$X\left( Z \right)$取log变成$\widehat {X}\left( Z\right)$ $\widehat {X}\left( Z\right)=log X\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的反转换$\widehat {x}\left[n\right]=\begin{cases} \log A & \mbox{if }n=0 \\ - \sum_{k=1}^{m_i} \cfrac{{a_k}^n}{n} + \sum_{k=1}^{P_i} \cfrac{{c_k}^n}{n} & \mbox{if }n>0 \\ \sum_{k=1}^{m_0} \cfrac{{b_k}^{-n}}{n} - \sum_{k=1}^{P_0} \cfrac{{d_k}^{-n}}{n} & \mbox{if }n<0 \end{cases}$

注意事项1.$\widehat {x}\left[ n\right]$总是IIR(无限脉冲响应)

2.对于FIR(有限脉冲响应)的情况, $c_k=0, d_k=0$

### 利用Z转换与微分

$Z \cdot \widehat {X}'\left( Z\right)=Z \cdot \cfrac{{X}'\left( Z\right)}{{X}\left( Z\right)}$ $Z {X}'\left( Z\right) = Z \widehat {X}'\left( Z\right) \cdot {X}\left( Z\right)$ 对其做Z的反转换$nx[n]= \sum_{k=- \infty }^{\infty } k \widehat {x}\left[ k\right] x[n-k]$ 故$x[n]= \sum_{k=- \infty }^{\infty } \frac{k}{n} \widehat {x}\left[ k\right] x[n-k] \quad for \ n \neq 0$

分别对于xn的四种不同的状况做延伸1.对于xn是因果(causal)和最小相位(minimum phase) i.e. $x[n]=\widehat {x}\left[ n\right]=0, n<0$ 对于$x[n]= \sum_{k=- \infty }^{\infty } \frac{k}{n} \widehat {x}\left[ k\right] x[n-k] \quad for \ n \neq 0$ 可得出$x[n]= \sum_{k=0}^{\infty } \frac{k}{n} \widehat {x}\left[ k\right] x[n-k] \quad for \ n>0$ 故$x[n]= \widehat {x}\left[ n\right]x[0] + \sum_{k=0}^{n-1} \frac{k}{n} \widehat {x}\left[ k\right] x[n-k]$ 2.对于xn是最小相位(minimum phase)

$\widehat {x}\left[ n\right]=\begin{cases} 0 & \mbox{if }n<0 \\ \cfrac{x[n]}{x[0]} - \sum_{k=0}^{n-1} \cfrac{k}{n} \widehat {x}\left[ k\right] \cfrac{x[n-k]}{x[0]} & \mbox{if }n>0 \\ \log A & \mbox{if }n=0 \end{cases}$ 3.对于xn是反因果(anti-causal)且最大相位(maximum phase) i.e. $x[n]=\widehat {x}\left[ n\right]=0, n>0$ $\begin{align} x[n] & = \sum_{k=n}^{0} \cfrac{k}{n}\widehat {x}\left[ k\right] x[n-k] \quad for \ n<0 \\ & = \widehat {x}\left[ n\right] x[0] + \sum_{k=n+1}^{0} \cfrac{k}{n} \widehat {x}\left[ k\right] x[n-k] \\ \end{align}$ 4.对于xn是最大相位(maximum phase)

$\widehat {x}\left[ n\right]= \begin{cases} 0 & \mbox{if }n>0 \\ \cfrac{x[n]}{x[0]} - \sum_{k=n+1}^{0} \cfrac{k}{n} \widehat {x}\left[ k\right] \cfrac{x[n-k]}{x[0]} & \mbox{if }n<0 \\ \log A & \mbox{if }n=0 \end{cases}$

## 特性

- 复数倒频谱至少以$\frac{1}{n}$的速度衰退$|\widehat {x}\left[ n\right]|=c |\frac{{\alpha}^n}{n}| \quad - \infty<n<\infty$ 其中 $\alpha =max(a_k,b_k,c_k,d_k)$ 2. 如果X(Z)没有在单位圆以外的零点和极点, 则$\widehat {x}\left[ n\right]=0 \quad for \ all \ n<0$ 因为$b_k,d_k=0$ 3. 如果X(Z)没有在单位圆以内的零点和极点, 则$\widehat {x}\left[ n\right]=0 \quad for \ all \ n>0$ 因为$a_k,c_k=0$ 4. 如果xn是有限长度, 则$\widehat {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}^{\infty}}(-1)^{k+1}{\frac{\alpha^{k}}{k}}Z^{-kN_p}$$

将${\hat{P}}(Z)$转换到时域

$${\hat{p}}[n]={\sum_{k=1}^{\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) 任意波传播的多路径分析

## 梅尔频率倒频谱

梅尔频率倒频谱是倒频谱的一种应用，梅尔频率倒频谱常应用在声音讯号处理，对于声音讯号处理比倒频谱更接近人耳对声音的分析特性，而梅尔频率倒频谱与倒频谱的差别在于:

- 梅尔频率倒频谱的频带分析是根据人耳听觉特性所设计，人耳对于频率的分辨能力，是由频率的"比值"决定，也就是说，人耳对200赫兹和300赫兹之间的差别与2000赫兹和3000赫兹之间的差别是相同的。
- 梅尔频率倒频谱是针对讯号的能量取对数，而倒频谱是针对讯号原始在频谱上的值取对数。
- 梅尔频率倒频谱是使用离散余弦转换，倒频谱是用离散傅立叶变换。
- 梅尔频率倒频谱系数足够描述语音的特征。

梅尔频率倒频谱系数(MFCCs)的推导步骤：

- 将信号做傅立叶变换
- 频谱上的值取绝对值再平方成为能量，在乘上频谱上对应的梅尔频率倒频谱三角重叠窗(window)的系数。
- 对每个梅尔频率取对数。
- 作离散余弦转换。
- 求得梅尔频率倒频谱系数。

## 梅尔频率倒频谱应用

- 梅尔频率倒频谱系数常利用在辨认语音技术上，例如辨认电话中说话的人的身份。
- 利用每种乐风、或乐器在梅尔频域上有不同特性来分析音乐的种类与类型，并且可以加以分类。

## 噪声敏感性

梅尔频率倒频谱系数很容易被外来的噪声所破坏，因此有些研究结果指出，在求梅尔频率倒频谱系数时，在作离散余弦转换前，提升适当的能量(大约2或3倍)，以减少噪声在低能量成分的影响。

## 梅尔频率倒频谱优点

相较于原始的倒频谱

- 有绝对值平方

## 卷积

倒频谱领域上的一项重要的特性为二信号卷积之产生，其产生之程序为二倒频谱值(cepstra)之相加：

$$x_1 * x_2 \rightarrow x'_1 + x'_2$$

## 微分倒频谱(differential cepstrum)

### 定义

```
<math>\widehat {x}_d(n) = Z^{-1}\frac{X'(Z)}{X(Z)} </math> 或 <math>\widehat {x}_d[n] = \int_{-\frac{1}{2}}^{\frac{1}{2}}\frac{X'(F)}{X(F)} e^{i2 \pi F}dF</math><br>
```

$(\frac{d}{dZ}\widehat {X}_d(Z) = \frac{d}{dZ}logX(Z) = \frac{X'(Z)}{X(Z)})$ If $x(n) = x_1(n)*x_2(n)$ $X(Z) = X_1(Z)X_2(Z)$ $X'(Z) = X_1'(Z)X_2(Z)+X_1(Z)X_2'(Z)$ $\frac{X'(Z)}{X(Z)} = \frac{X_1'(Z)}{X_1(Z)})+\frac{X_2'(Z)}{X_2(Z)})$ $\therefore \widehat {x}_d(n) = \widehat {x}_{1d}(n)+\widehat {x}_{2d}(n)$ 优点:

(a)没有模糊的相位(b)可以处理延迟问题

### 特性

(1)微分倒频谱在shift和scaling时，结果不改变。

ex: $y[n] = AX[n-r]$ $\Rightarrow \widehat {y}_d(n) = \begin{cases} \widehat {x}_d(n),n \neq 1 \\ <pre> -r + \widehat {x}_d(1),n = 1 \end{cases}$

(proof):

$Y(z) = Az^{-r}X(z)$ $Y(z) = Az^{-r}X'(z)-rAz^{-r-1}X(z)$ $\frac{Y'(z)}{Y(z)} = \frac{X'(z)}{X(z)} - rz^{-1}$ (2)复数倒频谱$\widehat {C}[n]$ 与 微分倒频谱 $\widehat {x}_d[n]$和原讯号xn有关$\widehat {C}(n) = \frac{-\widehat {x}_d(n+1)}{n}, n \neq 0$ diff cepstrum $-(n-1)x(n-1) = \sum_{k = -\infty}^{\infty} \widehat {x}_d(n)x(n-k)$ recursive formula $\Rightarrow$复数频谱做得到的事情, 微分倒频谱也做得到(3)如果xn是最小相位(minimum phase),则$\widehat {x}_d[n]=0$,当$n \leq 0$ minimum phase 意思为 no poles 或 zeros 在单位圆外(4)如果xn是最大相位(maximum phase),则$\widehat {x}_d[n]=0$,当$n \geq 2$ maximum phase 意思为 no poles 或 zeros 在单位圆内(5)如果x(n)为有限区间,则$\widehat {x}_d[n]$为无限区间

- 复数倒频谱的衰减率反比于n
- 微分倒频谱的衰减率下降

$\therefore \widehat {x}_d(n+1) = n\widehat {c}(n) \varpropto n \frac{1}{n} = 1$

## 范例

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

$x[n] \quad \stackrel{Z transform}{\longrightarrow} \quad X(Z)\quad \stackrel{log}{\longrightarrow} \quad \widehat {X}(Z) \quad \stackrel{Z^{-1}}{\longrightarrow}\quad \widehat {x}[n]$

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

- $\widehat x[0]=1$ ,otherwise 0 , Find its inverse cepstrum.

$\widehat x[n] \quad \stackrel{Z transform}{\longrightarrow} \quad \widehat X(Z)\quad \stackrel{exp}{\longrightarrow} \quad {X}(Z) \quad \stackrel{Z^{-1}}{\longrightarrow}\quad {x}[n]$

step 1. Z transform: $\widehat X[n]=Z^{-1}$ step 2. exp: $e(\frac{1}{z}) = \sum_{n=0}^\infty \frac{\frac{1}{z^n}}{n!}$ step 3. reverse Z transform: $x[n] = \begin{cases} \frac{1}{n!},n \geq 0 \\ 0, otherwise\\\end{cases}$

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

$x[n] \quad \stackrel{Z transform}{\longrightarrow} \quad X(Z)\quad \stackrel{log}{\longrightarrow} \quad \widehat {X}(Z) \quad \stackrel{Z^{-1}}{\longrightarrow}\quad \widehat {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: $\widehat {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: $\widehat {h}[n] = \begin{cases} log(-2),n= 0 \\ \displaystyle{-\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 }\\ \displaystyle{\frac{{(\frac{1}{2})}^{-n}}{n}},n<0\\\end{cases}$

## 参考文献

- 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.
- D. G. Childers, D. P. Skinner, R. C. Kemerait, " The Cepstrum: A Guide to Processing," Proceedings of the IEEE, Vol. 65, No. 10, October 1977, pp. 1428-1443.
- Jian-Jiun Ding, Advanced Digital Signal Processing class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2008
- Jian-Jiun Ding, Advanced Digital Signal Processing class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2024