Introduction
对一段音频信号进行处理,例如添加一个音效,有时候我们要考虑这样的处理对音频的影响是什么?是提高了低频,还是提高高频?或者频率做周期变化?
当然,你完全可以用耳朵来听,但是这么做难度太大了,毕竟并不是每个人都有一副金耳朵。即使你的耳朵非常好使,能够分辨最为细微的差别,但是你如何用人类能够理解的方式描述这种区别呢?毕竟并不是每个人都有莫言般的描述能力。
综上,我们需要一种可以量化的方式来描述”处理“的过程。毋庸置疑,我们只能借助数学来描述这种复杂过程。这篇文章中,我们将以多个视角来观察”处理“过程。
一些知识铺垫
下面介绍一些基础的信号处理知识,并不会很详细,都是点到为止,目的是将整体内容串起来,形成自我闭环。更多深入的介绍,大家有兴趣的话找本 dsp 的书来看吧。
信号表示
首先,如何表示一个音频信号(我们只讨论离散信号)?很简单,我们用一个数组来描述它,一个长度为 N 的信号,可以表示为:
x
[
0
]
,
x
[
1
]
,
.
.
.
,
x
[
N
−
1
]
(1)
x[0],x[1],...,x[N-1] \tag{1}
x[0],x[1],...,x[N−1](1)
我们称 x[n] 为离散时间序列
单位脉冲信号
它的数学表示是这样的:
δ
(
n
)
=
{
1
n = 0
0
otherwise
(2)
\delta(n)= \begin{cases} 1& \text{n = 0} \\ 0& \text{otherwise} \end{cases} \tag{2}
δ(n)={10n = 0otherwise(2)
也就是长这样:
δ
(
n
)
=
[
1
,
0
,
0
,
0
,
.
.
.
,
0
]
\delta(n) = [1,0,0,0,...,0]
δ(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
]
(3)
y[n] = b_0x[n] + b_1x[n-1] + ... + b_nx[n-N]- a_1y[n-1] - ... a_My[n-M] \tag{3}
y[n]=b0x[n]+b1x[n−1]+...+bnx[n−N]−a1y[n−1]−...aMy[n−M](3)
其中 x 为输入序列,y为输出序列。
b
0
,
.
.
.
,
b
N
b_0,...,b_N
b0,...,bN和
a
0
,
.
.
.
,
a
N
a_0,...,a_N
a0,...,aN 为系数。
脉冲响应
将单位脉冲信号输入到一个 LTI 滤波器中得到的结果,我们称之为脉冲响应,通常表示为 h ( n ) h(n) h(n)
卷积
信号
x
(
n
)
x(n)
x(n)经过 LTI 滤波器后得到
y
(
n
)
y(n)
y(n),这个过程可以表示为脉冲响应
h
(
n
)
h(n)
h(n)与
x
(
n
)
x(n)
x(n)的卷积:
y
(
n
)
=
(
h
∗
x
)
(
n
)
(4)
y(n) = (h*x)(n) \tag{4}
y(n)=(h∗x)(n)(4)
Z 变换
接下来我们介绍下 Z 变换,它非常重要。
Z变换可将离散时间序列变换为在复频域的表达式。额…好吧,前面那句话是百度百科抄的,对于那些从没有接触过 Z 变换的同学而言,这句话简直了。我们还是通过具体的例子,从感性的角度来理解下它。
Z变换的公式很简单的,如下:
X
(
z
)
=
∑
n
=
0
∞
x
[
n
]
z
−
n
(5)
X(z) = \sum_{n=0}^{\infty}x[n]z^{-n} \tag{5}
X(z)=n=0∑∞x[n]z−n(5)
举个例子,离散时间序列为:x[0]=1, x[1]=2, x[2]=3,那么:
X
(
z
)
=
1
+
2
z
−
1
+
3
z
−
2
=
1
+
2
z
−
1
+
3
(
z
−
1
)
2
X(z) = 1 + 2z^{-1} + 3z^{-2} = 1 + 2z^{-1} + 3(z^{-1})^2
X(z)=1+2z−1+3z−2=1+2z−1+3(z−1)2
通过 Z变换,我们将一组序列(N个值)转换为一个关于z的表达式(一个值)。
Z变换可以看成是线性变换,就比如上面的例子,可以改写为两个向量相乘:
[
1
z
−
1
z
−
2
]
[
1
2
3
]
=
[
1
+
2
z
−
1
+
3
(
z
−
1
)
2
]
\begin{bmatrix} 1 & z^{-1} & z^{-2} \end{bmatrix} \begin{bmatrix} 1\\ 2\\ 3 \end{bmatrix} = \begin{bmatrix} 1 + 2z^{-1} + 3(z^{-1})^2 \end{bmatrix}
[1z−1z−2]⎣⎡123⎦⎤=[1+2z−1+3(z−1)2]
关于 Z变换有一个重要的性质,卷积。两个信号在时域上做卷积,等于在 z变换域中相乘
x
0
[
n
]
∗
x
1
[
n
]
=
X
0
(
z
)
X
1
(
z
)
(6)
x_0[n]*x_1[n] = X_0(z)X_1(z) \tag{6}
x0[n]∗x1[n]=X0(z)X1(z)(6)
在前面提到,信号经过滤波器可以表示为卷积过程(4),那么对
y
(
n
)
y(n)
y(n)进行Z变化得到:
Y
(
z
)
=
H
(
z
)
X
(
z
)
(7)
Y(z)=H(z)X(z) \tag{7}
Y(z)=H(z)X(z)(7)
从这个公式的角度来看,
Y
(
z
)
Y(z)
Y(z)其实是
X
(
z
)
X(z)
X(z)进行了一次线性变换得到的,其中
H
(
z
)
H(z)
H(z) 被我们成为 传递函数(Transfer function; TF)
H
(
z
)
=
Y
(
z
)
X
(
z
)
(8)
H(z)=\frac{Y(z)}{X(z)} \tag{8}
H(z)=X(z)Y(z)(8)
上面的公式同时表明,脉冲响应
h
(
n
)
h(n)
h(n)的Z变化 等于 转换方程
H
(
z
)
=
Y
(
z
)
/
X
(
z
)
H(z)=Y(z)/X(z)
H(z)=Y(z)/X(z)。这也是求
h
(
n
)
h(n)
h(n) 的一种方法,通过计算
H
(
z
)
H(z)
H(z),然后做 z变换的逆变换,求得
h
(
n
)
h(n)
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
)
(12)
Y(z) = b_0X(z) + b_1z^{-1}X(z) + ... + b_nz^{-N}X(z)- a_1z^{-1}Y(Z) - ... a_Mz^{-M}Y(z) \tag{12}
Y(z)=b0X(z)+b1z−1X(z)+...+bnz−NX(z)−a1z−1Y(Z)−...aMz−MY(z)(12)
合并同类项得:
(
1
+
a
1
z
−
1
+
⋯
+
a
M
z
−
M
)
Y
(
z
)
=
(
b
0
+
b
1
z
−
1
+
⋯
+
b
n
z
−
N
)
X
(
z
)
(1+a_1z^{-1}+\dots+a_Mz^{-M})Y(z) = (b_0+b_1z^{-1}+\dots+b_nz^{-N})X(z)
(1+a1z−1+⋯+aMz−M)Y(z)=(b0+b1z−1+⋯+bnz−N)X(z)
如果 M ≠ N M\ne N M=N,我们添加一些 0 系数进去,让左右两边的个数一直。
可得传递方程
H
(
z
)
H(z)
H(z):
H
(
z
)
=
Y
(
z
)
X
(
z
)
=
b
0
+
b
1
z
−
1
+
⋯
+
b
n
z
−
N
1
+
a
1
z
−
1
+
⋯
+
a
N
z
−
N
(13)
H(z) = \frac{Y(z)}{X(z)} = \frac{b_0+b_1z^{-1}+\dots+b_nz^{-N}}{1+a_1z^{-1}+\dots+a_Nz^{-N}} \tag{13}
H(z)=X(z)Y(z)=1+a1z−1+⋯+aNz−Nb0+b1z−1+⋯+bnz−N(13)
多角度描述一个滤波器
差分方程 Difference equation
如公式(3),这就是差分方程。让我们举两个例子来感受下它。
y
[
n
]
=
x
[
n
]
+
0.5
y
[
n
−
1
]
(14)
y[n] = x[n] + 0.5y[n-1] \tag{14}
y[n]=x[n]+0.5y[n−1](14)
和
y
[
n
]
=
x
[
n
]
+
2
y
[
n
−
1
]
(15)
y[n] = x[n] + 2y[n-1] \tag{15}
y[n]=x[n]+2y[n−1](15)
将一个脉冲信号输入至 (14) 中,可得脉冲响应 h ( n ) = 1 , 0.5 , 0.25 , … h(n) = 1,0.5,0.25,\dots h(n)=1,0.5,0.25,…, h ( n ) h(n) h(n)的值逐渐接近 0,我们这样的滤波器是稳定的。
将一个脉冲信号输入至 (15) 中,可得脉冲响应 h ( n ) = 1 , 2 , 4 , 8 … h(n) = 1,2,4,8\dots h(n)=1,2,4,8…, h ( n ) h(n) h(n)的值逐渐无穷大,我们称这样的滤波器是不稳定的。
块状图 Block diagram
有时候用图形的形式来展示滤波器更加符合人类的认知过程,以公式(3)为例,它的块砖图如下:
好吧,它看起比较复杂,我们看一个简单的例子,延迟信号:
y
[
n
]
=
x
[
n
]
+
g
x
[
n
−
N
]
y[n] = x[n] + gx[n-N]
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
)
(16)
y(n) = x(n) + 0.5^3x(n-3) - 0.9^5y(n-5) \tag{16}
y(n)=x(n)+0.53x(n−3)−0.95y(n−5)(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) h(n) 的方法:对 h ( n ) h(n) h(n) 进行 z变换可以得到传递函数 H ( z ) H(z) H(z),因此对 H ( z ) H(z) H(z) 做逆 z变化得到 h ( n ) h(n) h(n)。
H ( z ) H(z) H(z)可以通过公式(13)得到,相当的容易。逆 z变换我还没学会,所以没法和大家解释了。但是我们可以来验证 h ( n ) h(n) h(n)的z变换就是公式(13)的结果。
对公式(16)做z变换,得到传递方程:
H
(
z
)
=
Y
(
z
)
X
(
z
)
=
1
+
0.
5
3
z
−
3
1
+
0.
9
5
z
−
5
H(z) = \frac{Y(z)}{X(z)} = \frac{1+0.5^3z^{-3}}{1+0.9^5z^{-5}}
H(z)=X(z)Y(z)=1+0.95z−51+0.53z−3
取 z = 2 z=2 z=2,可得 H ( z ) = 1 + 0. 5 3 2 − 3 1 + 0. 9 5 2 − 5 = 0.9972 H(z) = \frac{1+0.5^32^{-3}}{1+0.9^52^{-5}} = 0.9972 H(z)=1+0.952−51+0.532−3=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 ) h(n) h(n) 做z变换可以得到,也可用公式(13)得到。在分析 LTI 滤波器时,传递方程是一种重要的方法。
极点、零点及增益 Poles zeros gain
我们对公式(13)上下同时乘上
z
N
z^N
zN,得到一组多项式,并因式分解
H
(
z
)
=
b
0
z
N
+
b
1
z
N
−
1
+
⋯
+
b
n
z
N
+
a
1
z
N
−
1
+
⋯
+
a
N
=
g
(
z
−
q
1
)
(
z
−
q
2
)
…
(
z
−
q
N
)
(
z
−
p
1
)
(
z
−
p
2
)
…
(
z
−
p
N
)
(17)
\begin{aligned} H(z) &= \frac{b_0z^{N}+b_1z^{N-1}+\dots+b_n}{z^N+a_1z^{N-1}+\dots+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{aligned} \tag{17}
H(z)=zN+a1zN−1+⋯+aNb0zN+b1zN−1+⋯+bn=g(z−p1)(z−p2)…(z−pN)(z−q1)(z−q2)…(z−qN)(17)
我们称 g g g 为 增益(Gain)
我们称 p 1 , … , p N p_1,\dots,p_N p1,…,pN 为极点(Poles),当 z ∈ p z\in{p} z∈p, H ( z ) = ∞ H(z) = \infty H(z)=∞
我们称 q 1 , … , q N q_1,\dots,q_N q1,…,qN 为 零点(Zeros),当 z ∈ q z\in{q} z∈q, H ( z ) = 0 H(z)=0 H(z)=0
我们将所有极点和零点画在一个复平面上,如果所有极点都在单位圆内,那么滤波器是稳定的;如果所有零点在单位圆内,那么我们称这个滤波器是最小相位滤波器。
频谱响应 Frequency Response
频谱响应指的是,信号经过滤波器后,在频域发生的变化,定义为:输出信号的频谱除以输入信号的频谱
一个信号经过傅里叶变化(DTFT)便可以得到其频谱信息:
D
T
F
T
ω
(
x
)
=
∑
n
=
0
∞
x
(
n
)
e
−
j
ω
T
n
DTFT_{\omega(x)} = \sum_{n=0}^{\infty}x(n)e^{-j\omega Tn}
DTFTω(x)=n=0∑∞x(n)e−jωTn
你看DTFT是不是很像 z变换?公式(5)中
z
=
e
j
ω
T
z=e^{j\omega T}
z=ejωT,那么就完美对应上了 DTFT。
又有
Y
(
z
)
=
H
(
z
)
X
(
z
)
Y(z)=H(z)X(z)
Y(z)=H(z)X(z),因此得到:
Y
(
e
j
ω
T
)
=
H
(
e
j
ω
T
)
X
(
e
j
ω
T
)
Y(e^{j\omega T}) = H(e^{j\omega T})X(e^{j\omega T})
Y(ejωT)=H(ejωT)X(ejωT)
那么频谱响应其实就是
H
(
e
j
ω
T
)
=
Y
(
e
j
ω
T
)
X
(
e
j
ω
T
)
H(e^{j\omega T}) = \frac{Y(e^{j\omega T}) }{X(e^{j\omega T})}
H(ejωT)=X(ejωT)Y(ejωT)
其中 T = f / f s T=f/f_s T=f/fs,其中 f f f为信号频率其范围在 ( 0 , f s / 2 ) (0, f_s/2) (0,fs/2),因此 ω T ∈ ( 0 , π ) \omega T \in (0, \pi) ωT∈(0,π)
频谱响应又分 幅度响应(Amplitude response) 和 相位响应(Phase response)
幅度响应其实就是 ∣ H ( e j ω T ) ∣ \vert H(e^{j\omega T}) \vert ∣H(ejωT)∣,取个绝对值就ok了。幅度响应能够说明滤波器对每个频率幅度的影响。
相位响应其实就是 ∠ H ( e j ω T ) \angle H(e^{j\omega T}) ∠H(ejωT),它说明了滤波器对每个频率相位的影响
实际的例子
千言万语不如一个实际的例子,假设我们有这样一个滤波器:
y
[
n
]
=
4
x
[
n
]
−
4
x
[
n
−
1
]
+
x
[
n
−
2
]
−
y
[
n
−
1
]
−
y
[
n
−
2
]
/
2
y[n] = 4x[n]-4x[n-1] + x[n-2] - y[n-1] - y[n-2]/2
y[n]=4x[n]−4x[n−1]+x[n−2]−y[n−1]−y[n−2]/2
做z变换的:
Y
(
z
)
=
4
X
(
z
)
−
4
X
(
z
)
z
−
1
+
X
(
z
)
z
−
2
−
Y
(
z
)
z
−
1
−
Y
(
z
)
z
−
2
/
2
Y(z) = 4X(z) - 4X(z)z^{-1} + X(z)z^{-2} - Y(z)z^{-1} - Y(z)z^{-2}/2
Y(z)=4X(z)−4X(z)z−1+X(z)z−2−Y(z)z−1−Y(z)z−2/2
计算传递方程:
H
(
z
)
=
Y
(
z
)
X
(
z
)
=
4
z
2
−
4
z
+
1
z
2
+
z
+
1
/
2
H(z) = \frac{Y(z)}{X(z)} = \frac{4z^2-4z+1}{z^2+z+1/2}
H(z)=X(z)Y(z)=z2+z+1/24z2−4z+1
做因式分解
H
(
z
)
=
4
(
z
−
1
/
2
)
(
z
−
1
/
2
)
(
z
−
−
1
+
j
2
)
(
z
−
−
1
−
j
2
)
H(z) = \frac{4(z-1/2)(z-1/2)}{(z-\frac{-1+j}{2})(z-\frac{-1-j}{2})}
H(z)=(z−2−1+j)(z−2−1−j)4(z−1/2)(z−1/2)
得到两个极点: ( − 1 + j ) (-1+j) (−1+j) 和 ( − 1 − j ) (-1-j) (−1−j)
得到两个零点,它们都在 1 2 \frac{1}{2} 21,将它们画在复平面上,可以发现所有极点都在单位圆内,因此这个滤波器是稳定的;所有零点都在单位圆内,因此它是最小相位的
当 z = e j ω T z=e^{j\omega T} z=ejωT时,且 ω T ∈ ( 0 , π ) \omega T\in (0, \pi) ωT∈(0,π),计算其幅度响应和相位响应,直接上代码计算:
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
π
0.785\pi
0.785π(在 44.1kHz 采样率下为 17.3kHz)增益最大。
对于特别高和特别低的频率信息,这个滤波器在相位上只有很小的影响。但是在
0.5
π
0.5\pi
0.5π(在 44.1kHz 采样率下为 11.025kHz)处,其相位信息偏移了 116度
总结
本文介绍了多种方式来认识一个滤波器,
- 差分方程
- 块状图
- 脉冲响应
- 传递方程
- 极点、零点和增益
- 频谱响应
并列举相关例子,让大家有更为深入的理解和体会。
转载请注明出处 https://blog.csdn.net/weiwei9363/article/details/104950007