详情
本书介绍了使用Python语言实现数字信号处理的方法，内容共有11章，以Python代码为示例由浅入深地向读者介绍了数字信号处理的相关知识及其应用。书中涉及周期信号及其频谱、波形的谐波结构、非周期信号以及频谱图、噪声、自相关函数、离散余弦变换和离散傅里叶变换、滤波、卷积、微分与积分、调制采样等数字信号处理相关技术。
书名：Python数字信号处理应用
ISBN：978-7-115-46952-6
本书由人民邮电出版社发行数字版。版权所有，侵权必究。
您购买的人民邮电出版社电子书仅供您个人使用，未经授权，不得以任何方式复制和传播本书内容。
我们愿意相信读者具有这样的良知和觉悟，与我们共同保护知识产权。
如果购买者有侵权行为，我们可能对该用户实施包括但不限于关闭该帐号等维权措施，并可能追究法律责任。
• 著 [美] 艾伦•唐尼（Allen B. Downey）
译 缪 文
责任编辑 陈冀康
执行编辑 武晓燕
• 人民邮电出版社出版发行 北京市丰台区成寿寺路11号
邮编 100164 电子邮件 315@ptpress.com.cn
• 读者服务热线：(010)81055410
反盗版热线：(010)81055315
Copyright ©2016 by O’Reilly Media, Inc.
Simplified Chinese Edition, jointly published by O’Reilly Media, Inc. and Posts & Telecom Press, 2017. Authorized translation of the English edition, 2017 O’Reilly Media, Inc., the owner of all rights to publish and sell the same.
All rights reserved including the rights of reproduction in whole or in part in any form.
本书中文简体版由O’Reilly Media, Inc.授权人民邮电出版社出版。未经出版者书面许可，对本书的任何部分不得以任何方式复制或抄袭。
版权所有，侵权必究。
数字信号处理（DSP）是面向电子信息学科的专业基础知识，也是多门新兴学科的理论基础。作为广泛应用的脚本语言，Python在DSP领域中也很常用。
本书介绍了如何通过Python语言实现数字信号处理的应用。全书共有 11 章，以Python代码为示例由浅入深地向读者介绍了数字信号处理的相关知识及其应用。书中涉及周期信号及其频谱、波形的谐波结构、非周期信号及频谱图、噪声、自相关函数、离散余弦变换和离散傅里叶变换、滤波、卷积、微分与积分、调制采样等数字信号处理相关技术。每一章都是从示例开始，引导读者通过编程的方式来准确地理解概念。除此之外，本书每章还提供了练习和代码示例来帮助读者理解这些知识。
本书适合对数字信号处理感兴趣且有一定Python基础的读者阅读，也适合电子和通信相关专业的学生阅读。
信号处理是我最喜欢的主题之一。在科学和工程技术的很多领域，信号处理都发挥了很大的作用。如果你对信号处理的基础概念有所了解的话，就应该知道，它为我们洞察世界（特别是为我们所听到的东西）提供了工具。
但是，大部分人很可能没有学习信号处理的机会，除非恰好学习过电子信息工程或者机械工程。这方面的大部分书籍的问题（以及使用这些书籍的课程的问题）在 于，它们以向量这样的数学抽象概念为开始，整本书堆砌了各种公式、充满了理论，但在应用和实践方面非常欠缺。
本书的初衷在于，如果读者知道如何编程，那就能利用这些编程技巧学习其他东西，并从中获得快乐。
通过这种基于编程的方式，本书能准确传达最重要的概念。仅在学完第1章之后，读者就能分析录音数据以及其他各种信号，并生成新的音频数据。在每章中，本书都会介绍一种新的技术，并提供相关的真实信号的实践教程。在每一步中，读者都会先了解如何应用它，然后再理解其原理。
这样的方式更有实践意义，也希望读者能体味到其中的乐趣。
本书中的例子和辅助代码都是基于Python语言的。读者必须拥有Python的基础，并对面向对象编程比较熟悉——至少对使用并非自己定义的对象比较熟悉。
如果读者还不太了解Python，不妨从我的另一本书《像计算机科学家一样思考Python》(《Think Python》)开始，这本Python入门书是为从未接触过编程的人写的，你也可以从Mark Lutz所著的《Learning Python》开始，这本书更适合于有一些编程经验的人。
我使用NumPy和SciPy的频度很高。如果你对它们已经很熟悉了，那再好不过。不熟悉的话也没关系，我会解释用到的函数和数据结构。
我希望读者对一些基础数学概念有所了解。需要读者做的计算并不多，知道积分和微分的概念就够了。本书还使用了线性代数，在遇到的时候都会给出解释。
本书中的代码音频数据都可以从GitHub仓库（repository）获取：https://github.com/ AllenDowney/ThinkDSP。Git是一个版本控制工具，能让读者追踪项目文档的进 度——如果读者不熟悉Git或GitHub的话。在Git控制中写的文档集合被称为“仓库”。GitHub为Git仓库提供了存储功能，并提供了一个方便使用的网站界面。
有多种方式可以使用GitHub仓库的主页的代码。
$ jupyter notebook
所有的代码无需更改就可以在Python 2和Python 3环境下运行。
$ conda install jupyter
本书使用的是Continuum Analytics的Anaconda，这是一个免费的Python发布版，包含了运行这些代码所需要的全部包（并且还有更多其他的功能）。我认为，Anaconda安装很容易。在默认情况下Anaconda是以用户的形式安装的，而且不会影响到系统层面，所以读者在安装的过程中不需要提供管理员权限。而且它还同时支持Python 2和Python 3。 读者可以从http://continuum.io/downloads下载Anaconda。
如果读者不想使用Anaconda, 那么就需要下面的包。
虽然这些包很常用，但是它们并没有被Python安装版收入，而且在某些情况下安装它们可能非常困难。如果你在安装过程中碰到这种情况，我推荐安装Anaconda或者其他的包含这些包的Python发布版。
大部分的练习都基于Python脚本，不过也有部分使用的是 Jupyter Notebook。如果读者在这之前没有使用过Jupyter，可以从http://jupyter.org/获取其信息。
使用Jupyter Notebook的方式有3种。
祝大家好运，并享受学习过程！
在线Safari图书（http://safaribooksonline.com/）是一个按需服务的数字图书馆，以图书和视频的形式提供来自世界顶级的技术作者和商业作者的专业内容。
技术专家、软件开发者、网页设计者和商业创新专家将在线Safari图书作为他们研究、解决问题、学习和培训认证的首选资源。
在线Safari图书为企业、政府机构、学校和个人提供丰富的产品组合和定价模式。订阅者可以通过全文搜索数据库访问出版商提供的成千上万的图书、培训视频及正式出版前的手稿。这些出版商包括但不限于：O’Reilly Media、Prentice Hall Professional、Addison-Wesley Professional、Microsoft Press、Sams、Que、Peachpit Press、Focal Press、Cisco Press、JohnWiley & Sons、Syngress、Morgan Kaufmann、IBM Redbooks、Packt、Adobe Press、FTPress、Apress、Manning、New Riders、McGraw-Hill、Jones & Bartlett和Course Technology。更多信息请访问我们的网站。
请将对本书的评论和问题发送给出版商。
美国：
O’Reilly Media, Inc.
1005 Gravenstein Highway North
Sebastopol, CA 95472
中国：
北京市西城区西直门南大街2号成铭大夏C座807室（100035）
奥莱利技术咨询（北京）有限公司
我们为本书制作了一个网页，列出了勘误表、代码示例及各种附加信息。读者可以通过这个网址访问：http://bit.ly/think-dsp。
评论或技术问题可发电子邮件到bookquestions@oreilly.com。
更多信息请访问http://www.oreilly.com网站。
我们的Facebook: http://facebook.com/oreilly。
我们的Twitter: http://twitter.com/oreillymedia。
我们的YouTube: http://www.youtube.com/oreillymedia。
读者如果有意见或校正信息，请发送邮件到downey@allendowney.com。如基于您的反馈进行校正，您将被列于贡献者列表中（除非您要求匿名）。
如果您在提交时可以列出来错误所在的句子，这将对我们非常有帮助。列出页码和章节题目也不错，但这样处理起来要麻烦一些。谢谢！
特别感谢技术审校Eric Peters、Bruce Levens和John Vincent。他们给出了很多建议、说明和修正。
同时感谢Freesound, 我在本书中使用的很多声音样本都来自于此。本书的GitHub仓库中也包含了一些他们的波形文件，文件名保持原样，因此找到他们的来源应该也很容易。
遗憾的是，大部分的Freesound用户都没有提供他们的真名，所以我只能在致谢的时候使用他们的用户名。本书中使用的样本由Freesound用户iluppai、wcfl10、 thirsk、docquesting、kleeb、landup、zippi1、themusicalnomad、bcjordan、rockwehrmann、marcgascon7和jcveliz贡献。感谢你们！
Allen B.Downey是欧林工程学院（Olin College of Engineering）计算机科学学院的教授，曾任教于卫斯理学院（WellesLey College）、科尔比学院（Colby College）以及加州大学伯克利分校（U. C. Berkeley）。他在加州大学伯克利分校获得了计算机科学的博士学位，在麻省理工大学（MIT）获得了硕士和学士学位。
本书的封面动物是滑嘴犀鹃（smooth-billed ani, Crotophaga ani），它是属于杜鹃科的一种大型鸟类，主要栖于佛罗里达、巴哈马、加勒比群岛以及中美洲和南美洲的部分地区。
滑嘴犀鹃全身羽毛为黑色，尾长且喙为齿状。它们在地面觅食，可食用的东西包括白蚁、昆虫，甚至还有小蜥蜴和青蛙。这种鸟喜欢宽阔的、混有平地和矮灌木的栖息地。由于人类的定居以及对森林的砍伐，滑嘴犀鹃的居住环境受到了较大影响，它们也适应性地开始出入于农牧场，以畜牧业产生的昆虫为食。
滑嘴犀鹃非常特别，为群居生物。多对配偶聚在一起，轮流在树的高处搭建椭圆型鸟巢，以此安放鸟蛋，并抚育幼鸟。每个雌性会产下4~7个蛋，鸟巢中最多时会有29个鸟蛋。
O’Reilly封面上出现的很多动物都面临濒危，它们对世界非常重要。如果想要贡献自己的一点力，请访问animals.oreilly.com。
信号代表随着时间变化的量。这个定义是非常抽象的，所以我们从具象的例 子——声音开始。声音源自于空气压力的改变。声音信号代表的是空气压力随着时间的变化。
传声器是测量上述变化并产生表示所测声音的电信号的设备。扬声器是通过输入的电信号产生声音的设备。传声器和扬声器都被称为换能器（transducer），因为它们将信号从一种形式转化成另一种形式，也就是变换了能量的形式。
这本书重点关注信号处理，包括信号的合成、转换和分析。本书将着重于声音信号，但是其方法也同样适用于电信号、机械振动和其他各个领域的信号。
这套方法也同样适用于随空间而不是时间变化的信号，比如远足路线上的海拔变化。还适用于不止一个维度的信号，比如图像——读者可以将它想象为在二维空间中变化的信号。又或者电影，它是二维空间中随着时间变化的信号。
本书从简单的一维声音信号开始。
本章中的代码在这本书资料库的chap01.ipynb中。同时你也可以通过访问http://tinyurl.com/thinkdsp01获取。
先从周期信号开始，周期信号是在一段时间之后重复出现的信号。比如，你在敲钟时候，钟会振动从而产生声音。把这段声音录制下来，然后绘制出转换后的信号，具体如图1-1所示。
这个信号和三角函数相似，也就是说其形状和正弦三角函数的形状一样。
图1-1中的信号是周期性的。本书所选择的时段显示了3个完整的重复，这也被称为循环。每个循环的时长被称为周期，图1.1中信号的周期大约2.3ms。
图1-1 钟声录音片段
信号的频率是每秒钟内周期的数量，即周期的倒数。频率的单位是每秒钟循环数，或者称为赫兹（Hertz），英文符号为“Hz”。
图1-1所示的信号的频率大约为439 Hz，比交响音乐对音标准的440 Hz稍低。这个音符的音名叫作A，或者更确切地说是A4。如果你对“科学音调记号法”不熟悉，这里的数字后缀表示的是音符所在的八度。A4意思是中央C上面的那个A，而A5则比其高出一个八度。详见http://en.wikipedia.org/wiki/Sientific_pitch_notation。
音叉产生的是正弦信号，因为线形的变化是一种简单的谐波运动形式。大部分乐器产生的是周期信号，但是这些信号的形状不是正弦式的。例如，图1-2所示的是小提琴演奏鲍凯利尼E大调第五弦乐五重奏第三乐章的录音片段。
图1-2所示的信号是周期性的，但是信号的形状更加复杂。周期信号的形状被称为波形。大部分乐器产生的波形比正弦信号的复杂。波形的形状决定了音乐的音色，也就是我们对声音品性的感受。对于比正弦信号更复杂的复杂波形，人们主观上认为它更加富有内涵、更温暖而有趣。
图1-2 小提琴录音片段
本书中最重要的主题是频谱分析，在这个意义下，任何信号都可以表示成一系列不同频率的正弦信号的叠加和。
本书中最重要的数学概念是离散傅里叶变换（Discrete Fourier Transform, DFT）, DFT也就是将信号转换成频谱。频谱是指相加产生信号的正弦波的集合。
本书中最重要的算法是快速傅里叶变换（Fast Fourier Transform, FFT）, 它是计算离散傅里叶变换的一种高效方式。
例如，图1-3显示的是图1-2中小提琴录音的频谱。x轴表示的是合成这个信号的频率范围，y轴表示各个频率元素的强度，或者说是振幅。
其中频率最低的元素被称为基频。图1-2所示信号的基频在440 Hz左右（实际上要低一点，或者说“降调”）。
图1-3 小提琴录音片段的频谱
此信号中，基频具有最高的幅度，所以它也是主频。通常情况下，感知到的声音的音高是由其基频决定的，即使它不是主频。
频谱中的其他峰值的频率还有880、1 320、1 760和2 200，它们是基频的整数倍。这些频率元素被称为谐波，因为它们在乐理概念上跟基频和谐：
1 320 Hz接近E6，也就是A5的纯五度。如果读者对“纯五度”这样的音程概念不熟悉，可以参考https://en.wikipedia.org/wiki/Interval_（music）；
1 760 是A6的频率，比基频高出两个八度；
这些谐波共同构成了A大调和弦，虽然它们并不在同一个八度内。其中的一些音只是近似，因为西方音乐已经按平均律作了调整（参见http://en.wikipedia.org/wiki/ Equal_temperament）。
在给定谐波和其各自的振幅之后，读者就可以通过加和正弦波来重建信号了。这我们将在后文中讲解。
我写的Python模块thinkdsp.py中包含信号和频谱[1]分析的类和函数。读者可以在本书的仓库中找到这个模块（参见前言中“代码示例的使用”）。
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是相位差，按弧度定义。相位差定义了信号周期的开始。比如，offset=0的正弦信号从sin 0开始，也就是从0开始。offset=π/2则从sinπ/2开始，也就是从1开始。
各信号都有__add__方法。所以读者可以对它们使用+操作：
mix = sin_sig + cos_sig其结果为SumSignal，表示了两个或者多个信号的叠加和。
Signal实际上是数学函数的Python表示。大部分的信号（Signal）在全部的t值下做了定义——从负无穷到正无穷。
读者如果不对Signal进行求值，那它就没什么作用。在这个语境下，“求值”的意思是，根据一系列的时间点ts，来计算出对应的信号值ys。我在Wave对象中使用封装了表示ts和ys的NumPy数组。
Wave表示的是信号在一系列时间点下求出的值。每个时间点被称为帧（这是借用了电影和视频的概念）。测度本身被称为样本，虽然“帧”和“样本”在一些情况下可以互换使用。
Signal提供了make_wave, 后者返回一个新的Wave对象：
wave = mix.make_wave(duration=0.5, start=0, framerate=11025)duration参数指的是Wave的长度，单位为秒。start意为开始的时间，单位也是秒。framerate是每秒帧数，它是一个整数，意思也就是每秒内的样本数。
11 025帧每秒是在音频文件格式中最常使用的几个帧率之一，这些音频文件包括波形音频文件（Waveform Audio File, WAV）和MP3。
这个例子对信号从t=0到t=0.5进行求值，共取得5 513个间距相等的帧（5 513是11 025的一半）。帧间的时间差被称为时间步长，值为1/11 025 s，约等于91μs。
Wave提供了plot方法，使用方法为pyplot。读者可以这样绘制波形：
wave.plot()
pyplot.show()pyplot是matplotlib的一部分。很多Python发布版都含有它，否则读者就需要自行安装了。
对于freq=440，在0.5秒之内有220个周期，所以其绘制结果看起来会像一大块色块。要放大显示几个周期的波形，读者可以使用segment，它复制了Wave的一段，然后返回一个新波形：
period = mix.period
segment = wave.segment(start=0, duration=period*3)period是Signal的属性，它返回的是信号的周期，单位为秒。
start和duration的单位都是秒。这个例子从mix中复制了前3个周期。其结果是一个Wave对象。
如果我们绘制segment，它的结果如图1-4所示。这个信号包含两个频率元素，所以它比音叉的的信号更为复杂，但是比小提琴的简单。
图1-4 两个三角函数混合信号的片段
thinkdsp提供的read_wave可以读取WAV文件并返回一个Wave对象：
violin_wave = thinkdsp.read_wave('input.wav')Wave提供了write函数，它能写WAV文件：
wave.write(filename='output.wav')读者能通过媒体播放器播放WAV文件，从而听到这个Wave对象。在Unix系统上我使用的是aplay，它是一个简单鲁棒的播放器，很多Linux发行版都有它。
thinkdsp同时还提供play_wave, 它能以子进程运行媒体播放器：
thinkdsp.play_wave(filename='output.wav', player='aplay')一般情况下，我使用的是aplay，读者当然也可以使用其他播放器。
Wave提供了make_spectrum，它返回的是Spectrum:
spectrum = wave.make_spectrum()而Spectrum提供了plot:
spectrum.plot()
thinkplot.show()我所写的thinkplot模块提供了pyplot中一些函数的封装。它包含于本书的Git仓库中（参见前言的“代码示例的使用”）。
Spectrum提供了3种修改频谱的方法。
下面的例子将所有高于600的频率衰减了99%：
spectrum.low_pass(cutoff=600, factor=0.01)低通滤波器会移除明亮的高频声音信号，所以其结果的声音比较压抑而且昏暗。读者可以将Spectrum转化回Wave，然后播放来感受一下：
wave = spectrum.make_wave()
wave.play('temp.wav')play方法将波形写入一个文件并播放它。如果读者使用Jupyter Notebook，就可以使用make_audio，它会创建一个音频部件用于播放声音。
thinkdsp.py里面并没有什么复杂的东西。其提供的大部分函数只是对NumPy和SciPy函数的稀薄封装（thin wrappers）。
thinkdsp的初始库包括Signal、Wave和Spectrum。给定Signal，读者就能创建一个Wave。给定一个Wave，读者就能创建一个Spectrum，反之亦然。其关系如图1-5所示。
Wave对象包括3个特性：包含信号参数的NumPy数组ys；信号开始采样和取值的时间点数组ts；每单位时间的采样数framerate。时间的单位通常是秒，但也可以不是。在我所列举的一些例子中，时间的单位都是天。
图1-5 thinkdsp中各类之间的关系
Wave还提供了3个只读属性：start、end和duration。如果读者修改了ts，这些属性也会相应改变。
读者可以通过直接修改ts和ys来改变波形。例如：
wave.ys *= 2
wave.ts += 1第一行将波形扩大了2倍，让它听起来更响亮。第二行将波形按时间移动了，让它晚一秒开始。
此外Wave还提供了执行通用操作的方法。例如，与上面相同的两个变换可以写成以下形式：
wave.scale(2)
wave.shift(1)读者可以在http://greenteapress.com/thinkdsp.html阅读关于这些方法和其他内容的文档。
Signal是一个父类，它向各种信号提供通用的函数，比如make_wave。子类继承了这些方法并提供evaluate，也就是在给定的时间序列内对信号取值。
例如，Sinusoid是Signal的子类，定义如下：
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__init__的参数有：
freq
频率，其含义为每秒周期数，单位是Hz。
amp
幅度。幅度的单位比较随意，通常设定为1.0，对应为传声器的最大输入或给扬声器的最大输出。
offset
其含义为信号周期的起始。offset的单位为弧度。
func
它是一个Python函数，用来对指定时间点的信号求值。通常它不是np.sin就是np.cos，对应的分别是正弦信号和余弦信号。
跟很多初始化方法一样，它也是把参数存起来以备未来使用。
Signal提供了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是每秒帧数（采样数）。
n代表采样的数量，而ts是采样时间的NumPy数组。
为计算ys，make_wave需要引用evaluate，它是由Sinusoid提供的:
def evaluate(self, ts):
phases = PI2 * self.freq * ts + self.offset
ys = self.amp * self.func(phases)
return ys下面我们来逐步地解析这个函数。
1．self.freq是频率，代表每秒周期数，而ts的各个元素都是以秒计的，所以它们的乘积是从起始时间开始的周期数。
2．PI2是一个常数，其值为2π。乘上PI2之后就把周期转换成了相位。读者可以把相位理解为弧度形式的“从起始时间开始的周期数”，而每个周期的弧度是2π。
3．self.offset是t=0时刻的相位。其作用是把信号在时域上向左或者向右移动一定距离。
4．如果self.func是np.sin或者np.cos, 其值便是处于–1~1。
5．乘上self.amp之后产生的信号范围为–self.amp到self.amp。
在数学意义上，evaluate的形式如下：
其中A是幅度，f是频率，t是时间，φ0是相位偏移。这里看起来像是我写了大段的代码，而仅仅为了对一个简单的表达式求值，但是正如我们将要看到的，这段代码提供了处理所有信号的构架，而不仅是三级函数信号的。
读者在开始下面的练习之前，应该先按照前言的“代码示例的使用”说明来下载本书的代码。
这些练习的答案在chap01.soln.ipynb中。
练习1-1
如果读者已经有了Jupyter，请下载chap01.ipynb，仔细阅读代码，然后运行其中的例子。读者可以在http://tinyurl.com/thinkdsp01阅读该笔记。
练习1-2
通过访问http://freesound.org/下载音高标定好的声音样本，这些样本可以是音乐、讲话或者其他声音形式的。选择一段长约半秒的音高稳定的片段，计算并绘制出你所选择的片段的频谱。读者能据此看出声音的音色和频谱中的谐音结构之间的什么联系呢？
使用high_pass、low_pass和band_stop来过滤掉部分谐音。然后把频谱转换成回声波再听一下。读者所听到的声音和所采取的频谱操作有什么关系呢？
练习1-3
创建SinSignal和CosSignal对象，然后把它们加和起来合成一个复合信号。对这个信号求值得到一个波形，然后听一下。计算其频谱并绘制出来。当在其中加入一些与基频不是倍数关系的频率元素时，听到的会是什么呢？
练习1-4
写一个名为stretch的函数，它通过Wave和伸展因素两个参数来修改ts和framerate，从而对波形进行修改，使它加快或变慢。提示：只需两行代码就可以实现。
[1] 频谱”的复数形式通常写作“spectra”，但是我更喜欢用标准的英语复数形式（原文此处作者用了spectrums）。如果读者对“spectra”更熟悉，希望我的用法看起来不会太奇怪。
本章中我将介绍多种波形。通过观察它们的频谱来理解它们的谐波结构，也就是构成它们的三角函数集合的结构。
本章还要介绍数字信号处理中最重要的现象之一：混叠。同时解释一下Spectrum类的工作原理。
本章的对应代码是chap02.ipynb, 它位于本书的仓库中（参见前言中的“代码示例的使用”）。读者在http://tinyurl.com/thinkdsp02也可以查看到。
一个正弦波仅包含一个频率元素，所以其频谱只有一个尖峰。对于更加复杂的波形，例如图1-2所示的小提琴的录音，其离散傅里叶变换的结果中有很多尖峰。在本章中，我们将研究波形和其对应的频谱之间的关系。
我们从三角波开始，三角波就像是正弦波的直线版本。图2-1显示的三角波的频率是200 Hz。
读者可使用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 ysTriangleSignal从Sinusoid继承了__init__, 所以它们使用的参数是一样的，都有freq、amp和offest。
图2-1 200Hz的三角波片段
唯一的差别在于evalute。如我们在前面看到过的，ts是我们想对信号求值的采样时间的序列。
产生三角波的方法很多。其细节并不重要，但我们需要知道evaluate的工作原理。
1．cycles是自起始时间的周期数。np.modf是将周期的数量分解成小数部分和整数部分，前者被存于frac，后者被忽略了[1]。
2．frac序列介于0~1，其频率是给定的。减去0.5之后其值介于-0.5~0.5。对其取绝对值就得到了在0~0.5穿梭的波形。
3．unbias将波形下移使得其中心位于0，然后normalize将波形归一化到给定的振幅amp。
以下是产生图2-1的代码：
signal = thinkdsp.TriangleSignal(200)
signal.plot()然后我们可以使用这里的Signal来创建Wave，从而使用所创建的Wave来生成Spectrum：
wave = signal.make_wave(duration=0.5, framerate=10000)
spectrum = wave.make_spectrum()
spectrum.plot()图2-2显示的是其结果的两个视图。下图中的右图经过了缩放，使得谐波显示得更清楚。正如所预测的那样，最高峰位于基频200 Hz，频率是基频200的倍数的谐波频率峰也被显示在图上。
图2-2 频率为200 Hz的三角波信号的频谱，以两种不同的纵坐标尺度显示。右侧图中切去了基频的部分，从而更清楚地显示谐波
但令人奇怪的是，在偶数倍频率，比如400 Hz、800 Hz等的位置没有峰值。三角波的谐波全都是基频的奇数倍，比如600 Hz、1 000 Hz、1 400 Hz等。
这个频谱的另一个特征在于谐波的幅度与频率的关系。各谐波的幅度随着频率的平方等比例地衰减。比如，最前面的两个谐波之间的频率（分别为200 Hz和600 Hz）比率为3，幅度比为9。接下来的两个谐波（分别为600 Hz和1 000 Hz）频率比为1.7，对应的幅度比大约是1.72 = 2.9。这种关系被称为谐波结构。
thinkdsp还提供了SquareSignal，用来处理方波。其类的定义如下：
class SquareSignal(Sinusoid):
def evaluate(self, ts):
cycles = self.freq * ts + self.offset / PI2
frac, _ = np.modf(cycles)
ys = self.amp * np.sign(unbias(frac))
return ys和SquareSignal一样，它从Sinusoid继承了__init__，所以它们的参数是一样的。
同时evaluate的方式也类似。同样，cycles是从起始时间开始的周期数目，而frac是小数部分，其值在每个周期介于0~1。
Unbias将frac移至-0.5~0.5，然后np.sign将负值投射到-1，而正值投射到1。乘以amp之后就得到取值于-amp到+amp之间的方波。
图2-3显示的是频率100 Hz的方波的3个周期，图2-4显示的是其频谱。
和三角波一样，方波也只含有奇数谐波，这也就是其频谱中在300 Hz、500 Hz和700 Hz这样的值处有尖峰的原因。但相比来说，它的谐波衰减得要慢一些。也就是，其幅度随频率（而不是频率的平方）等比例地降低。
读者可以通过本章结尾位置的练习来探索其他的波形以及对应的谐波结构。
图2-3 100 Hz频率的方波的一部分
图2-4 100 Hz频率方波的频谱
必须承认，在前面的章节中，为了避免展示给读者容易困惑的内容，我精心挑选了例子来展示。但是现在的内容开始让人困惑了。
图2-5显示的是一个频率1 100 Hz的三角波的频谱，采样率为10 000帧每秒。同样，右图显示的是放大之后的谐波。
图2-5 一个频率为1 100 Hz的三角波的频谱，采样率为10 000帧每秒。右图显示的是放大之后的谐波
这个波的谐波应该位于3 300 Hz、5 500 Hz、7 700 Hz和9 900 Hz。在上图中，和预期的一样，1 100 Hz和3 300 Hz处有尖峰，但是第三个尖峰在4 500 Hz，而不是5 500 Hz。第4个尖峰在2 300 Hz，而不是7 700 Hz。如果读者细看的话，应对位于9 900 Hz处的尖峰现在实际上在100 Hz处。这是什么问题呢?
问题在于，如果对信号在离散的时间点上进行取值，就会失去在采样点之间的信息。对于低频元素，这不是问题，因为会在每个周期中采样多次。
但是如果对频率为5 000 Hz的信号每秒采样10 000次，那在每个周期之中就仅有两个采样。这仅仅是恰好够用，但如果频率更高的话，就不够了。
为了演示，我们来产生频率4 500 Hz和5 500 Hz的余弦信号，然后对它们每秒采样10 000次：
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所示。在绘制这些信号时我使用的是灰色线，采样用的是垂直线表示，以便于比较这两个波形。这样问题应该清楚了：虽然两个信号是不同的，但是波形是一样的！
图2-6 频率为4 500 Hz和5 500 Hz的余弦信号，每秒采样10 000次。信号是不同的，但采样结果是一样的
当我们对一个5 500 Hz的信号进行10 000帧每秒的采样时，其结果与4 500 Hz的无法区分。同样，7 700 Hz的信号也无法与2 300 Hz的信号区分，9 900 Hz的信号无法与100 Hz的区分。
这种效应被称为混叠，也就是采样高频信号后，其结果和采样低频信号的一样。
在这个例子中，我们能测定的最高频率是5 000 Hz，也就是采样率的一半。高于5 000 Hz的信号被折叠到5 000 Hz以下，这也就是这个频率有时被称为“折叠频率”的原因。另外一些时候也被称为奈奎斯特频率，参见http://en.wikipedia.org/wiki/ Nyquist_ frequency。
如果混叠频率低于0，折叠的模式一样。比如，1 100 Hz的三角波信号的第5个谐波位于12 100 Hz。以5 000 Hz折叠，就会出现在-2 100 Hz，但是它在0 Hz处被再次折叠，从而回到2 100 Hz。实际上，读者能在图2-4的2 100 Hz处观察到一个小峰值，而下一个在4 300 Hz。
我们已经多次讲过make_spectrum的Wave方式。这里有一个实现（稍后再关注细节）：
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是一个Wave对象。n是波形的采样数，而d是采样频率的倒数，也就是样本间的时间差。
np.fft是NumPy模块，提供快速傅里叶变换（FFT）的相关函数，这是一种计算离散傅里叶变化的高效算法。
make_spectrum使用了rfft, 其意义是“傅里叶变换实数部分（real FFT）”，因为Wave只有实数值，而没有复数值。之后我们将了解完整的FFT，它能处理复数信号（见7.9节）。我将rfft的结果赋值给hs，这是一个NumPy复数数组，表示的是波形中各个频率元素的振幅和相位差。
我将rfftfreq的结果赋值给fs，这是一个包含对应hs各元素的频率的数组。
为了便于理解hs的值，可以参考以下两种理解复数的方式。
。你可以将x和y理解成笛卡尔坐标系下的数值。
的乘积
，其中A是量值，
是弧度角，也被称为“幅角”。读者可以在极角坐标系下理解A和。hs中的每一个值对应一个频率元素。其量值与对应的元素的振幅成比例，其角度为相位差。
Spectrum类提供了两个只读属性：amps和angles。它们返回的是NumPy数组，代表的分别是hs的振幅和角度。在我绘制Spectrum对象时，通常绘制的是fs对应的amps。有时绘制fs对应的angles也是有用的。
虽然查看hs的实数和虚数部分总是令人神往，但是读者几乎完全不需要这样做。我建议将DFT想象为振幅和角度的向量，只不过恰好以复数的形式来表示了。
读者可以通过直接使用hs来修改Spetrum，例如：
spectrum.hs *= 2
spectrum.hs[spectrum.fs > cutoff] = 0第一行将hs的元素乘以2，这将所有匀速的振幅都变成原来的两倍。第二行将hs元素中频率超过截至频率的所有元素都置零。
不过Spectrum还提供了如下的操作方式：
spectrum.scale(2)
spectrum.low_pass(cutoff)读者可以在http://greenteapress.com/thinkdsp.html阅读以上所述的这些方式的相关文档。
现在，读者应该对Signal、Wave和Spectrum这些类的工作方式有了更好的了解了，但我还没有解释快速傅里叶变换的原理。这将在之后几章讲解。
这些练习的答案在chap02soln.ipynb。
练习2-1
如果读者使用Jupyter, 请下载chap02.ipynb。然后试试里面的例子。读者可以在http://tinyurl.com/thinkdsp02阅读相关笔记。
练习2-2
锯齿信号的波形从-1线性地增加到1，然后下降到-1，如此重复。详见http://en.wikipedia.org/wiki/Sawtooth_wave。
写一个名为SawtoothSignal的类，以扩展Signal并提供对锯齿波信号求值的evaluate方式。
计算锯齿信号的频谱。其谐波结构跟三角波和方波谐波结构相比有什么特点？
练习2-3
创建一个1 100 Hz的方波型号，然后创建一个Wave来采样它，采样帧率为10 000帧每秒。如果读者绘制其频谱，可以看到大部分谐波都混叠了。读者在听这个波形时，能听到混叠的谐波吗？
练习2-4
如果你有一个频谱对象Spectrum, 打印出spectrum.fs的前几个值。你会看到它们从0开始。所以spectrum.hs[0]是频率为0元素的量值。不过这意味着什么呢？
尝试下面的试验：
1．创建一个频率440 Hz的三角波，然后创建一个时长0.01 s的Wave，绘制其波形。
2．创建一个Spectrum对象并打印spectrum.hs[0]。这个元素的振幅和相位是多少？
3．设定spectrum.hs[0] = 100。这对波形的影响是什么？提示：Spectrum提供了make_wave方式，它根据Spectrum来计算对应的Wave。
练习2-5
写一个以Spectrum为参数的函数，修改它，将hs的各个元素除以对应的频率fs。提示：由于除以0是无法定义的，读者可以设定spectrum.hs[0] = 0。
使用方波、三角波和锯齿波验证你的函数。
1．计算Spectrum并绘制它。
2．使用你的函数修改Spectrum，再绘制它。
3．使用Spectrum.make_wave从修改过的Spectrum绘制一个波形，试听它。这个操作对信号的影响是什么？
练习2-6
三角波和方波只有奇数谐波，锯齿波既有奇数的也有偶数的谐波。方波和锯齿波的谐波按1/f比例降低，三角波的则按1/（f2）比例降低。你能找到一个波形，含有奇数和偶数谐波，并且按1/（f2）比例降低吗？
提示：读者有两种方法来达成目的。既可以通过叠加三角波来做到，也可以以一个类似于你想要的波形开始，然后逐步修改。
[1] 下划线作为变量是一种习惯做法，其意思是“我并不打算使用这个值”。
