【音频处理】如何“认识”一个滤波器？

Introduction

对一段音频信号进行处理，例如添加一个音效，有时候我们要考虑这样的处理对音频的影响是什么？是提高了低频，还是提高高频？或者频率做周期变化？

当然，你完全可以用耳朵来听，但是这么做难度太大了，毕竟并不是每个人都有一副金耳朵。即使你的耳朵非常好使，能够分辨最为细微的差别，但是你如何用人类能够理解的方式描述这种区别呢？毕竟并不是每个人都有莫言般的描述能力。

综上，我们需要一种可以量化的方式来描述”处理“的过程。毋庸置疑，我们只能借助数学来描述这种复杂过程。这篇文章中，我们将以多个视角来观察”处理“过程。

一些知识铺垫

下面介绍一些基础的信号处理知识，并不会很详细，都是点到为止，目的是将整体内容串起来，形成自我闭环。更多深入的介绍，大家有兴趣的话找本 dsp 的书来看吧。

信号表示

首先，如何表示一个音频信号（我们只讨论离散信号）？很简单，我们用一个数组来描述它，一个长度为 N 的信号，可以表示为：
$$x[0],x[1],...,x[N-1]\text{\hspace{1em}(1)}$$
我们称 x[n] 为离散时间序列

单位脉冲信号

它的数学表示是这样的：
$$\delta (n)=\{\begin{array}{ll}1& \text{n = 0}\\ 0& \text{otherwise}\end{array}\text{\hspace{1em}(2)}$$
也就是长这样：$\delta (n)=[1,0,0,0,...,0]$

线性时不变滤波器

所谓滤波器，其实就是对信号”操作“：一个信号，经过滤波器，变成了另一个信号。我们这里关注的是线性时不变(Linearity and time-invariance; LTI)滤波器。LTI滤波器应用广泛，是滤波器家族中最为壮大的一簇。我们可以差分方程来表示一个 LTI 滤波器：
$$y[n]=b_{0}x[n]+b_{1}x[n-1]+...+b_{n}x[n-N]-a_{1}y[n-1]-...a_{M}y[n-M]\text{\hspace{1em}(3)}$$
其中 x 为输入序列，y为输出序列。$b_0,...,b_N$和$a_0,...,a_N$ 为系数。

脉冲响应

将单位脉冲信号输入到一个 LTI 滤波器中得到的结果，我们称之为脉冲响应，通常表示为$h(n)$

卷积

信号$x(n)$经过 LTI 滤波器后得到 $y(n)$，这个过程可以表示为脉冲响应$h(n)$与$x(n)$的卷积：
$$y(n)=(h\ast x)(n)\text{\hspace{1em}(4)}$$

Z 变换

接下来我们介绍下 Z 变换，它非常重要。

Z变换可将离散时间序列变换为在复频域的表达式。额…好吧，前面那句话是百度百科抄的，对于那些从没有接触过 Z 变换的同学而言，这句话简直了。我们还是通过具体的例子，从感性的角度来理解下它。

Z变换的公式很简单的，如下：
$$X(z)=\sum _{n=0}^{\infty }x[n]z^{-n}\text{\hspace{1em}(5)}$$

举个例子，离散时间序列为：x[0]=1, x[1]=2, x[2]=3，那么：
$$X(z)=1+2z^{-1}+3z^{-2}=1+2z^{-1}+3(z^{-1}{)}^2$$

通过 Z变换，我们将一组序列（N个值）转换为一个关于z的表达式(一个值)。
Z变换可以看成是线性变换，就比如上面的例子，可以改写为两个向量相乘：
$$\left[\begin{array}{ccc}1& z^{-1}& z^{-2}\end{array}\right]\left[\begin{array}{c}1\\ 2\\ 3\end{array}\right]=\left[\begin{array}{c}1+2z^{-1}+3(z^{-1}{)}^2\end{array}\right]$$

关于 Z变换有一个重要的性质，卷积。两个信号在时域上做卷积，等于在 z变换域中相乘
$$x_0[n]\ast x_1[n]=X_0(z)X_1(z)\text{\hspace{1em}(6)}$$

在前面提到，信号经过滤波器可以表示为卷积过程(4)，那么对$y(n)$进行Z变化得到：
$$Y(z)=H(z)X(z)\text{\hspace{1em}(7)}$$
从这个公式的角度来看，$Y(z)$其实是$X(z)$进行了一次线性变换得到的，其中 $H(z)$ 被我们成为 传递函数(Transfer function; TF)
$$H(z)=\frac{Y(z)}{X(z)}\text{\hspace{1em}(8)}$$
上面的公式同时表明，脉冲响应$h(n)$的Z变化 等于 转换方程 $H(z)=Y(z)/X(z)$。这也是求 $h(n)$ 的一种方法，通过计算 $H(z)$，然后做 z变换的逆变换，求得 $h(n)$

关于 Z 变换还有一些性质我们需要掌握：
Z { x 1 [ n ] + x 2 [ n ] } = Z { x 1 [ n ] } + Z { x 2 [ n ] } = X 1 ( z ) + X 2 ( z ) (9) Z\{x_1[n] + x_2[n]\} = Z\{x_1[n]\} + Z\{x_2[n]\} = X_1(z) + X_2(z) \tag{9} Z{x1​[n]+x2​[n]}=Z{x1​[n]}+Z{x2​[n]}=X1​(z)+X2​(z)(9)

Z { a x [ n ] } = a Z { x [ n ] } = a X ( z ) (10) Z\{ax[n]\} = aZ\{x[n]\} = aX(z) \tag{10} Z{ax[n]}=aZ{x[n]}=aX(z)(10)

Z { x [ n − 1 ] } = ∑ n = − ∞ ∞ x [ n − 1 ] z − n = z − 1 ∑ n = − ∞ ∞ x [ n − 1 ] z − ( n − 1 ) = z − 1 Z { x [ n ] } = z − 1 X ( z ) (11) \begin{aligned} Z\{x[n-1]\} &= \sum_{n=-\infty}^{\infty} x[n-1]z^{-n} \\ &= z^{-1}\sum_{n=-\infty}^{\infty}x[n-1]z^{-(n-1)} \\ &= z^{-1}Z\{x[n]\} = z^{-1}X(z) \end{aligned} \tag{11} Z{x[n−1]}​=n=−∞∑∞​x[n−1]z−n=z−1n=−∞∑∞​x[n−1]z−(n−1)=z−1Z{x[n]}=z−1X(z)​(11)

因此对于一个 LTI 滤波器(3)，我们对其进行 Z变换得到：
$$Y(z)=b_{0}X(z)+b_1{z}^{-1}X(z)+...+b_n{z}^{-N}X(z)-a_1{z}^{-1}Y(Z)-...a_M{z}^{-M}Y(z)\text{\hspace{1em}(12)}$$

合并同类项得:
$$(1+a_1{z}^{-1}+\cdots +a_M{z}^{-M})Y(z)=(b_0+b_1{z}^{-1}+\cdots +b_n{z}^{-N})X(z)$$

如果 $M\ne N$，我们添加一些 0 系数进去，让左右两边的个数一直。

可得传递方程 $H(z)$:
$$H(z)=\frac{Y(z)}{X(z)}=\frac{b_0+b_1{z}^{-1}+\cdots +b_n{z}^{-N}}{1+a_1{z}^{-1}+\cdots +a_N{z}^{-N}}\text{\hspace{1em}(13)}$$

多角度描述一个滤波器

差分方程 Difference equation

如公式(3)，这就是差分方程。让我们举两个例子来感受下它。

$$y[n]=x[n]+0.5y[n-1]\text{\hspace{1em}(14)}$$
和
$$y[n]=x[n]+2y[n-1]\text{\hspace{1em}(15)}$$

将一个脉冲信号输入至 (14) 中，可得脉冲响应$h(n)=1,0.5,0.25,\dots$，$h(n)$的值逐渐接近 0，我们这样的滤波器是稳定的。

将一个脉冲信号输入至 (15) 中，可得脉冲响应$h(n)=1,2,4,8\dots$，$h(n)$的值逐渐无穷大,我们称这样的滤波器是不稳定的。

块状图 Block diagram

有时候用图形的形式来展示滤波器更加符合人类的认知过程，以公式(3)为例，它的块砖图如下：

好吧，它看起比较复杂，我们看一个简单的例子，延迟信号：
$$y[n]=x[n]+gx[n-N]$$
它的块状图如下：

脉冲响应 Impulse response

正如前面提到的那样，将一个脉冲信号（公式(3)）输入至滤波器中，然后观察输出，就能得到脉冲响应。

例如有个滤波器为：
$$y(n)=x(n)+0.5^{3}x(n-3)-0.9^{5}y(n-5)\text{\hspace{1em}(16)}$$

你可以用下面这段 MATLAB 代码来画出该滤波器的脉冲响应

g1 = 0.5^3;  B = [1 0 0 g1];      % Feedforward coeffs
g2 = 0.9^5;  A = [1 0 0 0 0 g2];  % Feedback coefficients

h = filter(B,A,[1,zeros(1,50)]);  % Impulse response
% h = impz(B,A,50); % alternative in octave-forge or MSPTB

% Matlab-compatible plot:
clf; figure(1); stem([0:50],h,'-k'); axis([0 50 -0.8 1.1]);
ylabel('Amplitude'); xlabel('Time (samples)'); grid;

当然，我们还有更为精确的计算 $h(n)$ 的方法：对$h(n)$ 进行 z变换可以得到传递函数 $H(z)$，因此对 $H(z)$ 做逆 z变化得到 $h(n)$。

$H(z)$可以通过公式(13)得到，相当的容易。逆 z变换我还没学会，所以没法和大家解释了。但是我们可以来验证$h(n)$的z变换就是公式(13)的结果。

对公式(16)做z变换，得到传递方程：
$$H(z)=\frac{Y(z)}{X(z)}=\frac{1+0.5^3{z}^{-3}}{1+0.9^5{z}^{-5}}$$

取 $z=2$，可得 $H(z)=\frac{1+0.5^3{2}^{-3}}{1+0.9^5{2}^{-5}}=0.9972$

然后我们在代码中计算 H(z) 的值，也是得到 0.9972，完全一致。

g1 = 0.5^3;  B = [1 0 0 g1];      % Feedforward coeffs
g2 = 0.9^5;  A = [1 0 0 0 0 g2];  % Feedback coefficients

h = filter(B,A,[1,zeros(1,50)]);  % Impulse response
z_value = 2;                      % z = 2
z = zeros(1,length(h)) + z_value; 
Z = z.^-(0:length(h)-1);
HZ = h*Z';                        %计算 H(z)，结果为 0.9972

传递方程 Transfer function

传递方程已经在上面说过了，对 $h(n)$ 做z变换可以得到，也可用公式(13)得到。在分析 LTI 滤波器时，传递方程是一种重要的方法。

极点、零点及增益 Poles zeros gain

我们对公式(13)上下同时乘上 $z^N$，得到一组多项式，并因式分解
$$\begin{array}{rl}H(z)& =\frac{b_0{z}^N+b_1{z}^{N-1}+\cdots +b_n}{z^N+a_1{z}^{N-1}+\cdots +a_N}\\ & =g\frac{(z-q_1)(z-q_2)\dots (z-q_N)}{(z-p_1)(z-p_2)\dots (z-p_N)}\end{array}\text{\hspace{1em}(17)}$$

我们称 $g$ 为 增益(Gain)

我们称 $p_1,\dots ,p_N$ 为极点(Poles)，当 $z\in p$，$H(z)=\infty$

我们称 $q_1,\dots ,q_N$ 为 零点(Zeros)，当 $z\in q$，$H(z)=0$

我们将所有极点和零点画在一个复平面上，如果所有极点都在单位圆内，那么滤波器是稳定的；如果所有零点在单位圆内，那么我们称这个滤波器是最小相位滤波器。

频谱响应 Frequency Response

频谱响应指的是，信号经过滤波器后，在频域发生的变化，定义为：输出信号的频谱除以输入信号的频谱

一个信号经过傅里叶变化(DTFT)便可以得到其频谱信息：
$$DTF{T}_{\omega (x)}=\sum _{n=0}^{\infty }x(n)e^{-j\omega Tn}$$

你看DTFT是不是很像 z变换？公式(5)中 $z=e^{j\omega T}$，那么就完美对应上了 DTFT。
又有 $Y(z)=H(z)X(z)$，因此得到：
$$Y(e^{j\omega T})=H(e^{j\omega T})X(e^{j\omega T})$$

那么频谱响应其实就是
$$H(e^{j\omega T})=\frac{Y(e^{j\omega T})}{X(e^{j\omega T})}$$

其中 $T=f/f_s$，其中$f$为信号频率其范围在 $(0,f_s/2)$，因此 $\omega T\in (0,\pi )$

频谱响应又分 幅度响应（Amplitude response） 和 相位响应（Phase response）

幅度响应其实就是 $|H(e^{j\omega T})|$，取个绝对值就ok了。幅度响应能够说明滤波器对每个频率幅度的影响。

相位响应其实就是 $\angle H(e^{j\omega T})$，它说明了滤波器对每个频率相位的影响

实际的例子

千言万语不如一个实际的例子，假设我们有这样一个滤波器：
$$y[n]=4x[n]-4x[n-1]+x[n-2]-y[n-1]-y[n-2]/2$$

做z变换的：
$$Y(z)=4X(z)-4X(z)z^{-1}+X(z)z^{-2}-Y(z)z^{-1}-Y(z)z^{-2}/2$$
计算传递方程：
$$H(z)=\frac{Y(z)}{X(z)}=\frac{4z^2-4z+1}{z^2+z+1/2}$$

做因式分解
$$H(z)=\frac{4(z-1/2)(z-1/2)}{(z-\frac{-1+j}{2})(z-\frac{-1-j}{2})}$$

得到两个极点: $(-1+j)$ 和 $(-1-j)$

得到两个零点，它们都在 $\frac{1}{2}$，将它们画在复平面上，可以发现所有极点都在单位圆内，因此这个滤波器是稳定的；所有零点都在单位圆内，因此它是最小相位的

当$z=e^{j\omega T}$时，且$\omega T\in (0,\pi )$，计算其幅度响应和相位响应，直接上代码计算：

w = 0:0.001:pi;
E = exp(1j*w);

Hz = (4*(E.^2) - 4*E + 1) ./ ((E.^2) + E + 0.5);

amplitude_Hz = abs(Hz);
phase_Hz = angle(Hz);

figure(1);
plot(20*log10(amplitude_Hz));

figure(2);
plot(phase_Hz*180/pi);

其结果大概长这样：

从图可以看出，这个滤波器抑制信号中的低频信息，增益高频信息。在 $0.785\pi$（在 44.1kHz 采样率下为 17.3kHz）增益最大。
对于特别高和特别低的频率信息，这个滤波器在相位上只有很小的影响。但是在 $0.5\pi$（在 44.1kHz 采样率下为 11.025kHz）处，其相位信息偏移了 116度

总结

本文介绍了多种方式来认识一个滤波器，

• 差分方程
• 块状图
• 脉冲响应
• 传递方程
• 极点、零点和增益
• 频谱响应

并列举相关例子，让大家有更为深入的理解和体会。

转载请注明出处 https://blog.csdn.net/weiwei9363/article/details/104950007