🥥Python数字信号处理

关键词

Python | 数字信号 | 离散傅里叶变换 | 光谱分析 | 循环巻积 | 等效噪声带宽 | 有限脉中响应滤波器 | 数学 | scipy | Matplotlib | Numpy

梗概

介绍

涵盖了用 Python 代码说明的信号处理中的基本概念，侧重于信号处理的核心、基本原理；对应的代码使用了科学 Python 工具链的核心功能；对于那些希望将信号处理代码迁移到 Python 的人来说，本文说明了可以简化这种转换的关键信号和绘图模块。对于那些已经熟悉科学 Python 工具链的人，本文阐述了信号处理中的基本概念，并为进一步的信号处理概念提供了途径

取样定理

以下是定理的通常陈述：如果函数$x(t)$不包含高于 B 赫兹的频率，则完全通过在一系列间隔为$1 /(2 B)$秒的点处给出其纵坐标来确定。

因为$x(t)$是从实线到实线的函数，任意两个连续样本之间有无限多的点，因此采样是数据的大量减少，因为它只需要有限数量的点来完全表征函数。 我们之前在傅立叶级数展开式中已经看到了将函数简化为离散数字集的概念，其中 (对于周期性$x(t)$) 我们有，

$a_{n}=\frac{1}{T} \int_{0}^{T} x(t) \exp \left(-j \omega_{n} t\right) d t,\qquad(1)$

相应的重建为：

$x(t)=\sum_{n} a_{n} \exp \left(j \omega_{n} t\right),\qquad(2)$

但是在这里我们通过对整个函数$x(t)$进行积分来生成离散点$a_n$，而不仅仅是在单个点上对其进行采样。 这意味着我们收集有关整个函数的信息来计算单个离散点$a_n$，而通过采样我们只是孤立地获取单个点。

另一方面，假设我们有一组样本$\left[x_{1}, x_{2}, \ldots, x_{N}\right]$，然后我们被要求重建函数。 我们会做什么？ 也许最自然的做法是在每个点之间画一条直线，就像在线性插值中一样。下面代码在单个周期内采集正弦样本并在样本之间画一条线。

from __future__ import division

 fig,ax = subplots()
 f = 1.0 # Hz, signal frequency
 fs = 5.0 # Hz, sampling rate (ie. >= 2*f)
 t = arange(-1,1+1/fs,1/fs) # sample interval, symmetric
                            # for convenience later
 x = sin(2*pi*f*t)
 ax.plot(t,x,’o-’)
 ax.set_xlabel(’Time’,fontsize=18)
 ax.set_ylabel(’Amplitude’,fontsize=18)

上述代码，第一行确保我们使用浮点除法而不是 Python 的默认整数除法，其中 1/4=0 而不是 0.25。

第 3 行设置图形和相应的轴。 这遵循使用 Matplotlib 的现代约定。 fig 是绑定到图形的对象， ax 是在该图形中绘制的轴。 我们希望这些分开的原因是更复杂的图可能在同一图中嵌入了多个轴。 因为 Matplotlib 可以在 GUI 窗口、文件中或嵌入其他应用程序（等等）中渲染绘图，所以假设您的绘图命令正在对“当前活动”图形进行操作并不是一个好习惯。 使用 fig 和 ax 避免了这个陷阱。

在第 6 行，arange 函数创建了一个 Numpy 数字数组，从 -1 开始，直到 1+1/fs，步长为 1/fs。 回想一下 Python 中的内置 range 命令产生一个半开区间（即不包括右端点），arange 也是如此。 下一行计算我们刚刚创建的点数组的正弦函数并返回一个 Numpy 数组。 请注意，虽然 Python 本身在内置 math 模块中带有正弦函数，但这个正弦是 Numpy 版本。不同之处在于 Numpy 版本可以摄取一个 Numpy 数组并产生相应的输出，而无需 math 模块的任何额外循环。

ax.plot 命令将绘图附加到轴 ax。 接下来的两行将文本标签放在 x 轴和 y 轴上。 这些标签功能还为文本格式提供了许多可能性，包括潜在的 LATEX 字体。

图1 显示，即使函数是最曲线的$(t=1 /(4 f)和t=3 /(4 f))$，我们拥有与其他任何地方相同的点密度。这是因为采样定理并没有指定我们应该在哪里采样，只要我们以周期性间隔采样。这意味着在正弦的上下斜率上， 哪些是线性的，哪些需要更少的样本来表征，我们将样本密度取为靠近弯曲峰。代码清单1放大到第一个峰来说明这一点

 fig,ax = subplots()
 ax.plot(t,x,’o-’)
 ax.axis(xmin = 1/(4*f)-1/fs*3,
				 xmax = 1/(4*f)+1/fs*3,
				 ymin = 0,
				 ymax = 1.1 )
 ax.set_xlabel(’Time’,fontsize=18)
 ax.set_ylabel(’Amplitude’,fontsize=18)

上述代码，创建了一个新图（例如线条、标记）。 标记规范与 MATLAB 相同，因此“o-”表示使用“o”标记并用实线连接它们。 下一行通过设置轴限制来放大绘图。 请注意，我们再次使用关键字函数样式，它清楚地表明了我们在轴上所做的更改。此代码对应图示2。

为了说明这一点，我们可以构造分段线性插值并使用 numpy.piecewise 比较近似的质量。代码清单2使用 hstack 函数分段构建 Numpy 所需的参数，该函数将 Numpy 数组水平“堆叠”为更大的 Numpy 数组。 Python 列表 tp 具有分段近似的区间，而 apprx 列表包含每个区间的相应线性近似。 logical_and 是必要的，因为我们想要一个元素明智的逻辑运算。 最后， tp 和 apprx 被分段输入，将它们打包在一起以供以后用作函数。

interval=[] # piecewise domains
apprx = [] # line on domains
 # build up points *evenly* inside of intervals
tp = hstack([linspace(t[i],t[i+1],20,False) for i in range(len(t)-1)])
 # construct arguments for piecewise
for i in range(len(t)-1):
	 interval.append(logical_and(t[i] <= tp,tp < t[i+1]))
	 apprx.append((x[i+1]-x[i])/(t[i+1]-t[i])*(tp[interval[-1]]-t[i]) + x[i])
 x_hat = piecewise(tp,interval,apprx) # piecewise linear approximation

完成所有设置后，我们可以检查插值中的平方误差。 代码清单3绘制了带有我们刚刚构建的线性插值的填充误差的正弦曲线。

 fig,ax1=subplots()
 # fill in the difference between the interpolant and the sine
 ax1.fill_between(tp,x_hat,sin(2*pi*f*tp),facecolor=’red’)
 ax1.set_xlabel(’Time’,fontsize=18)
 ax1.set_ylabel(’Amplitude’,fontsize=18)
 ax2 = ax1.twinx() # create clone of ax1
 sqe = (x_hat-sin(2*pi*f*tp))**2 #compute squared-error
 ax2.plot(tp, sqe,’r’)
 ax2.axis(xmin=-1,ymax= sqe.max() )
 ax2.set_ylabel(’Squared error’, color=’r’,fontsize=18)
 ax1.set_title(’Errors with Piecewise Linear Interpolant’,fontsize=18)

上述代码，使用 fill_between 函数用红色填充线性插值和正弦函数之间的凸面区域（facecolor=‘red’）。 下面的 twinx() 调用创建了当前轴（即 ax2）的副本，右侧带有 y 轴标签。 ax2.set_ylabel 调用在右侧附加 y 轴标签，其中颜色将字体颜色更改为红色。 拥有单独的轴变量（即 ax1 和 ax2）允许我们在代码中分别引用它们。 接下来，我们使用这个新创建的轴以红色 (‘r’) 绘制平方误差 (SE)。此代码对应图示3

离散时间傅立叶变换

介绍光谱分析

有限脉冲响应滤波器

🏈指点迷津 | Brief