网站建设资源平台,钢结构人才网,wordpress 微博,免费网络推广软件使用python设计滤波器 文章目录 使用python设计滤波器完整滤波器设计代码#xff08;未经完整验证#xff0c;博主还在不断完善中#xff09;关键原理与代码对应说明1. 滤波器类型选择2. 阶数估算原理3. 性能分析技术4. 设计参数调整指南 习惯了python后#xff0c;matlab逐…使用python设计滤波器 文章目录 使用python设计滤波器完整滤波器设计代码未经完整验证博主还在不断完善中关键原理与代码对应说明1. 滤波器类型选择2. 阶数估算原理3. 性能分析技术4. 设计参数调整指南 习惯了python后matlab逐渐成为了牛夫人尤其是漂亮国对龙国制裁后我作为有骨气其的码农虽然有和谐版的matlab可用但是终究是放不下码农的尊严不到万不得一已经很少用matlab了。绝大部分仿真工作都已经移步到python但最近要进行滤波器设计为了不使用matlab的fdatool便尝试着用python设计滤波器。
闲话少说代码说话 完整滤波器设计代码未经完整验证博主还在不断完善中
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from matplotlib.ticker import EngFormatter# 设计参数
fs 10000 # 采样率 10kHz
nyq fs / 2 # 奈奎斯特频率
passband 1500 # 通带截止频率 (Hz)
stopband 2000 # 阻带截止频率 (Hz)
pass_ripple 1 # 通带波动 (dB)
stop_atten 40 # 阻带衰减 (dB)# 计算归一化频率
wp passband / nyq
ws stopband / nyqdef design_iir_filter():设计椭圆IIR滤波器# 计算最小阶数和自然频率order, wn signal.ellipord(wp, ws, pass_ripple, stop_atten)# 设计椭圆滤波器b, a signal.ellip(order, pass_ripple, stop_atten, wn, btypelow, analogFalse)return b, a, orderdef design_fir_filter():设计FIR滤波器凯塞窗方法# 计算过渡带宽度width abs(stopband - passband) / nyq# 估算凯塞窗参数A stop_attenbeta 0.1102*(A - 8.7) if A 50 else 0.5842*(A - 21)**0.4 0.07886*(A - 21)# 计算所需阶数numtaps int((A - 8) / (2.285 * 2 * np.pi * width)) 1# 设计FIR滤波器taps signal.firwin(numtaps, wnpassband/nyq, window(kaiser, beta), pass_zerolowpass)return taps, numtapsdef analyze_filter(b, aNone):分析滤波器性能fig, (ax_mag, ax_phase, ax_grp, ax_zp) plt.subplots(4, 1, figsize(10, 12))# 幅频响应w, h signal.freqz(b, a, worN8000, fsfs)ax_mag.plot(w, 20 * np.log10(np.abs(h)), colorblue)ax_mag.set_title(幅频响应)ax_mag.set_ylabel(幅度 (dB))ax_mag.grid(True, whichboth, linestyle--)ax_mag.axvline(passband, colorg, linestyle--, alpha0.7)ax_mag.axvline(stopband, colorr, linestyle--, alpha0.7)ax_mag.set_ylim([-stop_atten-20, 5])# 相频响应angles np.unwrap(np.angle(h))ax_phase.plot(w, angles, colorgreen)ax_phase.set_title(相频响应)ax_phase.set_ylabel(相位 (弧度))ax_phase.grid(True)# 群延迟grp_delay -np.diff(angles) / np.diff(w)ax_grp.plot(w[:-1], grp_delay, colorred)ax_grp.set_title(群延迟)ax_grp.set_ylabel(采样点数)ax_grp.grid(True)# 零极点图if a is not None: # IIR滤波器z, p, k signal.tf2zpk(b, a)ax_zp.scatter(np.real(z), np.imag(z), markero, facecolorsnone, edgecolorsb, s80)else: # FIR滤波器z np.roots(b)p np.zeros(len(z))ax_zp.scatter(np.real(z), np.imag(z), markero, facecolorsnone, edgecolorsb, s80)ax_zp.scatter(np.real(p), np.imag(p), markerx, colorr, s80)unit_circle plt.Circle((0,0), radius1, fillFalse, colorgray, linestyle--)ax_zp.add_patch(unit_circle)ax_zp.set_title(零极点图)ax_zp.set_xlabel(实部)ax_zp.set_ylabel(虚部)ax_zp.grid(True)ax_zp.axis(equal)ax_zp.axhline(0, colorblack, linewidth0.5)ax_zp.axvline(0, colorblack, linewidth0.5)plt.tight_layout()return fig# 主程序
if __name__ __main__:# 设计IIR滤波器b_iir, a_iir, iir_order design_iir_filter()print(fIIR滤波器阶数: {iir_order})# 设计FIR滤波器fir_taps, fir_order design_fir_filter()print(fFIR滤波器阶数: {fir_order})# 分析IIR滤波器fig_iir analyze_filter(b_iir, a_iir)fig_iir.suptitle(椭圆IIR滤波器分析, fontsize16)# 分析FIR滤波器fig_fir analyze_filter(fir_taps)fig_fir.suptitle(FIR滤波器分析, fontsize16)# 应用滤波器示例t np.linspace(0, 1, fs) # 1秒时长sig np.sin(2*np.pi*1000*t) 0.5*np.sin(2*np.pi*3000*t)# IIR滤波filtered_iir signal.lfilter(b_iir, a_iir, sig)# FIR滤波filtered_fir signal.lfilter(fir_taps, 1.0, sig)# 绘制结果plt.figure(figsize(10, 6))plt.plot(t[:500], sig[:500], b-, label原始信号)plt.plot(t[:500], filtered_iir[:500], r-, labelIIR滤波后)plt.plot(t[:500], filtered_fir[:500], g-, labelFIR滤波后)plt.title(时域滤波效果对比 (前500点))plt.xlabel(时间 (秒))plt.ylabel(幅度)plt.legend()plt.grid(True)plt.tight_layout()plt.show()关键原理与代码对应说明
1. 滤波器类型选择 IIR滤波器椭圆滤波器Elliptic提供最陡峭过渡带 order, wn signal.ellipord(wp, ws, pass_ripple, stop_atten)
b, a signal.ellip(order, pass_ripple, stop_atten, wn, btypelow)ellipord计算最小阶数和自然频率椭圆滤波器在通带/阻带均有等波纹特性 FIR滤波器凯塞窗Kaiser提供参数化控制 beta 0.1102*(A - 8.7) # 窗形参数计算
taps signal.firwin(numtaps, cutoff, window(kaiser, beta))凯塞窗通过beta控制主瓣宽度和旁瓣衰减平衡
2. 阶数估算原理 IIR阶数估算 N K ( ω s ) K ( 1 − ω p 2 ) K ( ω p ) K ( 1 − ω s 2 ) N \frac{K(\omega_s)K(\sqrt{1-\omega_p^2})}{K(\omega_p)K(\sqrt{1-\omega_s^2})} NK(ωp)K(1−ωs2 )K(ωs)K(1−ωp2 ) 其中K为第一类完全椭圆积分 FIR阶数估算 N ≈ A − 8 2.285 ⋅ Δ ω N \approx \frac{A - 8}{2.285 \cdot \Delta\omega} N≈2.285⋅ΔωA−8 Δ ω \Delta\omega Δω为过渡带宽度A为阻带衰减(dB)
3. 性能分析技术 幅频响应signal.freqz计算复数频率响应 w, h signal.freqz(b, a, worN8000, fsfs)
20*np.log10(np.abs(h)) # 转换为dB相位特性 angles np.unwrap(np.angle(h)) # 解卷绕相位群延迟相位导数计算 -np.diff(angles)/np.diff(w) # 群延迟 -dφ/dω稳定性分析 IIR极点是否在单位圆内FIR恒稳定全零点系统
4. 设计参数调整指南 过渡带陡度 IIR增加阶数FIR增加窗长度 阻带衰减 IIR增加stop_atten参数FIR增大凯塞窗的beta值 计算效率 IIR阶数低但非线性相位FIR高阶但线性相位 此方案应该能替代MATLAB FDATool的核心功能提供从设计到分析的完整工作流。可根据具体应用调整滤波器类型低通/高通/带通和设计参数。但是比起FDATool的可视化设计还是有明显的不足需要不断完善。 研究学习不易点赞易。 工作生活不易收藏易点收藏不迷茫
