第十章:线性时不变系统

这章中,我们使用音乐作为例子来介绍一些信号与系统的理论知识, 解释卷积定理的证明,以及线性时不变系统的特性(马上我们就会看到线性时不变系统的定义)。

这章的代码 chap41.ipynb 可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp10 查看。

10.1 信号与系统

在信号处理领域, 系统(system) 代表的是接受一个信号为输入并产生一个信号输出的一种抽象概念, 如电子放大器,它是一个把输入信号放大后输出的电路。再比如,在你听音乐的时候,整个房间也可以看做是一个系统, 它将声音从它产生的地方(输入),然后传到你的耳朵里(输出)形成与原来不完全相同的声音。

线性时不变系统(linear, time-invariant system, LTI) 就是有如下两个性质的一类系统:

  1. 线性:如果你同时输入两个不同的信号,那么得到的输出和这两个信号各自的输出之和相同。 用数学语言描述为,如果输入 x1 产生输出 y1 ,输入 x2 产生输出 y2,那么 输入 ax1+bx2 产生的输出为 ay1+by2 。这里 ab 为常数。
  2. 时不变性:系统的作用效果不随时间的变化而变化,仅依赖于系统本身的状态。因此,如果输入变化仅是在时间上平移, 那么产生的输出变化也仅仅是相同的时间平移,其他都是一样的。

很多物理系统都具有这两个性质,或是近似的具有这两个性质。

  • 仅包含电阻,电容和电感的电路系统是LTI。这里我们只把这些元件考虑为理想模型。
  • 包含弹簧,块和阻尼装置的机械系统同样是LTI。这里假定弹簧和阻尼器是线性的(力与位移成正比,力与速度成正比)
  • 本书中的讲的最多的声音的传输介质(包括空气,水和固体等)也可以使用LTI作为近似模型。

LTI可以由线性微信方程描述,而它的解可以用复指数函数来表示, 见 http://en.wikipedia.org/wiki/Linear_differential_equation

我们可以运用这个结论来计算输入信号作用与一个LTI后的输出:

  1. 将输入信号表示为不同频率的复指数函数之和。
  2. 对每个分量计算相应的输出
  3. 将所有的输出加起来就得到了最终的输出

这看上去又很熟悉,它与 8.6 计算卷积使用的算法以及 9.3 计算微分的算法 是一个道理。这种方法我们成为 频谱分解 ,因为我们将输入信号“分解”成了频谱中的各个成分然后分别进行处理。

为了运用这个方法,我们需要找到每个频率成分作用于系统后产生的效果,以此来描述这个系统。对于,机械系统来说,有一个简单 有效的方法来完成这件事,就是触发一下并记录输出结果。

更准确的说,我们把触发一下这个操作叫做对系统输入一个 脉冲(impulse) , 而产生的输出我们称为 脉冲响应(impulse response) 。也就是说,对于线性系统,脉冲响应可以完全的表征系统的特性, 我们可以从脉冲的DFT结果看出其中的原因。下面我们生成了一个脉冲信号:

impulse = np.zeros(8)
impulse[0] = 1
impulse_spectrum = np.fft.fft(impulse)

这个信号的值是:

[ 1.  0.  0.  0.  0.  0.  0.  0.]

而它的频谱的值是:

[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j  1.+0.j  1.+0.j  1.+0.j  1.+0.j]

频谱中所有频率成分都是1。也就是说脉冲信号是包含所有频率成分并且它们具有相同幅值的信号。 注意不要与白噪声的频谱混淆了,白噪声的所有频率成分只是平均功率一样。

基于上面的分析,由于系统是线性时不变的,当我们对系统输入一个脉冲信号时其实相当于输入了所有的频率成分也就得到了 所有频率成分的响应。

译者注

从时域上看,任何信号均可以分解为不同时移脉冲信号乘以幅值之和,由于LTI的叠加性,任何信号的输出都可以由脉冲响应来 出,实际上输出就是输入信号和脉冲响应的卷积。

从频域上看,复指数信号是LTI的特性信号,即输入复指数信号后的输出是同频率的复指数信号,只有幅值相位的改变。因此只要 知道了系统对所有频率的复指数信号的响应就可以完全的确定系统的特性。而脉冲信号正是一个包括所有频率的信号。

10.2 窗和滤波器

为了展示使用脉冲响应来描述系统是可行的,我先从一个简单例子开始:一个包含两个元素的移动平均窗。我们可以把它当作是一个 对输入信号进行平滑的系统。

这个例子中,我们事先知道了系统使用的窗函数,因此我们可以直接计算它的滤波响应,但是实际情况,我们通常事先是不知道窗函数或滤波器 的形式的,之后我们会分析这种情况。

下面的代码生成了这个移动平均窗,并计算了它的DPF得到了滤波响应:

window_array = np.array([0.5, 0.5, 0, 0, 0, 0, 0, 0,])
window = thinkdsp.Wave(window_array, framerate=8)
filtr = window.make_spectrum(full=True)

结果如 图10.1 ,它是一个低通滤波器,滤波响应形状近似高斯曲线。

DFT of a 2-element moving average window.

图10.1: 两个元素的移动平均窗的DFT结果

现在,假定我们不知道窗函数或者滤波响应,我们就可以通过输入脉冲信号得到的脉冲响应,来描述这个系统。

这个例子中,我们将脉冲信号的频谱乘以滤波响应来得到了脉冲响应的频谱,然后由频谱生成了波形:

product = impulse_spectrum * filtr
filtered = product.make_wave()

由于 impulse_spectrum 的值都是1,相乘之后结果与滤波器是相同的,同样,最后得到的波形和窗函数也是一样的。

这个例子说明两点:

  • 由于脉冲信号的频谱都是1,因此脉冲响应的DFT与系统的滤波响应是相同的。
  • 脉冲响应与描述系统的卷积窗也是相同的。

10.3 声学响应

由上一节的结论我们知道,为了描述一个房间或开放空间的声学响应特性,可以简单的生成一个声音脉冲,可以用气球爆炸或开枪的声音 作为近似,然后得到的声音就是脉冲响应的近似。

我们以枪声作为例子生成脉冲响应来买描述一个房间的声学特性,然后使用它来模拟这个房间中录制的小提琴声音。

你可以在 代码库 中的 chap10.ipynb 上(或http://tinyurl.com/thinkdsp10)运行这个例子,然后听一听。

下面的代码是枪声:

response = thinkdsp.read_wave('180961__kleeb__gunshots.wav')
response = response.segment(start=0.26, duration=5.0)
response.normalize()
response.plot()

为了去除开枪前的空数据,我选择了从0.26s开始的一段。 图10.2 左图展示了枪声的波形。

Waveform of a gunshot.

图10.2: 枪声波形

接下来,我们计算 response 的DFT:

transfer = response.make_spectrum()
transfer.plot()

结果如 图10.2 右图。这个频谱就代表了这个房间的声学响应特性。频谱中每个频率成分都包含了一个表示幅值和初始相位的复数, 由于它包含了系统从输入到输出的所有信息,因此又被称为 传递函数(transfer function)

现在,我们可以模拟在这个房间内的小提琴声音,下面的代码是我们在 1.1 中使用的小提琴录音:

violin = thinkdsp.read_wave('92002__jcveliz__violin-origional.wav')
violin.truncate(len(response))
violin.normalize()

它和枪声的采样率都是44100Hz,它们的长度也大致相同,我调整了小提琴波形的长度使它们的长度一致。

然后,计算出小提琴波形的DFT:

spectrum = violin.make_spectrum()

现在,我们知道了信号的每个频率成分的幅值和初始相位,也知道了系统的传递函数,它们的乘积就是最终输入信号的DFT, 以此我们就可以反过来算出输出信号的波形:

output = (spectrum * transfer).make_wave()
output.normalize()
output.plot()
The waveform of the violin recording before and after convolution

图10.3: 经过系统前后的小提琴波形图

图10.3 中上图为输入信号波形,下图为输出信号波形。它们的波形是不一样的,并且你可以听出这个差别。 可以在 chap10.ipynb 中听一听,从中也许你可以感觉出这个房间的样子,对于我们来说,它听起来像是一个 狭长的房间,有硬质的地板和顶,像是一个靶场。

其实在这个例子中,为了避免带来一些困惑,我忽略了一件事,就是我们使用的原始的小提琴声音也是在某个房间录制的, 也就是说它已经经过了一次系统的传递,所以我们计算出的声音其实是经过了两个系统传递后的声音。如果想要真正的 模拟在另一个房间中的声音,应该先对录音进行一次逆向的传递函数处理。

10.4 系统和卷积

上一节中对于脉冲响应和传递函数的理解是这样的:

  • 脉冲信号是有幅值均为1的所有频率组成的
  • 脉冲响应是系统对所有频率的响应之和
  • 脉冲响应的DFT,即传递函数,它表征了系统对所有频率的响应。
  • 任何信号经过系统的输出,均可以表示为输入信号的所有频率成分经过系统后响应之和。

这种理解也许对你来说有点头疼,我们再来看看另一种理解:卷积。 由卷积定理,频域的乘积等于时域的卷积,可以得到,系统的输出等于输入与脉冲响应的卷积。

这个理解方式的关键在于:

  • 输入信号的采样值可以表示为一系列幅值不同时移不同的脉冲信号
  • 由于系统是时不变的,每个脉冲信号都会产生相应的不同幅值不同时移的脉冲响应
  • 由于系统是线性的,系统的输出应该等于这些不同的脉冲响应之和

我们一步步的来演示这个计算过程:首先假设我们开了两枪, t=0 时刻开了一枪,幅值为1, t=1 时刻开了一枪,幅值为0.5 。

我们可以通过将原来的脉冲响应经过时移和缩放后相加来得到最后的输出,下面这个函数可以计算出 信号经过时移和缩放后的波形:

def shifted_scaled(wave, shift, factor):
    res = wave.copy()
    res.shift(shift)
    res.scale(factor)
    return res

其中 shift 表示时移的秒数, factor 为缩放的因子。

下面我们使用它来计算输出:

shift = 1
factor = 0.5
gun2 = response + shifted_scaled(response, shift, factor)

结果如 图10.4 ,你可以在 chap10.ipynb 中听听它的声音。 当然,它听上去就是连续的两个枪声,前面的一声要大些,后面一声要小些。

Sum of a wave and a shifted, scaled copy

图10.4: 波形经过时移和缩放后的和

现在我们用100个枪声来替代之前的两个枪声,并且他们以每秒441Hz的速度开枪,那么 输出可以像下面这样来计算:

dt = 1 / 441
total = 0
for k in range(100):
    total += shifted_scaled(response, k*dt, 1.0)

由于每秒中有441次枪声,因此这次你已经分辨不出单独的一声枪响了,你感觉到的声音像是441Hz的周期信号。 如果你播放这个声音,你会发现它听起来像是在车库里按喇叭的声音。

这说明了我们可以把一个波形理解为一系列的不同幅值的脉冲。 同样,我们以一个441Hz的锯齿波作为例子:

signal = thinkdsp.SawtoothSignal(freq=441)
wave = signal.make_wave(duration=0.1,
                        framerate=response.framerate)

现在,我计算出了组成这个信号的一系列的脉冲产生的脉冲响应之和:

total = 0
for t, y in zip(wave.ts, wave.ys):
    total += shifted_scaled(response, t, y)

结果应该听起来像是在靶场中播放这个锯齿波的声音。你可以在 chap10.ipynb 中播放它。

Diagram of the sum of scaled and shifted copies of g

图10.5: g 时移和缩放后求和的过程

图10.5 展示了整个计算的过程,其中 f 是锯齿波, g 是脉冲响应, h 是计算结果。

例如:

\[h[2] = f[0]g[2] + f[1]g[1] + f[2]g[0]\]

更一般的来说:

h[n] = sumlimits_{m = 0}^{N - 1} {f[m]g[n - m]}

8.2 我们就见过这个式子了,它就是 fg 的卷积。 说明如果系统的脉冲响应为 g ,那么系统输入 f 产生的输出应该为 fg 的卷积。

概况起来,我们有两种形式来描述系统对输入信号的影响效果:

  1. 把输入理解为一系列的脉冲信号,那么输出就是脉冲响应进行时移和缩放后的和, 也就是输入与脉冲响应的卷积。
  2. 脉冲响应的DFT是系统的传递函数,它包含了系统对每个频率成分的影响效果,而输入可以理解为 不同频率的分量之和,因此将输入的DFT乘以传递函数就可以得到输出的DFT。

这两种描述无疑是等价的,这是由卷积定理所得到的:时域的卷积等价于频域的乘积。

这里我们也可以懂得为什么卷积的形式中 g 是反向的,这在我们学习平滑的时候提到过。原因就是 卷积的定义是在研究LTI的响应的时候自然得到的。

10.5 卷积定理的证明

现在我们是时候来解释卷积定理的证明过程了:

\[DFT(f*g) = DFT(f)DFT(g)\]

式中 fg 是长度均为 N 的两个向量。

证明过程分为两步:

  1. 我会先从 f 为复指数信号这个特例开始,说明它与 g 的卷积相当于对 f 乘以一个标量
  2. 然后我们把 f 推广到一般的信号,它可以表示为不同频率的复指数信号之和,然后通过乘以一个标量 来计算各个频率成分的卷积,再把结果加起来

为了证明卷积定理,我们先来从一些基本的式子开始。首先, g 的DFT,写作G,等于:

\[DFT(g)[k] = G[k] = \sum\limits_n {g(n){e^{ - 2\pi ink/N}}}\]

上式中, k 表示 0~N-1 的频率成分, n 表示 0~N-1 的采样时间。 这个结果是包含 N 个复数的向量。

f 的IDFT 写作 Ff 等于:

\[IDFT(F)[n] = f[n] = \sum\limits_k {F[k]{e^{2\pi ink/N}}}\]

卷积的定义为:

\[(f*g)[n] = \sum\limits_m {f[m]g[n - m]}\]

式中 m 的取值范围也为 0~N-1 ,由于卷积是满足交换律的,因此上式等于:

\[(f*g)[n] = \sum\limits_m {f[n - m]g[m]}\]

现在我们考虑 f 为复指数的特殊情况,我们将频率为 k 的复指数写作 \({e_k}\) ,那么:

\[f[n] = {e_k}[n] = {e^{2\pi ink/N}}\]

这里 k 表示频率, n 表示时间。 将 \({e_k}\) 代入之前卷积的第二个定义中得到:

\[({e_k}*g)[n] = \sum\limits_m {{e^{2\pi i(n - m)k/N}}g[m]}\]

把上式分解后得到:

\[({e_k}*g)[n] = {e^{2\pi ink/N}}\sum\limits_m {{e^{ - 2\pi imk/N}}g[m]}\]

可以看出,这个式子的第一项为 \({e_k}\) 本身, 第二项为 G[k] ,因此又可以写成:

\[({e_k}*g)[n] = {e_k}[n]G[k]\]

这个式子表明,复指数 \({e_k}\)g 的卷积等于 \({e_k}\) 乘以 G[k] 。 用数学语言来描述就是,与 g 进行卷积这个算子的特征函数是 \({e_k}\) , 对应的特征值为 G[k] (见 9.3

现在我们进行证明的第二步。当输入信号 f 不是复指数信号的时候,可以通过DFT将它表示为多个 复指数之和的形式,写作 F[k] 。其中 k0~N-1 的频率值, F[k] 就是包含了不同的频率 的复指数的幅值和初始相位。

根据之前的证明,对于输入信号的每个频率成分 F[k] ,输出信号的每个频率成分的系数应该为 F[k]G[k] 。 由于系统是线性的,输出就可以表示为:

\[(f*g)[n] = \sum\limits_k {F[k]G[k]{e_k}[n]}\]

代入 \({e_k}\) 的定义得到:

\[(f*g)[n] = \sum\limits_k {F[k]G[k]{e^{2\pi ink/N}}}\]

上式的右面部分实际上就是 FG 的乘积 FG 的IDFT,因此:

\[f*g = IDFT(FG)\]

由于 F=DFT(f) 以及 G=DFT(g) ,因此:

\[f*g = IDFT(DFT(f)DFT(g))\]

最后就得到了卷积定理的公式:

\[DFT(f*g) = DFT(f)DFT(g)\]

证明完毕。

10.6 练习

下面练习的答案可以参考文件 chap10soln.ipynb

练习110.4 中我把卷积描述为信号经过时移和缩放后的和。 但是,在 10.3 中我们把信号的DFT乘以传递函数,这个操作实际对应的 应该是 循环卷积 ,因为DFT假定信号是周期的。其结果会导致输出波形的前面有一些额外的影响。

幸运的是,这个问题有一种标准的解决方案,在信号的尾部加入足够多的0,然后再进行DFT,就可以避免 这个问题。

chap10.ipynb 中的例子进行补0后再计算,以消除这个问题对结果的影响。

练习2 OpenAIR ( http://www.openairlib.net )中提供了很多在线的声学脉冲响应数据。 浏览并下载一个你感兴趣的脉冲响应数据。找一段与你下载的脉冲响应同采样率的录音, 模拟这个录音通过这个脉冲响应的系统后产生的输出。使用两种方式来计算:一是通过输入与脉冲响应的卷积; 二是通过输入的DFT与脉冲响应的滤波响应相乘。