采样率变换和多速率filter

前言

多速率滤波器技术大约诞生于20世纪70年代，由于其能降低采样率，进而降低计算复杂度，传输速率和存储空间，在信息压缩和通信领域得到了广泛的应用。距【4】介绍，最早的应用是80年代的数字音乐领域，通过过采样+降采样的方法令反混叠滤波器的指标要求大大降低，从而节约了成本。然后在48k vs 44.1k的采样率互换问题上，多速率滤波器也体现了巨大的运算优越性。而令多速率滤波技术走上神坛的则是正交滤波器组和子带分解，子带分解给音视频压缩提供了强大的引擎，同时也在宽带通信、宽带语音信号处理领域得到了大量的应用。说到这里的第一个问题是子带分解和DFT啥关系，仍然引用【4】中的一副图，非常形象的找出了二者的异同。

据说在QMF的帮助下，临频混叠问题相对于DFT会得到完美的解决。

sampling和aliase

采样和混叠是数字信号处理永远也绕不过去的话题，引自【5】的一张混叠示意图，当采样率不满足Nyquist定理，混叠会在频域表现出这样的状态。

上图所表示的带限信号截止频率为$B$，如果$f_s/2>B$，如（b）的情况，经过采样后的频域延拓信息没有发生混叠，$f_s/2<B$，如（c）的情况，经过采样后的频域延拓信息发生混叠，从理论上无法恢复信号。根据傅里叶变换理论，而对于任何一个带限内的频率分量，我们都可以视为一个正弦波（频率是连续的），表示为$x_a(t)=Acos(2\pi ft+\theta )$，如果以$f_s=1/T_s$的采样率进行采样，则会得到：$$x_a(nT)\equiv x(n)=Acos(2\pi fnT+\theta )=Acos(2\pi fn/f_s+\theta )$$
时间分量被$n$指数分量代替，模拟到数字的鸿沟就此被搭桥。而根据上面的推导，可以得出一个归一化频率$F=f/f_s$，由于奈奎斯特定理限制，$-1/2\le F\le 1/2$。如果根据上式，令角频率$\Omega =2\pi f$，则有数字角频率$\omega =2\pi f/f_s=2\pi F$。由于奈奎斯特定理限制，$-\pi \le \omega \le \pi$。

decimation vs lpf

decimation

最最简单的抽取就是隔一抽取，很显然这样就会令采样率变成原来的一半，根据上一节的图例，原本清晰奈奎斯特关系（b）变成了混叠的情况（c），但有的时候非抽取不可呢？答案就是先上一个滤波器，截止频率比新采样率的一半还低，问题就解决了。一般地，假设$M$为整数抽取器的取值，如果输入序列为$x(n)$，为了避免混淆，此处定义$$x_1=\{\begin{array}{ll}x(n),& n=0,\pm M,\pm 2M,...\\ 0,& \text{其他}\end{array}$$那么输出序列：$$y(n)=x(Mn)=x_1(Mn)$$表示每到第$M$个进行抽取一个采样点，采样率就变成原来的$1/M$。$z$变换表示$$\begin{array}{rl}Y(z)=\sum _{n=-\infty }^{\infty }y(n)z^{-n}& =\sum _{n=-\infty }^{\infty }x(Mn)z^{-n}\\ & =\sum _{n=-\infty }^{\infty }{x}_1(Mn)z^{-n}\\ & =\sum _{n=-\infty }^{\infty }{x}_1(n)z^{-n/M}\\ & =X_1(z^{1/M})\end{array}$$考虑到(渊源欠奉)$$p(n)=\frac{1}{M}\sum _{k=0}^{M-1}{e}^{j2\pi kn/M}=\frac{1}{M}\sum _{k=0}^{M-1}{W}^{-kn},W=e^{-j2\pi /M}$$的性质可以将$x_1(n)$表示为$$x_1(n)=x(n)p(n)$$这样$$\begin{array}{rl}{X}_1(z)& =\sum _{n=-\infty }^{\infty }x(n)p(n)z^{-n}\\ & =\sum _{n=-\infty }^{\infty }x(n)[\frac{1}{M}\sum _{k=0}^{M-1}{e}^{j2\pi kn/M}]z^{-n}\\ & =\frac{1}{M}\sum _{k=0}^{M-1}\sum _{n=-\infty }^{\infty }x(n)e^{j2\pi kn/M}{z}^{-n}\\ & =\frac{1}{M}\sum _{k=0}^{M-1}\sum _{n=-\infty }^{\infty }x(n)[e^{-j2\pi k/M}z{]}^{-n}\\ & =\frac{1}{M}\sum _{k=0}^{M-1}\sum _{n=-\infty }^{\infty }x(n)[W^{k}z{]}^{-n}\\ & =\frac{1}{M}\sum _{k=0}^{M-1}X(z{W}^k)\end{array}$$
$$\begin{array}{rl}Y(z)& =X_1(z^{1/M})\\ & =\frac{1}{M}\sum _{k=0}^{M-1}X(z^{1/M}{W}^k)=\frac{1}{M}\sum _{k=0}^{M-1}X(e^{j(w-2\pi k)/M})\end{array}$$
这个公式阐述了几个特性：

1. 拉伸了$X(e^{jw})$在频率轴的宽度，假设M=2，则在归一化频率轴上向外扩大2倍；
2. 拉伸的响应函数基础上，建立$M-1$个copy，每个copy平移$2\pi$，那么如果频率不能小于$\pi /M$，下一步平移相加后就产生了混叠；
3. 将1.和2.的结果相加，然后除以$M$；

对于带限信号$X(e^{jw})$，可以利用下面的图【4】来理解，即无限延拓频谱（离散信号）再次沿着频率轴正向延拓$M$次，得到了抽取后的频谱（密集恐惧）。这种情况很容易产生混叠，所以需要增加一个滤波器，对抽取的信息进行低通滤波。

lpf

假设增加一个理想低通滤波器【5】$$H_M(w)=\{\begin{array}{ll}1,& |w|\le \pi /M\\ 0,& \text{其他}\end{array}$$经过滤波器的信号$$v(n)=\sum _{k=0}^{\infty }h(n)x(n-k)$$$$y(m)=v(mM)=\sum _{k=0}^{\infty }h(n)x(mM-k)$$结合前文中的p(n)得到$$\hat{v}(n))=v(n)p(n)$$这样$$y(m)=\hat{v}(mM)=v(mM)p(mM)=v(mM)$$$z$变换$$Y(z)=\sum _{m=-\infty }^{\infty }y(m)z^{-m}=\sum _{m=-\infty }^{\infty }\hat{v}(mM)z^{-m}=\sum _{m=-\infty }^{\infty }\hat{v}(m)z^{-m/M}$$
利用p(n)的变换$$\begin{array}{rl}Y(z)& =\sum _{m=-\infty }^{\infty }v(m)[\frac{1}{M}\sum _{k=0}^{M-1}{e}^{j2\pi mk/M}]z^{-m/M}\\ & =\frac{1}{M}\sum _{k=0}^{M-1}\sum _{m=-\infty }^{\infty }v(m)[e^{-j2\pi mk/M}{z}^{1/M}{]}^{-m}\\ & =\frac{1}{M}\sum _{k=0}^{M-1}V(e^{-j2\pi mk/M}{z}^{1/M})\\ & =\frac{1}{M}\sum _{k=0}^{M-1}{H}_M(e^{-j2\pi mk/M}{z}^{1/M})X(e^{-j2\pi mk/M}{z}^{1/M})\end{array}$$相应的单位圆上的傅里叶变换$$Y(jw)=\frac{1}{M}\sum _{k=0}^{M-1}{H}_M(e^{j(w-2\pi k)/M})X(e^{j(w-2\pi k)/M})$$由于开始滤波器的假设条件可知，除了首项其他都为0，即$$Y(jw)=\frac{1}{M}{H}_M(e^{j(w)/M})X(e^{j(w)/M})$$至此抽取的时频关系大抵如此。

interpolation vs lpf

而插值的情况稍微复杂一些， 假设原始采样信号为$x(n)$，在两个连续值之间插入$I-1$个零得到新的序列$y(m)$，可以表示为$$y(m)=\{\begin{array}{ll}x(m/I),& 0,\pm I,\pm 2I,...\\ 0,& \text{其他}\end{array}$$
$z$变换如下$$Y(z)=\sum _{n=-\infty }^{\infty }y(n)z^{-n}=\sum _{k=-\infty }^{\infty }y(kI)z^{-kI}=\sum _{k=-\infty }^{\infty }x(k)z^{-kI}=X(z^I)$$

利用DTFT变换$$y(m)={\int }_{-\pi }^{\pi }Y(\omega )e^{j\omega }\mathrm{d}\omega$$
$$\begin{gathered} Y(\omega )=\sum _{m=-\infty }^{\infty }v(m)e^{-j\omega m} \\ =\sum _{m=-\infty }^{\infty }x(m)e^{-j\omega mI}=X(Iw) \end{gathered}$$而通常用$e^{jw}$表示$$Y(e^{jw})=X(e^{jwI})$$
根据上述推导，再回顾假设原来$x(n)$的归一化角频率为${\omega }_x=2\pi f/f_x$，新的采样率是原来的$I$倍，${\omega }_y=2\pi f/(f_x\ast I)={\omega }_x/I$，新采样后的在归一化的频谱上被压缩了，从上图的©变成了(b)，但问题是不仅如此，镜像频率也进入到区间（如果$I$=2,相当于原来$f_s$处的频响变成$f_s/2$）,从频域的角度，仍然需要一个低通滤波器，在插值完毕，滤掉$B/I$之外的频率响应（镜像），那么再恢复回时域就完成了采样率的提高。

多相滤波器结构

基本推导

任何系统函数的$Z$变换都可以如下表示[8]$$\begin{array}{rlrl}H(z)=& ...+h(0)& +h(M)z^{-M}& +...\\ & ...+h(1)z^{-1}& +h(M+1)z^{-(M+1)}& +...\\ & ...+h(2)z^{-2}& +h(M+2)z^{-(M+2)}& +...\\ & ...& & \\ & ...+h(M-1)z^{-(M-1)}& +h(2M-1)z^{-(2M-1)}& +...\end{array}$$
对第$i$行提取因子$z^{-(i-1)}$,$$\begin{array}{rlrl}H(z)=& ...+h(0)& +h(M)z^{-M}& +...\\ & ...+z^{-1}[h(1)& +h(M+1)z^{-M}& +...]\\ & ...+z^{-2}[h(2)& +h(M+2)z^{-M}& +...]\\ & ...& & \\ & ...+z^{-(M-1)}[h(M-1)& +h(2M-1)z^{-M)}& +...]\end{array}$$令$$E_i(z)=\sum _{n=-\infty }^{\infty }h(nM+i)z^{-n}$$系统函数可以改写为$$H(z)=\sum _{i=0}^{M-1}{z}^{-i}{E}_i(z^M)$$
网络拓扑结构如下【4】：

假设$M=3$，系统输入的$Z$变换为$X(z)$，系统输出的$Z$变换为$Y(z)$，结合上式可表达为： Y ( z ) = H ( z ) X ( z ) = P 0 ( z 3 ) X ( z ) + z − 1 { P 1 ( z 3 ) X ( z ) + z − 1 [ P 0 ( z 3 ) X ( z ) ] } \begin{aligned} Y(z)&=H(z)X(z)\\ &=P_0(z^3)X(z)+z^{-1}\{P_1(z^3)X(z)+z^{-1}[P_0(z^3)X(z)]\} \end{aligned} Y(z)​=H(z)X(z)=P0​(z3)X(z)+z−1{P1​(z3)X(z)+z−1[P0​(z3)X(z)]}​
上面这个式子表示的的是频域上已经分为三个滤波器了，而且由于$z^3$的关系，每个滤波器实质上已经是隔三抽取的降采样版本。

何为多相

polyphase两个词根表示多个相位，相位又从何谈起呢？这个主要是从傅里叶变换的频率轴开始的。下面翻译一下【4】中一个DFT Filter Bank的例子，帮助理解多相和滤波器组的概念。回顾一下这个矩阵的样子。

$$\begin{gathered} \mathbf{W}=\left[\begin{array}{ccccc}{e}^{-j0}& e^{-j0}& e^{-j0}& \cdots & e^{-j0}\\ e^{-j0}& e^{-j2\pi /N}& e^{-j4\pi /N}& \cdots & e^{-j2(N-1)\pi /N}\\ e^{-j0}& e^{-j4\pi /N}& e^{-j8\pi /N}& \cdots & e^{-j4(N-1)\pi /N}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ e^{-j0}& e^{-j2(N-1)\pi /N}& e^{-j4(N-1)\pi /N}& \cdots & e^{-j2(N-1)(N-1)\pi /N}\end{array}\right] \\ =\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& e^{-j2\pi /N}& e^{-j4\pi /N}& \cdots & e^{-j2(N-1)\pi /N}\\ 1& e^{-j4\pi /N}& e^{-j8\pi /N}& \cdots & e^{-j4(N-1)\pi /N}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& e^{-j2(N-1)\pi /N}& e^{-j4(N-1)\pi /N}& \cdots & e^{-j2(N-1)(N-1)\pi /N}\end{array}\right] \end{gathered}$$而$${\mathbf{W}}^{\ast }=\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& e^{j2\pi /N}& e^{j4\pi /N}& \cdots & e^{j2(N-1)\pi /N}\\ 1& e^{j4\pi /N}& e^{j8\pi /N}& \cdots & e^{j4(N-1)\pi /N}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& e^{j2(N-1)\pi /N}& e^{j4(N-1)\pi /N}& \cdots & e^{j2(N-1)(N-1)\pi /N}\end{array}\right]$$
首先利用DFT 矩阵$W^{\ast }$，构建序列$x(n)$的延时链，令$s_i(n)=x(n-i)$，摘图

另外，由于$W=e^{-j2\pi /M}$，则令$$x_k(n)=\sum _{i=0}^{M-1}{s}_i(n)W^{-ki}=\sum _{i=0}^{M-1}{s}_i(n)e^{ki\ast j2\pi /M}$$这里要先跳出stft分帧的思想，每个时间指数$n$，上面的框架图会求取$M$个信号$x_k(n)$，而输入$s_i(n)$则是信号$x(n)$的$i$个延时（长度仍然为$M$），而上面的公式（共轭）是典型的IDFT关系（除了没有定标因子$1/M$）。这个框图右侧信号的$Z$变换有什么特点呢？$$X_k(z)=\sum _{i=0}^{M-1}{S}_i(z)W^{-ki}=\sum _{i=0}^{M-1}{z}^{-i}{W}^{-ki}X(z)=\sum _{i=0}^{M-1}(z{W}^k{)}^{-i}X(z)$$联想一下$X_k(z)=H_k(z)X(z)$，那么$H_k(z)={\sum }_{i=0}^{M-1}(z{W}^k{)}^{-i}=H_0(z{W}^k)$，这里的$$H_0(z)=1+z^{-1}+...+z^{-(M-1)}$$过程有点绕，但结果已经很清晰了，上述等价一个分析滤波器组（analysis bank）$H_k(z)$。这个滤波器频域幅度相应可以单用sinc函数来表达$|H_0(e^{jw})|=|sin(Mw/2)/sin(w/2)|$，而第$k$个滤波器的相应为：$$H_k(e^{jw})=H_0(e^{jw}{e}^{-j2k\pi /M})=H_0(e^{j(w-2k\pi /M)})$$这相当于$H_0(e^{jw})$移（相）位版本。同样借用文中的图例如下

上面这个也被称为归一化的DFT 滤波器组（Uniform DFT banks）。这个例子主要是为了理解所以多相的概念，即多个移动相位的滤波器组成的一个组合，但这个例子并不是多相展开，切忌混淆。

参考文档

1.从“零”开始学习一下DCT
2.The Filter Bank Summation (FBS) Interpretation of the Short Time Fourier Transform (STFT)
3.Multirate Filter Banks from CCRMA
4.Multirate Systems and Filter Banks, P.Vaidyanathan
5.Understanding Digital Signal Processing 3rd Edition
6.多相滤波
7.多相滤波器
8.数字信号处理 第四版
9.一窥数字傅里叶变换矩阵