《ThinkDSP》¶
前言¶
信号处理是我很感兴趣的一个学科,它被广泛的应用在很多领域。 如果理解了信号处理中的一些基本概念,可以帮助我们更好的理解我们生活中所看到 和所听到的东西。
但是,如果你不是学相关专业的,你很难学懂信号处理这个学科。 因为大多数的相关的书籍(包括学校使用的教材)都是自底向上来讲解的, 先从最基本的数学抽象开始讲起(例如向量)然后逐渐深入,通常缺少理论与实际的联系以及具体的应用示例。
阅读本书之前,你需要有些编程的基础,我们需要使用这项技能来完成一些有趣的东西。 使用编程式的讲解方式,可以很直观的展示大多数的概念,比如在第一章结束后,你就可以通过程序来 分析声音以及其他的信号,生成新的声音了。每一章中,我们都会介绍一项新的技术,以及它在处理真实 信号中的应用。我们会先学习怎样使用这个技术,然后再学习它的工作原理。
我相信这样的学习形式是更实用也更有趣的(希望你们是这样认为的)。
本书的读者¶
本书使用Python作为编程语言,因此需要读者有一些Python的编程基础并且对面向对象的程序设计有所了解。
如果你对Python还不太熟的话,你可以先看看我写的另一本书,ThinkPython,它是一本Python的入门书籍, 适用于那些没有编程基础的人。如果你有编程经验但是对Python不熟悉,那么你也可以看Mark Lutz的Learning Python。
本书中还会用到Numpy和Scipy两个Python的扩展库,如果你对此不熟悉也没有关系,我会在使用它们的函数和类的时候, 进行简单的解释。
本书的代码¶
本书的代码和书中使用到的声音文件都托管在Github上: https://github.com/AllenDowney/ThinkDSP 。 Git是一个版本管理工具,使用它你可以很好的管理和跟踪整个项目中的所有文件,被Git管理项目我们称之为代码库。 Github是一个托管代码库的网站,它提供了很方便的用户界面来管理你的代码库。
本书代码库的Github主页上提供了几种方式来使用其中的代码:
- 你可以把我的代码库Fork到自己的账号下。(如果你没有Github账户,你可以免费注册一个)。 Fork后你就有了一个自己的代码库了,你可以使用它来跟踪管理你在学习过程中编写的代码。 然后,你需要使用clone命令来把代码库中的代码复制到你的本地电脑中。
- 如果你没有Github账号,你也可以直接clone我的代码库。但是这样你就不能把你写的代码上传到Github上了,不过你 依然可以在本地跟踪和管理你的代码。
- 如果你不想使用Git,你还可以把代码库下载成一个zip文件。(Github页面的右下角有一个下载的按钮)
本书中的所有代码均可以工作与Python2和Python3。
我在开发这些代码的时候,使用了Continuum Analytics公司开发的Anaconda,它是一个免费的Python发行版,其中包含了 我们将要用到的大部分Python包。Anaconda的安装非常简单,默认情况下它可以在用户级权限下进行安装,不需要管理员权限。 它同时支持Python2和Python3。你可以在 http://continuum.io/downloads 下载Anaconda。
如果你不想使用Anaconda,你需要自己手动安装以下的Python包:
- Numpy:用于数值计算, http://www.numpy.org/
- Scipy: 用于科学计算, http://www.scipy.org/
- matplotlib:用于作图, http://matplotlib.org/
虽然这些程序包都是很常用的,但是它们不会包含在基本的Python的安装包中,因此你需要手动安装。 在某些环境下,它们有可能不太好安装,如果你在安装的时候遇到问题,我建议你使用Anaconda或其他的Python发行版, 一般这些发行版的Python都默认包含了这些程序包。
本书的练习有的需要使用Jupyter notebook,如果你对此不熟悉,可以参考 http://jupyter.org 。 使用Jupyter有三种方式:
在电脑中运行Jupyter
如果安装了Anaconda,其中默认包含了Jupyter,你可以使用命令行来启动Jupyter服务:
$ jupyter notebook
如果没有安装Jupyter,可以使用下面的命令行进行安装:
$ conda install jupyter
当你启动了Jupyter服务后,它会在浏览器中自动的打开一个新的页面也就是Jupyter notebook的主页。
在Binder中运行Jupyter
Binder是一个运行Jupyter的网络服务,通过 http://mybinder.org/repo/AllenDowney/ThinkDSP 你可以直接打开 本书的Jupyter主页,并且你可以在编辑和运行里面的代码,但是你编写和改动的代码不会被保存起来,如果你把页面关闭了或者 长时间没有操作(1个小时),那么这些代码会消失。
在nbviewer中查看
本书还会提供在noviewer中的链接,在nbviewer中你仅仅可以静态展示代码和结果。你可以通过nbviewer的链接来阅读本书的代码, 也可以播放示例中的声音,但是你不能进行改动也不能运行它们,交互式的控件也使用不了。
Good Luck,and have fun!
Contributor List¶
If you have a suggestion or correction, please send email to downey@allendowney.com. If I make a change based on your feedback, I will add you to the contributor list (unless you ask to be omitted).
If you include at least part of the sentence the error appears in, that makes it easy for me to search. Page and section numbers are fine, too, but not as easy to work with. Thanks!
- Before I started writing, my thoughts about this book benefited from conversations with Boulos Harb at Google and Aurelio Ramos, formerly at Harmonix Music Systems.
- During the Fall 2013 semester, Nathan Lintz and Ian Daniher worked with me on an independent study project and helped me with the first draft of this book.
- On Reddit’s DSP forum, the anonymous user RamjetSoundwave helped me fix a problem with my implementation of Brownian Noise. And andodli found a typo.
- In Spring 2015 I had the pleasure of teaching this material along with Prof. Oscar Mur-Miranda and Prof. Siddhartan Govindasamy. Both made many suggestions and corrections.
- Giuseppe Masetti sent a number of very helpful suggestions.
- Kim Cofer copyedited the entire book, found many errors, and made many helpful suggestions.
Other people who found typos and errors include Silas Gyger and Abe Raher.
Special thanks to the technical reviewers, Eric Peters, Bruce Levens, and John Vincent, for many helpful suggestions, clarifications, and corrections.
Also thanks to Freesound, which is the source of many of the sound samples I use in this book, and to the Freesound users who contributed those samples. I include some of their wave files in the GitHub repository for this book, using the original file names, so it should be easy to find their sources.
Unfortunately, most Freesound users don’t make their real names available, so I can only thank them by their user names. Samples used in this book were contributed by Freesound users: iluppai, wcfl10, thirsk, docquesting, kleeb, landup, zippi1, themusicalnomad, bcjordan, rockwehrmann, marcgascon7, jcveliz. Thank you all!
第一章:声音和信号¶
信号 代表一个随时间变化的量。这个定义很抽象,我们用声音信号来作为一个更具体的例子。 声音是空气压力的变化产生的,声音信号则代表着空气压力随时间的变化。
麦克风能够测量这种变化并转换成可以表示声音的电信号。扬声器可以接收这种电信号并产生相应的声音。 像麦克风和扬声器这类的设备统称为 信号变换器,他们能够将一种信号变换为另一种信号。
本书会介绍信号处理的基础内容,包括信号的合成,变换和分析等。本书主要关注声音信号,但是同样的方法 完全可以运用于处理其他领域的信号,比如电信号,机械振动等。
这些方法同样也可以运用于其他随空间而非时间变化的信号,比如沿运动路径的高度变化信号。他们也可以运用于处理 多维的信号,例如图像或电影。图像可以看做是随二维空间变化的信号,电影则可以看做是随二维空间和时间同时 变化的信号。
但是,我们会从最简单的一维信号 声音信号 开始。
这章的代码 chap01.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp01 查看。
1.1 周期信号¶
我们先从 周期信号 开始。周期信号是按一定的周期重复出现的信号。例如当我们敲钟的时候,它会振动并发出声音, 如果把这个声音信号记录下来,可以画出如 图1.1 。
这个信号的形状与正弦曲线很相似,这就意味着它可以用三角函数来大致描述出来。
可以看出这个信号是周期的,在这个图中,我截取了三个重复循环( cycles )的信号,每个循环所经历的时间, 我们通常称作周期( peroid ), 这个信号的周期大概是 2.3 ms。
信号的频率( frequence )是指信号每秒重复循环的次数,是周期的倒数。频率的单位是赫兹(Hz)。
这个信号的频率大约是439Hz,略小于国际标准音高440Hz,用音名A表示,更准确说应该是A4,其中数字后缀表示所属八度音程, A4与中央C属于同一个八度,A5则比A4高一个八度。如果不熟悉标准音高记法(scientific pitch notation)可以参考: http://en.wikipedia.org/wiki/Scientific_pitch_notation 。
调音用的音叉会产生一个正弦信号,因为它尖端的振动是一个简谐运动。大部分乐器所产生的声音都是周期的,但是这些信号的形状 并不是简单的正弦曲线。例如,图1.2 显示了小提琴演奏的鲍凯利尼(Boccherini)的第5号E大调弦乐五重奏的第三乐章中 截取的一部分录音。
从图中我们可以看出,这个信号是周期的,但是信号的形状要比正弦波要复杂一些。我们把信号的形状称为信号的波形( waveform )。 大多数乐器产生的波形都比正弦信号复杂,而不同的波形实际上决定了不同乐器的音色,也就影响着我们对声音的心理感受。相比于正弦波, 人们通常会感觉这些复杂的信号更加的丰富,更加的温暖或者更加的有趣。
1.2 频谱分解¶
本书中一个重要的主题就是频谱分解。信号可以进行频谱分解是由于信号可以表示为不同频率的正弦信号的叠加结果, 这些不同频率的正弦信号的集合就是 频谱 。信号经过 离散傅立叶变换 就可以得到它的频谱, 离散傅立叶变换(简称DFT),也是本书最重要的一个数学概念之一。
一般情况下,我们通过一个叫做 快速傅立叶变换(FFT) 的算法来计算信号的DFT,FFT也是本书最重要的算法之一。
例如,图1.3 显示了 图1.2 中小提琴声音的频谱。其中x轴代表了组成这个信号的正弦波频率的范围, y轴显示了它们的强度(也叫做幅度)。
我们把频谱中最小的频率称作 基频 。这个信号的基频大概是440Hz(稍低一些)。
这个信号基频的幅度也是最大的,所以它也是这个信号的 主导频率 。通常我们对音高的感受取决于声音的 基频(即使它不是主导频率)。
这个频谱中其他幅度比较大的频率是880Hz,1320Hz,1760Hz和2200Hz,它们都是基频的整数倍,这些频率分量 称为信号的 谐波(harmonics) 。之所以叫做谐波,是因为它们和基频一起出现我们会觉得声音听起来比较和谐。
译者注
对于声音信号来说,谐波的另一个名称是泛音。实际上,对于任何周期性的信号,进行傅里叶变换后都会得到谐波。 因此,在翻译的时候,统一将harmonics翻译为谐波。
- 880Hz是A5的频率,它比基频高一个 八度 。一个八度的音程在频率上是两倍的关系。
- 1320Hz大约是E6的频率,它比A5高一个纯五度。如果你不熟悉音程(如纯五度)的概念,可以参考 https://en.wikipedia.org/wiki/Interval_(music).
- 1760Hz是A6的频率,它比基频高两个八度。
- 2200Hz 大约是C#7,它比A6高一个大三度。
这些谐波的音如果在同一个八度中都是三和弦的组成部分。其中有些音的频率是近似的,这是由于组成音阶的音的频率在 平均律下有一些小的调整。详见 http://en.wikipedia.org/wiki/Equal_temperament。
给定这些谐波的频率以及相应的幅度,我们就可以通过直接叠加的方法重构出整个信号。接下来我们会学习这个方法。
1.3 信号¶
我编写了一个Python模块 thinkdsp.py
,其中包含了我们用于处理信号和频谱的类和函数。
你可以在本书的 代码库 中找到这些代码。
thinkdsp
提供一个 Signal
类来表示信号,它是其他信号类的父类,包括正(余)弦信号 Sinusoid
等。
使用 thinkdsp
可以像下面这样来构造一个正弦信号或是余弦信号:
cos_sig = thinkdsp.CosSignal(freq=440, amp=1.0, offset=0)
sin_sig = thinkdsp.SinSignal(freq=880, amp=0.5, offset=0)
freq
参数代表频率(Hz), amp 参数代表幅度(无量纲,1.0代表可以录制和播放的最大的幅度),
offset
参数代表初始相位(弧度),它决定了信号在 0 时刻的值。例如,一个 offset=0
的正弦信号
在 0 时刻的值为 sin(0)
也就是 0 ,而 offset=pi/2
的正弦信号在 0 时刻的值为 sin(pi/2)
也就是 1 。
Signal
类中定义了 __add__
方法,因此我们可以使用 +
来对两个信号进行叠加:
mix = sin_sig + cos_sig
两个信号相加的结果为 SumSingal
,他可以代表任意两个信号的叠加结果。
Singal
实际上是由一个数学函数所定义的,大部分 Signal
定义了时间从正无穷到负无穷上信号的取值。
我们可以在一个时间区间上对一个信号进行求值( evaluate ),也就是输入一个时间的序列 - ts
,计算信号在这些时间点上
相应的值 - ys
。我们使用 Numpy
的数组来表示 ts
和 ys
,并且将他们包装成一个叫做 Wave
的类中。
Wave
表示信号在某个时间序列上的值,其中每个时间点叫做帧( frame ),也可以叫做采样( sample ),这两种叫法经常
可以混用。
Signal
类提供了一个 make_wave
返回一个新的 Wave
对象:
wave = mix.make_wave(duration=0.5, start=0, framerate=11025)
其中, duration
参数代表了需要求值的时间区间的长度(单位为秒), start
参数代表开始时间(单位为秒),
framerate
是一个整数,表示每秒的帧数,也叫做 采样率 (单位为Hz)
这个例子中,11025Hz的采样率经常用于音频信号的采集中,如 WAV 和 MP3 。该例中我们求得了信号在 t=0
到 t=0.5
区间内
等间隔的5513个采样值,我们把两次采样之间的时间间隔称为时间步长( timestep ),为采样率的倒数,这里的时间步长为 1/11025
大约为 91 微妙
Wave
提供了 plot
方法来画出波形图(使用 pyplot
模块实现):
wave.plot()
pyplot.show()
pyplot
是 matplotlib
库的一部分,是很常用的作图模块,他被包含在大部分的Python发行版中,当然你也可以手动安装它:
pip install matplotlib
这个440Hz的信号在0.5s内有220个周期,因此上面的代码画出的图形看起来会像是一条很粗的实线。
我们可以使用 segment
方法来截取一个更小的时间范围的波形:
period = mix.period
segment = wave.segment(start=0, duration=period*3)
period
是信号的一个属性,他返回信号的周期值(单位为秒)。
start
和 duration
的单位也是秒。这段代码截取了 mix
信号一开始的三个周期,其结果 segment
也是一个 Wave
对象。
我们画出 segment
的波形图,如 图1.4 。这个信号包含了两个不同频率的成分,因此它的波形看起来会比音叉发出的正弦信号要复杂一些,
但是比小提琴发出的声音信号要简单一些。
1.4 读写波形数据¶
thinkdsp
提供 read_wave
函数从 WAV 文件中读取数据并返回一个 Wave
对象:
violin_wave = thinkdsp.read_wave('input.wav')
Wave
对象提供了 write_wave
方法将数据写入到 WAV 文件中:
wave.write(filename='output.wav')
你可以用任意的媒体播放器来播放这些 WAV 文件。 在 UNIX 系统中,我通常使用 aplay
,这是一个简单而稳定的播放器,
多数的Linux发行版中都包含这个程序
thinkdsp
也提供了一个直接播放声音的函数 play_wave
,它会在一个子进程中运行播放器来播放音频数据:
thinkdsp.play_wave(filename='output.wav', player='aplay')
上面的代码中使用了默认的播放器 aplay
,当然你也可以通过 player
来指定其他的播放器。
译者注
我没有在Windows下找到播放器来替代 aplay
,如果有读者知道请
1.5 频谱¶
Wave
中提供了 make_spectrum
来生成频谱 Spectrum
spectrum = wave.make_spectrum()
Spectrum
同样也提供了 plot
方法用于作图:
spectrum.plot()
thinkplot.show()
我在 thinkplot
模块中包装了一些常用的 pyplot
方法,这个模块也包含在本书的 代码库 中。
Spectrum
提供了三个方法来对频谱进行变化:
low_pass
会对频谱应用一个低通滤波器,也就是说,高于给定截止频率的分量会被衰减(幅度减小), 衰减的程度由factor
指定,通常为一个 [0, 1] 的数,默认为 0 (完全衰减)。high_pass
会对频谱应用一个高通滤波器,也就是说,低于给定截止频率的分量会被衰减。band_stop
会对频谱应用一个带通滤波器,也就是说,在给定截止频率区间以外的分量会被衰减。
以下的代码将频谱的600Hz以上的频率成分衰减了99%:
spectrum.low_pass(cutoff=600, factor=0.01)
低通滤波器去除了声音中的明亮的高频声音,使声音变得比较低沉。你可以通过将频谱转换为波形后来播放它:
wave = spectrum.make_wave()
wave.play('temp.wav')
play
方法会将波形数据写入文件并且进行播放。如果使用 Jupyter notebooks
,你可以用 make_audio
生成一个可以直接播放的网页音频部件。
1.6 波形对象¶
其实 thinkdsp.py
中并没有什么复杂的东西,它提供的大多数方法仅仅是对 Numpy
和 Scipy
的包装。
其中主要有三个类: Signal
, Wave
和 Spectrum
。
给定一个 Signal
可以生成一个 Wave
,
给定一个 Wave
可以生成一个 Spectrum
,反之亦然。 图1.5 展示了这些关系。
Wave
包含三个属性: ys
是包含信号值的Numpy数组; ts
是对应的时间数组; framerate
是采样率。
其中单位时间通常是秒,但是有些例子中也会使用其他的单位时间,例如天。
Wave
还包含三个只读属性: start
, end
和 duration
,
这些属性由 ts
所决定,改变ts后这些属性会相应的改变。
我们可以通过直接改变 ts
以及 ys
来改变波形,例如:
wave.ys *= 2
wave.ts += 1
第一行代码将信号放大了两倍,使其音量变的更大。第二行代码将波形右移了一个单位时间,使其声音晚一秒钟才开始。
Wave
也提供了很多方法来进行更常规的操作,例如以下两个变换与之前的代码效果一样:
wave.scale(2)
wave.shift(1)
这些方法的文档在 http://greenteapress.com/thinkdsp.html 中。
1.7 信号对象¶
Signal
是所有信号的父类,其中提供了信号的基础方法,如 make_wave
。子类信号通过继承 Signal
并重写
evaluate
方法来实现。 evaluate
方法用于计算信号在任意时刻的值。
例如, Sinusoid
子类的定义如下:
class Sinusoid(Signal):
def __init__(self, freq=440, amp=1.0, offset=0, func=np.sin):
Signal.__init__(self)
self.freq = freq
self.amp = amp
self.offset = offset
self.func = func
其中构造参数包括:
- freq:信号的频率(Hz)
- amp:信号的幅度,通常单位为1
- offset:信号的初始相位,单位为弧度
- func:用于计算给定时间点的信号值的函数。可以为
np.sin
或np.cos
,对应为正弦信号和余弦信号。
Singal
类中的 make_wave
方法的代码如下:
def make_wave(self, duration=1, start=0, framerate=11025):
n = round(duration * framerate)
ts = start + np.arange(n) / framerate
ys = self.evaluate(ts)
return Wave(ys, ts, framerate=framerate)
其中, start
duration
为开始时间和持续时间(单位为秒),framerate
是采样率(单位为Hz)。
n
是采样点的总数, ts
用Numpy数组表示的采样时间
make_wave
会调用 evaluate
方法来计算信号在每个采样点的值 ys
,
例如: Sinusoid
中的 evaluate
是这样的:
def evaluate(self, ts):
phases = PI2 * self.freq * ts + self.offset
ys = self.amp * self.func(phases)
return ys
让我们详细解释一下这个函数:
self.freq
是频率,ts
是采样时间序列,因此他们的乘积为采样的cycle
PI2
是常数 \(2\pi\) ,把cycle
与 \(2\pi\) 相乘 就得到了相位( phase )。我们如果将波形循环一周的长度视作360°,即 \(2\pi\) , 那么相位就是信号在一周内所处的位置。self.offset
是初始相位,也就是t=0
时刻信号的相位。它实际上代表了波形的左右平移。- 如果
self.func
是np.sin
或np.cos
则计算的值会在[-1,1]的范围内。 - 乘以
self.amp
使得最终的结果范围为[-self.amp, self.amp]。
evaluate
用数字公式表示为:
其中 \(A\) 是幅度,\(f\) 是频率,\(t\) 是时间,\({\varphi _0}\) 是相位。 看起来好像我们用了很多代码来描述了一个简单的公式,实际上,我们得到了一个通用的框架来描述 所有类似的信号,而不仅仅是正余弦信号。
译者注
在 thinkdsp.py
中除了正余弦信号外,有很多信号都继承自 Sinusoid
,
包括三角信号 TriangleSignal
,方波信号 SquareSignal
,
锯齿信号 SawtoothSignal
等。这些信号的特征是都具有频率,幅度和初始相位的属性,
而唯一的区别就是它们的 evaluate
方法。
1.8 练习¶
在开始下面的练习之前,你可以从本书的 代码库 中下载本书的源码。
下面练习的答案可以参考文件 chap01soln.ipynb
练习1 如果你安装了 Jupyter
,使用它来打开 chap1.ipynb
,
阅读并且运行上面的代码示例。 如果没有 Jupyter
,可以在
http://tinyurl.com/thinkdsp01 浏览和运行它。
练习2 在 http://freesound.org 上下载一段清楚的声音,可以是音乐,语音或其他的声音,
使用代码截取其中音高固定的半秒声音,并计算并画出这段声音频谱,观察一下这个声音的音色
和它的频谱之间有什么样的关系。然后使用 high_pass
, low_pass
, band_stop
来
滤除其中的一些谐波分量,把他们在反过来转换为波形对象并播放,听一听与原来的声音有什么区别。
练习3 产生一些正弦信号和余弦信号,并将他们相加合成一个复合信号。 然后生成并画出信号的波形以及频谱。播放这个声音听一听,看看如果频率分量不是基频的整数倍的时候, 声音是怎么样的。
练习4 编写一个 stretch
函数,接收一个 Wave
对象以及一个伸缩因子,通过改变 ts
和 framerate
来让波形变快或变慢。提示:仅需要两行代码就能实现这个功能。
第二章:谐波¶
在这章中我们会介绍几种新的波形以及他们的频谱, 理解频谱的谐波结构( harmonic structure ),也就是构成整个频谱的正弦信号的集合。
另外,还会介绍在数字信号处理中的另一个重要的概念:混叠( aliasing )。
然后我会解释一下 Spectrum
是如何工作的。
这章的代码 chap02.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp02 查看。
2.1 三角波¶
正弦信号仅包含一个频率分量,因此它的频谱只有一个峰值。大多数复杂的信号,如小提琴声,它们的DFT会包含多个峰值。 这一小节,我们来研究这些波形和他们的频谱之间的关系。
我们先从三角波开始,图2.1 展示了一个200Hz的三角波,它的波形看起来像是把正弦信号拉直了。
你可以使用 thinkdsp.TriangleSignal
来生成一个三角波:
class TriangleSignal(Sinusoid):
def evaluate(self, ts):
cycles = self.freq * ts + self.offset / PI2
frac, _ = np.modf(cycles)
ys = np.abs(frac - 0.5)
ys = normalize(unbias(ys), self.amp)
return ys
TriangleSignal
继承自 Sinusoid
,
因此它同样包含 freq
, amp
, offset
三个属性。
不同在于它复写了 evaluate
。其中 ts
依然是采样点的时间,我们来看看这个 evaluate
是如何产生
三角波的:
cycles
表示从采样点的循环数。np.modf
把它的小数部分提取出来放到frac
中,它的整数部分 被丢弃了 [1]frac
现在是一个给定频率的0~1变化的斜坡信号,将它减去0.5会使其范围变到-0.5~0.5。然后取绝对值后, 它的值就成了从0增加到0.5,再从0.5减小到0unbias
会把整个信号向下移动,使其相对于y轴居中,然后normalize
将信号的幅度放大到amp
。
产生 图2.1 的代码如下:
signal = thinkdsp.TriangleSignal(200)
signal.plot()
接下来,我们可以用这个信号来产生波形对象,然后再生成它的频谱:
wave = signal.make_wave(duration=0.5, framerate=10000)
spectrum = wave.make_spectrum()
spectrum.plot()
图2.2 显示了频谱图,右面的图在Y轴上进行了放大,这样可以更清晰的显示谐波结构。 像我们期望的那样,基频200Hz的幅值是最大的,其它谐波频率分量的峰值出现在200Hz的整数倍频率上。
有一个奇怪的现象是:谐波里面没有基频的偶数倍的频率(400Hz,800Hz等), 而只有奇数倍的频率(600Hz,1000Hz,1400Hz等)。
这些频率成分的另一个特性是,随着谐波频率的增加,幅度的减弱与频率的平方大致呈比例的关系。 例如,600Hz的谐波是基频200Hz的3倍,他的幅度和基频的比例大约是9倍的关系(3的平方)。 而1000Hz的谐波是600Hz的1.7倍左右,他们的幅度之比大概是 \({1.7^2} = 2.9\) 。 我们把这种关系就称为信号的 谐波结构 。
2.2 方波¶
thinkdsp
还提供了 SquareSignal
类来表示方波信号,这个类的定义如下:
def evaluate(self, ts):
cycles = self.freq * ts + self.offset / PI2
frac, _ = np.modf(cycles)
ys = self.amp * np.sign(unbias(frac))
return ys
类似 TriangleSignal
和 SquareSignal
这样继承自 Sinusoid
的类,
他们的共同点是都具有相同的构造参数:频率,幅度,初始相位。
SquareSignal
的 evaluate
方法也具有类似的结构。其中 ts
依然是采样时间序列,
frac
是它的小数部分,它的值从0到1周期的变化。
unbias
将 frac
调整到-0.5~0.5,然后 np.sign
将结果的负值映射到-1,正值映射到1。
最后,乘以 amp
将信号的幅值调整到 -amp~amp 。
图2.3 显示了100Hz方波信号的三个周期, 图2.4 显示了它的频谱。
和三角波一样,方波同样只包含奇数倍的谐波频率,它们的峰值在300Hz,500Hz,700Hz等。 但是它们幅度的减弱要比三角波慢一些(不是平方的关系)。
在本章后面的练习题中,你可以看到一些其他的波形和它们的谐波结构。
2.3 混叠¶
坦白说,之前介绍的几个波形,都是我刻意挑选的,避免了比较的复杂的波形和频谱给大家带来困惑。 但是,接下来我会介绍一些比较复杂的情况。
图2.5 显示了一个1100Hz的三角信号在10KHz采样率下的频谱。右图是左图的放大后的图像,这样可以看的更清楚。
这个信号的谐波应该在3300Hz,5500Hz,7700Hz和9900Hz。图中可以看到我们期望的1100Hz和3300Hz的频率, 但是第三个峰值的频率是在4500Hz而不是5500Hz,第四个峰值的频率是在2300Hz而不是7700Hz,下一个峰值的频率 是100Hz而不是9900Hz,这是怎么回事呢?
造成这个情况的原因是,在计算整个信号的波形的过程中,实际上是在采样点在对信号进行了离散化的处理,因此 在连续信号的各个采样点之间会丢失掉一些信息。对于低频的信号丢失的信息不多,因为同样的采样率下,频率低 的信号在一个周期内可以有更多的采样点。
但是如果你用10000Hz的采样率来采集5000Hz的信号,一个信号周期内就仅有两个采样点了。实际上两个采样点是足够的, 但是如果信号的频率再高一点,一个周期内采样点小于两个,那么就会产生问题了。
为了解释这个现象,让我们来看两个余弦信号(4500Hz和5500Hz),我们使用10000Hz的采样率来计算他们的波形:
framerate = 10000
signal = thinkdsp.CosSignal(4500)
duration = signal.period*5
segment = signal.make_wave(duration, framerate=framerate)
segment.plot()
signal = thinkdsp.CosSignal(5500)
segment = signal.make_wave(duration, framerate=framerate)
segment.plot()
图2.6 中灰色的线是信号本身,而蓝色的竖线是采样后的信号。对比这两个图,可以发现, 两个不同的信号却产生了相同的采样值。
事实上,当我们用10000Hz采样率对5500Hz信号进行采样的时候,其结果与4500Hz的信号是相同的。 正是因为这样,7700Hz的信号和2300Hz的信号,9900Hz的信号和100Hz的信号在采样后也是相同的。
信号采样后产生的这种现象,我们就称为 混叠(aliasing) ,简单来说,就是高频信号采样后的与 某些特定的低频信号是无法区分出来的。
在这个例子中(10000Hz采样率),我们最高可以采集的频率为5000Hz,也就是采样率的一半,高于5000Hz的 频率成分会被折叠到5000Hz以内,因此我们把这个频率叫做折叠频率(floding frequence), 又称为 奈奎斯特频率(Nyquist frequency) 。参见 http://en.wikipedia.org/wiki/Nyquist_frequency 。
我们可以这样来计算折叠后的频率:如果信号的频率大于采样率,通过对信号频率与采样率相除求余,来得到在0到采样率之间 的频率,然后如果这个频率大于折叠频率,则用采样率减去这个频率,最后就得到了折叠后的结果。 例如,之前波形的第五个 谐波频率是12100Hz,求余后为2100Hz,就是折叠后的频率了。你也可以从 图2.4 上看到这个2100Hz的频率。同样,也可以 看到4300Hz的频率(14300Hz,折叠后为4300Hz)。
2.4 频谱的计算¶
在之前的章节中,我们多次使用了 make_spectrum
,它的代码(省略了一些细节)是这样的:
from np.fft import rfft, rfftfreq
# class Wave:
def make_spectrum(self):
n = len(self.ys)
d = 1 / self.framerate
hs = rfft(self.ys)
fs = rfftfreq(n, d)
return Spectrum(hs, fs, self.framerate)
self
参数代表的是波形对象本身, n
是波形的采样点数目, d
是采样率的倒数,也就是采样时间步长。
np.fft
是Numpy提供的FFT方法(一种高效的计算DFT的算法)。
make_spectrum
使用了 rfft
,它的意思是“实数FFT”,如果信号是实数而不是复数,我们就可以使用它。
之后,我们会看到“完整FFT”,它可以处理复信号(见 7.9 )。
rfft
的结果 hs
是一个复数的Numpy数组,
它表示了各个频率分量的复数幅值(幅度和初始相位另一种表示形式)。
rfftfreq
的结果 fs
包含了与 hs
对应的频率值。
对于 hs
中的复数,我们可以这样理解:
- 复数是实部和虚部的和,通常写成: \(x + iy\) ,其中 \(i\) 是单位虚数, 也就是 \(\sqrt { - 1}\) 。 我们可以把复数的x和y看做是复数在复平面下的坐标(以实轴为横坐标,虚轴为纵坐标的直角坐标系)
- 复数也可以表示为幅值和复指数的形式,写成: \(A{e^{i\varphi }}\) ,其中 \(A\) 为模, \(\varphi\) 为幅角。我们可以把它看做是复数在极坐标下的表示。
译者注
\(x + iy\) 的极坐标表示为: \(A\cos (\varphi ) + A\sin (\varphi )i\) , 根据欧拉公式 \({e^{ix}} = \cos (x) + isin(x)\) ,可以得出 \(A{e^{i\varphi }}\)
hs
中的每个复数就代表了该频率分量的复数幅值:它的模值就是该频率的幅度,它的幅角就是该频率的初始相位。
Spectrum
类中提供了两个只读的属性: amps
和 angles
,
用来得到这些幅值和初始相位(它们都被放在Numpy数组中)。
我们在画频谱图的时候,一般会画出相对于 fs
下的 amps
,或者相对于 fs
的 angles
。
在实际使用中,我们几乎不会直接去关注 hs
的实部和虚部。
我们也可以直接通过改变 hs
的值来改变频谱,例如:
spectrum.hs *= 2
spectrum.hs[spectrum.fs > cutoff] = 0
第一行代码将 hs
中的元素乘了2, 相当于将所有谐波的幅值增加了2倍。
第二行代码将大于 cutoff
的频率分量的幅值设置到了0。
Spectrum
类中提供了简单的方法来完成这两个操作:
spectrum.scale(2)
spectrum.low_pass(cutoff)
你可以在 http://greenteapress.com/thinkdsp.html 上查看这些方法的文档说明。
至此,你应该对 signal
, Wave
和 Spectrum
这几个类的工作方式有了比较清晰的了解,
但我们还没有解释FFT的原理,接下来的几章我们会慢慢的介绍。
2.5 练习¶
下面练习的答案可以参考文件 chap02soln.ipynb
。
练习1 使用 Jupyter
打开 chap2.ipynb
,阅读并且运行上面的代码示例。
或者在 http://tinyurl.com/thinkdsp02 浏览和运行它。
练习2 锯齿波的波形是周期性的从-1到1线性变化然后立即下降到-1再循环的一种信号,
详见 http://en.wikipedia.org/wiki/Sawtooth_wave 。编写一个锯齿波的类,继承自 Signal
,
复写它的 evaluate
方法来生成锯齿信号。
计算出锯齿信号的频谱,看看它的谐波结构和方波,三角波有什么区别。
练习3 生成一个1100Hz的方波信号以及在10000Hz采样率下的波形,画出频谱图,可以看见大部分的谐波 被混叠了。那么当你听这个声音的时候,你能听到这些混叠后的频率吗?
练习4 生成任一个信号的频谱,把它的 fs
打印出来,你可以看到他们第一个值是0。
也就是说第一个频率分量是0Hz,这代表什么意义呢?试着做如下实验:
- 生成一个440Hz的三角信号,生成并画出它0.01s的波形。
- 生成频谱并打印出
hs[0]
,看看他的幅值和初始相位是多少? - 设置
hs[0]=100
,看看这样做会在波形上产生什么样的影响。 提示:频谱对象有一个make_wave
方法可以生成对应的波形。
练习5 写一个函数,接受一个频谱对象作为参数,将 hs
中的每个元素除以 fs
中
对应的频率。提示:因为除以0是不可行的,你应该设置 hs[0]=0
。用方波,三角波和锯齿波
来测试这个函数。
- 画出原始的频谱图
- 调用你写的函数,并画出计算得到的频谱图
- 使用
make_wave
生成这个频谱的波形,并听听看产生了什么效果。
练习6 三角波和方波只有奇次谐波,锯齿波有奇次谐波也有偶次谐波。方波和锯齿波的谐波 按照 \(\frac{1}{f}\) 的规律衰减,而三角波的谐波是以 \(\frac{1}{{{f^2}}}\) 的规律衰减。你能找到一个波形它包含奇次和偶次的谐波,并且它们是以 \(\frac{1}{{{f^2}}}\) 的规律衰减吗? 提示:有两个方法可以完成这个任务,你可以通过将不同频率的正弦波相加来构造这个波形, 或者你也可以从一个类似的信号开始,然后改变它来生成想要的波形。
[1] | 使用下划线代表一个变量的时候,表示之后不会使用它,这是一个编码惯例。 |
第三章:非周期信号¶
到目前为止,我们学习到的信号都是周期的,它们会无限的重复并循环下去。 这也意味着,它们的频谱不会随着时间的变化而变化。 本章我们会学习非周期的信号,它们的频谱是随时间而变化。 几乎现实中的所有声音信号都是这样的非周期信号。
本章还会介绍声谱图(spectrograms),用于非周期信号的频谱的可视化。
这章的代码 chap03.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp03 查看。
3.1 线性啁啾声
我们先从啁啾声( chirp )开始,这是一种频率随时间变化的信号。 thinkdsp
中提供了一个 Chirp
类,
用来生成频率在一定范围内线性变化的正弦啁啾信号。
以下的代码生成了一个频率从220Hz变化到880Hz(从A3到A5升了两个八度)的 Chirp
信号:
signal = thinkdsp.Chirp(start=220, end=880)
wave = signal.make_wave()
图3.1 的三个图分别显示了这个信号开始,中间和结束的三段波形图。图中可以清晰的看出频率的变化。
我们先来看看 Chirp
是怎么实现的,下面是这个类的代码:
class Chirp(Signal):
def __init__(self, start=440, end=880, amp=1.0):
self.start = start
self.end = end
self.amp = amp
def evaluate(self, ts):
freqs = np.linspace(self.start, self.end, len(ts)-1)
return self._evaluate(ts, freqs)
def _evaluate(self, ts, freqs):
dts = np.diff(ts)
dphis = PI2 * freqs * dts
phases = np.cumsum(dphis)
phases = np.insert(phases, 0, 0)
ys = self.amp * np.cos(phases)
return ys
构造函数 __init__
中 start
和 end
分别表示开始和结束的频率(单位为Hz)。 amp
表示幅值。
evaluate
方法实现了信号的计算,其中 ts
表示采样点的时间序列。简单起见,我们假设采样的时间间隔是固定的。
假设 ts
的长度为 n, 就有 n-1 个时间段, 我们可以使用 np.linspace
求出对应的频率,也就是从 start
到 end
的等间隔和 n-1 个值(Numpy数组)
_evaluate
私有方法完成了接下来的数学运算 [1] 。其中 np.diff
计算了 ts
中相邻采样点的时间间隔,
结果保存在 dfs
中。如果 ts
中的元素是等距的,那么 dts
中的值应该都是一样的。
接下来我们来计算在每个时间间隔内的相位变化。在 1.7 中,我们看到当频率是常量的时候,相位 \(\varphi\) 是相对时间线性变化的:
当频率随时间变化的时候,相位在短时间 \(\Delta t\) 内的变化量是:
因为 freqs
实际上就是 \(f(t)\) ,而 dts
就是 \(\Delta t\) ,因此上式可以写成如下的Python代码:
dphis = PI2 * freqs * dts
现在 dphis
中包含了相位的变化量,我们可以通过累加来得到各个时间点的相位:
phases = np.cumsum(dphis)
phases = np.insert(phases, 0, 0)
np.cumsum
方法计算出了累加值,可以看出来这个值的第一个元素不为0,因此我们需要使用 np.insert
在前面添加一个0值。
最后, 我们使用 np.cos
计算出了整个信号的值(记住相位是用弧度表示的)。
实际上,如果用微积分来表示,当 \(\Delta t\) 足够小的时候:
两边同时除以 \(dt\) 得到:
也就是说,相位的微分就是频率。反过来,相位应该是频率的积分。所以,我们可以使用 np.cumsum
累加来得到相位,
因为累加实际上就是积分的近似计算方法。
译者注
译者觉得这里其实没必要使用累加来计算,因为线性变化的频率很容易可以求得解析解。考虑到相位的微分就是频率,
如果希望频率线性变化,那么 \(f(t) = at + b\),
那么相位就是频率的不定积分 \(\varphi = \frac{a}{2}{t^2} + bt + c\) ,
考虑到单位转换后,最后的信号随时间变化公式为:\(s = \cos (2\pi (\frac{a}{2}t + b)t + c)\)
已知 ts
和 start
end
,很容易可以得到 a 和 b 的值,最终 evaluate
的代码如下:
def evaluate(self, ts):
k = (self.end - self.start) / (ts[-1] - ts[0]) / 2
freqs = k*(ts-ts[0])+self.start
phases = 2 * np.pi * freqs * ts
ys = self.amp * np.cos(phases)
return ys
有兴趣的读者可以尝试使用这个方式来构造信号,看和书中提供的方法产生的信号是否一致。
3.2 指数啁啾声¶
当你听这个啁啾声的时候,你会发现一开始音高上升的很快,然后会慢下来。这个啁啾声跨越了两个八度, 跨越第一个八度只用了 1/3s 时间,而第二个八度用了 2/3s。
造成这个现象的原因是我们感受到的音高取决于频率的对数,也就是说我们听到的两个声音的音高间隔 取决于它们之间的频率比值,而不是差值。用音乐的术语来说,两个音高之间的间隔,被称为音程( interval )
例如,一个八度指的是频率之比为2的两个音高之间的音程。因此从220Hz到440Hz为一个八度, 从440Hz到880Hz又是一个八度。虽然他们之间的频率差更大,但是他们的音程是一样的。
因此,如果频率是线性升高的,那么听起来音高是按对数升高的。
如果我们想得到音高按线性变化的信号,那么信号的频率就得按指数变化。这种信号我们成为指数啁啾声。代码如下:
class ExpoChirp(Chirp):
def evaluate(self, ts):
start, end = np.log10(self.start), np.log10(self.end)
freqs = np.logspace(start, end, len(ts)-1)
return self._evaluate(ts, freqs)
这里我们使用了 np.logspace
来替代 np.linspace
,它可以产生按指数变化的序列值。
其他的代码与之前的 Chirp
是一样的,我们使用它来生成一个指数啁啾声:
signal = thinkdsp.ExpoChirp(start=220, end=880)
wave = signal.make_wave(duration=1)
你可以在 chap03.ipynb
中听一听这些信号的区别。
3.3 啁啾声的频谱¶
啁啾声的频谱图是怎样的呢?这里我们构造了一个1s内八度的信号,并且计算出了它的频谱:
signal = thinkdsp.Chirp(start=220, end=440)
wave = signal.make_wave(duration=1)
spectrum = wave.make_spectrum()
图3.2 展示了这个频谱图。可以看到,这个信号包含从220Hz到440Hz的所有频率成分。 还可以注意到,在220Hz到440Hz区间内,频谱图大概是平的,这就表明频率在时间上是均匀变化的。 那么我们可以猜测指数啁啾声的频谱是什么样子的吗?
实际上,从频谱图中,我们可以得到信号的频率成分的信息,但是却掩盖了频率随时间变化的信息。 我们不能从频谱中看出信号的频率是随时间变大还是变小了。
3.4 声谱图¶
为了展示信号频率随时间变化的关系,我们可以把信号分段后分别计算频谱,然后画出每段的频谱图。 这种方法我们成为 短时傅立叶变换(STFT) 。
我们常用声谱图( spectrogram )来可视化STFT的结果。声谱图的x轴是时间,y轴是频率。 声谱图中的每列显示了一小段时间内信号的频谱,使用灰度值(或颜色亮度)来表示幅值大小。
我们以 Chirp
信号作为例子来计算声谱图:
signal = thinkdsp.Chirp(start=220, end=440)
wave = signal.make_wave(duration=1, framerate=11025)
spectrogram = wave.make_spectrogram(seg_length=512)
spectrogram.plot(high=700)
Wave
类提供了 make_spectrogram
来生成声谱图。其中 seg_length
表示每段包含的采样点数。
这里使用了512,通常情况使用2的n次方的值可以提升FFT的效率。图3.3 为生成的声谱图。
图中,x轴的时间范围从0s到1s,y轴频率范围从0Hz到700Hz。因为信号频率成分比较低,为了更清除的展示, 我把整个声谱图的上部分裁剪了,实际上完整的频率范围是0~5512.5Hz,即采样率的一半。
声谱图清楚的展示了信号频率随时间的变化情况。但是,我们也可以注意到,图中每列的峰值都模糊的覆盖了2-3个色块, 实际上这反应了声谱图的频率分辨率是有限的。
3.5 Gabor limit¶
声谱图在时域上的分辨率是分段的时间长度,在图中对应的是每个色块的宽度。上例中,每段是512个采样点, 在11025Hz的采样率下,大概是0.046s。
而声谱图在频域上的分辨率是频谱上两个相邻频率的间隔,在图中对应的是每个色块的高度。 上例中,对每段512个采样点进行频谱计算后,可以得到分布在0~5512.5Hz的256个频率分量, 也就是说频率分辨率为 5512.5/256 大约为 21.6Hz。
更普遍的来说,如果 n 是分段的长度,那么频谱应该包含 n/2 个分量。 在采样率为 r 的情况下,最大的频率分量应该为 r/2 。 因此频率分辨率为: \(\frac{{r/2}}{{n/2}}\),即 r/n 。 另一方面时间分辨率为分段的长度,即 n/r 。
我们通常希望时间分辨率的值越小越好,这样我们才能够反应出频率的快速变化。 同时,我们也希望频率分辨率的值越小越好,这样才能更准确的描述频率的分布情况。 然而,事实上这两点并不能同时满足,因为时间分辨率 n/r 正好是频率分辨率 r/n 的倒数,也就是说如果一个值变小,那么另一个值就会相应的变大。
例如,如果我们把分段的长度变长两倍,那么同时我们得到的频率分辨率就会减小一半。 即使我们提高采样率也无济于事,因为虽然采样点多了,但是同时获得的频谱的范围也相应的增大了。
这是在进行时频分析时的一个基本原理,称为为 Gabor limit 。
译者注
我没有找到 Gabor limit 的准确翻译,实际上它就是时频分析的不确定性原理, 可以参考 https://en.wikipedia.org/wiki/Uncertainty_principle#Signal_processing 。
3.6 频谱泄露¶
在介绍 make_spectrogram
的工作原理前,我想先介绍一下窗的概念。
这里,我们先研究一个叫做频谱泄露的问题。
计算离散傅立叶变换(DFT)的时候,我们是把有限长度的信号看成是周期信号来处理的, 也就是说,DFT假定进行变换的有限长度信号是一个无限长度的周期信号的一个完整周期。 但是,这个假设通常是错的,并且会产生一些问题。
因为DFT计算的时候是把信号的开始拼接到信号的末尾来构成无限循环的周期信号的。 因此一个普遍的问题是在这个信号的开始和结束的值并不相等,使得最终扩展的周期信号不连续。 这种不连续会造成频谱中包含一些本来不属于信号本身的频率分量。
我们用一个440Hz的正弦信号作为例子,理论上它的频谱应该只有一个440Hz的频率分量:
signal = thinkdsp.SinSignal(freq=440)
如果我们选择一个整数倍周期的时间段来计算频谱,这样它的首尾相连就是连续的, 那么结果没有什么问题:
duration = signal.period * 30
wave = signal.make_wave(duration)
spectrum = wave.make_spectrum()
就像我们期望的那样,这个频谱只有一个440Hz的峰值, 图3.4 的左图展示了这个结果。
但是如果我们选择一个非整数倍周期的长度来计算频谱,就有问题了。
例如 duration = signal.period * 30.25
,这个信号从0开始,以1结束。
图3.4 的中图展示了这个信号的频谱图。可以看到除了440Hz的峰值之外, 频谱中还有一些其他的频率分量,它们分布在240~640Hz之间。 这种现象,我们就称为 频谱泄露(spectral leakage) ,因为整个信号的能量从基频 泄露了一部分到其他的频率上。
在这个例子就是由于我们选了一个非周期的时间段,使得DFT扩展后的周期信号不连续,从而导致了频谱泄露。
3.7 窗函数¶
我们可以通过把不连续的首尾连接处做平滑处理来减小泄露的产生, 加窗(windowing) 就是进行平滑的一种方法。
窗函数就是用来将非周期的信号段转变为周期的信号段的函数。 图3.5 的上图显示了 一个首尾不连续的信号段的波形。
图3.5 的中图展示了一个汉明窗(Hamming window)的形状,这是一个常用的窗函数。 每种窗函数都有各自不同的应用场景,没有哪一个是完美的,而汉明窗是一个比较好的通用窗函数。
图3.5 的下图展示了信号加窗后(信号与窗函数相乘)的波形。这样处理之后,在窗函数值的1的时候, 信号是不变的,而在窗函数值为0的时候信号就被衰减了。由于窗函数的中间部分值比较大,而两端慢慢的 减小到0,因此信号的两端也被逐渐的衰减到0,最终的结果就是加窗后的信号首尾连接处被平滑了。
图3.4 的右图展示了这个加窗后的信号的频谱图,可以看到,加窗后在很大程度上的减小了频率泄露, 但是并没有完全消除。
Wave
类中提供了加窗的方法 window
,代码如下:
#class Wave:
def window(self, window):
self.ys *= window
Numpy提供了一些常用的窗函数,其中 hamming
可以产生一个给定长度的汉明窗,
下面的代码对波形应用了一个汉明窗:
window = np.hamming(len(wave))
wave.window(window)
Numpy提供的窗函数包括: bartlett
, blackman
, hanning
和 kaiser
。
你可以在本章后面的练习中对这些窗函数进行试验。
3.8 声谱图的实现¶
现在我们来看看 make_spectrogram
是怎么实现的,代码如下:
#class Wave:
def make_spectrogram(self, seg_length):
window = np.hamming(seg_length)
i, j = 0, seg_length
step = seg_length / 2
spec_map = {}
while j < len(self.ys):
segment = self.slice(i, j)
segment.window(window)
t = (segment.start + segment.end) / 2
spec_map[t] = segment.make_spectrum()
i += step
j += step
return Spectrogram(spec_map, seg_length)
这段代码是本书中最长的一段代码,如果你能弄懂了它,那么其他的代码就肯定没什么问题了。
首先 self
参数时 Wave
对象本身, seg_length
是每段的长度。
window
是与每段长度相同的汉明窗。
i
和 j
用来指示循环中处理的每段的开始和结束的位置。 step
是窗口每次移动的步长,
这里 step
是 seg_length
的一半,也就是说每次处理的段都是有一半是重叠的,
图3.6 展示了这些重叠的窗口。
spec_map
是一个字典 dictionary
,用于存储时间到频谱的映射关系。
循环的内部,我们通过 i
j
从波形中截取一段,然后应用窗函数,
然后进行频谱计算并把计算存放到 spec_map
中。其中, t
是每段的中间时间点。
接下来就是对 i
j
进行递增,并在 j
到达波形结尾后结束循环。
最后,通过 spec_map
和 seg_length
构造了 Spectrogram
对象。
Spectrogram
类的定义如下:
class Spectrogram(object):
def __init__(self, spec_map, seg_length):
self.spec_map = spec_map
self.seg_length = seg_length
Spectrogram
的构造函数就是简单的保存了参数,这个类提供的 plot
方法
使用伪彩色画出了时间和频率的变化关系图,也就是声谱图。
3.9 练习¶
下面练习的答案可以参考文件 chap03soln.ipynb
。
练习1 使用 Jupyter
打开 chap3.ipynb
,阅读并且运行上面的代码示例。
或者在 http://tinyurl.com/thinkdsp03 浏览和运行它。
在频谱泄露的那个例子中,试着将汉明窗替换成Numpy中的其他窗函数,看看结果会怎么样。
参考 http://docs.scipy.org/doc/numpy/reference/routines.window.html 。
练习2 编写一个锯齿啁啾声类 SawtoothChirp
继承自 Chirp
,
并复写 evaluate
方法生成一个频率线性增大或减小的锯齿波信号。
提示:可以结合 SawtoothSignal
和 Chirp
两个类的 evaluate
方法。
凭你自己的想象,在纸上大致画出这个信号的声谱图,然后再用代码画出来。 可以看出,频率混叠的影响是很明显的,如果你仔细的听应该可以听出来。
练习3 生成一个频率范围为2500到3000Hz的 SawtoothChirp
,并生成在20kHz采样率下的1s的波形,
在纸上画出你认为的频谱图,然后用代码画出频谱图,看看你画的对不对。
练习4 在音乐术语中,滑音指的是一种从一个音高渐变到另一个音高,与啁啾声类似。 找一段滑音的录音,画出它前几秒的声谱图。建议:George Gershwin 的蓝色狂想曲(Rhapsody)就是从一段著名的 滑音开始的,你可以从这里下载到它:http://archive.org/details/rhapblue11924 。
练习5 长号的滑奏是通过滑动滑管来发出滑音的。当滑管滑动的时候,管子的长度在变化,发出的声音的音高与长度成 反比例的变化。假设演奏者匀速的移动滑管,发出的声音频率随时间应该是怎么样变化的呢?
编写一个长号滑奏的类 TromboneGliss
,继承自 Chirp
并复写 evaluate
,
生成模拟长号从C3到F3然后再从F3到C3的滑音。C3大概是262Hz,F3是349Hz。
画出这个声音的声谱图,看看长号的滑音更接近线性啁啾声还是指数啁啾声?
练习6 生成或者找一些元音的录音,画出他们的声谱图,看看你能识别出不同元音的声谱吗?
[1] | 方法名前面加下划线表示这个方法是私有的,不应该在外部进行调用 |
第四章:噪声¶
通常我们所说的噪声(noise)指的是那些听起来混乱不舒服的声音。而在信号处理中,噪声还有另外两层含义:
- 噪声是指那些我们不期望出现的信号。例如,两个信号重合在一起相互干扰, 那么对于其中的一路信号来说,另一路信号就是噪声。
- 噪声还指那些频谱中包含几乎所有的频谱成分,不具有前几章那些周期信号的谐波结构的信号。
本章主要关注噪声的第二种含义。
这章的代码 chap04.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp04 查看。
4.1 不相关噪声¶
了解噪声的一种最简单的方法就是生成一段噪声。不相关的均匀噪声(UU noise)就是一种很容易生成的噪声。 均匀(Uniform)的意思是信号的值是随机的并且服从均匀分布,也就是说在一定范围内信号出现任意值的概率 都是一样的。不相关(Uncorrelated)的意思是信号在不同时间的值是独立不相关的;也就是说, 某一时刻的信号值不能够提供其他任何时候的信号值的信息。
下面这个类就是UU噪声:
class UncorrelatedUniformNoise(_Noise):
def evaluate(self, ts):
ys = np.random.uniform(-self.amp, self.amp, len(ts))
return ys
UncorrelatedUniformNoise
继承自 _Noise
,而 _Noise
又继承自 Signal
。
evaluate
方法传入参数 ts
,也就是信号求值的时间序列,然后使用 np.random.uniform
生成
服从在[-amp,amp]范围均匀分布的值。
下面的代码生成了在11025Hz采样率下的0.5s的UU噪声:
signal = thinkdsp.UncorrelatedUniformNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
播放这个波形,你会发现它听起来像是我们切换广播频道时的那种静态噪声。 图4.1 显示了这个噪声的波形,就像之前所说的,它的值是随机的。
接下来我们来看看它的频谱:
spectrum = wave.make_spectrum()
spectrum.plot_power()
spectrum
的 plot_power
方法与 plot
方法类似,区别在与它画的是
功率(Power)而不是幅值,功率是幅值的平方。在研究噪声的情况下,我们使用功率
比使用幅值更方便。
图4.2 展示了这个频谱图。与波形图一样,它的频谱图看起来也是随机的。 实际上,它确实是随机的,但是在这里,我们得更加精确的定义一下“随机”这个词。 对于噪声信号来说,我们需要关注它的三种性质:
- 分布(Distribution):随机信号的分布值的是他们可能产生的值以及产生这些值的概率。 例如,在均匀噪声中,信号的值可以在[-1,1]的范围中,产生这些值的概率都是一样的。 而高斯噪声( Gaussian Noise )则不同,它的值可以从负无穷到正无穷,但是越接近0 的值出现的概率越高,它的概率分布图是一个高斯分布曲线(也叫正态分布曲线或钟型曲线)。
- 相关性(Correlation):表示信号的不同值之间是独立的或是有依赖关系。在UU噪声中, 所有的值都是独立的。而布朗噪声( Brownian noise )则不同,它的每个值都是前一时刻 信号的值再加上一个随机的“步值(step)”。因此,如果某一个时刻信号的值比较大, 我们可以预测它接下来依然会保持比较大的值,反之亦然。
- 功率和频率之间的关系:在UU噪声中,所有频率成分的功率看起来是均匀分布的,也就是说信号的功率 被平均分布到的所有的频率成分上。而粉红噪声( Pink Noise )则不同, 它的功率与频率成反比例关系,也就是说对于频率为 f 的成分来说,它的功率为 1/f
4.2 频谱积分¶
对于UU噪声来说,我们可以通过对频谱进行积分( integrated spectrum )来更清晰的观察到频率和 功率之间的关系。频谱积分就是计算频率从0到f上的功率的累加和。
spectrum
类提供了一个方法来产生积分后的频谱:
def make_integrated_spectrum(self):
cs = np.cumsum(self.power)
cs /= cs[-1]
return IntegratedSpectrum(cs, self.fs)
self.power
是包含每个频率成分的功率的Numpy数组, np.cumsum
用于计算累加和。
除以最后一个元素的值可以将积分后的值进行归一化。最后返回了 IntegratedSpectrum
对象,
这个类的定义如下:
class IntegratedSpectrum(object):
def __init__(self, cs, fs):
self.cs = cs
self.fs = fs
和 spectrum
一样, IntegratedSpectrum
也提供 plot_power
方法用来计算并作图:
integ = spectrum.make_integrated_spectrum()
integ.plot_power()
如 图4.3 所示,积分后的频谱是一条直线,表示所有频率成分的平均功率是一个常值。 像这样的噪声,我们称为 白噪声(white noise) ,之所以叫白噪声,是因为如果灯发出所有频率的可见光 并且功率都一样,那么灯光的颜色应该是白色的。
4.3 布朗噪声¶
UU噪声是不相关的,每个时刻的值互相都是独立的。而布朗噪声就是相关的,它的值是上一时刻的值加上一个随机的步长。 之所以叫布朗噪声,是由于它和布朗运动很类似。悬浮在液体上微小粒子由于受到分子的撞击产生的无规则运动,被称为 布朗运动,在数学上,可以用 随机游走(random walk) 来描述,也就是每步运动的距离服从一定的随机分布。
以一维的随机游走来说,每次微粒向上或者向下移动一个随机的距离,那么微粒的位置应该等于之前移动的距离之和。 我们可以用这个思路来生成一个布朗噪声:先生成一个不相关的随机步长,然后将他们累加,如下面的代码:
class BrownianNoise(_Noise):
def evaluate(self, ts):
dys = np.random.uniform(-1, 1, len(ts))
ys = np.cumsum(dys)
ys = normalize(unbias(ys), self.amp)
return ys
evaluate
使用了 np.random.uniform
来生成不相关的均匀分布的步长,
然后用 np.cumsum
来计算累加值。
然后用 unbias
来调整信号的偏移量,使其均值为0,最后用 normalize
调整到给定的最大幅值。
下面的代码生成了一个布朗噪声,并画出了波形图:
signal = thinkdsp.BrownianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
wave.plot()
如 图4.4 所示,布朗噪声的波形在上上下下的变化,但是连续的两个值之间有明显的相关性, 当值比较大的时候,下一时刻的值也会比较大,反之亦然。
在线性刻度下,布朗噪声的频谱图的几乎所有的频率成分都是在低频,高频分量几乎不可见,如 图4.5 左图。
为了更清楚的显示频谱图,我们将它画在对数刻度坐标中,代码如下:
spectrum = wave.make_spectrum()
spectrum.plot_power(linewidth=1, alpha=0.5)
thinkplot.config(xscale='log', yscale='log')
如 图4.5 ,可以看出频率和功率的关系虽然也呈现噪声的特点,但是又有一定的线性规律。
Spectrum
类中还提供了 ectimate_slope
方法使用 Scipy
中的最小二乘法来对功率谱进行拟合:
#class Spectrum
def estimate_slope(self):
x = np.log(self.fs[1:])
y = np.log(self.power[1:])
t = scipy.stats.linregress(x,y)
return t
由于 log0 是未定义的,因此我们的代码去掉了频率为0的分量。
estimate_slope
使用了 scipy.stats.linregress
来计算最小二乘,返回值 t
中包含了拟合直线的
斜率和截距,相关系数,P值以及标准差。这里,我们只关心斜率和截距。
对于布朗噪声来说,斜率是 -2 (见 第九章 ),我们可以把频率和功率的关系写成如下形式:
这里 P 代表功率, f 代表频率, k 代表直线的截距(不重要)。将等式两边求幂后得到:
这里的 K 等于 \({e^k}\) (依然不重要)。我们可以看到布朗噪声的一个很重要的特性就是 它的功率与 \(1/{f^2}\) 成正比例关系。
与白噪声类似,由于呈现这种频率-功率关系的灯光是红色的,我们也把布朗噪声称为 红噪声(red noise)。
4.4 粉红噪声¶
对于红噪声来说,频率和功率的关系是:
这里,频率平方的关系其实可以更一般化的表示为指数 \(\beta\) :
当 \(\beta = 0\) 的时候,所有频率分量的功率都是一个常数,也就是白噪声。 当 \(\beta = 2\) 的时候,是红噪声。
当 \(\beta\) 在0~2以内的时候,叫做 粉红噪声(pink noise) 。
生成粉红噪声的一种简单的方式是在UU噪声的基础上应用一个低通滤波器。
thinkdsp
中提供一个粉红噪声的类 PinkNoise
class PinkNoise(_Noise):
def __init__(self, amp=1.0, beta=1.0):
self.amp = amp
self.beta = beta
def make_wave(self, duration=1, start=0, framerate=11025):
signal = UncorrelatedUniformNoise()
wave = signal.make_wave(duration, start, framerate)
spectrum = wave.make_spectrum()
spectrum.pink_filter(beta=self.beta)
wave2 = spectrum.make_wave()
wave2.unbias()
wave2.normalize(self.amp)
return wave2
其中, amp
是噪声的幅值, beta
为之前所说的频率的指数。
make_wave
用于生成波形。其中 duration
是需要生成的长度。
start
是波形的开始时间, framerate
是采样率。
由于要与所有信号类的 make_wave
方法保持一致,因此这些参数都必须写出来,
虽然对于随机的噪声来说,开始时间和采样率都是没用的。
make_wave
方法首先生成了一个白噪声以及它的频谱,然后应用了一个滤波器 pink_filter
。
最后进行 unbias
和 normalize
操作。
Spectrum
中提供的 pink_filter
方法如下:
def pink_filter(self, beta=1.0):
denom = self.fs ** (beta/2.0)
denom[0] = 1
self.hs /= denom
它将所有频率成分的幅值除以 \({f^{\beta /2}}\) 。由于功率是幅度的平方,这个操作使得 功率变化为原来的 \(1/{f^\beta }\) 。当然,为了避免与0相除,我们将第0个元素设置为了1, 这意味着信号的直流分量的幅值不会改变。
图4.6 展示了粉红噪声的波形,与布朗噪声有点类似,它有一些保持当前的运动趋势的特点, (也就是相邻两个值之间的相关性),但是又增加了一些随机性。下一章我们会再回过头来准确的 解释我所说的“相关性”和“随机性”。
最后,图4.7 展示了白噪声,粉红噪声和红噪声的频谱(在对数刻度下)。图中,我们可以很明显的 看出频率和功率在对数刻度下的线性关系,也就是斜率 \(\beta\) 。
4.5 高斯噪声¶
UU噪声中所有频率成分的功率都是一个固定的常数,因此UU噪声是白噪声。 但是通常我们提到“白噪声”的时候,指的并不是UU噪声, 而是不相关的高斯噪声(uncorrelated Gaussian Noise),简称UG噪声。
thinkdsp
中提供了一个UG噪声的类:
class UncorrelatedGaussianNoise(_Noise):
def evaluate(self, ts):
ys = np.random.normal(0, self.amp, len(ts))
return ys
np.random.normal
计算出了均值( \(\mu\) )为0,
方差( \(\sigma\) )为 self.amp
的高斯分布。
在这个分布下,信号值的范围是正无穷到负无穷之间,
但是99%的值都分布在 \((\mu - 3\sigma ,\mu + 3\sigma )\) 区间。
UG噪声在很多方面都与UU噪声相似,它的频谱的所有分量都有常值的平均功率, 因此属于白噪声。一个有趣的特性是,UG噪声的频谱也是UG噪声,更准确的说是, UG噪声的频谱的实部和虚部均服从高斯分布。
为了验证这个说法,我们来生成UG噪声的频谱,然后生成一个可以直观的判断 数据是否服从高斯分布的“正态概率图”:
signal = thinkdsp.UncorrelatedGaussianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
spectrum = wave.make_spectrum()
thinkstats2.NormalProbabilityPlot(spectrum.real)
thinkstats2.NormalProbabilityPlot(spectrum.imag)
NormalProbabilityPlot
方法包含在 thinkstats2
中,你可以在本书的代码库中
找到这个文件。
关于正态概率图,可以参考 《ThinkStats》中的第五章, http://thinkstats2.com 。
在正态概率图中,如果数据与直线拟合的很好,则说明数据是服从高斯分布的。 图4.8 中灰线代表数据的拟合结果,蓝线是数据。可以看出,除了很少的点外, 数据和直线是拟合的很好的,说明了UG噪声的频谱也是UG噪声。
实际上,根据中心极限定理,只要信号的分布具有有限的均值和标准差并且采样点足够多, 则任何不相关的噪声的频谱都是服从高斯分布的。
4.6 练习¶
下面练习的答案可以参考文件 chap04soln.ipynb
。
练习1 网站 A Soft Murmur 中可以找到 很多自然界的混合噪声,包括雨声,波浪声,风声等。下载一些声音文件,计算它们的频谱, 看看他们的功率谱是否像白噪声,粉红噪声或是布朗噪声?这些频谱相对于时间是如何变化的?
练习2 像这样的噪声,它的频率成分通常都是随时间变化的。长期来看,所有频率成分的功率
应该是一样的,但是在一段短时间内,功率又是随机的。为了估计长时间的平均功率,我们可以把信号
分段,并分别计算每段时间内的功率谱,最后计算这些分段的功率的平均值,
这个算法叫做 Bartlett’s method ,参见 http://en.wikipedia.org/wiki/Bartlett’s_method 。
试着实现这个算法,并用它来计算一个噪声信号的功率谱。 提示:可以参考 make_spectrogram
的实现。
练习3 在 http://www.coindesk.com 上下载比特币的每日价格数据(csv文件),使用代码将数据 加载进来,并计算它的功率谱。看看它是否与白噪声,粉红噪声或是布朗噪声类似。
练习4 盖革-米勒计数器(Geiger counter)一种专门探测电离辐射强度的记数仪器。每当离子激发探测器,
就会输出一个电流脉冲。在单位时间输出脉冲的个数可以用一个不相关的泊松噪声(UP噪声)来描述,也就是说在单位时间
内检测到的离子的个数服从泊松分布。
编写一个 UncorrelatedPoissonNoise
类,继承自 _Noise
并复写 evaluate
方法。
使用 np.random.poisson
来生成服从泊松分布的随机信号值。 该类提供一个 lam
属性表示
单位时间发生次数的均值,你可以使用 amp
来作为 lam
使用。例如,当采样率是10kHz的时候,
如果 amp
是0.001,意味着没秒钟大约会出现10次。
生成1s的UP噪声并听一听,你会发现对于较小的 amp
,比如0.001,它会听起来像是盖革-米勒计数器。
而对于较大的 amp
,它会听起来像是白噪声。计算并画出它的功率谱,看看是否像白噪声。
练习5 本章中介绍的计算粉红噪声的方法相对来说比较简单,但是效率比较低。 Voss-McCartney 算法是一种更高效的算法。研究并实现这个算法,计算出结果的功率谱,看看是否具有期望的 频率与功率的关系。
第五章:自相关¶
上一章中,我们说白噪声是不相关的,它任意时刻的值相互之间都是独立的。 而布朗噪声是相关的,它某个时刻的值依赖于上一个时刻。 这章中,我们引入 自相关函数(autocorrelation function) 来精确的定义这些概念,这在信号分析中是很有用的一种方法。
这章的代码 chap05.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp05 查看。
5.1 相关性¶
通常来说,如果一个变量与另一个变量相关,那么意味着如果知道一个变量的值就可以得到另一个变量的一些信息。 我们有好几种方式来度量相关性, 其中最常见的一种是皮尔森相关系数(Pearson product-moment correlation coefficient), 记为: \(\rho\) 。对于两个均包含N个测量值的变量 x 和 y ,
上式中的 \({\mu _x}\) 和 \({\mu _y}\) 是 x 和 y 的均值, \({\sigma _x}\) 和 \({\sigma _x}\) 是它们的标准差。
皮尔森相关系数总是在[-1,1]的范围内,如果是正值,表示两个变量正相关,也就是说一个变量如果比较大, 那么另一个也趋向于变大,相反,如果是负值,那么表示两个变量负相关,也就是说一个变量如果比较大, 那么另一个会趋向于变小。
\(\rho\) 的大小表示两个变量之间的相关程度,越大就越相关。如果 \(\rho\) 是1或-1, 则说明它们是完全相关的,如果我们知道其中一个的值,就能准确预测出另一个值。如果 \(\rho\) 接近0,说明它们之间的相关性可能很弱,也就是知道其中一个的值,对于我们预测另一个的值来说没有什么意义。
这里,我用了“可能很弱”,是因为如果它们之间是非线性的关系,那么相关系数并不能很好的表征这种非线性关系。 在统计学中,非线性相关性也很重要,但是在信号处理中没有那么常见,因此这里我们先不考虑它。
Python中有很好几种方法来计算相关性。其中 np.corrcoef
可以计算多个变量两两之间的相关系数,结果用
一个 相关矩阵(correlation matrix) 来表示。
这里以两个变量为例,首先我构造了一个生成不同初始相位的正弦信号的函数,如下:
def make_sine(offset):
signal = thinkdsp.SinSignal(freq=440, offset=offset)
wave = signal.make_wave(duration=0.5, framerate=10000)
return wave
然后我们生成了两个不同相位的正弦信号:
wave1 = make_sine(offset=0)
wave2 = make_sine(offset=1)
图5.1 显示了这两个信号的前几个周期的波形。可以看到,一个信号的值比较大的时间,另一个的值也相对比较大, 感觉它们之间应该是相关的。
它们之间的相关系数计算如下:
>>> corr_matrix = np.corrcoef(wave1.ys, wave2.ys)
[[ 1. 0.54]
[ 0.54 1. ]]
结果是一个相关矩阵,每个位置的值都代表的是对应的行号和列号两个信号之间的相关系数。
例如,第一行第一列元素的值代表了 wave1
和自身的相关系数,同样第二行第二列的值代表了
wave2
和自身的相关系数。
其中非对角元素的值是我们所关心的,也就是 wave1
和 wave2
之间的相关系数0.54,
表明了它们之间有一定的相关性。
当它们之间的相位差增大时,它们的相关系数会减小,直到相位差为180°,这时相关系数为-1。然后, 随着相位差继续增大,相关系数值又会增大,直到360°,这个时候它们的相关系数为1.
图5.2 显示了它们之间的相关系数随相位差变化的情况。这个曲线的形状你应该很熟悉,它是余弦曲线。
thinkdsp
中提供一个简单的方法来计算波形之间的相关性:
>>> wave1.corr(wave2)
0.54
5.2 序列相关性¶
信号一般情况下都是对某个量的在一段时间内的连续测量值。 例如,声音信号就是在一段时间内我们对声压相对应的电压(电流)的测量值。
像这样的测量值都具有序列相关性,它表示一段测量序列中的不同时刻的测量值之间的相关性。 我们可以通过对信号进行延迟,再于原始信号计算相关系数来得到序列相关性:
def serial_corr(wave, lag=1):
n = len(wave)
y1 = wave.ys[lag:]
y2 = wave.ys[:n-lag]
corr = np.corrcoef(y1, y2, ddof=0)[0, 1]
return corr
serial_corr
接受一个波形对象和 lag
(一个整数,表示延迟的大小)作为参数,
然后计算出延迟后的波形和原始波形间的相关系数。
我们用上一章中的UG噪声信号来对这个函数进行测试:
signal = thinkdsp.UncorrelatedGaussianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)
这个代码运行的结果是0.006,这表示信号的序列相关性很小。你可以试试用不同的 lag
值,
你会发现它们同样都很小。这是因为UG噪声是不相关的。
而对于布朗噪声来说,它的值是上一时刻的值加一个随机的步长,因此应该具有较大的序列相关性:
signal = thinkdsp.BrownianNoise()
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)
我可以保证上面这段代码运行的结果要大于0.999。
对于粉红噪声来说,它介于UG噪声和布朗噪声之间,相关系数也应该位于中间:
signal = thinkdsp.PinkNoise(beta=1)
wave = signal.make_wave(duration=0.5, framerate=11025)
serial_corr(wave)
当 \(\beta = 1\) 的时候,它的序列相关系数为0.851; 当 \(\beta = 0\) 的时候,即UG噪声; 当 \(\beta = 2\) 的时候,即布朗噪声。 如 图5.3 所示,粉红噪声的序列相关系数的范围为[0,1]。
5.3 自相关¶
上一小节中我们计算了相邻两个时刻的信号值之间的相关性(也就是 lag=1
)。
实际上,我们可以很容易的计算出不同 lag
的相关系数。
实际上,你可以把 serial_corr
理解为一个从 lag
到相关系数的映射。
我们可以循环的计算出不同 lag
的相关系数:
def autocorr(wave):
lags = range(len(wave.ys)//2)
corrs = [serial_corr(wave, lag) for lag in lags]
return lags, corrs
autocorr
接收一个波形对象作为参数,并返回一个序对形式的自相关函数:
lags
是从0到波形长度一半的整数; corrs
是相对应的序列相关系数。
图5.4 显示了三个不同 \(\beta\) 的粉红噪声的自相关函数的曲线图。
可以看出,对于较小的 \(\beta\) ,信号基本是不相关的,它的自相关函数很快的下降到0附近。
对于较大的 \(\beta\) ,信号是强相关的,它的自相关函数比较大,并且下降比较慢。
当 \(\beta=1.7\) 的时候,即使很大的 lag
值也有较强的相关性。这种情况,我们称为
长期相关(long-range dependence) ,因为信号的值依赖于长时间之前的值。
5.4 周期信号的自相关性¶
粉红噪声的自相关性虽然有趣,但是在实际应用中比较有限。 相比而言,周期信号的自相关性更有用一些。
这里我以一个从 freesound.org 上下载的声音作为例子。
这个文件可以在 代码库 中找到,是一个人唱出的啁啾声的录音。
你可以在 chap05.ipynb
中播放它。
图5.5 展示了这段声音的声谱图,图中可以清晰的看出基频和其他的谐波成分。 这个啁啾声从大概500Hz开始下降到大概300Hz,大致是从C5到E4。
我可以使用频谱来估计某一个时刻的音高,但是实际上,这种方法不是太好。 我们来看看这是为什么,这里我截取了一小段声音波形并画出了它的频谱:
duration = 0.01
segment = wave.segment(start=0.2, duration=duration)
spectrum = segment.make_spectrum()
spectrum.plot(high=1000)
这段波形从0.2s开始持续了0.01s的时间。频谱见 图5.6 。图中可以清晰的看到400Hz的峰值。 但是我们不能就此断定音高就是400Hz。由于波形的长度是0.1s, 因此频率分辨率是100Hz(参见 3.5 )。 这就意味着估计的音高总是会有50Hz的偏差,这对于音乐来说就差的太远了,因为从350Hz到450Hz 已经跨越了5个半音,是一个很大的音程了。
当然,我们可以使用一个较长的分段来提高频率的分辨率,但是这样一来,音高随时间的变化就会变得模糊; 也就是说峰值会分布在这段波形的开始和结束的音高上,见 3.3 。
运用自相关性,我们可以更精确的估算音高。如果信号是周期的,那么当 lag
等于周期的时候,信号的
自相关函数应该是最大的。
下面的代码,我画出了录音的不同的两段:
def plot_shifted(wave, offset=0.001, start=0.2):
thinkplot.preplot(2)
segment1 = wave.segment(start=start, duration=0.01)
segment1.plot(linewidth=2, alpha=0.8)
segment2 = wave.segment(start=start-offset, duration=0.01)
segment2.shift(offset)
segment2.plot(linewidth=2, alpha=0.4)
corr = segment1.corr(segment2)
text = r'$\rho =$ %.2g' % corr
thinkplot.text(segment1.start+0.0005, -0.8, text)
thinkplot.config(xlabel='Time (s)')
其中一段声音从0.2s开始,另一段从这之后0.0023s开始。 如 图5.7 所示,这两段信号很相似,他们的相关系数是0.99。 这个结果表明信号的周期大约是0.0023s,对应的频率是435Hz。
上面,我是通过试错的方法找出了周期。我们也可以用自相关函数来自动的计算:
lags, corrs = autocorr(segment)
thinkplot.plot(lags, corrs)
如 图5.8 ,这段信号的自相关函数从0.2s开始,在 lag=101 处达到峰值。 对应的频率可以像下面这样计算:
period = lag / segment.framerate
frequency = 1 / period
这样估算出来的基频是437Hz。我们可以用同样的方法计算出 lag
为100和102时的
频率值为432Hz和441Hz,也就是说这里估算的频率的精度小于10Hz,这与使用频谱估计的
100Hz的精度相比,好太多了。
在音乐中,我们可以分辨30分的音高,大约是一个半音的1/3。
5.5 点积¶
在本章开始引入了皮尔森相关系数:
然后我们使用 \(\rho\) 来定义了序列相关性和自相关性。 这些定义实际上都是在统计学上使用的,在信号处理中,这些定义还稍微有些不同。
信号处理中,我们通常处理的是无偏(均值为0)归一化(标准差为1)的信号,在这种情况下, \(\rho\) 的定义简化为:
进一步简化为:
这个定义不是标准的,它的值的范围不是[-1,1],但是它有一些别的有用的特性。
如果把 x 和 y 看做是向量,那么这个公式就是 点积(dot product) 的公式。 参见 http://en.wikipedia.org/wiki/Dot_product 。
也就是说点积表征了信号之间的相似度,如果他们都是归一化的,那么:
这里, \(\theta\) 就是两个向量之间的夹角,这也就解释了 图5.2 是一个余弦曲线的原因。
5.6 使用Numpy¶
Numpy中提供了计算互相关函数的方法 correlate
,我们可以用它来计算上一小节的所说的自相关函数:
corrs2 = np.correlate(segment.ys, segment.ys, mode='same')
mode
参数表示用于计算的 lag
范围,当使用 same
的时候,范围为 -N/2 到 N/2 (N为信号长度)。
如 图5.9 所示,自相关函数是对称的,这是因为两个相同信号进行相关运算的时候,
正的lag值与负的lag值产生的影响是一样的。我们选择后半部分来和 autocorr
进行对比:
N = len(corrs2)
half = corrs2[N//2:]
对比 图5.8 与 图5.9 的图像,你会发现,用 np.correlate
计算的自相关函数随着lag
值的增加会更小一些,这是因为 np.correlate
使用的是相关性的非标准定义,随着lag的增大,
两个信号之间的重合部分变小了,因此相关系数也变小了。
我们可以通过除以长度来校正这个问题:
lengths = range(N, N//2, -1)
half /= lengths
经过这样调整后, autocorr
和 np.correlate
的值就基本上相同的。虽然还是有1%~2%的差别,
但那其实已经不那么重要了。造成这个差别的原因是,是由于 autocorr
对于每个lag值先进行了归一化
操作,而 np.correlate
没有。
现在你应该对自相关有所了解了,并且知道如何使用它来估计信号的基频,以及计算它的两种方式。
5.7 练习¶
下面练习的答案可以参考文件 chap05soln.ipynb
。
练习1 chap05.ipynb
中有一个交互式计算自相关函数的单元,可以直观的计算不同的lag值的自相关系数。
使用这个方法来估计本章中的人声啁啾声在不同时间的音高。
练习2 chap05.ipynb
中有一段使用自相关来估计周期信号基频的代码。
将这段代码封装成 estimate_fundamental
函数,并用这个函数计算一段录音的音高变化曲线。
把计算结果和声谱图画到一起,看看效果怎么样?
练习3 使用上一章练习中的比特币价格数据,计算它的自相关函数。看看它是否会快速的减小? 分析数据是否有明显的周期性?
练习4 本书的 代码库 中有一个 saxophone.ipynb
文件,里面讨论了自相关性,音高感知
以及一个叫 “基频缺失” 的现象。阅读并运行这个文件,试着使用不同的声音来重新运行里面的实例。
观看视频 What is up with Noises? (The Science and
Mathematics of Sound, Frequency, and Pitch) ,
它演示了基频缺失现象并解释了音高感知的原理。
第六章:离散余弦变换¶
本章的主题是 离散余弦变换(Discrete Cosine Transform, DCT) , 它在MP3和其他格式的音频文件压缩中经常使用,在图像和视频的压缩中也会用到。
DCT与频谱分解时用到的DFT很类似。学习了DCT之后可以更容易的理解DFT。
我们分几个步骤来学习DCT:
- 先从合成信号开始:给定一组不同频率的余弦信号以及它们的幅值,构造一个合成的信号波形。
- 使用Numpy来重写合成的函数,这样可以提升计算效率并且为下一步打好基础。
- 分解信号:给定一个信号和每个频率分量,找到每个频率的幅值。我们先介绍一个简单直接的方法, 但是比较慢。
- 最后,我们使用一些线性代数的知识来改善计算效率。如果你有线性代数的基础是最好的, 不过没有也没关系,我会把使用到的知识都简单介绍一下。
这章的代码 chap06.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp06 查看。
6.1 合成¶
给定一组余弦信号的频率和幅值,我们可以将它们相加来合成一个信号。
在 thinkdsp
模块中, 可以使用 synthesis
函数来完成这个操作:
def synthesize1(amps, fs, ts):
components = [thinkdsp.CosSignal(freq, amp)
for amp, freq in zip(amps, fs)]
signal = thinkdsp.SumSignal(*components)
ys = signal.evaluate(ts)
return ys
其中 amps
是幅值, fs
是对应的频率值, ts
是用于计算波形的时间序列。
components
是一组由给定幅值和频率构造的余弦信号, SumSignal
计算出来这些信号的和。
最后,使用 evaluate
计算出了 ts
的信号波形。
可以像下面这样来测试这个函数:
amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = [100, 200, 300, 400]
framerate = 11025
ts = np.linspace(0, 1, framerate)
ys = synthesize1(amps, fs, ts)
wave = thinkdsp.Wave(ys, framerate)
这个例子中,合成信号的频率成分为100Hz的基频和三个谐波频率(100Hz相当于G♯2)。 然后使用合成的信号在11025Hz采样率下生产了波形对象。
虽然像这样合成信号是非常简单的,但是对我们 分解信号 没有什么帮助。 分解信号就是合成的逆操作:给定一个信号,将它分解成不同频率和幅值的余弦信号。
6.2 使用Numpy数组的合成¶
合成信号的另一个方式是这样:
def synthesize2(amps, fs, ts):
args = np.outer(ts, fs)
M = np.cos(PI2 * args)
ys = np.dot(M, amps)
return ys
这段代码看上去和之前的完全不同,但是它们实际上完成了相同的事情, 让我们看看它是怎么工作的:
np.outer
计算了ts
与fs
的 外积(outer product) 。 外积的结果是一个n行m列的矩阵(n为ts的长度,m为fs的长度),其中的每个元素 都是对应位置的时间和频率的乘积 ft 。- 将
args
乘以 \(2\pi\) 并应用np.cos
计算出了 \(\cos (2\pi ft)\) 。也就是说 M 中的每一列都是一个频率分量 在ts
时间序列下的值。 np.dot
计算出了M
与amps
的点积,也就是将amps
的每个值与M
中每行对应的值相乘并求和。这样就求出了每个时间的合成信号的值。
图6.1 展示了这个计算的过程。 矩阵 M 的每一行对应了0~1s的时间序列, 其中第n个元素对应的时间用 \({t_n}\) 表示;每一列对应了100~400Hz的频率序列, 其中第k个元素对应的频率用 \({k_n}\) 表示。
我使用了a,b,c,d来表示每一行的四个元素,这里可以得出: \(a = \cos [2\pi (100){t_n}]\)
点积的结果 ys
的每个值对应了 M 的一行,记为 e ,它可以写成:
这样,我们就可以计算出 ys
的所有元素的值,也就是四个频率分量在各个不同时间点下的值与相应
的幅值乘积的和。这就是我们想要的合成结果。
我们可以比较一下两个版本的合成结果是否一致:
ys1 = synthesize1(amps, fs, ts)
ys2 = synthesize2(amps, fs, ts)
max(abs(ys1 - ys2))
可以得出 ys1
和 ys2
最大的差值为 1e-13,这仅仅是浮点数的精度误差。
我们可以用线性代数的形式来描述这个计算过程,如下:
这里 a 表示幅值的向量, t 表示时间向量, f 表示频率向量。 ⊗ 表示外积。
6.3 分解¶
现在我们开始来解决分解的问题。假设给你一个信号并且告诉你这个信号是一组余弦信号合成的。
怎么求出这些余弦信号的幅值呢?换句话说,就是给定 ys
ts
和 fs
,求 amps
。
其实,第一步与合成信号是一样,我们先计算出 \(M = \cos (2\pi t \otimes f)\\\) ,
然后就是找到 a 使得 \(y = Ma\) ,也就是解这个线性方程。Numpy提供了 np.solve
来解决这个问题,代码如下:
def analyze1(ys, fs, ts):
args = np.outer(ts, fs)
M = np.cos(PI2 * args)
amps = np.linalg.solve(M, ys)
return amps
前两行和合成函数是一样,用 ts
和 ys
计算出了 M
,然后使用 np.linalg.solve
计算出了 amps
。
但是这里其实有一个问题,通常解一个线性方程组,需要方程的数量(矩阵的行数)和
未知数的数目(矩阵的列数)相同,也就是说 M
应该是方阵。
但是这个例子中,我们只有四个频率,却有11025个采样时间点,也就是方程的数量比未知数的数量多。
通常情况下,如果 ys
包含的元素超过4个,我们不会只计算4个频率成分。但是在这个例子中,
我们已知了只有这4个频率成分,所有我们可以仅使用 ys
的4个值就可以算出 amps
。
简单起见,我们使用了 ys
的前四个采样值来运行 analyze1
,如下:
n = len(fs)
amps2 = analyze1(ys[:n], fs, ts[:n])
计算结果确为:
[ 0.6 0.25 0.1 0.05 ]
这个计算方法虽然简单但是很慢,因为解线性方程组的时间复杂度为 \(O({n^3})\) , n 为 M 的列数。 接下来我们来优化它。
6.4 正交矩阵¶
求解线性方程组的一个方法是对矩阵求逆。对于方阵 M 来说,逆矩阵表示为 \({M^{ - 1}}\) 使得 \({M^{ - 1}}M = I\) 。 \(I\) 表示单位矩阵,它的对角线元素均为1,其他元素均为0 。
因此,为了求解方程 \(y = Ma\) ,我们将等式两边同时左乘 \({M^{ - 1}}\) 得到:
将右式的 \({M^{ - 1}}M\) 用 \(I\) 代替,得到:
因为任何向量左乘单位向量都等于其本身,所以:
也就是说,如果我们求出 \({M^{ - 1}}\) ,就可以使用一个简单的点积来计算出 a 。 这样的时间复杂度仅为 \(O({n^2})\) ,比 \(O({n^3})\) 要好很多。
但是矩阵求逆也比较慢,但是有一些特殊情况比较快的计算出结果。例如如果 M 是正交矩阵, 那么它的逆就是它的转置,记为 \({M^T}\) 。在Numpy中,转置是一个常量时间复杂度的操作, 因为它并没有实际去改变矩阵的元素的位置,而是生成了一个视图,在我们存取元素的时候,实际上 存取的是相应位置转置后的值。
再次强调一遍,正交矩阵的转置与逆相等, \({M^T} = {M^{ - 1}}\) , 也就意味着 \({M^T}M = I\) 。我们可以使用这个特性来判断一个矩阵是否是正交的。
我们来看看 synthesize2
中的矩阵是不是正交的。在之前的例子中, M 有11025行,
我们需要使用一个更小的版本:
def test1():
amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
time_unit = 0.001
ts = np.arange(N) / N * time_unit
max_freq = N / time_unit / 2
fs = np.arange(N) / N * max_freq
ys = synthesize2(amps, fs, ts)
amps
是与之前一样的幅值向量。我们有4个频率成分,因此需要4个采样点就可以了,以保证 M 为方阵。
ts
是从0到1个单位时间的采样时间点序列,这里我随便选了一个0.001s作为单位时间,实际上无论选择什么
样的单位时间,计算结果都是一样的。
我们在单位时间内采样了N个点,奈奎斯特频率应为 N/单位时间/2 ,这里为2000Hz。 fs
为0~2000Hz的
频率序列。矩阵 M 的计算结果为:
[[ 1. 1. 1. 1. ]
[ 1. 0.707 0. -0.707]
[ 1. 0. -1. -0. ]
[ 1. -0.707 -0. 0.707]]
很容易发现,这个矩阵是对称的,也就是说在 (j,k) 位置上的元素与 (k,j) 位置上的元素是相等的, 也就是它的转置与自身相等: \({M^T} = M\) 。
但不幸的是,它不是正交的,因为 \({M^T}M\) 不是单位矩阵:
[[ 4. 1. -0. 1.]
[ 1. 2. 1. -0.]
[-0. 1. 2. 1.]
[ 1. -0. 1. 2.]]
6.5 DCT-IV¶
但是我们有好几种选择方法,可以使得选择 ts
和 fs
后的 M 是正交的,
因此也就形成了不同的DCT的计算方法。
其中一种方式是将 ts
和 fs
移动半个单位,被称作DCT-IV,“IV”表示罗马数字4,
因为它是8个DCT方法中的第四个。
更新后的代码如下:
def test2():
amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
ts = (0.5 + np.arange(N)) / N
fs = (0.5 + np.arange(N)) / 2
ys = synthesize2(amps, fs, ts)
与之前的代码比较后,你会发现两个改动,一个在 ts
和 fs
上加了0.5的偏移,
二是去掉了 time_unit
从而简化了fs。
结果, M 等于:
[[ 0.981 0.831 0.556 0.195]
[ 0.831 -0.195 -0.981 -0.556]
[ 0.556 -0.981 0.195 0.831]
[ 0.195 -0.556 0.831 -0.981]]
\({M^T}M\) 等于:
[[ 2. 0. 0. 0.]
[ 0. 2. -0. 0.]
[ 0. -0. 2. -0.]
[ 0. 0. -0. 2.]]
由于浮点数精度的原因,其中一些非对焦元素的值被显示为了-0,可以把它当做0来看。 这个矩阵非常接近 \(2I\) 了,也就意味着 M 几乎是正交的,只是多了2倍的因子。 这样已经足够满足我们的需求了。
由于 M 是对称且正交的,因此, M 的逆为 M/2 ,我们可以由此把代码改写为:
def analyze2(ys, fs, ts):
args = np.outer(ts, fs)
M = np.cos(PI2 * args)
amps = np.dot(M, ys) / 2
return amps
我们把 np.linalg.solve
换成了与 M/2 进行点积。
结合 test2
和 analyze2
,我们就实现了DCT-IV的计算:
def dct_iv(ys):
N = len(ys)
ts = (0.5 + np.arange(N)) / N
fs = (0.5 + np.arange(N)) / 2
args = np.outer(ts, fs)
M = np.cos(PI2 * args)
amps = np.dot(M, ys) / 2
return amps
ys
是波形数据,我们不用给定 ts
和 fs
, dct_iv
会通过 ys
的长度来自行计算 ts
和 fs
。
如果之前的推导正确,这个函数应该可以用来对给定的 ys
进行分解求出 amps
。
我们来对它进行一下测试:
amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
ts = (0.5 + np.arange(N)) / N
fs = (0.5 + np.arange(N)) / 2
ys = synthesize2(amps, fs, ts)
amps2 = dct_iv(ys)
max(abs(amps - amps2))
首先,我们预设了一组幅值 amps
,并计算出了 ts
和 fs
,
然后使用 synthesize2
计算出了合成信号,然后用 dct_iv
对合成信号进行分解,
得到 amps2
。最后,将 amps
与 amps2
比较,他们的最大差值仅为1e-16
(由浮点数精度导致),说明计算正确。
6.6 逆DCT¶
最后,我们注意到 analyze2
和 synthesize2
的代码几乎是一样的,唯一的区别在于
analyze2
的结果除了2。因此,我们可以这样来计算DCT的逆运算:
def inverse_dct_iv(amps):
return dct_iv(amps) * 2
inverse_dct_iv
其实就是信号的合成:它把输入幅值向量并输出合成信号 ys
。
我们可以像下面这样测试 inverse_dct_iv
amps = [0.6, 0.25, 0.1, 0.05]
ys = inverse_dct_iv(amps)
amps2 = dct_iv(ys)
max(abs(amps - amps2))
同样,结果最大的差值是1e-16。
6.7 DCT类¶
thinkdsp
中提供了一个 Dct
类对DCT进行封装(类似于 Spectrum
类对FFT进行的封装)。
我们可以通过波形的 make_dct
来生成Dct对象:
signal = thinkdsp.TriangleSignal(freq=400)
wave = signal.make_wave(duration=1.0, framerate=10000)
dct = wave.make_dct()
dct.plot()
这个一个400Hz的三角波的DCT结果,见 图6.2 。DCT的结果可以是正值也可以是负值, 负值代表的是负的余弦,也就相当于移动了180°的相位。
make_dct
内部用的是DCT-II方法来计算DCT,在 scipy.fftpack
中提供这个方法:
import scipy.fftpack
# class Wave:
def make_dct(self):
N = len(self.ys)
hs = scipy.fftpack.dct(self.ys, type=2)
fs = (0.5 + np.arange(N)) / 2
return Dct(hs, fs, self.framerate)
dct
的结果保存在 hs
中,对应的频率值为 fs
(见 6.5 DCT-IV )。
然后使用它们和采样率来生成一个 Dct
对象。
Dct
类也提供了 make_wave
方法来计算逆DCT,我们来测试一下:
wave2 = dct.make_wave()
max(abs(wave.ys-wave2.ys))
wave
和 wave2
的 ys
差别也大概为1e-16(同样是浮点数误差)。
make_wave
使用了 scipy.fftpack.idct
# class Dct
def make_wave(self):
n = len(self.hs)
ys = scipy.fftpack.idct(self.hs, type=2) / 2 / n
return Wave(ys, framerate=self.framerate)
逆DCT默认不会结果进行归一化,因此我们需要除以 2N 。
6.8 练习¶
下面练习的答案可以参考文件 chap06soln.ipynb
。
练习1 之前我说 analyze1
的时间复杂度是 \(O({n^3})\) ,
而 analyze2
的时间复杂度为 \(O({n^2})\) 。使用不同长度的信号作为输入,
运行这两个函数并计时,看看我这个说法对不对。提示:可以使用魔法命令 %timeit
来计时。
如果你将运行时间和输入数据的长度画在一个对数坐标下,应该可以得到一条直线,对于
analyze1
来说直线的斜率为3,而 analyze2
斜率为2。
用同样的方法测试 dct_iv
和 scipy.fftpack.dct
。
练习2 DCT主要应用在音频和图像的压缩中。简单来说,基于DCT的压缩原理是:
- 把长段的数据分段。
- 计算每段的DCT。
- 把幅值很小的频率成分去掉,保存剩下的频率和幅值。
- 解压的时候,将频率和幅值进行逆DCT运算。
把这个算法实现一下,并应用到一段音乐或语音上,看看多少频率成分被去除后,解压后的 声音能够感觉到与原始声音有区别。
为了让这个方法更实用,我们需要存储稀疏矩阵(大部分的元素为0)。 Numpy提供了几种方法,参见 http://docs.scipy.org/doc/scipy/reference/sparse.html 。
练习3 本书的 代码库 中有一个 phase.ipynb
的文件,讨论了相位对于声音的感受的影响。
阅读并运行里面的代码,并找一个其他的录音进行试验。你能找到相位结构与声音感觉之间的关系吗?
第七章:离散傅里叶变换¶
从第一章开始,我们就都在使用DFT来分析信号,但是我们还没有解释 它的工作原理,这章我们就来学习DFT的原理。
通过上一章,你应该理解了DCT的原理,这样理解DFT就很容易了。 它们之间唯一的区别在于DCT使用的是余弦函数,而DFT使用的是复指数函数。 像上一章一样,我会从复指数信号开始,一步一步的解释DFT:
- 首先是合成问题:给定一组频率和幅值,构造它们的合成信号。 合成问题其实就是DFT的逆运算。
- 用矩阵形式重写合成算法。
- 分解问题:给定一个信号,找到构成它的频率分量的幅值和初始相位,即DFT
- 最后我们使用线性代数的知识来提高计算DFT的效率。
这章的代码 chap07.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp07 查看。
7.1 复指数¶
一直以来,推动数学发展的动力很多都来自于对公式的推广,就是将公式的适用范围扩大。 例如,我们通常定义的阶乘只能应用在正整数上。你一定很好奇,如何对一个小数进行阶乘,比如3.5。 由于之前的定义对小数不适用了,我们需要找另外的方法来计算小数的阶乘并且需要同样适用于整数。
1730年欧拉发现了一种在实数和复数域上扩展的阶乘,被称为伽玛函数。 (参见 http://en.wikipedia.org/wiki/Gamma_function)
欧拉还发现了很多类似的扩展数学定义的函数,复指数函数就是其中之一。
通常指数的定义是重复的乘以自身, 例如 \({\varphi ^3} = \varphi \cdot \varphi \cdot \varphi\) 。 但是定义不能应用在非整数指数上。
然后我们发现,指数可以表示成幂级数的形式:
这个定义就可以适用于任何实数,虚数以及复数了。如果我们使用一个纯虚数 \(i\varphi\) 代入得到:
对上式的各项重新排列后,我们可以发现它等价于:
这个公式被称为欧拉公式,详见 http://en.wikipedia.org/wiki/Euler’s_formula 它意味着 \({e^{i\varphi }}\) 是一个模为1的复数,对应于复平面中单位圆上的一点, 而参数 \(\varphi\) 就是这个点的向量与x轴(实轴)之间的夹角。
如果指数是一个复数的话,会得到:
这里 A 表示幅值, \({e^{i\varphi }}\) 是单位复指数,表示相角。
Numpy提供了 exp
方法来计算复指数:
>>> phi = 1.5
>>> z = np.exp(1j * phi)
>>> z
(0.0707+0.997j)
Python中使用 j
来表示单位虚数(就是上面使用的 i )。一个以 j
结尾的数,表示了一个纯虚数,
例如, 1j
其实就是 i 。
当 np.exp
的参数是虚数或者复数的时候,其结果是也是一个复数类型( np.complex128
) ,
实际上就是两个分别代表实部和虚部的浮点数。这个例子中,结果为 0.0707+0.997j
。
复数有两个属性, real
和 imag
,分别代表它的实部和虚部:
>>> z.real
0.0707
>>> z.imag
0.997
可以使用Python内置的 abs
或是Numpy中的 np.absolute
来计算复数的模值:
>>> abs(z)
1.0
>>> np.absolute(z)
1.0
使用 np.angle
可以计算复数的相角:
>>> np.angle(z)
1.5
这个例子证明了 \({e^{i\varphi }}\) 确实是一个模为1,相角为1.5弧度的复数。
7.2 复信号¶
如果 \(\varphi (t)\) 变成时间的函数,那么 \({e^{i\varphi (t)}}\) 也变成了时间的函数。 准确的说,可以写成:
这个式子描述了一个随时间变化的量,因此它是一个信号,准确的说叫做 复指数信号 。
这个例子中,信号的频率是一定的,如果把 \(\varphi (t) = 2\pi ft\) ,那么:
把初始相位 \({\varphi _0}\) 加入后得到:
thinkdsp
中实现了复信号 - ComplexSinusoid
class ComplexSinusoid(Sinusoid):
def evaluate(self, ts):
phases = PI2 * self.freq * ts + self.offset
ys = self.amp * np.exp(1j * phases)
return ys
ComplexSinusoid
类继承了父类 Sinusoid
的构造函数,并复写了 evaluate
以计算信号的值,
它与 Sinusoid
相比只是用 np.exp
代替了 np.sin
。
一个复信号的值是复数的Numpy数组:
>>> signal = thinkdsp.ComplexSinusoid(freq=1, amp=0.6, offset=1)
>>> wave = signal.make_wave(duration=1, framerate=4)
>>> wave.ys
[ 0.324+0.505j -0.505+0.324j -0.324-0.505j 0.505-0.324j]
这个信号的频率是1Hz,幅值是0.6,初始相位是1弧度。例子中计算了从0~1s的四个采样点的值, 结果都是复数。
7.3 合成¶
与实的正弦信号一样,我们可以将不同频率和幅值的复信号加起来生成一个复合信号, 也就是复信号的合成问题:给定不同频率和幅值的复信号,计算合成信号。
最简单的方法就是生成 ComplexSinusoid
并把他们加起来:
def synthesize1(amps, fs, ts):
components = [thinkdsp.ComplexSinusoid(freq, amp)
for amp, freq in zip(amps, fs)]
signal = thinkdsp.SumSignal(*components)
ys = signal.evaluate(ts)
return ys
这个函数与 6.1 中的 synthesize1
函数几乎一样,
唯一的区别是把 CosSignal
换成了 ComplexSinusoid
。
下面是使用这个函数进行合成的例子:
amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = [100, 200, 300, 400]
framerate = 11025
ts = np.linspace(0, 1, framerate)
ys = synthesize1(amps, fs, ts)
结果为:
[ 1.000 +0.000e+00j 0.995 +9.093e-02j 0.979 +1.803e-01j ...,
0.979 -1.803e-01j 0.995 -9.093e-02j 1.000 -5.081e-15j]
看起来复信号就是一组复数值的序列,我们应该怎样来理解它呢?对于实数的信号, 我们可以很直观的理解为随时间变化的一个量,例如声音信号就是随时间变化的空气压力值。 然后,对于复信号,我们在现实中却找不到与之对应的物理量。
究竟什么是复信号?对于这个问题,我们还没有找到很直观的答案,但是有两种观点:
- 复信号是对信号进行计算和分析时的一种数学抽象,在现实中不能直接对应于任何实际的信号。
- 你可以把复信号想象为彼此独立的两个分别代表实部和虚部的实数信号。
按照第二种观点,我们可以把之前的复信号分解为实部和虚部两个信号:
n = 500
thinkplot.plot(ts[:n], ys[:n].real, label='real')
thinkplot.plot(ts[:n], ys[:n].imag, label='imag')
如 图7.1 ,实部是余弦信号的合成结果,虚部是正弦信号的合成结果。 虽然他们的波形看起来不一样,但是他们包含同样的频率成分,并且听起来也是一样的, 因为他们的频率成分仅仅相差了90°,而对于常值的相位差,我们人耳是分辨不出区别的。
7.4 使用矩阵进行合成¶
我们可以使用类似 6.2 的方法将之前的合成函数改写成矩阵的形式:
PI2 = np.pi * 2
def synthesize2(amps, fs, ts):
args = np.outer(ts, fs)
M = np.exp(1j * PI2 * args)
ys = np.dot(M, amps)
return ys
这里的 amps
是以Numpy数组保存的幅值序列。
fs
是不同的频率成分, ts
是计算波形的采样时间序列。
args
计算了 ts
和 fs
的点积,参见 6.2 。
M
的每一列都是在某个频率的复指数信号在 ts
下采样的值。
把 M
与 amps
进行点积的结果就是不同频率的复指数信号在 ts
下的采样值之和,即我们想要的合成结果。
我们再来测试一下:
>>> ys = synthesize2(amps, fs, ts)
>>> ys
[ 1.000 +0.000e+00j 0.995 +9.093e-02j 0.979 +1.803e-01j ...,
0.979 -1.803e-01j 0.995 -9.093e-02j 1.000 -5.081e-15j]
这与上一节的结果是一样。
这个例子中,幅值是实数。实际上,幅值也可以是复数。复数的幅值是什么效果呢?记得我们之前所说的, 一个复数可以是一个实数加一个虚数,也可以是一个实数乘以一个复指数,这样一来,复数的幅值就可以 表示成:
也就是说,复数 \(A{e^{i{\varphi _0}}}\) 的幅值表示的是信号幅值为 A , 初始相位为 \({\varphi _0}\)
为了验证这一结论,我们在之前的例子中加入了复数的幅值:
phi = 1.5
amps2 = amps * np.exp(1j * phi)
ys2 = synthesize2(amps2, fs, ts)
thinkplot.plot(ts[:n], ys.real[:n])
thinkplot.plot(ts[:n], ys2.real[:n])
amps
乘以 np.exp(1j * phi)
就得到了模为 amps
相角为 phi
的复数。
信号的波形如 图7.2 所示,因为加入了 \({\varphi _0} = 1.5\) 的初始相位,相当与
对每个频率分量都移动了大约半个周期。由于每个频率分量的周期都不一样,导致了每个分量移动的
时间长度也不一样,于是合成的信号波形也不一样了。
现在,我们有了一个可以处理复数幅值的合成函数,接下我们该学习信号的分解了。
7.5 分解¶
分解是合成的逆运算:给定一个信号的采样序列 y ,以及组成信号的不同分量的频率, 我们需要计算出这些分量的复数幅值 a 。
与 6.3 一样,我们可以通过求解线性方程 \(Ma = y\) 来得到 a
def analyze1(ys, fs, ts):
args = np.outer(ts, fs)
M = np.exp(1j * PI2 * args)
amps = np.linalg.solve(M, ys)
return amps
analyze1
输入复数信号序列 ys
,以及一系列的频率 fs
和采样时间序列 ts
,
然后使用 np.linalg.solve
计算并返回了结果 amps
。
使用上一节合成的信号,我们可以测试这个 analyze1
运行的结果是否正确。
由于求解线性方程需要 M
为方阵,因此我们把 ys
fs
和 ts
都设置为相同的长度,
都等于 fs
的长度,代码如下:
>>> n = len(fs)
>>> amps2 = analyze1(ys[:n], fs, ts[:n])
>>> amps2
[ 0.60+0.j 0.25-0.j 0.10+0.j 0.05-0.j]
可以看到,结果与之前的幅值是一致的(由于浮点数精度误差导致了有很小的虚数部分)。
7.6 高效的分解¶
不幸的是,解线性方程的过程很慢。对于DCT来说,我们通过选择 fs
和 ts
使得
M
是正交矩阵来提高计算的速度,也就是 M
的逆与它的转置相同,利用矩阵乘法来
代替了解方程的过程。
这里对于DFT,我们也可以使用同样的手法,但是又稍有不同。由于 M
是复数矩阵,
因此我们希望它是 酉矩阵(unitary) 而非正交矩阵,即 M
的逆为它的共轭转置,
就是将 M
转置并将所有元素的虚数部分取相反数(复共轭)。
参见 http://en.wikipedia.org/wiki/Unitary_matrix 。
Numpy中提供了 conj
和 transpose
函数来完成复共轭和转置两个操作,
下面的代码计算了 N=4 的情况下的 M
N = 4
ts = np.arange(N) / N
fs = np.arange(N)
args = np.outer(ts, fs)
M = np.exp(1j * PI2 * args)
如果 M 是酉矩阵,那么根据酉矩阵的定义,我们可以得出 \({M^ * }M = I\) , 这里的 \({M^ * }\) 表示 M 的共轭转置。我们可以利用这个性质来测试 M 是否是酉矩阵:
MstarM = M.conj().transpose().dot(M)
不考虑浮点数精度误差的话,结果等于 4I ,也就是说 M 是有一个因子 N 的酉矩阵。 这与在DCT类似,那里的因子是2。
使用这个性质,我们将之前的算法改写为:
def analyze2(ys, fs, ts):
args = np.outer(ts, fs)
M = np.exp(1j * PI2 * args)
amps = M.conj().transpose().dot(ys) / N
return amps
然后来测试一下:
N = 4
amps = np.array([0.6, 0.25, 0.1, 0.05])
fs = np.arange(N)
ts = np.arange(N) / N
ys = synthesize2(amps, fs, ts)
amps3 = analyze2(ys, fs, ts)
结果又一次证明了算法的正确性:
[ 0.60+0.j 0.25+0.j 0.10-0.j 0.05-0.j]
7.7 DFT¶
由于 analyze2
需要在特定的 ts
和 fs
下才能工作,
使用起来很不方便,于是我想把 ts
fs
的选择过程放到函数里面,
重写为仅仅输入 ys
就可以自己计算 ts
和 fs
并得到结果。
首先,我定义了一个计算矩阵 M 的函数:
def synthesis_matrix(N):
ts = np.arange(N) / N
fs = np.arange(N)
args = np.outer(ts, fs)
M = np.exp(1j * PI2 * args)
return M
然后实现了这个 analyze3
def analyze3(ys):
N = len(ys)
M = synthesis_matrix(N)
amps = M.conj().transpose().dot(ys) / N
return amps
analyze3
和DFT已经很接近了,唯一的区别是通常DFT的定义中不会除以 N
def dft(ys):
N = len(ys)
M = synthesis_matrix(N)
amps = M.conj().transpose().dot(ys)
return amps
现在,我们可以用这个函数与Numpy中的 np.fft
进行对比测试:
>>> dft(ys)
[ 2.4+0.j 1.0+0.j 0.4-0.j 0.2-0.j]
结果等于 amps*N
。下面是使用 np.fft
计算的结果:
>>> np.fft.fft(ys)
[ 2.4+0.j 1.0+0.j 0.4-0.j 0.2-0.j]
它们的结果是相同的。
逆DFT基本上也是一样的道理,只是不是对 M 进行共轭转置。 而且,由于之前没有除以 N ,这里就需要把除以 N 加上:
def idft(ys):
N = len(ys)
M = synthesis_matrix(N)
amps = M.dot(ys) / N
return amps
最后,我们来把这两个函数联合起来测试一下,确保 dft(idft(amps))
的结果为 amps
>>> ys = idft(amps)
>>> dft(ys)
[ 0.60+0.j 0.25+0.j 0.10-0.j 0.05-0.j]
如果回到以前让我来定义DFT的话,我宁愿把除以 N 加在DFT的定义中,而逆DFT中不加, 那样与我们解释的原理更一致些。或者两个定义里面都除以 \(\sqrt N\) , 这样DFT和IDFT的形式就更一致了。当时可惜那是不可能的,因此这里我沿用了 最常规的DFT定义,虽然在形式上有点奇怪,但是并不影响我们使用。
7.8 DFT的周期性¶
上一节中我们将DFT表示成了矩阵乘法:先计算矩阵 M 然后计算它的酉矩阵 \({M^ * }\) ,
然后将它与信号 y
相乘。这一系列的步骤,如果展开来看的话,可以写成:
这里的 k 是频率分量的序号0~N-1, n 是采样时间序列的序号,同样也是0~N-1. DFT(y)[k] 就是对y进行DFT结果的第 k 个值。
通常我们只需计算k为0~N-1的值就可以了,因为之后的值都是循环重复的,即第 k 个值 与第 k+N 个值是相同的,同样与 k+2N 和 k-N 的值也是相同的。
为了证明这一点,我们可以把 k+N 代入上式中:
于是可以分解为:
上式中最后一项,指数均为 \({2\pi }\) 的整数倍,因此它恒等于1,于是:
可以看到,这个式子就等于 DFT(y)[k] 。因此,我们说DFT是周期性的,它的周期是 N 。 在本章后面的练习中会让你实现快速傅立叶变换算法(FFT),那里就要用到这个性质。
值得一提的是,将DFT写成上面的和式的形式,可以帮助我们加深对它的理解。 我们回过头来看 6.2 中的框图,你会发现合成矩阵 M 中的每一列 都是不同频率的信号在一个采样时间序列下的值。而分解矩阵就是 M 的共轭转置, 它的每行代表了一个信号在一个采样时间序列下的值。
因此,和式正好是信号 y
与一个复指数信号的互相关,也就是说DFT的每一个元素代表了
信号 y
与不同频率的复指数信号的相关性。
7.9 实信号的DFT¶
thinkdsp
中的频谱类使用了 np.fft.rfft
来进行DFT计算,它实际上做的就是实信号的DFT,
也就是说它仅可以处理值为实数的信号。前面我们讲到的DFT都是针对复数信号的(完整DFT),
那么当我们使用完整DFT来计算实信号的DFT的时候,会怎么样呢?我们来看一个例子:
signal = thinkdsp.SawtoothSignal(freq=500)
wave = signal.make_wave(duration=0.1, framerate=10000)
hs = dft(wave.ys)
amps = np.absolute(hs)
这里我们产生了一个500Hz的锯齿波信号,采样率为10kHz。 hs
计算了波形的完整DFT,
结果为 amps
,它包含了每个频率成分的幅值。但是这些幅值对应的频率究竟是多少呢?
从 dft
的实现上,我们可以看到:
fs = np.arange(N)
由于 dft
并不知道实际的采样率是多少,因此 fs
并不能代表实际的频率值。
实际上,dft
中假设波形是在1个单位时间内产生的,因此采样率是 N/单位时间 ,
为了得到频率值,我们需要将单位时间转换为秒:
fs = np.arange(N) * framerate / N
这样转换之后,我们得到了从0到10kHz(采样率)的频率值,我们将它画成频谱图:
thinkplot.plot(fs, amps)
thinkplot.config(xlabel='frequency (Hz)',
ylabel='amplitude')
如 图7.3 ,左半部分就是我们希望的结果:基频在500Hz,然后谐波分量的幅值 以 1/f 的速度逐渐衰减。
但是图的右半部分却很奇怪,从5000Hz开始,谐波的幅值开始逐渐增大, 直到9500Hz达到峰值。这个为什么呢?
答案是混叠。在 2.3 混叠 中,我们说过10000Hz的采样率,折叠频率是5000Hz,
这种情况下,5500Hz的频率和4500Hz的频率在采样后是没有区别的。当我们计算在5500Hz
频率下的DFT的时候,就会得到与4500Hz下相同的结果。同样的道理,6000Hz和4000Hz等
也是一样的。
也就是说,实信号的DFT是相对于折叠频率对称的。因此,我们仅计算一半的值就可以了,
这样可以节省一半的时间, np.fft.rfft
就是这样做的。
译者注
译者认为用混叠来解释对称的现象并不合适。在实际使用DFT的时候,我们通常并不会 选择频率范围为0~采样率,而是选择 -折叠频率 ~ +折叠频率。引入了负频率之后, 容易证明,正负两个复为相反数的频率的信号之和为实数的条件是,他们的幅值(复数) 应该是共轭的,就是模相等,相角互为相反数。因此,DFT后幅值是偶对称,而初始相位 是奇对称的。
7.10 练习¶
下面练习的答案可以参考文件 chap07soln.ipynb
。
练习1 chap07.ipynb
中包含了一些额外的实例和说明,阅读并运行这些代码。
练习2 本章中我解释了使用矩阵乘法来进行DFT以及IDFT的方法。 这种算法的时间复杂度是 \(O({n^2})\) , n 为波形的长度。虽然在某些情况下, 这个速度以及足够快了,但是其实还有一个更快的算法,就是快速傅立叶变换(FFT)算法。 FFT的时间复杂度是 \(O(n\log n)\) 。FFT算法的关键是Danielson-Lanczos定理:
上式中, DFT(y)[n] 表示 y 的DFT的第n个元素值, e 表示 y 中偶数位置上的 元素, o 表示 y 中奇数位置上的元素。 这个定理表明,DFT可以分解为奇偶两个部分的DFT计算(分治),这样就可以递归的来计算 整个DFT:
- 给定波形
y
,将其分成奇偶两个部分e
和o
。 - 递归的计算
e
和o
的DFT。 - 根据上式把两项相加得到最后的DFT结果。
递归的结束条件可以是波形的长度为1的时候,这个时候 DFT(y)=y 。 也可以选择在长度很小的时候结束,通过矩阵乘法来计算最后这个DFT结果。
提示:我建议你从一个非递归的版本开始一步步的来实现这个算法,在上面的第2步,
用 dft
( 7.7 DFT ) 或 np.fft.fft
来代替递归调用。然后,确保
第三步的结果是正确的。然后把递归的结束条件加上,并将第2步替换成递归调用。
提示:记住DFT是周期性的,你会发现 np.tile
是很有用的。
关于FFT的更多知识,可以参考 https://en.wikipedia.org/wiki/Fast_Fourier_transform 。
第八章:滤波和卷积¶
本章我会介绍信号处理中的一个重要的定理-卷积定理(Convolution Theorem)。 不过在这之前,我们需要先了解一下卷积这个概念。 这里,我们先以一个简单的平滑作为例子开始我们的学习。
这章的代码 chap08.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp08 查看。
8.1 平滑¶
平滑是一种试图消除信号短期的变化的操作,信号经过平滑操作后,可以更好的 揭示长期的变化趋势。例如,如果把股票的每日价格变化画出来,它会看起来像噪声一样, 但是经过平滑之后,我们就可以比较清楚的看到他是上升的趋势还是下降的趋势了。
一个最普通的平滑操作就是移动平均法。它将信号之前 n 个时刻的值的平均值作为 新的信号的第 n 个值。
例如, 图8.1 展示了Fackbook的股票从2012年5月17日到2015年11月8日的每日收盘价格曲线。灰色的曲线代表原始数据, 而蓝色的曲线则是进行了30日移动平均之后的数据。可见,平滑后,信号的短期变化被很大程度的消除了, 这样一来我们也就更容易的看出数据的长期变化趋势。
平滑操作也能用于声音信号。我们先从一个440Hz的方波信号为例,如 2.2 中所述,方波的谐波衰减的 比较慢,它包含很多高频的成分。
我们先生成一个方波信号和它的两段波形:
signal = thinkdsp.SquareSignal(freq=440)
wave = signal.make_wave(duration=1, framerate=44100)
segment = wave.segment(duration=0.01)
wave
是1s的波形, segment
是一个较短的波形(0.01s),我会用它来作图。
为了计算信号移动平均值,我会使用一个类似 3.7 的窗口。在那一章,我们使用了汉明窗来避免 因为信号前后不连续引起的频谱泄露。这里,我们使用窗函数来计算一段信号的加权和。
上例中,我们可以通过一个包含相等且和为1的11个元素的窗,来计算移动平均:
window = np.ones(11)
window /= sum(window)
然后将这个窗应用到波形上:
ys = segment.ys
N = len(ys)
padded = thinkdsp.zero_pad(window, N)
prod = padded * ys
sum(prod)
padded
将窗函数尾部添加0来保证其长度与 segment.ys
一致。这种补0的方法,称为 padding 。
prod
是波形数据与窗函数的乘积, 他们的和其实就是窗口内前11个元素的平均值。
我们需要把窗口滚动到下一个位置。这个例子中,前11个元素都是-1,因此平均值也是-1。
为了计算移动平均的下一个值,我们需要将窗口向右移动一个位置,然后再进行同样的计算:
rolled = np.roll(rolled, 1)
prod = rolled * ys
sum(prod)
使用同样的方法,我们可以计算剩余的所有元素,下面的代码将之前的代码都放到了一个循环中, 这样它就可以循环的处理整段信号并把结果放入了数组中:
def smooth(ys, window):
N = len(ys)
smoothed = np.zeros(N)
padded = thinkdsp.zero_pad(window, N)
rolled = padded
for i in range(N):
smoothed[i] = sum(rolled * ys)
rolled = np.roll(rolled, 1)
return smoothed
smooth
就是用来保存结果的数组, padded
就是补0后的窗函数, rolled
是 padded
的一份拷贝,它在每次循环结束后都会右移一个元素。在循环内,我们将 ys
和 rolled
相乘并求和,
然后依次放入 smooth
数组中。
结果见 图8.2 ,其中灰线是原始信号,蓝线是平滑后的信号。可见,平滑后的信号在原始信号突变的时候 是以斜坡的方式变化的,也就是信号的变化没有那么“尖锐”了。
8.2 卷积¶
上一节中,我们运用窗函数对信号逐段进行的操作就叫做 卷积(convolution)
卷积是一个经常使用的操作,Numpy中提供了一个更加简单和快速的实现:
convolved = np.convolve(ys, window, mode='valid')
smooth2 = thinkdsp.Wave(convolved, framerate=wave.framerate)
np.convolve
计算了波形和窗函数的卷积。参数 valid
表示只计算窗函数和波形
完全重叠没有交叉的部分,因此它从窗函数与波形数据左对齐开始,直到窗口移动到与波形右对齐。
计算结果与 图8.2 是完全相同的。
严格的来说,之前的代码和它还是有一定的区别,之前的代码计算的是互相关函数:
这里 f 代表了长度为 N 的波形数据, g 是窗函数, ⋆ 是互相关的算符。 当计算第n个值的时候,实际上我们需要将 g 右移 n 个位置,这也是 g[n+m] 的意义。
而卷积的定义是:
\(*\) 是卷积算符,这个式子与之前那个的区别在于 g 的 m 是负的, 也就是 g 的值是反过来的,计算的时候也应该反过来。由于上一个例子中, 窗函数是对称的,因此它们的结果才会一样。如果使用其他的窗函数,就必须要注意这一点。
为什么要像这样来定义卷积呢?有两个原因:
- 这个定义在很多场合下都是很自然而然的得出来的,尤其是在分析信号处理系统的时候, 我们会在第十章的时候学习。
- 这个定义是卷积定理的基础,我们马上就会学习。
之后我们还会学习到循环卷积。
8.3 频域¶
在时域上,平滑操作会使信号的变化没有那么剧烈,那么在频域上,有什么变化呢? 让我们先来看看原始信号的频谱:
spectrum = wave.make_spectrum()
spectrum.plot(color=GRAY)
再来看看平滑后的信号频谱:
convolved = np.convolve(wave.ys, window, mode='same')
smooth = thinkdsp.Wave(convolved, framerate=wave.framerate)
spectrum2 = smooth.make_spectrum()
spectrum2.plot()
same
模式表示计算结果需要和输入的长度一致。虽然,这样会产生一些不太好的值,
不过不会应该我们之后的分析。
如 图8.3 所示,基频的幅值几乎没有变化,而接下来几个谐波分量的幅值都有所衰减。 再往高频的分量就几乎被消除了。因此,平滑操作具有低通滤波的效果,见 1.5 频谱 以及 4.4 粉红噪声 。
为了知道各个频率成分衰减的程度,我们计算出了两个频谱分量的比例:
amps = spectrum.amps
amps2 = spectrum2.amps
ratio = amps2 / amps
ratio[amps<560] = 0
thinkplot.plot(ratio)
ratio
就是平滑前和平滑后的幅值比例。 当 amps
很小的时候,这个比例会很大并且没什么意义,
因此我们简单的把它设置为了0 。
如 图8.4 ,在低频的时候 ratio
是比较大的,接近1,然后逐渐减小直到截止频率4000Hz。
但是,在4000Hz之后, ratio
又开始在0~0.2之间变动,这是怎么回事呢?
8.4 卷积定理¶
答案就是卷积定理,它的数据公式如下:
式中, f 是波形数据, g 是窗函数,也就是说根据卷积定理, f 与 g 的卷积后的信号的DFT, 与它们分别进行DFT然后再相乘的结果是一样的。
由于波形是随时间变化的函数,因此卷积操作是在 时域(time domain) 上进行的。 而DFT后相乘是在 频域(frequency domain) 进行的,因为DFT的结果是频率的函数。 因此,我们又可以这样来描述卷积定理:
- 时域的卷积相当于频域的乘积
这也解释了 图8.4 的结果。当我们把波形和窗函数进行卷积,实际在频域上就是波形的频谱和窗函数频谱的 乘积。为了证明这一点,我们先计算出了窗函数的DFT:
padded = zero_pad(window, N)
dft_window = np.fft.rfft(padded)
thinkplot.plot(abs(dft_window))
padded
是补0后的窗函数, dft_window
是窗函数的DFT。
结果见 图8.5 ,平滑前后的幅值比例与窗函数的DFT结果 dft_window
是完全重合的。
用数学公式表示是:
在频域下,窗函数的DFT就称为 滤波器(filter) ,也就是说时域上一个窗函数的卷积,就对应了 频域下的一个滤波器。
8.5 高斯滤波器¶
上一节使用的移动平均窗是一个低通滤波器,但是滤波的效果并不是很好。因为它的DFT一开始衰减的很快, 但是后面开始上下震荡,这种情况被称为 旁瓣(sidelobes) 。由于移动平均窗实际上是一个方波 信号,它的频谱上包含高频的谐波成分并且幅值是按照 1/f 衰减的,这就比较慢了。
而高斯窗函数就比较好了。Scipy中提供了很多常用的卷积窗函数,其中就包括高斯窗:
gaussian = scipy.signal.gaussian(M=11, std=2)
gaussian /= sum(gaussian)
这里, M
是窗的长度, std
是高斯分布函数的标准差参数。 图8.6 展示了这个窗函数的形状,
它是高斯分布钟型曲线的离散化,图中还画出了移动平均窗的曲线,实际上就是一条平的直线,由于看起来像是
矩形,因此又称为矩形窗。
我用高斯窗函数重新进行了一次和上一节相同的计算,结果如 图8.7 ,可见通过高斯窗进行平滑后的 频谱衰减比例与高斯窗的DFT是一致的。
对于低通滤波的效果来说,高斯平滑要比移动平均要更好一些。图中,随着幅值的衰减,几乎没有旁瓣出现,因此它 可以更好的去除信号的高频成分。
我们之前提到过,高斯函数的DFT也是高斯函数,它的频谱是以 \({e^{ - {f^2}}}\) 衰减的,因此,它比 1/f 的矩形窗要好很多。
8.6 高效的卷积¶
FFT的一个重要的运用是,它与卷积定理结合起来可以提供一种高效的计算卷积,互相关和自相关函数的方法。
这里,我们重新写一次卷积定理的公式:
因此,我们得到一种计算卷积的方法:
IDFT
指的是逆DFT变换。我们之前使用的卷积计算方式,时间复杂度是 \(O({n^2})\) ;
而使用这种方式,时间复杂度仅为 \(O(n\log n)\) 。
为了验证这种方式的正确性,我们将分别使用两种方法来计算 图8.1 中Fackbook股票数据的卷积结果:
import pandas as pd
names = ['date', 'open', 'high', 'low', 'close', 'volume']
df = pd.read_csv('fb.csv', header=0, names=names)
ys = df.close.values[::-1]
上面的代码中使用了Pandas来从csv文件中读取数据,这个文件包含在 代码库 中。 Pandas是一个数据处理库,本书中涉及的不多,因此即使你不太熟悉也没关系。 如果有兴趣和话,可以阅读 Think Stats at http://thinkstats2.com 。
从csv读入的数据保存在 df
中,它是一个 DataFrame
对象(Pandas中的一个数据类),
close
是包含收盘数据的Numpy数组。
接下来,我将高斯窗应用到了这个数据上:
window = scipy.signal.gaussian(M=30, std=6)
window /= window.sum()
smoothed = np.convolve(ys, window, mode='valid')
下面的 fft_convolve
使用FFT来进行卷积:
from np.fft import fft, ifft
def fft_convolve(signal, window):
fft_signal = fft(signal)
fft_window = fft(window)
return ifft(fft_signal * fft_window)
我们将同样的高斯窗补0后,使用 fft_convolve
来进行同样的计算:
padded = zero_pad(window, N)
smoothed2 = fft_convolve(ys, padded)
smooth2
的开头包含 M-1 个无效的数据, M 是窗长度,我们可以像下面这样来去掉这些数据:
M = len(window)
smoothed2 = smoothed2[M-1:]
经过比较, fft_convolve
和 np.convolve
的结果是相同的。
8.7 高效的自相关¶
在 8.2 卷积 中我提到了互相关的定义,它和卷积之间仅有一个符号的差别。
上一节我们介绍了计算卷积的高效方法,同样,我们也可以高效的计算互相关和自相关。 我们继续使用Facebook的股票数据来计算它的自相关:
corrs = np.correlate(close, close, mode='same')
使用 same
模式使得结果与输入数据 close
的长度一致,相应的 lag
值是从 -N/2 到 N/2-1 。
计算结果见 图8.8 中的灰色曲线。除了 lag
为0的时候以外,没有其他的峰值,表明这个信号没有明显的
周期性。然而,自相关函数下降的比较慢,说明这个信号类似于 5.3 中的粉红噪声。
使用卷积来计算自相关之前,我们需要把信号补0到两倍信号的长度,这是因为FFT假定信号值周期的, 而进行自相关的信号使用这个假设是不成立的,因此我们把信号补0,最后再从结果中去掉无效的元素。
考虑到卷积和自相关的公式差一个负号,我们使用 np.flipud
把窗函数(本身)进行反向,再
调用 fft_convolve
计算卷积,最后的结果就是自相关了。 np.flipud
操作并没有真的对
内存进行反向的操作,只是产生了一个反向的视图,因此它的速度是很快的:
def fft_autocorr(signal):
N = len(signal)
signal = thinkdsp.zero_pad(signal, 2*N)
window = np.flipud(signal)
corrs = fft_convolve(signal, window)
corrs = np.roll(corrs, N//2+1)[:N]
return corrs
fft_convolve
的结果长度是2N,其中前 N/2 个元素和后 N/2 个元素是有效的,其余都是由于补0
导致的无效数据,于是我们将它右移并取出了前 N 个元素,也就对应了 lag
从 -N/2 到 N/2-1 的值。
结果如 图8.8 ,两种计算方法得到的结果是相同的。
需要注意的是, 图8.8 中自相关的值是很大的,我们需要将其进行归一化(使值的范围在-1到1之间), 见 5.6 。
计算互相关也可以使用同样的方法。首先将两个信号信号补0然后将其中一个反向,然后使用上述方法计算卷积,最后再 取出有效的元素。补0和截取有效数据的操作是比较繁琐的,因此Numpy中提供了一些函数来帮助我们完成这项操作。
8.8 练习¶
下面练习的答案可以参考文件 chap08soln.ipynb
。
练习1 阅读并运行 chap08.ipynb
中的代码。里面有一个交互性的实验来帮助你了解高斯窗中参数的变化
对于截止频率的影响。
当你增大高斯窗的 std
而不改变窗函数的长度 M 的时候,会出现什么问题?
练习2 之前我提到了高斯曲线的傅立叶变换同样是高斯曲线,对于DFT来说,这个结论也是近似正确的。
尝试一下当你改变 std
后,它的DFT有什么变化。
练习3 在第三章中,我们学到了汉明窗和其他的一些窗函数,用于抑制频谱泄露,并且我们通过他们的DFT结果来 了解了它们各自不同的抑制效果。
与之前使用的高斯窗一样,生成一个同样长度汉明窗,补0后画出它的DFT,看看这两种窗函数的低通滤波效果更好? 你会发现使用y轴对数坐标来画DFT是很有用的。
同样的,使用不同的窗函数和长度来进行这个实验。
第九章:微分和积分¶
这一章中,我们会继续分析时域中的窗函数和频域中的滤波的关系。 更具体的说,我们会分析有限差分(近似于微分)以及累加(近似于积分)两种操作的作用。
这章的代码 chap09.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp09 查看。
9.1 有限差分¶
在 8.1 平滑 中,我们对Fackbook的股价数据运用了一个平滑操作,发现在时域上平滑窗的作用 等价于在频域上的低通滤波。
在这一节中,我们来尝试计算它们每日价格差,你会发现在频域上,这相当于高通滤波的作用。
下面的代码从文件中读取数据保存到波形对象中并计算了它的频谱:
import pandas as pd
names = ['date', 'open', 'high', 'low', 'close', 'volume']
df = pd.read_csv('fb.csv', header=0, names=names)
ys = df.close.values[::-1]
close = thinkdsp.Wave(ys, framerate=1)
spectrum = wave.make_spectrum()
代码使用Pandas读取了数据,结果为 DataFrame
对象 df
,它的每一列代表了一天的开盘价,收盘价和最高最低价。
我选择了收盘价作为我们分析的数据,并把它保存到了波形对象中,这里的采样率为每天1次。
图9.1 展示了这些数据的波形和频谱图。直观上看,这个波形类似 4.3 中的布朗噪声。 它的频谱看起来近似一条带有噪声的直线,斜率大约是-1.9, 与布朗噪声一样,它的斜率是常值。
接下来,我们使用 np.diff
来计算它的每日价格变化情况:
diff = np.diff(ys)
change = thinkdsp.Wave(diff, framerate=1)
change_spectrum = change.make_spectrum()
图9.2 展示了结果的波形和频谱。可见,每日的变化曲线类似于白噪声,频谱近似一条平的直线,斜率大概是-0.06,接近于0, 这和白噪声是一样的。
9.2 频域¶
实际上,计算相邻元素的差值与差分窗 [1,-1] 的卷积是一样的。你有可能会觉得这两个元素的值反了,实际上没有, 记住卷积的操作会先对窗函数进行反向的操作。
我们计算出了这个窗的DFT来看看它在频域上的作用效果:
diff_window = np.array([1.0, -1.0])
padded = thinkdsp.zero_pad(diff_window, len(close))
diff_wave = thinkdsp.Wave(padded, framerate=close.framerate)
diff_filter = diff_wave.make_spectrum()
如 图9.3 左图所示,有限差分窗对应了一个高通滤波器,它的幅值随着频率的增大而增大, 在低频的时候是线性的,而随着频率增大,曲线慢慢变平。下一节中,我们会分析这个原因。
9.3 微分¶
上一节我们用到的差分窗其实是一阶微分的数值近似,因此滤波效果和微分也是近似的。
时域的微分在频域上也相当于一个滤波器,下面我们会用一些数学推导来说明这一点。
假设我们有一个频率为 f 的复指数信号:
它的一阶微分是:
这个式子又可以写成:
也就是说,对复指数信号进行微分等价于乘以 \(2\pi if\) , 这个一个模为 \(2\pi f\) 相角为 \(2\pi\) 的复数。
因此,我们可以构造如下的滤波器来代替微分运算:
deriv_filter = close.make_spectrum()
deriv_filter.hs = PI2 * 1j * deriv_filter.fs
为了和之前的数据长度和采样率保持一致,这里我先使用了 close
的频谱数据,
然后将它的 hs
替换为了 \(2\pi if\) 。 如 图9.3 左图所示,
它是一条直线。
如 7.4 所述,复指数信号乘以一个复数实际上产生了两种效果: 一是乘以一个幅值,这里是 \(2\pi f\) ;二是移动一个初始相位,这里是 \(2\pi\) 。
如果你了解特征函数的话,复指数函数实际上就是微分算子的特征函数, 对应的特征值是 \(2\pi if\) 。参见 http://en.wikipedia.org/wiki/Eigenfunction 。
如果你对此不熟悉,下面我会简单的解释一下它的含义:
- 算子是指从一个函数到另一个函数的映射关系。例如微分就是一个算子。
- 对于算子 A 和函数 g ,如果将 A 应用到 g 的结果 与 g 本身乘以一个标量值 \(\lambda\) 相等,即 \(Ag = \lambda g\) , 那么我们称 g 是 A 的特征函数。
- 相对应的,我们称标量 \(\lambda\) 为特征函数 g 的特征值。
- 一个给定的算子可以有多个特征函数,每个特征函数均对应一个特征值。
由于复指数信号是微分算子的特征函数,因此我们可以很容易的通过乘以一个复数来计算它的微分。
而对于包含多个频率成分的信号来说,要稍微复杂一些:
- 将信号表示成不同频率的复指数信号之和。
- 计算每个频率成分的微分(乘法)
- 将不同频率的微分结果累加
上面的过程看起来和 8.6 中卷积的算法是一样的,先计算DFT,然后运用一个滤波器, 然后再进行IDFT。
spectrum
类中提供了一个方法来计算差分滤波:
# class Spectrum:
def differentiate(self):
self.hs *= PI2 * 1j * self.fs
我们可以使用这个方法来计算Facebook数据的微分:
deriv_spectrum = close.make_spectrum()
deriv_spectrum.differentiate()
deriv = deriv_spectrum.make_wave()
图9.4 比较了使用差分 np.diff
和微分分别计算每日价格变化的曲线。
为了让结果显示的更清晰,这里我们截取了前50个元素来作图。
图中可以看出,微分的结果更接近于噪声,因为它的高频分量的幅值更大一些,见 图9.3 左图。 而且它的前面几个值比后面噪声的程度更大,这是由于使用DFT的微分是基于周期性假设的,它将信号 首尾相连而导致了在边界上的不连续。
总结一下,我们展示了:
- 计算相邻元素值之间的差分的方法,它又可以表示成对信号进行一个简单的卷积操作。最后的结果为一阶微分的近似值。
- 时域上的微分相当于频域上的一个滤波器。对于周期信号,它的结果是一阶微分,而对于非周期信号,它的结果近似于一阶微分。
使用DFT来计算微分是求解微分方程的频谱方法的基础,详见 http://en.wikipedia.org/wiki/Spectral_method 。 对于分析线性时不变系统,这个方法特别有用,我们会在第十章中进行介绍。
9.4 积分¶
上一节中,我们介绍了时域的微分相当于频域的滤波,它将每个频率分量乘以 \(2\pi if\) 。 而积分是微分的逆运算,实际上,它相当于把每个频率成分除以 \(2\pi if\) ,也是一个滤波器。
我们可以这样来计算这个滤波器:
integ_filter = close.make_spectrum()
integ_filter.hs = 1 / (PI2 * 1j * integ_filter.fs)
图9.3 右图是这个滤波器在对数Y轴坐标下的图像。
Spectrum
类提供了一个方法来计算积分滤波:
# class Spectrum:
def integrate(self):
self.hs /= PI2 * 1j * self.fs
为了确保这样计算时正确的,我们将它应用到之前的微分的频谱上:
integ_spectrum = deriv_spectrum.copy()
integ_spectrum.integrate()
需要注意的是在 f=0 的时候,我们会进行除0操作,那样会引起 NaN (表示不是一个数 Not a Number)。为了避免这个问题,我们对频率为0的分量 进行特殊的处理,让他的值简单的等于0,然后再生成波形:
integ_spectrum.hs[0] = 0
integ_wave = integ_spectrum.make_wave()
图9.5 中对比了积分运算的结果和原始的数据曲线。显然,积分计算的曲线相当于 原始曲线向下平移了一些。这是由于我们将 f=0 的频率幅值设置成了0,而这代表的是 信号的直流分量。不过,其实这也没有什么好奇怪的,因此微分操作后我们就完全丢失了直流 分量的信息,积分运算并不能对这个损失进行恢复。在一定程度上来说,计算结果中的 NaN 又表示了这个元素是未知的。
当然,如果我们知道这个所谓的“积分常数”,那么积分的结果就是一定的,并且可以保证这个积分滤波器 就是微分滤波器的逆运算。
9.5 累加¶
之前我们说差分操作是微分的近似,那么累加就是积分的近似。我通过一个锯齿信号来演示这个操作:
signal = thinkdsp.SawtoothSignal(freq=50)
in_wave = signal.make_wave(duration=0.1, framerate=44100)
图9.6 展示了结果的波形和频谱图。
Spectrum
类中提供了一个方法来计算波形的累加并返回一个新的波形:
# class Wave:
def cumsum(self):
ys = np.cumsum(self.ys)
ts = self.ts.copy()
return Wave(ys, ts, self.framerate)
我们使用这个方法来计算 in_wave
的累加结果:
out_wave = in_wave.cumsum()
out_wave.unbias()
图9.7 展示了结果的波形和频谱图。如果你认真的完成了第二章后面的练习的话,你就会发现这个波形 与抛物线信号很像。
抛物线信号的频谱和锯齿信号的频谱相比,幅值的衰减要快很多。第二章中,我们知道锯齿信号的幅值是按照 1/f 的规律衰减的,由于累加操作是积分的近似,而积分滤波的效果也相当于按 1/f 的规律衰减。 因此累加操作后,幅值就近似按照 \(1/{f^2}\) 衰减了,这和抛物线信号是一致的。
我们可以这样来计算累加操作所对应的滤波器:
cumsum_filter = diff_filter.copy()
cumsum_filter.hs = 1 / cumsum_filter.hs
由于 cumsum
是 diff
的逆运算,因此我们将 diff_filter
的所有值均设置为它的倒数。
图9.8 对比了累加和积分的频谱响应图。可见累加确实是积分的近似,只是在高频的时候积分滤波器衰减的
稍微要快一些。
为了确保这个滤波器确实和累加操作是一样的效果,我们计算了 in_wave
和 out_wave
的频谱比例并进行了对比:
in_spectrum = in_wave.make_spectrum()
out_spectrum = out_wave.make_spectrum()
ratio_spectrum = out_spectrum.ratio(in_spectrum, thresh=1)
上面的 ratio
方法如下:
def ratio(self, denom, thresh=1):
ratio_spectrum = self.copy()
ratio_spectrum.hs /= denom.hs
ratio_spectrum.hs[denom.amps < thresh] = np.nan
return ratio_spectrum
因为当 denom.amps
很小的时候,相除的结果没太大意义,我们将它直接设置为了 NaN 。
图9.9 中展示了它们的对比结果,可见,它们几乎是一样的。因此,我们可以确信差分滤波器的倒数 就是累加滤波器。
最后,我们来验证卷积定理同样适用于积分滤波器:
out_wave2 = (in_spectrum * cumsum_filter).make_wave()
结果 out_wave2
在浮点数精度误差内与我们通过 cumsum
计算出的 out_wave
是一样的,
所以卷积定理在这里也是适用的。
需要注意的是,这样的验证方法仅仅适用于周期信号。
9.6 噪声的积分¶
在 4.3 中,我们通过对白噪声进行累加来计算出了布朗噪声,现在我们知道了累加操作 在频域上的效果,因此现在我们可以更深入的看一下布朗噪声的频谱。
对于白噪声来说,它的所有频率成分都有同样的平均功率。那么累加操作后,由于每个频率的幅度都被除以了 f , 而功率是幅值的平方,因此对于功率来说,每个频率成分的功率是除以了 \({f^2}\) ,也就是说对于频率为 f 的分量来说,平均功率与 \(1/{f^2}\) 是成正比的,可以表示为:
上式中个的 K 是一个无关紧要的常数。对式子两边求对数后得到:
这也说明了为什么我们将布朗噪声的频谱画到对数坐标下会是一条斜率约为-2的直线。
在 9.1 中,我们画出了Fackbook收盘价的频谱图,并得出它的斜率约为-1.9,这和布朗噪声是基本一致的, 其实很多股票的价格都有类似的频谱。
而当我们使用 diff
计算差分后,我们实际上是在频域上应用了一个与频率 f 成正比的滤波器,也就意味着相应的功率变化
是 \({f^2}\) 。在对数坐标下,这个操作将功率谱的直线的斜率增加了2,因此计算结果的斜率约为0.1。(实际上要低一些,因为
diff
操作只是微分的近似)
9.7 练习¶
下面练习的答案可以参考文件 chap09soln.ipynb
。
练习1 阅读并运行 chap09.ipynb
中的代码。
在 9.5 中我提到对于非周期信号是不适用的。试着使用Facebook数据代替锯齿信号进行同样的计算,
看看会有什么问题。
练习2 这个练习的目的是探索对信号进行差分和微分的效果。生成一个三角波并画出波形,对他进行 diff
操作,
并画出结果的波形。计算三角波的频谱并运用微分滤波器然后画出频谱图,再将这个结果反过来计算出波形并画出波形图。
看看这两个波形有什么区别。
练习3 这个练习的目的是探索对信号进行累加和积分的效果。生成一个方波并画出波形,对他进行 cumsum
操作,
并画出结果的波形。计算方波的频谱并运用积分滤波器然后画出频谱图,再将这个结果反过来计算出波形并画出波形图。
看看这两个波形有什么区别。
练习4 这个练习的目的是探索进行两次积分后的效果。生成一个锯齿波,计算出它的频谱,然后进行两次积分滤波操作。 画出结果的波形和频谱。看看这个波形在数学上是什么形式的曲线?为什么它会像是正弦信号呢?
练习5 这个练习的目的是探索二阶差分和二阶微分的效果。生成一个立方信号 CubicSignal
(在 thinkdsp
中实现了这个类)
通过计算两次差分 diff
运算来进行二阶差分,看看结果是什么样的?通过计算两次微分操作来计算二阶微分,看看结果是否一样?
画出相应的二阶差分和二阶微分的滤波器响应图进行对比。
提示:为了是滤波器在相同的刻度下,需要使波形的采样率等于1。
第十章:线性时不变系统¶
这章中,我们使用音乐作为例子来介绍一些信号与系统的理论知识, 解释卷积定理的证明,以及线性时不变系统的特性(马上我们就会看到线性时不变系统的定义)。
这章的代码 chap41.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp10 查看。
10.1 信号与系统¶
在信号处理领域, 系统(system) 代表的是接受一个信号为输入并产生一个信号输出的一种抽象概念, 如电子放大器,它是一个把输入信号放大后输出的电路。再比如,在你听音乐的时候,整个房间也可以看做是一个系统, 它将声音从它产生的地方(输入),然后传到你的耳朵里(输出)形成与原来不完全相同的声音。
而 线性时不变系统(linear, time-invariant system, LTI) 就是有如下两个性质的一类系统:
- 线性:如果你同时输入两个不同的信号,那么得到的输出和这两个信号各自的输出之和相同。 用数学语言描述为,如果输入 x1 产生输出 y1 ,输入 x2 产生输出 y2,那么 输入 ax1+bx2 产生的输出为 ay1+by2 。这里 a 和 b 为常数。
- 时不变性:系统的作用效果不随时间的变化而变化,仅依赖于系统本身的状态。因此,如果输入变化仅是在时间上平移, 那么产生的输出变化也仅仅是相同的时间平移,其他都是一样的。
很多物理系统都具有这两个性质,或是近似的具有这两个性质。
- 仅包含电阻,电容和电感的电路系统是LTI。这里我们只把这些元件考虑为理想模型。
- 包含弹簧,块和阻尼装置的机械系统同样是LTI。这里假定弹簧和阻尼器是线性的(力与位移成正比,力与速度成正比)
- 本书中的讲的最多的声音的传输介质(包括空气,水和固体等)也可以使用LTI作为近似模型。
LTI可以由线性微信方程描述,而它的解可以用复指数函数来表示, 见 http://en.wikipedia.org/wiki/Linear_differential_equation 。
我们可以运用这个结论来计算输入信号作用与一个LTI后的输出:
- 将输入信号表示为不同频率的复指数函数之和。
- 对每个分量计算相应的输出
- 将所有的输出加起来就得到了最终的输出
这看上去又很熟悉,它与 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 ,它是一个低通滤波器,滤波响应形状近似高斯曲线。
现在,假定我们不知道窗函数或者滤波响应,我们就可以通过输入脉冲信号得到的脉冲响应,来描述这个系统。
这个例子中,我们将脉冲信号的频谱乘以滤波响应来得到了脉冲响应的频谱,然后由频谱生成了波形:
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 左图展示了枪声的波形。
接下来,我们计算 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()
图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
中听听它的声音。
当然,它听上去就是连续的两个枪声,前面的一声要大些,后面一声要小些。
现在我们用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
中播放它。
图10.5 展示了整个计算的过程,其中 f 是锯齿波, g 是脉冲响应, h 是计算结果。
例如:
更一般的来说:
h[n] = sumlimits_{m = 0}^{N - 1} {f[m]g[n - m]}
在 8.2 我们就见过这个式子了,它就是 f 与 g 的卷积。 说明如果系统的脉冲响应为 g ,那么系统输入 f 产生的输出应该为 f 与 g 的卷积。
概况起来,我们有两种形式来描述系统对输入信号的影响效果:
- 把输入理解为一系列的脉冲信号,那么输出就是脉冲响应进行时移和缩放后的和, 也就是输入与脉冲响应的卷积。
- 脉冲响应的DFT是系统的传递函数,它包含了系统对每个频率成分的影响效果,而输入可以理解为 不同频率的分量之和,因此将输入的DFT乘以传递函数就可以得到输出的DFT。
这两种描述无疑是等价的,这是由卷积定理所得到的:时域的卷积等价于频域的乘积。
这里我们也可以懂得为什么卷积的形式中 g 是反向的,这在我们学习平滑的时候提到过。原因就是 卷积的定义是在研究LTI的响应的时候自然得到的。
10.5 卷积定理的证明¶
现在我们是时候来解释卷积定理的证明过程了:
式中 f 和 g 是长度均为 N 的两个向量。
证明过程分为两步:
- 我会先从 f 为复指数信号这个特例开始,说明它与 g 的卷积相当于对 f 乘以一个标量
- 然后我们把 f 推广到一般的信号,它可以表示为不同频率的复指数信号之和,然后通过乘以一个标量 来计算各个频率成分的卷积,再把结果加起来
为了证明卷积定理,我们先来从一些基本的式子开始。首先, g 的DFT,写作G,等于:
上式中, k 表示 0~N-1 的频率成分, n 表示 0~N-1 的采样时间。 这个结果是包含 N 个复数的向量。
而 f 的IDFT 写作 F , f 等于:
卷积的定义为:
式中 m 的取值范围也为 0~N-1 ,由于卷积是满足交换律的,因此上式等于:
现在我们考虑 f 为复指数的特殊情况,我们将频率为 k 的复指数写作 \({e_k}\) ,那么:
这里 k 表示频率, n 表示时间。 将 \({e_k}\) 代入之前卷积的第二个定义中得到:
把上式分解后得到:
可以看出,这个式子的第一项为 \({e_k}\) 本身, 第二项为 G[k] ,因此又可以写成:
这个式子表明,复指数 \({e_k}\) 与 g 的卷积等于 \({e_k}\) 乘以 G[k] 。 用数学语言来描述就是,与 g 进行卷积这个算子的特征函数是 \({e_k}\) , 对应的特征值为 G[k] (见 9.3 )
现在我们进行证明的第二步。当输入信号 f 不是复指数信号的时候,可以通过DFT将它表示为多个 复指数之和的形式,写作 F[k] 。其中 k 为 0~N-1 的频率值, F[k] 就是包含了不同的频率 的复指数的幅值和初始相位。
根据之前的证明,对于输入信号的每个频率成分 F[k] ,输出信号的每个频率成分的系数应该为 F[k]G[k] 。 由于系统是线性的,输出就可以表示为:
代入 \({e_k}\) 的定义得到:
上式的右面部分实际上就是 F 和 G 的乘积 FG 的IDFT,因此:
由于 F=DFT(f) 以及 G=DFT(g) ,因此:
最后就得到了卷积定理的公式:
证明完毕。
10.6 练习¶
下面练习的答案可以参考文件 chap10soln.ipynb
。
练习1 在 10.4 中我把卷积描述为信号经过时移和缩放后的和。 但是,在 10.3 中我们把信号的DFT乘以传递函数,这个操作实际对应的 应该是 循环卷积 ,因为DFT假定信号是周期的。其结果会导致输出波形的前面有一些额外的影响。
幸运的是,这个问题有一种标准的解决方案,在信号的尾部加入足够多的0,然后再进行DFT,就可以避免 这个问题。
将 chap10.ipynb
中的例子进行补0后再计算,以消除这个问题对结果的影响。
练习2 OpenAIR ( http://www.openairlib.net )中提供了很多在线的声学脉冲响应数据。 浏览并下载一个你感兴趣的脉冲响应数据。找一段与你下载的脉冲响应同采样率的录音, 模拟这个录音通过这个脉冲响应的系统后产生的输出。使用两种方式来计算:一是通过输入与脉冲响应的卷积; 二是通过输入的DFT与脉冲响应的滤波响应相乘。
第十一章:调制和采样¶
在 2.3 中我们知道,当信号的采样率为10000Hz时, 5500Hz的频率与4500Hz的频率是区分不出来的,这时的折叠频率为5000Hz, 也就是采样率的一般。当时,我并没有解释原因。
这章我们来研究采样以及采样定理,解释混叠以及折叠频率的原因。
我会先从脉冲的卷积开始来解释调幅(AM)的过程,这对理解采样定理是很有用的。
这章的代码 chap11.ipynb
可以在本书的 代码库 中找到,你也可以在 http://tinyurl.com/thinkdsp11 查看。
1.1 脉冲卷积¶
在 10.4 中我们得知信号与一系列脉冲的卷积的结果就是信号进行时移和缩放后的和。 这里,我以一个蜂鸣声作为例子:
filename = '253887__themusicalnomad__positive-beeps.wav'
wave = thinkdsp.read_wave(filename)
wave.normalize()
然后构造包含四个脉冲的波形:
imp_sig = thinkdsp.Impulses([0.005, 0.3, 0.6, 0.9],
amps=[1, 0.5, 0.25, 0.1])
impulses = imp_sig.make_wave(start=0, duration=1.0,
framerate=wave.framerate)
然后对他们进行卷积:
convolved = wave.convolve(impulses)
结果如 图11.1 ,输入信号为左上图,脉冲序列为左下图,卷积结果为右图。
你可以在 chap11.ipynb
中播放这个声音,它听起来就是四个逐渐衰减的蜂鸣声。
这个例子演示了信号与脉冲序列卷积的结果,这对于之后的分析是很有用的。
11.2 幅度调制¶
幅度调制(调幅,AM)被广泛用于调幅广播和其他一些应用中。信号(包括语音,音乐等) 以与载波相乘的方式被调制到载波上,经过发射装置发送出去。一般载波都是适用于无线电广播的 高频信号,在美国,AM广播的频率为500~1600kHz,参见 https://en.wikipedia.org/wiki/AM_broadcasting 。
而在接收端,需要从广播信号中恢复出原始的信号,这个过程就是解调。令人惊讶的是,解调也是通过与同样的载波 相乘来计算出来的。为了说明工作原理,我将一个信号调制到了10kHz的载波上,信号如下:
filename = '105977__wcfl10__favorite-station.wav'
wave = thinkdsp.read_wave(filename)
wave.unbias()
wave.normalize()
载波如下:
carrier_sig = thinkdsp.CosSignal(freq=10000)
carrier_wave = carrier_sig.make_wave(duration=wave.duration,
framerate=wave.framerate)
我们可以直接使用乘法来进行调制的操作:
modulated = wave * carrier_wave
这个信号的声音很难听,你可以在 chap11.ipynb
中播放它。
如 图11.2 展示了调幅在频域的变化过程。顶部是原始信号的频谱,接下来是调制后的信号频谱,也就是与载波相乘之后的信号, 它包含了原始频谱的两个拷贝,分别向左右平移了10kHz。
我们知道时域的卷积等于频域的乘积,反过来其实也成立,时域的乘积等于频域的卷积。当我们将信号乘以载波的时候,在频域上, 就相当于把信号的频谱与载波的频谱进行了卷积。
由于载波是简单的余弦信号,它的DFT仅包含两个脉冲,10kHz以及-10kHz。由上一节可知,与脉冲序列的卷积等于信号平移和缩放后的和。 注意经过调制之后,频谱的幅值减小了,这个因为原始实现的能量变分散到了两个部分中。
我们再通过乘以载波来进行解调:
demodulated = modulated * carrier_wave
图11.2 的第3个图展示了结果的频谱,这个结果就是图2经过平移缩放后得到的。由于调制后的频谱具有两边对称的峰值, 每个峰值都被分成了左右平移10kHz,幅值为一半两个部分,两边的0kHz部分正好重合到了一起。而另外两个部分分别是 20kHz以及-20kHz。
解调后的信号听上去恢复的很好,虽然其中包含了一些原始信号没有的高频成分,但是这些频率太高了,大部分播放设备播放不出来, 而且我们的耳朵也听不到。当然,如果你的播放设备足够好,而你又有一个足够好的耳朵,你也许可以听到。
这种情况,我们其实可以使用一个低通滤波器将额外的高频成分滤掉:
demodulated_spectrum = demodulated.make_spectrum(full=True)
demodulated_spectrum.low_pass(10000)
filtered = demodulated_spectrum.make_wave()
这样结果就和原始信号更接近了,只是在解调和滤波后损失了一半的功率。这在实际应用中并不是一个问题, 因为实际上大部分的功率都损失到了广播信号的发射,传播和接收上,不过我们只要把解调后的信号再进行 放大就可以了。
11.3 采样¶
上一节我们介绍了幅度调制,它可以帮助我们理解采样。采样(sampling)就是在一系列的时间点上 (通常是等间隔的)对一个模拟信号进行测量得到一系列的数字测量值的过程。
例如我们之前使用WAV文件,就是使用模数转换器(ADC)对麦克风的输出进行采样后保存成的文件。 大多数音频的采样率为44.1kHz,这个标准的CD音质的采样率,而使用48kHz的采样率为标准DVD音质的采样率。
在48kHz采样率之下,折叠频率是24kHz,这比我们大部分能听到到最高的频率还要高。 参见 https://en.wikipedia.org/wiki/Hearing_range 。
大多数的声音波形,每个采样值为16位(bits),也就是有 \({2^{16}}\) 个等级,这也被称为位宽。 \({2^{16}}\) 的位宽对于声音来说就已经足够了,即使再增大,我们也感觉不到音质的提升了。 参见 https://en.wikipedia.org/wiki/Digital_audio 。
当然,对于音频信号以外的其他信号来说,有的需要更高的采样率以捕捉到更高的频率,有的需要更大的位宽以 提高波形的保真度。
为了演示采样的过程,我还是先从一个44.1kHz采样率的波形开始。然后中波形中用11kHz的采样率再采样一次。 这与直接从模拟信号上采样是不完全一样的,但效果其实差不多。
首先,我加载了一段单独的鼓声:
filename = '263868__kevcio__amen-break-a-160-bpm.wav'
wave = thinkdsp.read_wave(filename)
wave.normalize()
图11.3 上图显示了这个波形的频谱,下面是对这个波形进行采样的函数:
def sample(wave, factor=4):
ys = np.zeros(len(wave))
ys[::factor] = wave.ys[::factor]
return thinkdsp.Wave(ys, framerate=wave.framerate)
然后使用这个函数从原始信号每4个值中抽取一个值作为采样后的波形:
sampled = sample(wave, 4)
最后的结果与原始波形具有相同的采样率,但是其中大部分元素都是0。如果你播放这段采样后的波形, 它听起来不会很好听,它相比原始的信号引入了一些高频的成分。
图11.3 下图显示了采样后的频谱图。它是原始信号的频谱复制了四份(图中看上去是5份,但是其实边上的 两个其实一个被分成了两份)
与上一节类似,采样过程可以看做是信号与一系列的脉冲的乘积,因此产生了相同的效果。除了使用 sample
以外,我们还可以使用下面这个函数来进行采样,这种方法有时被叫做 脉冲序列(impulse train)
def make_impulses(wave, factor):
ys = np.zeros(len(wave))
ys[::factor] = 1
ts = np.arange(len(wave)) / wave.framerate
return thinkdsp.Wave(ys, ts, wave.framerate)
然后将信号与脉冲序列相乘得到采样信号:
impulses = make_impulses(wave, 4)
sampled = wave * impulses
结果与之前是相同的。由于时域的乘积等于频域的卷积,当我们乘以脉冲序列的时候,就相当于在频域上与脉冲序列的频谱 进行卷积。事实证明脉冲序列的DFT同样也是脉冲序列。
图11.4 中,上面两个图是11025Hz的脉冲序列的波形和频谱。可以看出它的DFT是包含4个脉冲的脉冲序列,这就是为什么 上例会得到原始频谱的四个复制。下面两个图是5512Hz的脉冲序列的波形和频谱。它的DFT是8个脉冲。一般来说,时域上 有更多从脉冲,在频谱上就会有更少的脉冲。
总结起来就是:
- 采样过程可以看做是信号与采样序列的乘积。
- 时域上乘以采样序列,相当于在频域上与一个采样序列的卷积
- 与采样序列的卷积后的频谱是原始频谱的周期延拓
11.4 混叠¶
在 11.2 中,我们对解调后的AM信号使用了低通滤波器来去掉多于的频谱部分。 而对于采样来说这样做并不是一个完美的方案。
图11.5 说明了为什么这样做会有什么问题。其中第一行的图是原始信号的频谱,它包含了超过10kHz的高频成分, 当我们对他进行采样,相当于频域上与一个脉冲序列进行卷积,第二行的图就是这个脉冲序列的频谱。卷积后得到了 采样信号的频谱(第三行),第四行的图是采样信号经过5512Hz截止频谱的低通滤波后的频谱。
虽然,将结果反算出波形后,与原始信号是相似的,但是其实有两个问题:
- 经过低通滤波后,原始信号中高于截止频谱5512Hz的频谱成分丢失了
- 并且低于5512Hz的频率成分中叠加其他左边原始频谱的右半部分。
也就是说,如果相邻的两个频谱复制重叠了,那么在采样后,我们就不能从中把原始信号恢复出来了。
但是,如果没有重叠的话,这个方法还是很好的,下面的例子我使用了一段贝斯的声音。
如 图11.6 所示,它的频谱(上图)中几乎看不到大于5512Hz的频谱成分,中图则是采样后的频谱, 下图是滤波后的频谱。由于滤波会导致能量的损失,因此滤波后的频谱的幅值变小了,但是频谱的形状与 原始信号是完全一样的,滤波后的波形听上去和原始声音也是完全一样的。
11.5 内插¶
上一节中我们使用的滤波器又叫做 矩形滤波器(brick wall filter) ,它对于高于截止频率的成分 是完全消除的。
图11.7 右图显示了这个滤波器的图形。由于频域上相乘等于时域上卷积,也就是与滤波器的卷积窗进行卷积。 它的卷积窗就是频谱的IDFT(右图)。这个函数有一个特殊的名字叫 辛格函数(sinc function) (实际上是 它的离散形式),参见 https://en.wikipedia.org/wiki/Sinc_function :
当应用这个滤波器的时候,就相当于与 sinc 函数进行卷积。
sinc 函数在0时刻的值是1,而当 x 为其他整数的时候,它的值均为0 。因此当与这个函数进行卷积的时候, 在采样点上得到的值严格的等于采样值,因此卷积的过程实际上就是对采样点进行内插。
图11.8 中使用一小段信号来展示了内插的过程。其中蓝线为原始信号,绿线是卷积中使用的一系列sinc函数经过时移和缩放的 曲线,它们的和等于原始信号。灰线是采样信号。 我再重新强调一次:
- 这些sinc函数的和等于原始信号
这是因为原始信号没有高于5512Hz的频率成分,在11025Hz的采样率下,我们可以完成的重建原始信号。当然,如果我们有原始信号的 频谱,也可以准确的重建原始信号。
上面的例子,我们先从一个44100Hz采样率的信号开始,然后以11025Hz的采样率进行重新采样,采样后的频谱是原始频谱的周期延拓, 周期是11025Hz。
如果原始信号没有高于5512Hz的频率成分,那么周期延拓后的频谱就不会重叠,在滤波后就不会丢失信息,也就可以完全的重建信号。
这就是所谓的香农采样定理,见 https://en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem 。
之前我们的讨论都是具有原始信号是在44.1kHz采样率下进行的。实际上,原始信号采样率无关紧要,它的采样率可以更高,甚至是连续的模拟信号, 这个定理都是成立的:当以采样率 f 对信号进行采样后,只要信号不包含高于 f/2 的频率成分,我们就可以准确的从采样信号中重建原始信号。 这种信号又叫做 有限带宽信号(bandwidth limited) 。
译者注
总结:时域中的连续信号经脉冲序列采样后会在频域产生周期延拓,周期等于采样频率; 只有采样频率足够高,大于原始信号频谱中最高频率的两倍,则采样后频谱就不会发生混叠。 为了恢复信号,可以把采样信号通过一个矩形滤波器,而矩形滤波器的窗函数就是sinc函数。 因此将采样后的信号与sinc信号进行卷积就可以重建信号,这个过程就是内插。
11.6 小结¶
恭喜你已经完成了本书的所有内容(除了本章的练习以外),在结束之前,我重新总结一下所有的内容:
- 我们从周期信号开始,介绍了
thinkdsp
模块中的几个关键类:Signal
,Wave
,Spectrum
。 - 学习了简单波形以及乐器的谐波结构,并且看到了混叠的效果。
- 学习了声谱图用来分析啁啾声以及其他的频谱随时间变化的信号。
- 生成并分析了噪声信号,对自然界的噪声信号特性有了一定的了解。
- 学习了自相关函数,并用它来进行音高估计,以及描述了噪声的另外一些特性。
- 学习了离散余弦变换,它在数据压缩中应该很广,然后又进一步学习了离散傅立叶变换。
- 使用复指数信号来合成复信号,并通过它的逆向过程得到了DFT。如果你做了那张的练习,就会实现快速傅立叶变换了,这个算法是20世纪最重要的算法之一。
- 从平滑开始,我们学习了卷积以及卷积定理,它描述了时域的操作(如平滑)在频域上的滤波作用。
- 探索了微分和积分操作的滤波作用,他们是求解微分方程的频谱方法的基础。我们还解释了前面章节中的一些遗留问题,例如白噪声和布朗噪声的关系。
- 学习了线性时不变系统,并使用卷积定理将LTI系统描述为它的脉冲响应。
- 演示了幅度调制的过程,这在广播通信中很重要。并且进一步学习了采样定理,这个定理是数字信号处理的关键。
至此,你应该具有了使用计算工具对信号和频谱进行分析的实践知识以及如何进行采样和滤波等工作的理论知识了。
希望你能够从中找到感兴趣的东西。谢谢!
11.7 练习¶
下面练习的答案可以参考文件 chap11soln.ipynb
。
练习1 阅读并且运行 chap11.ipynb
上面的代码示例,听一听示例中的音频。
练习2 Chris “Monty” Montgomery 有一个很好的视频, 叫做“D/A and A/D | Digital Show and Tell” 它演示了采样定理并且展示了很多关于采样的很好的知识, 观看这个视频 https://www.youtube.com/watch?v=cIQ9IXSUzuM 。
练习3 如果我们使用一个比较低的采样率进行采样,高于折叠频率的部分就会混叠,一旦发生混叠, 我们就不可能完全重建信号了,因为这些频率成分采样和低频的成分完全一样区别不出来了。 我们可以在采样之前对信号进行低通滤波,滤除这部分频率成分,这种滤波就叫 抗混叠滤波器(anti-aliasing filter) 。
回到之前鼓声的例子,在采样前对信号进行低通滤波,采样后在使用低通滤波去掉频谱周期延拓的部分。 看看最后恢复的信号和滤波后的信号是否是一样的。