用Python画图搞懂DFT和z变换一个11点滑动滤波器的可视化教程信号处理理论中的离散傅里叶变换DFT和z变换常常让初学者感到抽象难懂。本文将通过Python代码实现一个11点滑动滤波器的可视化分析帮助读者直观理解这些核心概念。我们将从零开始构建滤波器模型绘制三维系统函数图、零极点图和频率响应图让数学公式变成可视化的几何图形。1. 理论基础与Python环境准备1.1 理解滑动求和滤波器11点滑动求和滤波器是最简单的FIR滤波器之一其系统函数可以表示为H(z) 1 z⁻¹ z⁻² ... z⁻¹⁰这个表达式实际上是一个有限长度的几何级数。在时域中它相当于对输入信号的11个连续样本进行简单相加。这种滤波器常用于信号平滑处理能有效抑制高频噪声。1.2 DFT与z变换的关系DFT实际上是z变换在单位圆上的等间隔采样。具体来说z变换X(z) Σx[n]z⁻ⁿ定义在整个复平面DFTX[k] Σx[n]e^(-j2πkn/N)相当于在ze^(j2πk/N)处的采样1.3 Python环境配置我们需要以下Python库来实现可视化import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D建议使用Jupyter Notebook进行交互式实验可以实时观察图形变化。安装必要的库pip install numpy matplotlib2. 系统函数的三维可视化2.1 构建复平面网格为了绘制系统函数的三维图形我们需要在复平面上创建高密度采样点real_vals np.linspace(-1.5, 1.5, 500) imag_vals np.linspace(-1.5, 1.5, 500) real_grid, imag_grid np.meshgrid(real_vals, imag_vals) z real_grid 1j * imag_grid2.2 计算系统函数值对于每个复平面上的点z计算系统函数H(z)的值H_z np.zeros_like(z, dtypenp.complex128) for k in range(11): # 0到10 H_z z**(-k)特别处理z1处的奇异点mask (np.abs(z.real - 1) 1e-6) (np.abs(z.imag) 1e-6) H_z[mask] 11 0j2.3 绘制3D幅度响应使用matplotlib的3D功能绘制系统函数的幅度响应fig plt.figure(figsize(14, 10)) ax fig.add_subplot(111, projection3d) surf ax.plot_surface(real_grid, imag_grid, np.abs(H_z), cmaprainbow, alpha0.8) ax.set_xlabel(实部(Re)) ax.set_xlabel(虚部(Im)) ax.set_zlabel(|H(z)|) plt.title(11点滑动求和滤波器系统函数)这个三维图形清晰地展示了系统函数在整个复平面上的行为特别是单位圆附近的特性。3. 零极点分析与可视化3.1 解析零极点位置11点滑动求和滤波器的系统函数可以表示为H(z) (1 - z⁻¹¹)/(1 - z⁻¹) (z¹¹ - 1)/(z¹⁰(z - 1))由此可得零点z¹¹1的解即单位圆上的11等分点除去z1极点z0处的10阶极点以及z1处的一阶极点与零点相消3.2 Python实现零极点图绘制零极点分布的代码如下zeros np.exp(2j * np.pi * np.arange(1, 11) / 11) # 10个零点 poles np.zeros(10) # 10阶极点 plt.figure(figsize(8, 8)) ax plt.subplot(111) unit_circle plt.Circle((0,0), 1, fillFalse, colorblack, linestyle--) ax.add_patch(unit_circle) # 绘制零点蓝色圆圈和极点红色叉号 plt.plot(zeros.real, zeros.imag, bo, markersize8, label零点) plt.plot(poles.real, poles.imag, rx, markersize10, label极点) plt.axis(equal) plt.legend() plt.grid(True)零极点图直观展示了系统的稳定性特性。由于所有极点都在单位圆内原点这个滤波器是稳定的。4. 频率响应分析4.1 单位圆上的频率响应DTFT就是系统函数在单位圆上的取值。我们可以通过密集采样单位圆上的点来计算频率响应num_points 2000 frequencies np.linspace(-np.pi, np.pi, num_points) H_vals np.zeros(num_points, dtypenp.complex128) for i, omega in enumerate(frequencies): z np.exp(1j * omega) if np.abs(omega) 1e-10: # ω0处 H_vals[i] 11 else: H_vals[i] (1 - z**(-11)) / (1 - z**(-1))4.2 绘制幅频和相频特性plt.figure(figsize(12, 8)) plt.subplot(2, 1, 1) plt.plot(frequencies, np.abs(H_vals), b-) plt.title(幅频响应) plt.xlabel(频率(rad)) plt.ylabel(增益) plt.subplot(2, 1, 2) plt.plot(frequencies, np.angle(H_vals), r-) plt.title(相频响应) plt.xlabel(频率(rad)) plt.ylabel(相位(rad)) plt.tight_layout()从幅频响应曲线可以看出这是一个低通滤波器在ω0处增益最大随着频率增加增益逐渐减小。5. 综合分析与实际应用5.1 理解DFT与z变换的关系通过可视化我们可以直观看到DFT相当于在单位圆上等间隔采样频率响应DTFT是单位圆上连续的频率响应z变换提供了整个复平面的系统特性视图5.2 滤波器参数调整实验尝试修改滤波器长度观察图形变化def sliding_filter_response(N): frequencies np.linspace(-np.pi, np.pi, 1000) response np.zeros(1000, dtypenp.complex128) for i, omega in enumerate(frequencies): z np.exp(1j * omega) if np.abs(omega) 1e-10: response[i] N else: response[i] (1 - z**(-N)) / (1 - z**(-1)) return np.abs(response) for N in [5, 11, 21]: plt.plot(frequencies, sliding_filter_response(N), labelfN{N}) plt.legend() plt.title(不同长度滑动滤波器的幅频响应)可以看到随着滤波器长度的增加主瓣变窄频率分辨率提高旁瓣数量增加第一旁瓣相对幅度减小5.3 实际信号处理示例让我们用这个滤波器处理一个含噪声的信号t np.linspace(0, 1, 1000) signal np.sin(2*np.pi*5*t) 0.5*np.random.randn(1000) # 实现滑动平均滤波 filtered np.convolve(signal, np.ones(11)/11, modevalid) plt.figure(figsize(12, 4)) plt.plot(t, signal, b-, alpha0.5, label原始信号) plt.plot(t[5:-5], filtered, r-, label滤波后信号) plt.legend()这个简单的例子展示了滑动平均滤波器如何有效抑制高频噪声同时保留低频信号成分。
用Python画图搞懂DFT和z变换：一个11点滑动滤波器的可视化教程
用Python画图搞懂DFT和z变换一个11点滑动滤波器的可视化教程信号处理理论中的离散傅里叶变换DFT和z变换常常让初学者感到抽象难懂。本文将通过Python代码实现一个11点滑动滤波器的可视化分析帮助读者直观理解这些核心概念。我们将从零开始构建滤波器模型绘制三维系统函数图、零极点图和频率响应图让数学公式变成可视化的几何图形。1. 理论基础与Python环境准备1.1 理解滑动求和滤波器11点滑动求和滤波器是最简单的FIR滤波器之一其系统函数可以表示为H(z) 1 z⁻¹ z⁻² ... z⁻¹⁰这个表达式实际上是一个有限长度的几何级数。在时域中它相当于对输入信号的11个连续样本进行简单相加。这种滤波器常用于信号平滑处理能有效抑制高频噪声。1.2 DFT与z变换的关系DFT实际上是z变换在单位圆上的等间隔采样。具体来说z变换X(z) Σx[n]z⁻ⁿ定义在整个复平面DFTX[k] Σx[n]e^(-j2πkn/N)相当于在ze^(j2πk/N)处的采样1.3 Python环境配置我们需要以下Python库来实现可视化import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D建议使用Jupyter Notebook进行交互式实验可以实时观察图形变化。安装必要的库pip install numpy matplotlib2. 系统函数的三维可视化2.1 构建复平面网格为了绘制系统函数的三维图形我们需要在复平面上创建高密度采样点real_vals np.linspace(-1.5, 1.5, 500) imag_vals np.linspace(-1.5, 1.5, 500) real_grid, imag_grid np.meshgrid(real_vals, imag_vals) z real_grid 1j * imag_grid2.2 计算系统函数值对于每个复平面上的点z计算系统函数H(z)的值H_z np.zeros_like(z, dtypenp.complex128) for k in range(11): # 0到10 H_z z**(-k)特别处理z1处的奇异点mask (np.abs(z.real - 1) 1e-6) (np.abs(z.imag) 1e-6) H_z[mask] 11 0j2.3 绘制3D幅度响应使用matplotlib的3D功能绘制系统函数的幅度响应fig plt.figure(figsize(14, 10)) ax fig.add_subplot(111, projection3d) surf ax.plot_surface(real_grid, imag_grid, np.abs(H_z), cmaprainbow, alpha0.8) ax.set_xlabel(实部(Re)) ax.set_xlabel(虚部(Im)) ax.set_zlabel(|H(z)|) plt.title(11点滑动求和滤波器系统函数)这个三维图形清晰地展示了系统函数在整个复平面上的行为特别是单位圆附近的特性。3. 零极点分析与可视化3.1 解析零极点位置11点滑动求和滤波器的系统函数可以表示为H(z) (1 - z⁻¹¹)/(1 - z⁻¹) (z¹¹ - 1)/(z¹⁰(z - 1))由此可得零点z¹¹1的解即单位圆上的11等分点除去z1极点z0处的10阶极点以及z1处的一阶极点与零点相消3.2 Python实现零极点图绘制零极点分布的代码如下zeros np.exp(2j * np.pi * np.arange(1, 11) / 11) # 10个零点 poles np.zeros(10) # 10阶极点 plt.figure(figsize(8, 8)) ax plt.subplot(111) unit_circle plt.Circle((0,0), 1, fillFalse, colorblack, linestyle--) ax.add_patch(unit_circle) # 绘制零点蓝色圆圈和极点红色叉号 plt.plot(zeros.real, zeros.imag, bo, markersize8, label零点) plt.plot(poles.real, poles.imag, rx, markersize10, label极点) plt.axis(equal) plt.legend() plt.grid(True)零极点图直观展示了系统的稳定性特性。由于所有极点都在单位圆内原点这个滤波器是稳定的。4. 频率响应分析4.1 单位圆上的频率响应DTFT就是系统函数在单位圆上的取值。我们可以通过密集采样单位圆上的点来计算频率响应num_points 2000 frequencies np.linspace(-np.pi, np.pi, num_points) H_vals np.zeros(num_points, dtypenp.complex128) for i, omega in enumerate(frequencies): z np.exp(1j * omega) if np.abs(omega) 1e-10: # ω0处 H_vals[i] 11 else: H_vals[i] (1 - z**(-11)) / (1 - z**(-1))4.2 绘制幅频和相频特性plt.figure(figsize(12, 8)) plt.subplot(2, 1, 1) plt.plot(frequencies, np.abs(H_vals), b-) plt.title(幅频响应) plt.xlabel(频率(rad)) plt.ylabel(增益) plt.subplot(2, 1, 2) plt.plot(frequencies, np.angle(H_vals), r-) plt.title(相频响应) plt.xlabel(频率(rad)) plt.ylabel(相位(rad)) plt.tight_layout()从幅频响应曲线可以看出这是一个低通滤波器在ω0处增益最大随着频率增加增益逐渐减小。5. 综合分析与实际应用5.1 理解DFT与z变换的关系通过可视化我们可以直观看到DFT相当于在单位圆上等间隔采样频率响应DTFT是单位圆上连续的频率响应z变换提供了整个复平面的系统特性视图5.2 滤波器参数调整实验尝试修改滤波器长度观察图形变化def sliding_filter_response(N): frequencies np.linspace(-np.pi, np.pi, 1000) response np.zeros(1000, dtypenp.complex128) for i, omega in enumerate(frequencies): z np.exp(1j * omega) if np.abs(omega) 1e-10: response[i] N else: response[i] (1 - z**(-N)) / (1 - z**(-1)) return np.abs(response) for N in [5, 11, 21]: plt.plot(frequencies, sliding_filter_response(N), labelfN{N}) plt.legend() plt.title(不同长度滑动滤波器的幅频响应)可以看到随着滤波器长度的增加主瓣变窄频率分辨率提高旁瓣数量增加第一旁瓣相对幅度减小5.3 实际信号处理示例让我们用这个滤波器处理一个含噪声的信号t np.linspace(0, 1, 1000) signal np.sin(2*np.pi*5*t) 0.5*np.random.randn(1000) # 实现滑动平均滤波 filtered np.convolve(signal, np.ones(11)/11, modevalid) plt.figure(figsize(12, 4)) plt.plot(t, signal, b-, alpha0.5, label原始信号) plt.plot(t[5:-5], filtered, r-, label滤波后信号) plt.legend()这个简单的例子展示了滑动平均滤波器如何有效抑制高频噪声同时保留低频信号成分。
