[离散时间信号处理学习笔记] 14. 多采样率信号处理
多采样率信号处理一般是指利用增采样、减采样、压缩器和扩展器等方式来提高信号处理系统效率的技术（These multirate techniques refer in general to utilizing upsampling, downsampling, compressors, and expanders in a variety of ways to increase the efficiency of signal-processing systems. ）本文章主要讨论多采样率技术中的两个研究成果：滤波与压缩器/扩展器的互换；多相分解。
压缩器/扩展器的时域与频域表示
尽管上一篇文章中已经讨论过这部分内容，不过由于这部分是理解本文所必须的关键知识点，这里将在时域与频域展开更详细的分析。
压缩器
假设压缩器的压缩率为$M$，那么压缩器在时域上的表示为
$x_d[n] = x[nM]$
$x[n]$的采样频率为$T$，那么$x_d[n]$的采样频率为$T_d = MT$，按照离散序列与连续信号在频域上的关系，有
$\begin{align*}
X(e^{j\omega}) &= \frac{1}{T}\sum_{k=-\infty}^{\infty}X_c\left[ j\left(\frac{\omega}{T}-\frac{2\pi k}{T}\right)\right ]\\
X_d(e^{j\omega}) &= \frac{1}{MT}\sum_{r=-\infty}^{\infty}X_c\left[ j\left(\frac{\omega}{MT}-\frac{2\pi r}{MT}\right)\right ]
\end{align*}$
压缩前的序列频谱$X(e^{j\omega})$与压缩后的序列频谱$X_d(e^{j\omega})$之间有如下关系
$\begin{align*}
X_d(e^{j\omega}) &= \frac{1}{MT}\sum_{r=-\infty}^{\infty}X_c\left[ j\left(\frac{\omega}{MT}-\frac{2\pi r}{MT}\right)\right ] \\
& = \frac{1}{MT}\left\{\cdot\cdot\cdot+X_c\left[j\left(\frac{\omega}{MT}-\frac{-2\pi}{MT} \right ) \right ] +X_c\left[ j\left(\frac{\omega}{MT}-\frac{0}{MT}\right)\right ] + X_c\left[j\left(\frac{\omega}{MT}-\frac{2\pi}{MT} \right ) \right ]+\cdot\cdot\cdot \right \}\\
& = \frac{1}{MT}\left\{\cdot\cdot\cdot+X_c\left[j\left(\frac{\omega}{MT}-\frac{0}{MT} \right ) \right ]+\cdot\cdot\cdot +X_c\left[ j\left(\frac{\omega}{MT}-\frac{2(M-1)\pi}{MT}\right)\right ]\right.\\
&\quad\qquad\qquad\left.+ X_c\left[j\left(\frac{\omega}{MT}-\frac{2M\pi}{MT} \right ) \right ]+\cdot\cdot\cdot+ X_c\left[j\left(\frac{\omega}{MT}-\frac{2M\pi}{MT} -\frac{2(M-1)\pi}{MT}\right ) \right ]+\cdot\cdot\cdot \right \}\\
\end{align*}$
$\begin{align*}
\qquad\quad\ &= \frac{1}{MT}\left\{\cdot\cdot\cdot+\sum_{i=0}^{M-1}X_c\left[j\left(\frac{\omega}{MT}-\frac{2i\pi}{MT} \right ) \right ]+\sum_{i=0}^{M-1}X_c\left[j\left(\frac{\omega}{MT}-\frac{2i\pi}{MT}-\frac{2\pi}{T} \right ) \right ]+\cdot\cdot\cdot\right\}\\
&= \frac{1}{MT}\sum_{k=-\infty}^{\infty} \sum_{i=0}^{M-1}X_c\left[j\left(\frac{\omega}{MT}-\frac{2\pi i}{MT}-\frac{2\pi k}{T} \right ) \right ] \\
&=\frac{1}{M}\sum_{i=0}^{M-1}\left\{\frac{1}{T}\sum_{k=-\infty}^{\infty}X_c\left[j\left(\frac{\omega-2\pi i}{MT}-\frac{2\pi k}{T} \right ) \right ]\right\}\\
&=\frac{1}{M}\sum_{i=0}^{M-1}X(e^{j(\omega-2\pi i)/M})
\end{align*}$
如下图所示
扩展器
假设扩展器的扩展率为$L$，那么扩展器在时域上的表示为
$x_e[n] = \left\{\begin{matrix}
x[n/L], &n=0,\pm L,\pm 2L,\cdot\cdot\cdot \\
0, &else
\end{matrix}\right.$
扩展前的序列频谱$X(e^{j\omega})$与扩展后的序列频谱$X_e(e^{j\omega})$之间有如下关系
$\begin{align*}
X_e(e^{j\omega}) &= \sum_{n=-\infty}^{\infty}x_e[n]e^{-j\omega n}\\
&=\sum_{n=-\infty}^{\infty}x[n/L]e^{-j\omega n}\quad n=0,\pm L,\pm 2L,\cdot\cdot\cdot\\
&=\sum_{k=-\infty}^{\infty}x[k]e^{-j\omega kL}\quad letting\ n=kL\\
&=X(e^{j\omega L})
\end{align*}$
如下图所示
滤波与压缩器/扩展器的互换
滤波器与压缩器互换
如上一篇文章所描述的减采样就是一个滤波器与压缩器的级联系统。对于这种级联方式，以下两个系统是等价的
证明：
$\begin{align*}
Y_{R}(e^{j\omega})
&= \frac{1}{M}\sum_{i=0}^{M-1}X_b(e^{j(\omega-2\pi i)/M}) \quad compress\ X_b(e^{j\omega})\ with\ M\\
&= \frac{1}{M}\sum_{i=0}^{M-1}X(e^{j(\omega-2\pi i)/M})H(e^{j(\omega-2\pi i)/M \cdot M})\\
&= \frac{1}{M}\sum_{i=0}^{M-1}X(e^{j(\omega-2\pi i)/M})H(e^{j(\omega-2\pi i)})\\
&= H(e^{j\omega})\frac{1}{M}\sum_{i=0}^{M-1}X(e^{j(\omega-2\pi i)/M})\quad compress\ X(e^{j\omega})\ with\ M\\
&= H(e^{j\omega})X_a(e^{j\omega})\\
&= Y_L(e^{j\omega})
\end{align*}$
这个证明过程运用了前面一小节中的压缩器频谱运算，$Y_{R}(e^{j\omega}) = Y_{L}(e^{j\omega})$说明右边的系统跟左边的系统是完全等价的。这两个等价系统也很容易理解，观察下图
根据前面我们对扩展器的讨论，我们可以把$H(z^M)$理解为对$H(z)$做了因子为$M$的扩展。两个系统分别对应上图的第二第三行：
- 上图第二行先对$x[n]$进行因子为$M$的压缩得到$x_a[n]$，然后级联一个滤波器$H(z)$
- 上图第三行先用滤波器$H(z^M)$对$x[n]$进行滤波得到$x_b[n]$，然后对$x_b[n]$进行因子为$M$的压缩
滤波器与扩展器互换
如上一篇文章所描述的增采样就是一个滤波器与扩展器的级联系统。对于这种级联方式，以下两个系统是等价的
证明：
$\begin{align*}
Y_{L}(e^{j\omega})
&= X_a(e^{j\omega L}) \quad expand\ X_a(e^{j\omega})\ with\ L\\
&= X(e^{j\omega L})H(e^{j\omega L}) \\
&= X_b(e^{j\omega})H(e^{j\omega L})\\
&= Y_R(e^{j\omega})
\end{align*}$
这个证明过程运用了前面一小节中的扩展器频谱运算，$Y_{L}(e^{j\omega}) = Y_{R}(e^{j\omega})$说明右边的系统跟左边的系统是完全等价的。这两个等价系统也很容易理解，观察下图
同样，根据我们前面对扩展器的讨论，我们可以吧$H(z^L)$理解为对$H(z)$进行了因子为$L$的扩展。两个系统分别对应上图的第二第三行：
- 上图第二行先用滤波器$H(z)$对$x[n]$进行滤波得到$x_a[n]$，然后对$x_a[n]$进行因子为$L$的扩展
- 上图第三行先对$x[n]$进行因子为$L$的扩展得到$x_b[n]$，然后级联一个滤波器$H(z^L)$
多级抽取和内插
当抽取或内插率较大时，即$M$或者$L$会非常大，这种情况下对$M$或者$L$进行分解，如$M = M_1M_2$，把单级抽取转换为多级抽取将会使得系统更为高效，具体原因以后（书中第七章）讨论。如下是两级抽取器：
按照上面的步骤反过来，则能把一个一级抽取转换为二级抽取，该一级抽取的的滤波器的系统函数按照如下方式进行分解：
$H(z) = H_1(z)H_2(z^{M_1})$
原滤波器的脉冲响应则为分解出来的两个滤波器的脉冲响应的卷积，如下：
$\displaystyle{h[n] = h_1[n]* \sum_{k=-\infty}^{\infty}h_2[k]\delta[n-kM_1]}$
按照这种方式能把一级抽取器扩展到多级抽取。同理，上述理论也能应用到内插处理中。
多相分解
对一个序列进行多相分解，就是把该序列表示成M组子序列的叠加。如下图：
为了得到这M组子序列，可以采用以下的分解方法：
步骤解析如下
- 首先是对$h[n]$进行分解能得到$\color{red}{e_k[n]}$：
- 在z变换中，频域上的$H(z)\cdot z^k$就相当于时域上$h[n] \delta[n+k]$，即进行相位为$k$的移动，得到$h[n+k] $
- 对$h[n+k]$进行因子为$M$的压缩，得到$e_k[n]$
- 接下来是复原为原来的$h[n]$：
- 对$e_k[n]$进行因子为$M$的扩展，得到$h_k[n]$
- 分别对$M$个$h_k[n]$进行相应的移位，然后就能合并得到$h[n]$
下面是一个$M=4,k=3$的例子
该过程在频域上表示为：
$\displaystyle{ H(z)=\sum_{k=0}^{M-1}E_k(z^M)z^{-k} }$
即系统函数$H(z)$可以被分解成多个经过延迟的滤波器之和，因此就可以把该系统表示成如下并联结构：
抽取/内插滤波器的多相实现
多相分解能降低抽取以及内插实现的计算量。直观地解释就是在抽取的时候，实际上所需要的序列数量只有原序列数量的$\frac{1}{M}$；在内插的时候，对$x[n]$进行因子为$L$的扩展后，每$L$个样本才会有一个非零样本，而零样本是不需要计算的。详情请看下面的分析。
抽取滤波器
原抽取滤波器有如下形式
把抽取滤波器多相分解成如下图左边的形式，然后滤波器与压缩器进行互换就能得到右图的形式。
假设输入序列$x[n]$的采样率为一个单位时间，即一个单位时间采集一个样本。假设滤波器$h[n]$为有限长度滤波器（FIR），长度为N个点，那么原抽取滤波器在每个单位时间需要进行$N$次乘法以及$N-1$次加法；在采用多相分解实现后，对每个分量来说，每$M$个单位时间才采集一次，即每一个时间单位采集$\frac{1}{M}$次，并且滤波器变为$e_k[n]$，这些滤波器的长度为$\frac{N}{M}$，因此对于每一个分量来说，每个单位时间的运算量为$\frac{1}{M}\frac{N}{M}$次乘法、$\frac{1}{M}\left(\frac{N}{M}-1\right)$次加法。由于共有$M$个分量，因此整个系统在每个单位时间的运算量为乘法$\frac{N}{M}$次、加法$\left(\frac{N}{M}-1\right)+(M-1)$次。
内插滤波器
原内插滤波器有如下形式：
把内插滤波器多相分解为如下图左边的形式，然后滤波器与扩展器进行互换就能得到右图的形式。
假设输入序列$x[n]$的采样率为一个单位时间，即一个单位时间采集一个样本。假设滤波器$h[n]$为有限长度滤波器（FIR），长度为N个点，那么在对$x[n]$进行因子为$L$的增采样后，内插滤波器在每个单位时间需要进行$NL$次乘法以及$L(N-1)$次加法；在采用多相分解实现后，对每个分量来说，每一个单位时间采集一次，并且滤波器变为$e_k[n]$，这些滤波器的长度为$\frac{N}{L}$，因此对于每一个分量来说，每个单位时间的运算量为$\frac{N}{L}$次乘法、$\left(\frac{N}{L}-1\right)$次加法。由于共有$L$个分量，因此整个系统在每个单位时间的运算量为乘法$N$次、加法$N-L+(L-1)$次。
多采样滤波器组
这一小节主要描述了一个使用了多相抽取以及内插结构的系统。
理想双信道分解/合成系统
如下图所示为一个语音编码中常用的对音频信号进行双信道分解（analysis）以及合成（synthesis）系统的结构框图。对结构中各部分的解释如下：
- 系统分解部分目的是对音频进行高频以及低频的划分，两个信道的带宽分别为$|\omega|<\frac{\pi}{2}$以及$\frac{\pi}{2}<|\omega|<\pi$，这部分的分解分别由低通滤波器$h_0[n]$以及高通滤波器$h_1[n]$完成。通常取$h_1[n] = e^{j\pi n}h_0[n]$，这意味着只需要对低通滤波器频谱左移$\pi$即可得到高通滤波器频谱，$H_1(e^{j\omega}) = H_0(e^{j(\omega-\pi)})$。
- 由于对源信号并行进行了两组滤波，样本数增加了一倍。为了保持滤波后信号的总样本数不变，需要对滤波后的信号进行压缩，压缩倍率为2，也就是使得低频序列$v_0[n]$以及高频序列$v_1[n]$各占一半的样本数。对信号进行压缩意味着信号频谱的扩展，滤波后的信号频谱的带宽分别为$|\omega|<\frac{\pi}{2}$以及$\frac{\pi}{2}<|\omega|<\pi$，而压缩后会使得两个信道的信号频谱都扩展到了$|\omega|<\pi$。
- 然后可以对$v_0[n],v_1[n]$进行各种处理。
- 系统合成部分的目的就是把两个信道的音频恢复成为全带宽的音频，对于每一个信道来说，此时的样本数只有原样本数的一半，因此需要先把样本进行2倍扩展。
- 最后分别通过一个低通滤波器以及高通滤波器完成重构，再把得到的信号相加即可得到全带宽音频。
理想情况下会有如下频谱变化：
式子推导如下：
$\begin{align*}
Y(e^{j\omega})&=expand\bigg\{compress\bigg\{X(e^{j\omega})H_0(e^{j\omega})\bigg\}\bigg\}G_0(e^{j\omega})+expand\bigg\{compress\bigg\{X(e^{j\omega})H_1(e^{j\omega})\bigg\}\bigg\}G_1(e^{j\omega})\\
&=expand\left\{\frac{1}{2}\bigg[X(e^{j\omega/2})H_0(e^{j\omega/2})+X(e^{j(\omega-\pi)/2})H_0(e^{j(\omega-\pi)/2})\bigg]\right\}G_0(e^{j\omega})\\
&\quad + expand\ \left\{\frac{1}{2}\bigg[X(e^{j\omega/2})H_1(e^{j\omega/2})+X(e^{j(\omega-\pi)/2})H_1(e^{j(\omega-\pi)/2})\bigg]\right\}G_1(e^{j\omega})\\
&= \frac{1}{2}\bigg[X(e^{j\omega})H_0(e^{j\omega})+X(e^{j(\omega-\pi)})H_0(e^{j(\omega-\pi)})\bigg]G_0(e^{j\omega})\\
&\quad+\frac{1}{2}\bigg[X(e^{j\omega})H_1(e^{j\omega})+X(e^{j(\omega-\pi)})H_1(e^{j(\omega-\pi)})\bigg]G_1(e^{j\omega})\\
&=\frac{1}{2}\bigg[G_0(e^{j\omega})H_0(e^{j\omega})+G_1(e^{j\omega})H_1(e^{j\omega})\bigg]X(e^{j\omega})\\ &\quad+\frac{1}{2}\bigg[G_0(e^{j\omega})H_0(e^{j(\omega-\pi)})+G_1(e^{j\omega})H_1(e^{j(\omega-\pi)}) \bigg]X(e^{j(\omega-\pi)})
\end{align*}$
对于上述式子，如果是理想滤波器，分频的低通滤波器$H_0(e^{j\omega})$与重构的低通滤波器$G_0(e^{j\omega})$的取值范围是仅在$|\omega|<\frac{\pi}{2}$处不为0，因此$G_0(e^{j\omega})H_0^{e^{j(\omega-\pi)}}$的值为0；同理，理想滤波器的情况下，$G_1(e^{j\omega})H_1(e^{j(\omega-\pi)})$的值为0，所以上述式子的第二项在理想滤波器的情况下应为0。另外，第二项中的频谱包含$\omega-\pi$相关项，而$\omega-\pi$的相关项是由于信号的压缩带来的，因此这一项表征了压缩操作中潜在的混叠失真，而上述式子的第一项表征的是式子在理想情况下的分解与重构。
一种能消除混叠的重构系统
现实中不可能完全实现上述理想滤波器，不过我们也可以通过选取满足下述条件的滤波器来消除混叠：
$G_0(e^{j\omega})H_0(e^{j(\omega-\pi)})+G_1(e^{j\omega})H_1(e^{j(\omega-\pi)})=0$
该条件被称为混叠抵消条件，满足该式子的一组条件为
$\color{red}{\begin{align*}
h_1[n] = e^{j\pi n}h_0[n]&\Leftrightarrow H_1(e^{j\omega}) = H_0(e^{j(\omega-\pi)})\\
g_0[n] = 2h_0[n]&\Leftrightarrow G_0(e^{j\omega}) = 2H_0(e^{j\omega})\\
g_1[n] = -2h_1[n]&\Leftrightarrow G_1(e^{j\omega}) = -2H_0(e^{j(\omega-\pi)})
\end{align*}}$
把这些式子代入频谱函数得到
$Y(e^{j\omega}) = \Big[H_0^2(e^{j\omega})-H_0^2(e^{j(\omega-\pi)})\Big]X(e^{j\omega})$
因此完美重构（可能带M个样本延迟）要求
$\color{red}{H_0^2(e^{j\omega})-H_0^2(e^{j(\omega-\pi)}) = e^{-j\omega M}}$
也就是说，在设计分解/合成系统的时候，如果能找到一个符合上面式子的滤波器$H_0(e^{j\omega})$，就能使用前面的式子来设计$H_1(e^{j\omega}),G_0(e^{j\omega}),G_1(e^{j\omega})$，得到一个完美重构系统。不过这只是其中的一种分解/合成系统，这类系统的频谱并不能很好地对源信号进行分频，变换过程与前面描述的理想重构系统并不一样（书中稍有涉猎）。
对于这里所描述的这种系统，可以通过多相分解来大幅降低计算量。
因为系统的压缩以及扩展的倍率为2，因此把$h_0[n]$分成两个序列分量$e_{00}[n],e_{01}[n]$，把$g_0[n]$分成两个序列分量$f_{00}[n],f_{01}[n]$，
结合前面的条件，有
$\begin{align*}
e_{00}[n] &=h_0[2n]\\
e_{01}[n] &=h_0[2n+1]\\
e_{10}[n] &=h_1[2n]=e^{j2\pi n}h_0[2n]=e_{00}[n]\\
e_{11}[n] &=h_1[2n+1]=e^{j(2n+1)\pi}h_0[2n+1]=-e_{01}[n]\\
f_{00}[n] &=g_0[2n]=2h_0[2n]=2e_{00}[n]\\
f_{01}[n] &=g_0[2n+1]=2h_0[2n+1]=2e_{01}[n]\\
f_{10}[n] &=g_1[2n]=-2h_1[2n]=-2e_{00}[n]=-f_{00}[n]\\
f_{11}[n] &=g_1[2n+1]=-2h_1[2n+1]=2e_{01}[n]=f_{01}[n]
\end{align*}$
因此该系统可以优化为
