倒频谱

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

倒频谱的范例

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

• 字义：倒频谱（信号）是信号频谱取对数的傅里叶变换后的新频谱（信号），有时候会称频谱的倒频谱。
• 数学上：信号的倒频谱 = IFT ( log ( | FT (信号) | ) + j2πm )（m为实数）
• 算法：信号 -> 傅立叶变换 -> 取绝对值 -> 取对数 -> 相位展开 -> 逆傅立叶变换 -> 倒频谱

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

复数倒频谱与实数倒频谱

复数倒频谱

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

实数倒频谱

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

应用

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

倒频谱观念

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

倒滤波器

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

计算倒频谱的方法

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

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

利用Z变换的零点与极点

先对信号做Z变换, 并整理一下系数, 让他变成下面的形式
${\displaystyle 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)}}}$
其中${\displaystyle \left|a_{k}\right|,\left|b_{k}\right|,\left|c_{k}\right|,\left|d_{k}\right|\leq 1}$

分子:
第一项A是系数
第二项${\displaystyle Z^{r}}$是延迟
第三项是位于单位圆内的零点
第四项是位于单位圆外的零点

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

对${\displaystyle X\left(Z\right)}$取log变成${\displaystyle {\widehat {X}}\left(Z\right)}$
${\displaystyle {\widehat {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的逆变换
${\displaystyle {\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.${\displaystyle {\widehat {x}}\left[n\right]}$总是IIR(无限冲激响应)
2.对于FIR(有限冲激响应)的情况, ${\displaystyle c_{k}=0,d_{k}=0}$

利用Z变换与微分

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

分别对于x[n]的四种不同的状况做延伸
1.对于x[n]是因果(causal)和最小相位(minimum phase) i.e. ${\displaystyle x[n]={\widehat {x}}\left[n\right]=0,n<0}$
对于${\displaystyle x[n]=\sum _{k=-\infty }^{\infty }{\frac {k}{n}}{\widehat {x}}\left[k\right]x[n-k]\quad for\ n\neq 0}$
可得出
${\displaystyle x[n]=\sum _{k=0}^{\infty }{\frac {k}{n}}{\widehat {x}}\left[k\right]x[n-k]\quad for\ n>0}$
故
${\displaystyle 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.对于x[n]是最小相位(minimum phase)
${\displaystyle {\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.对于x[n]是反因果(anti-causal)且最大相位(maximum phase) i.e. ${\displaystyle x[n]={\widehat {x}}\left[n\right]=0,n>0}$
${\displaystyle {\begin{aligned}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{aligned}}}$
4.对于x[n]是最大相位(maximum phase)
${\displaystyle {\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}}}$

特性

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

同态解卷积的应用(Application of Homomorphic Deconvolution)

同态解卷积有非常多应用面，尤其是在声学工程和语音分析方面的实用性

(1) 回声的均衡化

${\displaystyle y[n]=x[n]+{\alpha }x[n-N_{p}]}$

其中${\displaystyle y[n]}$是接收到的信号，${\displaystyle x[n]}$是原始信号，${\displaystyle N_{p}}$是延迟的样本数，${\displaystyle {\alpha }}$是衰减系数

令${\displaystyle p[n]}$是冲激响应，描述原始信号与回声信号之间的关系

${\displaystyle p[n]={\delta }[n]+{\alpha }{\delta [n-N_{p}]}}$，其中${\displaystyle {\delta }}$是单位脉冲函数

${\displaystyle y[n]=x[n]+{\alpha }x[n-N_{p}]=x[n]*p[n]}$

系统函数 ${\displaystyle P(Z)=1+{\alpha }Z^{-N_{p}}}$

透过对系统函数进行对数变换，简化回声成分的分析和处理

${\displaystyle {\hat {P}}(Z)=log(1+{\alpha }Z^{-N_{p}})={\sum _{k=1}^{\infty }}(-1)^{k+1}{\frac {\alpha ^{k}}{k}}Z^{-kN_{p}}}$

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

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

(2) 声学工程

${\displaystyle y[n]=x[n]*h[n]}$，其中${\displaystyle y[n]}$是合成音乐，${\displaystyle x[n]}$是原始音乐，${\displaystyle h[n]}$是冲激响应(例如建筑物空间的影响)

(3) 语音分析

透过在complex cepstrum domain中进行滤波，分离这些成分，使得对语音信号的理解和处理更为精确。

${\displaystyle s[n]=g[n]*v[n]*p[n]}$，其中${\displaystyle s[n]}$是语音波，${\displaystyle g[n]}$是全局波形，${\displaystyle v[n]}$是声道脉冲，${\displaystyle p[n]}$是音高，*是卷积

(4) 地震信号分析

(5) 任意波传播的多路径分析

梅尔频率倒频谱

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

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

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

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

梅尔频率倒频谱应用

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

噪声敏感性

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

梅尔频率倒频谱优点

相较于原始的倒频谱

• 有绝对值平方

卷积

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

${\displaystyle x_{1}*x_{2}\rightarrow x'_{1}+x'_{2}}$

微分倒频谱(differential cepstrum)

定义

${\displaystyle {\widehat {x}}_{d}(n)=Z^{-1}{\frac {X'(Z)}{X(Z)}}}$ 或 ${\displaystyle {\widehat {x}}_{d}[n]=\int _{-{\frac {1}{2}}}^{\frac {1}{2}}{\frac {X'(F)}{X(F)}}e^{i2\pi F}dF}$

${\displaystyle ({\frac {d}{dZ}}{\widehat {X}}_{d}(Z)={\frac {d}{dZ}}logX(Z)={\frac {X'(Z)}{X(Z)}})}$
If ${\displaystyle x(n)=x_{1}(n)*x_{2}(n)}$
${\displaystyle X(Z)=X_{1}(Z)X_{2}(Z)}$
${\displaystyle X'(Z)=X_{1}'(Z)X_{2}(Z)+X_{1}(Z)X_{2}'(Z)}$
${\displaystyle {\frac {X'(Z)}{X(Z)}}={\frac {X_{1}'(Z)}{X_{1}(Z)}})+{\frac {X_{2}'(Z)}{X_{2}(Z)}})}$
${\displaystyle \therefore {\widehat {x}}_{d}(n)={\widehat {x}}_{1d}(n)+{\widehat {x}}_{2d}(n)}$
优点: (a)没有模糊的相位 (b)可以处理延迟问题

特性

(1)微分倒频谱在shift和scaling时，结果不改变。
ex: ${\displaystyle y[n]=AX[n-r]}$
${\displaystyle \Rightarrow {\widehat {y}}_{d}(n)={\begin{cases}{\widehat {x}}_{d}(n),n\neq 1\\-r+{\widehat {x}}_{d}(1),n=1\end{cases}}}$
(proof):
${\displaystyle Y(z)=Az^{-r}X(z)}$
${\displaystyle Y(z)=Az^{-r}X'(z)-rAz^{-r-1}X(z)}$
${\displaystyle {\frac {Y'(z)}{Y(z)}}={\frac {X'(z)}{X(z)}}-rz^{-1}}$
(2)复数倒频谱${\displaystyle {\widehat {C}}[n]}$ 与 微分倒频谱 ${\displaystyle {\widehat {x}}_{d}[n]}$和原信号x[n]有关
${\displaystyle {\widehat {C}}(n)={\frac {-{\widehat {x}}_{d}(n+1)}{n}},n\neq 0}$ diff cepstrum
${\displaystyle -(n-1)x(n-1)=\sum _{k=-\infty }^{\infty }{\widehat {x}}_{d}(n)x(n-k)}$ recursive formula
${\displaystyle \Rightarrow }$复数频谱做得到的事情, 微分倒频谱也做得到
(3)如果x[n]是最小相位(minimum phase),则${\displaystyle {\widehat {x}}_{d}[n]=0}$,当${\displaystyle n\leq 0}$
minimum phase 意思为 no poles 或 zeros 在单位圆外
(4)如果x[n]是最大相位(maximum phase),则${\displaystyle {\widehat {x}}_{d}[n]=0}$,当${\displaystyle n\geq 2}$
maximum phase 意思为 no poles 或 zeros 在单位圆内
(5)如果x(n)为有限区间,则${\displaystyle {\widehat {x}}_{d}[n]}$为无限区间

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

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

范例

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

${\displaystyle x[n]\quad {\stackrel {Ztransform}{\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: ${\displaystyle X(Z)=1+0.5Z^{-1},pole=-0.5}$
step 2. log: ${\displaystyle {\widehat {X}}(Z)=\sum _{k=1}^{m_{i}}log(1-(-0.5Z^{-1}))}$
step 3. reverse Z transform: ${\displaystyle {\widehat {x}}[n]=\sum _{n=0}^{N}-{\frac {-0.5^{n}}{n}},n>0}$

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

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

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

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

${\displaystyle x[n]\quad {\stackrel {Ztransform}{\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: ${\displaystyle 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: ${\displaystyle {\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: ${\displaystyle {\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}}}$

${\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 Processing（页面存档备份，存于互联网档案馆）," 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