书名：Python数字信号处理
ISBN：978-7-115-66924-7
本书由人民邮电出版社发行数字版。版权所有，侵权必究。
您购买的人民邮电出版社电子书仅供您个人使用，未经授权，不得以任何方式复制和传播本书内容。
我们愿意相信读者具有这样的良知和觉悟，与我们共同保护知识产权。
如果购买者有侵权行为，我们可能对该用户实施包括但不限于关闭该帐号等维权措施，并可能追究法律责任。
著 赵春江
责任编辑 张 涛
人民邮电出版社出版发行 北京市丰台区成寿寺路11号
邮编 100164 电子邮件 315@ptpress.com.cn
网址 http://www.ptpress.com.cn
读者服务热线：(010)81055410
反盗版热线：(010)81055315
本书从实用的角度出发，深入浅出地介绍数字信号处理的基本理论和方法，并通过大量实例详细讲解Python在数字信号处理中的应用。本书首先介绍Python的基础知识，以及本书用到的NumPy库、Matplotlib库和SciPy库；然后讲述数字信号处理的基本理论和方法，包括离散时间信号与系统的时域分析、信号与系统的频谱分析、离散傅里叶变换及其快速算法、数字滤波器的结构、IIR数字滤波器的设计和FIR数字滤波器的设计。
本书可作为高等学校数字信号处理相关课程的工具书，也可作为广大科研工程技术人员的参考用书。
多年来，数字信号处理已成为电子信息等相关学科的重要技术，并被应用于许多领域，如数字通信、医学影像、音视频系统、消费电子、机器人、遥感、金融等。数字信号处理是把信号用数字或符号表示，并利用数值计算的方法进行处理，以达到提取有用信息的目的。数字信号处理以算法为核心，内容繁多、概念抽象、公式复杂、理论性强。
Python是一种强大的通用编程语言，被用于网络开发、数据科学分析、软件原型创建等。Python之所以被广泛使用，是因为它具有丰富且易于理解的语法和众多的可用开源包。
本书将介绍Python在数字信号处理方面的应用，把抽象的理论通过实例具体化，内容主要包括Python及常用库（NumPy、Matplotlib、SciPy）的基础知识、离散时间信号与系统的时域分析、信号与系统的频谱分析、离散傅里叶变换及其快速算法、数字滤波器的结构、IIR数字滤波器的设计和FIR数字滤波器的设计。
本书由浅入深、循序渐进地介绍数字信号处理的基本理论和方法，并通过大量具有代表性的实例阐述Python在数字信号处理中的应用。这些实例内容新颖，本书尽量详尽地给出实例的完成步骤，将Python的使用方法与技巧呈现给读者，使读者在阅读时能够快速掌握书中所讲内容。通过对本书的阅读和学习，读者不仅可以加深对数字信号处理基本理论和方法的理解，还可以更好地掌握Python的用法。
为便于读者更好地阅读本书中的代码，区分它们的不同作用，本书的版式约定如下。
• 完整的、可运行的程序上下有两条线。示例如下。
software= ['Autocad', 'Photoshop', 'Indesign'] numbers = [45, 341] empty = [ ] print(software, numbers, empty)
• 可重复调用的函数左侧有竖双细线。示例如下。
import numpy as np
def stepseq(n1, n2, n0=0):
n = np.arange(n1, n2+1)
x = np.zeros(n2-n1+1)
x[n0-n1:] = 1
return x, n• 对Python内部函数的详细解释放在文本框中。示例如下。
| scipy.signal.unit_impulse(shape, idx=None, dtype=<class 'float'>)：用于生成单位脉冲序列。 参数如下。 • shape：int或int类型的元组，表示一维输出的样本数量。 …… |
在本书的出版过程中编者得到了所在单位的相关领导和同事（胡学友、胡国华、顾涓涓、张胜、倪敏生、李瑶、黄慧等）的大力支持，在此表示诚挚的感谢。
由于编者水平有限，书中难免会有不足之处，真诚欢迎各位读者批评指正。
注：书中的图为软件自动生成，为了保持原有样式，编者没有对图中英文进行翻译，在此特别说明。
编者
Python是一种编程语言，语法简洁，易于学习。Python程序可以在不同的计算机上运行。它可用于Web应用程序、科学应用程序、桌面应用程序、教育和通用软件应用程序等。
要运行Python程序，我们需要使用Python解释器，Python解释器用于执行Python代码（或称程序），如图1.1所示。一个程序可以是一个或多个Python文件，Python文件可以包含其他文件或模块。
图1.1 运行Python程序
如果我们不想在终端上工作，可以使用IDE（Integrated Development Environment，集成开发环境）。它是一个图形编辑器，我们可以在其中输入代码、处理多个文件、运行代码等。
本书中的所有代码都是在名为Spyder的IDE内编辑并运行的。Spyder是一个开源的跨平台的IDE。Spyder集成了NumPy、SciPy、Matplotlib与IPython，以及其他开源软件。与其他科学数值分析专用IDE（如MATLAB）相比，Spyder有下列特色：开放源代码、兼容于非自由软件许可协议。所以利用Spyder编写有关数字信号处理的代码非常合适。
本书所使用的Spyder的版本为5.2.2，其集成的Python的版本为3.7.9、NumPy的版本为1.19.3、SciPy的版本为1.7.3、Matplotlib的版本为3.5.1。
Python文件使用.py作为扩展名，如一个简单的程序，我们把它命名为hello.py，其内容为：
print('Hello world!')运行后系统会输出如下结果。
Hello world!
Python支持不同数据类型的变量，如整数、浮点数和字符串等。我们不需要事先指定变量的数据类型，就可以简单地将任何值分配给变量。Python将根据我们分配给变量的值来确定数据类型，如创建一个变量x，令x=3，则Python假定它的数据类型是整型。
变量名必须以字母（大写或小写）或下画线开头，不能以数字开头，并且区分大小写。
【例1-1】变量。
程序如下。
x = 3 #整数 f = 3.1415926 #浮点数 name = "Python" #字符串 print(x) print(f) print(name) combination = name + " " + name print(combination) sum = f + f print(sum)
运行结果如下。
3 3.1415926 Python Python Python 6.2831852
☐
列表很像字符串，它包含一系列的值。在字符串中，值是字符型的，而列表中的元素可以是任何类型的。值为列表的列表称为嵌套列表。不包含任何元素的列表称为空列表，我们可以使用空的方括号[ ]来创建一个空列表。
【例1-2】列表。
程序如下。
software= ['Autocad', 'Photoshop', 'Indesign'] numbers = [45, 341] empty = [ ] print(software, numbers, empty)
运行结果如下。
['Autocad', 'Photoshop', 'Indesign'] [45, 341] []
☐
加法、减法、乘法和除法运算符可以用于Python中的数值运算，使用两个乘法运算符表示指数运算，Python还支持将字符串与数字相乘以形成重复字符串。
【例1-3】运算符。
程序如下。
remainder = 14 % 3 #求余数 print(remainder) square = 3 ** 2 #求平方 cube = 3 ** 3 #求立方 print(square) print(cube) lol= "ha" * 10 #重复字符串 print(lol) even_numbers = [2,4,6,8] odd_numbers = [1,3,5,7] all_numbers = odd_numbers + even_numbers #列表相加 print(all_numbers)
运行结果如下。
2 9 27 hahahahahahahahahaha [1, 3, 5, 7, 2, 4, 6, 8]
☐
Python使用布尔逻辑来评估条件，当一个表达式被比较或评估时，会返回布尔值True或False。and和or布尔运算符用于构建复杂的布尔表达式。in运算符可用于检查指定对象是否存在于可迭代对象容器中。
条件判断if语句的使用涉及代码块，Python利用缩进来定义代码块，而不是括号，标准的Python缩进是4个空格。
【例1-4】条件判断。
程序如下。
name = "John"
age = 23
if name == "John" and age == 23: #与运算
print("Your name is John, and you are also 23 years old.")
if name == "John" or name == "Rick": #或运算
print("Your name is either John or Rick.")
if name in ["John", "Rick"]:
print("Your name is either John or Rick.")
x = 3
if x == 2:
print('two')
elif x == 3:
print('three')
elif x == 4:
print('four')
else:
print('something else')运行结果如下。
Your name is John, and you are also 23 years old. Your name is either John or Rick. Your name is either John or Rick. three
☐
Python中有两种循环，for循环和while循环。for循环遍历给定的序列，而对于while循环，只要满足某个布尔条件，它就会重复执行。break用于退出for循环或while循环，而continue则用于跳过当前语句块，并返回到for或while循环。与其他语言不同，Python还可以使用else循环，即当for或while循环的循环条件失败时，将执行else循环中的代码。
【例1-5】循环语句。
程序如下。
primes = [2, 3, 5, 7]
for prime in primes:
print(prime)
count = 0
while count < 5:
print(count)
count += 1运行结果如下。
2 3 5 7 0 1 2 3 4
☐
函数是一种将代码划分为有用语句块的便捷方式，它允许我们对代码进行排序，使其更具可读性、重用性。同时，函数也是定义接口的关键方法。
【例1-6】函数。
程序如下。
def my_function():
print("Hello From My Function!")
def my_function_with_args(username, greeting):
print("Hello, %s, From My Function! I wish you %s"%(username, greeting))
def sum_two_numbers(a, b):
return a + b
my_function()
my_function_with_args("John Doe", "a great year!")
x = sum_two_numbers(1, 2)运行结果如下。
Hello From My Function! Hello, John Doe, From My Function! I wish you a great year!
☐
程序中的模块是具有特定功能的语句。例如一个模块将负责绘图，另一个模块将负责数值计算。每个模块也都是一个不同的文件。Python模块本质上是扩展名为.py的Python文件，模块的名称将是文件的名称，Python模块可以定义和实现一组函数、类或变量。假设我们将创建两个Python文件，如下所示。
myprog/draw.py myprog/calc.py
calc.py将实现数值计算，并将使用draw.py中的函数draw_prog。使用import命令可导入其他模块，在该例中，calc.py可能如下。
import draw
def calc_value():
...
def main():
result = calc_value()
draw.draw_prog(result)#如果执行该文件，main()将执行
if __name__ == '__main__':
main()
我们也可以使用from命令将draw_prog直接导入主文件代码命令空间中：
from draw import draw_prog
def main():
result = calc_value()
draw_prog(result)NumPy是一个用于处理数组的Python库，它还具有处理线性代数、傅里叶变换和矩阵等领域的运算的能力。数组在科学运算中非常常用，因为数组在运算速度和所需资源方面比其他数据类型更有优势。而NumPy数组是Python列表的最佳替代品，它的主要优点就是速度快、易于使用。
由于数字信号处理主要以一维信号为研究对象，因此我们重点介绍一维数组的相关性质和运算。
【例1-7】NumPy数组。
程序如下。
import numpy as np height = [1.87, 1.87, 1.82, 1.91, 1.90, 1.85] weight = [81.65, 97.52, 95.25, 92.98, 86.18, 88.45] np_height = np.array(height) np_weight = np.array(weight) print(type(np_height)) bmi = np_weight / np_height ** 2 print(bmi) print(bmi[bmi > 25])
运行结果如下。
<class 'NumPy.ndarray'> [23.34925219 27.88755755 28.75558507 25.48723993 23.87257618 25.84368152] [27.88755755 28.75558507 25.48723993 25.84368152]
☐
【例1-8】创建各类特殊数组。
程序如下。
import numpy as np a = np.zeros(3) #创建元素为0的数组 print(a) b = np.ones(3) #创建元素为1的数组 print(b) c = np.empty(3) #创建一个空数组 print(c) d = np.arange(4) #创建一个包含一系列元素的数组 print(d) #创建指定了第一个元素、最后一个元素和步长的数组（不含最后那个元素） e = np.arange(2, 10, 2) print(e) f = np.linspace(0, 10, 5) #创建一个所含元素为指定间隔的线性分布的数组 print(f) g = np.random.random(4) #创建一个随机数组 print(g)
运行结果如下。
[0. 0. 0.] [1. 1. 1.] [1. 1. 1.] [0 1 2 3] [2 4 6 8] [0. 2.5 5. 7.5 10.] [0.94704006 0.24872617 0.02854942 0.04201419]
☐
【例1-9】数组的索引和分割。
程序如下。
import numpy as np a = np.arange(10) print(a) print(a[5]) #获取数组a中的第5个元素（从0开始算起） print(a[-3]) #获取数组a中的倒数第3个元素 #从数组a中分割出一个数组，其元素为数组a中的第1个至第5个元素（不含） print(a[1 : 5]) #从数组a中分割出一个数组，其元素为数组a中第7个元素（不含）之前的所有元素 print(a[ : 7]) #从数组a中分割出一个数组，其元素为数组a中第3个元素之后的所有元素 print(a[3 : ]) #从数组a中分割出一个数组，其元素为数组a中第1个元素至第7个元素（不含） #步长为2 print(a[1 : 7 : 2]) #从数组a中分割出一个数组，其元素为数组a中倒数第4个元素至倒数第2个元素（不含） print(a[-4 : -2]) #从数组a中分割出一个数组，其元素为数组a中所有元素，步长为3 print(a[ : : 3]) #从数组a中分割出一个数组，其元素为数组a中所有元素，但为逆序 #即对数组a进行翻转 print(a[ : : -1]) print(a[ : : -2]) #对数组a进行翻转，步长为2
运行结果如下。
[0 1 2 3 4 5 6 7 8 9] 5 7 [1 2 3 4] [0 1 2 3 4 5 6] [3 4 5 6 7 8 9] [1 3 5] [6 7] [0 3 6 9] [9 8 7 6 5 4 3 2 1 0] [9 7 5 3 1]
☐
【例1-10】数组的合并、拆分及元素的删除。
程序如下。
import numpy as np a = np.array([1, 2, 3]) b = np.array([4, 5, 6]) c = np.append(a, b) #合并数组 print(c) print(np.array_split(c, 3)) #拆分数组 print(np.delete(c, 3)) #删除指定元素
运行结果如下。
[1 2 3 4 5 6] [array([1, 2]), array([3, 4]), array([5, 6])] [1 2 3 5 6]
☐
为了进行统计推导，有必要将计算得到的数据可视化，而Matplotlib是Python的可视化解决方案之一。Matplotlib是一个功能非常强大的绘图库。Matplotlib中常用的模块是pyplot，它提供了类似于MATLAB的接口。
【例1-11】绘制折线图。
程序如下。
import matplotlib.pyplot as plt
plt.plot([1, 2, 3, 4], [1, 4, 9, 16]) #依次连接各坐标点，组成一条折线
plt.title('line chart') #图形的名称
plt.xlabel('X label') #横坐标名称
plt.ylabel('Y label') #纵坐标名称
plt.grid() #绘制网格运行结果如图1.2所示。
☐
图1.2 例1-11的绘制折线图
【例1-12】在同一个坐标系内绘制两个不同的图形。
程序如下。
import numpy as np import matplotlib.pyplot as plt x = np.arange(1, 5) y1 = x**2 y2 = x**3 plt.plot(x, y1, '--', label='$x^2$') plt.plot(x, y2, '-.', label='$x^3$') plt.legend() #显示图例
运行结果如图1.3所示。
☐
图1.3 例1-12的同一个坐标系内绘制两个不同的图形
【例1-13】绘制多幅图形。
程序如下。
import numpy as np import matplotlib.pyplot as plt x = np.arange(1, 5) y1 = x**2 y2 = x**3 plt.subplot(211) #多幅图形的分布为2行1列，当前为第1幅 plt.stem(x, y1) #绘制垂直线段图 plt.subplot(212) #多幅图形的分布为2行1列，当前为第2幅 plt.stem(x, y2) #绘制垂直线段图
运行结果如图1.4所示。
☐
图1.4 例1-13绘制的多幅图形
Matplotlib还有另一种绘图方法，这种方法要比前面介绍的方法更灵活。
【例1-14】axes.plot绘图方法。
程序如下。
import numpy as np import matplotlib.pyplot as plt t = np.linspace(0, 2*np.pi, 400) s1 = np.sin(2 * np.pi * 0.5 * t) s2 = np.sin(2 * np.pi * 0.25 * t) fig, ax = plt.subplots(2, 1) ax[0].plot(t, s1) ax[1].plot(t, s2)
运行结果如图1.5所示。
☐
图1.5 例1-14的axes.plot绘图方法
数学涉及大量的概念，这些概念非常重要，但也很复杂。然而，Python提供了成熟的SciPy库，为我们解决了这个问题。
SciPy是一个开源的Python库，用于解决科学和数学问题。它建立在NumPy库之上，允许用户使用各种高级命令来操作和可视化数据。在数字信号处理领域，主要用到的是SciPy中的fftpack模块和signal模块。
【例1-15】求
的根。
程序如下。
from scipy.optimize import root
from numpy import cos
def eqn(x):
return x + cos(x)
myroot = root(eqn, 0)
print(myroot.x)运行结果如下。
[-0.73908513]
☐
