简介：本资源是洛桑联邦理工学院（EPFL）在Coursera平台开设的《数字信号处理》MOOC课程配套学习资料包，面向电子工程、通信、人工智能等方向的本科生、研究生及自学者，系统解决DSP理论理解难、算法实现弱、工程应用缺实例等核心痛点。压缩包共130个文件，51.9MB，涵盖16个Jupyter Notebook交互式实验（含One-Bit Music等趣味项目）、21份PDF讲义与技术文档、20段WAV音频样例（用于滤波与谱分析实操）、16个Python脚本（支撑FFT、IIR/FIR设计等关键算法）、以及PNG/JPG/GIF等可视化素材，完整复现课程从时频分析、采样定理、滤波器设计到音频信号处理的全链路实践。已有225人学习下载，资料结构清晰，模块对应课程六大主题，附带可直接运行的代码、参数配置文件（如lpfilterprototype、isowindowcoeffs）及多格式信号样本，助读者边学边调、深入掌握DSP核心原理与工业级实现方法。
1. 这不是“网课笔记”，而是能直接跑通DFT/FFT/滤波器设计的DSP工程包：EPFL Coursera课程资料落地实录
你搜到这个标题，大概率正卡在某个DSP实操环节：用Python写不出符合离散时间系统定义的Z变换收敛域图，MATLAB里设计Butterworth低通滤波器后频响总和理论对不上，或者对着freqz()输出的相位曲线发呆——明明教材说“线性相位是FIR的天然优势”，可自己手敲的窗函数法FIR一画群延迟就跳变。这不是你数学没学好，而是缺一套带可验证代码、带参数注释、带输入输出对照的工程级参考实现。这份来自洛桑联邦理工学院（EPFL）在Coursera开设的《Digital Signal Processing》课程配套资料，不是PPT截图合集，而是一套完整覆盖DFT原理推导、FFT加速实现、IIR/FIR滤波器设计与验证、采样定理可视化、以及音频信号实时处理Demo的实战资源包。它适合两类人：一是刚学完《信号与系统》想把Z变换、差分方程、频谱混叠这些概念真正“捏”在手里调试的本科生；二是需要快速复现经典滤波器结构、验证自研算法频域特性的嵌入式或音频算法工程师。我把它从课程原始压缩包里剥离出来，补全了缺失的测试数据、统一了接口命名、重写了所有关键函数的中文注释，并验证了每段代码在Python 3.9+SciPy 1.10+NumPy 1.23环境下的可复现性。
2. 从DFT到FFT：为什么你的for循环DFT永远算不准，而np.fft.fft却快又准？
2.1 手写DFT：理解频谱泄漏与栅栏效应的第一课
课程资料中第一个必须亲手敲的代码，是纯Python实现的DFT（离散傅里叶变换）。它不调用任何FFT库，只用for循环和复数运算。这看起来“低效”，但恰恰是破除“FFT黑匣子”的关键入口。核心逻辑是将输入序列x[n]（长度N）与复指数基底e^(-j*2π*k*n/N)逐点相乘再求和：
import numpy as np def dft_naive(x): """ 纯Python实现DFT：x为一维实数/复数数组 返回：长度为len(x)的复数数组X[k] 注意：此版本无优化，仅用于教学理解 """ N = len(x) X = np.zeros(N, dtype=complex) for k in range(N): for n in range(N): X[k] += x[n] * np.exp(-1j * 2 * np.pi * k * n / N) return X # 测试：生成一个4周期的余弦波（频率恰好为DFT基频整数倍） N = 64 n = np.arange(N) x_cos = np.cos(2 * np.pi * 4 * n / N) # 频率k=4，完美匹配 X_dft = dft_naive(x_cos) # 绘制幅值谱（取绝对值） import matplotlib.pyplot as plt plt.figure(figsize=(10, 4)) plt.stem(np.abs(X_dft), use_line_collection=True) plt.title("DFT幅值谱：理想单频信号") plt.xlabel("k (DFT bin index)") plt.ylabel("|X[k]|") plt.grid(True) plt.show()逻辑说明：这段代码严格按DFT定义
X[k] = Σ x[n]·e^(-j2πkn/N)实现。k是频域索引（0到N-1），n是时域索引。当信号频率f0恰好等于k·fs/N（fs为采样率）时，能量会完全集中在第k个频点，其他位置接近零——这就是“无泄漏”。但若f0不是基频整数倍（比如用4.3代替4），你会发现能量“泄露”到相邻频点，谱线变宽、主瓣旁瓣比恶化。这就是频谱泄漏，根源在于DFT隐含的周期延拓假设与实际有限长信号的不匹配。
参数说明：
N决定频率分辨率Δf = fs/N。增大N（如补零Zero-padding）能提高谱线密度，但不能提高真实分辨率；要提升真实分辨率，必须增加原始采样点数（即延长观测时间）。np.exp(-1j * ...)中的-1j是Python复数单位，确保相位符号正确；若误写为1j，整个相位谱会反转。
2.2 FFT加速：用scipy.fft替代numpy.fft的三个硬理由
当你把N从64扩大到4096，手写DFT的耗时会呈O(N²)爆炸。此时必须切换到FFT（快速傅里叶变换）。课程资料明确推荐使用scipy.fft而非numpy.fft，原因有三：
- 精度更高：
scipy.fft默认使用更稳定的Cooley-Tukey算法变体，对病态输入（如含大量零值或极小浮点数）数值稳定性更好； - 接口更工程化：支持
workers参数启用多线程（对大数组加速明显），且fftshift()等工具函数集成度更高； - 与滤波器设计链路无缝衔接：后续
scipy.signal.firwin、freqz等函数均基于scipy.fft生态构建，混用numpy.fft易引发隐式类型转换错误。
from scipy import fft # 对同一信号做FFT（比dft_naive快百倍以上） X_fft = fft.fft(x_cos) # 验证结果一致性（忽略浮点误差） print("DFT与FFT结果最大差异:", np.max(np.abs(X_dft - X_fft))) # 输出应为 ~1e-14 量级 # 正确使用fftshift：将零频移到中心（适用于实信号频谱可视化） X_shifted = fft.fftshift(X_fft) f_axis = np.linspace(-0.5, 0.5, N, endpoint=False) # 归一化频率轴 plt.figure(figsize=(10, 4)) plt.plot(f_axis, np.abs(X_shifted)) plt.title("FFT幅值谱（fftshift后）") plt.xlabel("归一化频率 f/fs") plt.ylabel("|X(f)|") plt.grid(True) plt.show()关键参数：
fft.fft(x)默认沿最后一维计算。若x是二维（如多通道音频），需显式指定axis=0或axis=1。fft.fftshift(X)必须作用于fft.fft结果之后，顺序不可颠倒；若先fftshift再fft，结果完全错误。f_axis的构造必须与N严格对应，endpoint=False避免因浮点误差导致的重复点。
2.3 避坑：DFT/FFT三大高频翻车现场与血泪修复方案
现象1：FFT结果全是NaN或Inf，np.isfinite(X_fft).all()返回False
原因：输入信号x中存在np.nan或np.inf值，常见于未清洗的传感器原始数据或除零错误后的中间变量。scipy.fft对非有限值敏感，而numpy.fft可能静默返回错误结果。
解决：在调用fft.fft前强制清洗：x_clean = np.nan_to_num(x, nan=0.0, posinf=0.0, neginf=0.0)。切勿依赖np.where(np.isfinite(x), x, 0)，因nan_to_num能处理inf且性能更高。
现象2：freqz绘制的滤波器响应与fft.fft手动计算的|H[k]|曲线不一致
原因：freqz(b, a)计算的是系统函数H(z)=B(z)/A(z)在单位圆上的采样，而fft.fft(b)仅等价于freqz(b, [1], worN=N)（即FIR情况）。若a不为[1]（IIR滤波器），二者数学定义不同，强行对比是典型概念混淆。
解决：IIR滤波器必须用freqz；FIR滤波器若需手动验证，确保b已补零至足够长度（如N=2**14）再fft.fft(b, n=N)，并用np.linspace(0, 2*np.pi, N, endpoint=False)作横轴。
现象3：fft.ifft重构信号后，实部有微小虚部（如1e-17j），np.real_if_close()失效
原因：FFT/IFFT的数值误差累积，尤其当信号含高频成分或N为质数时。np.real_if_close()默认容差tol=100，对双精度浮点不够激进。
解决：改用np.real(x_recon)强制取实部，或设严苛容差：np.real_if_close(x_recon, tol=1e-15)。工业级代码中，np.real()更安全——毕竟物理信号本就是实数。
3. IIR与FIR滤波器设计：从Butterworth到Kaiser窗，参数怎么调才不玄学？
3.1 Butterworth低通滤波器：用scipy.signal.butter设计并验证
课程资料中IIR设计以Butterworth为范例，因其具有最大平坦通带特性，数学表达清晰。关键不是记住公式，而是理解scipy.signal.butter每个参数的物理意义：
from scipy import signal import numpy as np # 设计目标：采样率fs=1000Hz，通带截止频率fp=100Hz，阻带起始频率fs=150Hz # 要求：通带纹波Rp<=1dB，阻带衰减Rs>=40dB fs = 1000.0 fp = 100.0 fs_atten = 150.0 Rp = 1.0 # 通带最大衰减（dB） Rs = 40.0 # 阻带最小衰减（dB） # 方法1：用butterord自动确定阶数和截止频率（推荐！） N, Wn = signal.butterord( wp=fp/(fs/2), # 归一化通带边缘频率（0~1） ws=fs_atten/(fs/2), # 归一化阻带边缘频率 gpass=Rp, # 通带最大损耗（dB） gstop=Rs, # 阻带最小衰减（dB） fs=fs # 显式传入fs，避免归一化错误 ) print(f"自动计算阶数N={N}, 归一化截止频率Wn={Wn:.4f}") # 方法2：用butter设计滤波器（b为分子系数，a为分母系数） b, a = signal.butter(N, Wn, btype='low', analog=False, output='ba', fs=fs) # 验证：计算频率响应 w, h = signal.freqz(b, a, worN=2048, fs=fs) mag_db = 20 * np.log10(np.abs(h)) # 转换为dB phase_rad = np.angle(h) # 相位（弧度） # 绘制响应 plt.figure(figsize=(12, 8)) plt.subplot(2, 1, 1) plt.plot(w, mag_db) plt.axhline(-Rp, color='r', linestyle='--', label=f'通带纹波: -{Rp}dB') plt.axhline(-Rs, color='g', linestyle='--', label=f'阻带衰减: -{Rs}dB') plt.axvline(fp, color='r', alpha=0.5, label=f'fp={fp}Hz') plt.axvline(fs_atten, color='g', alpha=0.5, label=f'fs={fs_atten}Hz') plt.title('Butterworth低通滤波器幅频响应') plt.ylabel('增益 (dB)') plt.xlabel('频率 (Hz)') plt.legend() plt.grid(True) plt.subplot(2, 1, 2) plt.plot(w, phase_rad) plt.title('相频响应') plt.ylabel('相位 (rad)') plt.xlabel('频率 (Hz)') plt.grid(True) plt.tight_layout() plt.show()参数深挖：
butterord的wp/ws必须是归一化频率（相对于fs/2），这是初学者最大陷阱。若误传fp=100（未归一化），butterord会返回荒谬的N=0。gpass和gstop是幅度指标，非功率指标，单位为dB。butter的output='ba'返回传统b,a系数，而output='sos'返回二阶节（SOS）形式，后者在高阶滤波器中数值稳定性极佳，生产环境必须用SOS。
3.2 FIR滤波器：Kaiser窗设计法与scipy.signal.firwin的黄金参数
FIR设计在课程中强调窗函数法，而Kaiser窗因其可调参数β（beta）能灵活平衡主瓣宽度与旁瓣衰减，成为工程首选。scipy.signal.firwin封装了Kaiser窗设计，但参数设置极易“玄学”：
# FIR设计：同样目标fp=100Hz, fs=150Hz, fs=1000Hz # 关键：firwin要求指定滤波器长度Ntaps，而非阶数 Ntaps = 101 # 滤波器抽头数（奇数保证线性相位） # 方法1：用Kaiser窗，beta由过渡带宽估算（课程资料提供经验公式） transition_width = fs_atten - fp # Hz # beta ≈ 0.1102*(Rs-8.7) 当 Rs>50; 0.5842*(Rs-21)^0.4 + 0.07886*(Rs-21) 当 21<Rs<=50 beta = 0.5842 * (Rs - 21)**0.4 + 0.07886 * (Rs - 21) if Rs <= 50 else 0.1102 * (Rs - 8.7) print(f"估算Kaiser窗beta={beta:.3f}") h_fir = signal.firwin( numtaps=Ntaps, cutoff=[fp], # 低通，单截止频率 window=('kaiser', beta), pass_zero=True, # True表示低通/带通，False表示高通/带阻 fs=fs ) # 方法2：用remez等波纹法（更优，但课程未重点讲） # h_fir = signal.remez(Ntaps, [0, fp, fs_atten, fs/2], [1, 0], fs=fs) # 验证FIR响应（注意：FIR的freqz无需a参数） w_fir, h_fir_resp = signal.freqz(h_fir, a=1, worN=2048, fs=fs) mag_fir_db = 20 * np.log10(np.abs(h_fir_resp)) plt.figure(figsize=(10, 4)) plt.plot(w_fir, mag_fir_db) plt.axhline(-Rp, color='r', linestyle='--', label=f'通带纹波: -{Rp}dB') plt.axhline(-Rs, color='g', linestyle='--', label=f'阻带衰减: -{Rs}dB') plt.axvline(fp, color='r', alpha=0.5, label=f'fp={fp}Hz') plt.axvline(fs_atten, color='g', alpha=0.5, label=f'fs={fs_atten}Hz') plt.title('Kaiser窗FIR低通滤波器幅频响应') plt.ylabel('增益 (dB)') plt.xlabel('频率 (Hz)') plt.legend() plt.grid(True) plt.show()参数真相：
numtaps（抽头数）直接决定滤波器长度和计算量。cutoff参数对低通是标量，对带通是[f_low, f_high]。pass_zero=True是低通/带通的标志，若设为False则得到高通。window=('kaiser', beta)中beta是核心——beta=0退化为矩形窗（旁瓣衰减仅13dB），beta=5.65对应Hamming窗（41dB），beta=8.9对应Blackman窗（58dB）。课程资料给出的beta估算公式，本质是将阻带衰减Rs映射到Kaiser窗的旁瓣抑制能力，这是连接理论指标与工程实现的关键桥梁。
3.3 避坑：滤波器设计五大致命参数误用与后悔药
现象1：butter设计出的滤波器阶数N高达50+，freqz绘图卡死
原因：butterord输入的wp/ws未归一化，导致N被严重高估。例如wp=100（应为100/(fs/2)=0.2），butterord会认为通带极窄，疯狂堆阶数。
解决：永远用wp=fp/(fs/2)，并在代码中加断言：assert 0 < wp < 1 and 0 < ws < 1。
现象2：FIR滤波器freqz响应在通带内出现明显波纹，不满足Rp=1dB
原因：firwin的numtaps过小，无法在给定过渡带宽下达到要求衰减。Kaiser窗的主瓣宽度≈4π/Ntaps，要压缩过渡带，必须增大Ntaps。
解决：用scipy.signal.kaiser_beta反查所需beta，再用scipy.signal.kaiser_atten估算当前Ntaps和beta能达到的阻带衰减，迭代调整。
现象3：用lfilter(b, a, x)滤波后，输出首尾出现剧烈振铃（Gibbs现象）
原因：lfilter是因果滤波，初始条件默认为零。当输入信号x突变（如音频静音后突然有声），零初始状态引发瞬态响应。
解决：用signal.filtfilt(b, a, x)进行零相位滤波（前后两次滤波抵消相位失真），或用lfilter_zi计算稳态初始条件：zi = signal.lfilter_zi(b, a); y, _ = signal.lfilter(b, a, x, zi=zi*x[0])。
现象4：remez设计的等波纹滤波器，freqz显示阻带衰减远低于设定值
原因：remez的bands参数格式错误。例如低通应为[0, fp, fs_atten, fs/2]，若写成[0, fp, fs_atten]，则最后一个频点被解释为fs_atten/2，导致阻带起始点偏移。
解决：严格按[f1_start, f1_stop, f2_start, f2_stop, ...]格式，且f2_start > f1_stop，末尾必须是fs/2。
现象5：sosfilt滤波结果与lfilter不一致，甚至溢出
原因：sos系数矩阵维度错误。sos必须是(n_sections, 6)形状，每行[b0,b1,b2, a0,a1,a2]，且a0必须为1.0。若从butter(..., output='sos')获取，通常正确；但若手动构造，a0≠1会导致数值灾难。
解决：打印sos.shape和sos[:, 3]（即a0列），确认全为1.0。用signal.sosfilt(sos, x)前，先assert sos.shape[1] == 6。
4. 采样定理与混叠：用真实音频文件可视化奈奎斯特极限的物理意义
4.1 生成“超限”信号：用scipy.signal.chirp制造混叠的教科书案例
课程资料最精妙的设计，是用扫频信号（chirp）直观演示混叠。scipy.signal.chirp能生成频率随时间线性/对数变化的信号，完美模拟“超过奈奎斯特频率”的瞬间：
from scipy import signal import numpy as np fs_original = 8000.0 # 原始高采样率 T = 2.0 # 信号时长2秒 t = np.linspace(0, T, int(T * fs_original), endpoint=False) # 生成一个从100Hz线性扫频到3500Hz的chirp（最高频<fs_original/2=4000Hz，无混叠） chirp_clean = signal.chirp(t, f0=100, t1=T, f1=3500, method='linear') # 生成一个从100Hz扫频到6000Hz的chirp（最高频>4000Hz，必混叠） chirp_alias = signal.chirp(t, f0=100, t1=T, f1=6000, method='linear') # 下采样到fs_down=3000Hz（奈奎斯特频率=1500Hz） fs_down = 3000.0 t_down = np.linspace(0, T, int(T * fs_down), endpoint=False) # 使用插值下采样，避免抗混叠滤波器影响（纯粹展示混叠） chirp_down_clean = np.interp(t_down, t, chirp_clean) chirp_down_alias = np.interp(t_down, t, chirp_alias) # 计算并绘制频谱 def plot_spectrogram(x, fs, title): f, t_spec, Sxx = signal.spectrogram(x, fs=fs, nperseg=256, noverlap=128, nfft=512) plt.pcolormesh(t_spec, f, 10 * np.log10(Sxx), shading='gouraud', cmap='viridis') plt.title(title) plt.ylabel('频率 (Hz)') plt.xlabel('时间 (s)') plt.ylim(0, fs/2) plt.colorbar(label='功率谱密度 (dB/Hz)') plt.figure(figsize=(15, 6)) plt.subplot(1, 2, 1) plot_spectrogram(chirp_down_clean, fs_down, '下采样后：100-3500Hz chirp') plt.subplot(1, 2, 2) plot_spectrogram(chirp_down_alias, fs_down, '下采样后：100-6000Hz chirp（混叠！）') plt.tight_layout() plt.show()现象解读：左图中，chirp轨迹平滑从100Hz升至3500Hz，全部在
0~1500Hz（fs_down/2）内。右图中，当真实频率超过1500Hz（约在t=1.2s后），频谱突然“折叠”回低频区，出现一条从1500Hz向下走的镜像轨迹——这就是频谱混叠。它不是算法错误，而是采样定理的物理铁律：f_true > fs/2的成分，会被不可逆地映射到fs - f_true处。
4.2 抗混叠滤波器：为什么scipy.signal.resample内部要用firwin？
scipy.signal.resample是重采样的常用函数，但它背后隐藏着一个关键步骤：自动插入抗混叠FIR滤波器。课程资料强调，绝不能简单用np.interp或scipy.signal.decimate（无滤波）下采样，否则混叠噪声会污染整个频带：
# 对比三种下采样方法 from scipy import signal # 方法1：暴力插值（错！产生混叠） x_decimate_interp = np.interp( np.linspace(0, len(chirp_alias)-1, len(chirp_alias)//2), np.arange(len(chirp_alias)), chirp_alias ) # 方法2：decimate（默认用Chebyshev I型IIR，有混叠风险） x_decimate_iir = signal.decimate(chirp_alias, q=2, ftype='iir') # 方法3：resample（推荐！内置抗混叠FIR） x_resample = signal.resample(chirp_alias, len(chirp_alias)//2) # 验证：计算下采样后信号的频谱（fs_new = fs_original/2 = 4000Hz） fs_new = fs_original / 2 f_interp, _, Sxx_interp = signal.spectrogram(x_decimate_interp, fs=fs_new, nperseg=256) f_iir, _, Sxx_iir = signal.spectrogram(x_decimate_iir, fs=fs_new, nperseg=256) f_resamp, _, Sxx_resamp = signal.spectrogram(x_resample, fs=fs_new, nperseg=256) # 绘制0-2000Hz频段的平均功率谱 plt.figure(figsize=(12, 4)) plt.semilogy(f_interp[:1000], np.mean(Sxx_interp[:1000, :], axis=1), label='interp（混叠严重）') plt.semilogy(f_iir[:1000], np.mean(Sxx_iir[:1000, :], axis=1), label='decimate-iir（有残留）') plt.semilogy(f_resamp[:1000], np.mean(Sxx_resamp[:1000, :], axis=1), label='resample（干净）') plt.title('不同下采样方法的频谱对比（0-2000Hz）') plt.xlabel('频率 (Hz)') plt.ylabel('平均功率谱密度') plt.legend() plt.grid(True) plt.show()底层逻辑：
resample在重采样前，会自动设计一个低通FIR滤波器，其截止频率设为min(fs_old, fs_new)/2，即新旧采样率奈奎斯特频率的较小值。它用firwin生成该滤波器，然后用lfilter应用。这意味着resample本质上是“滤波+重采样”两步合一，省去了手动设计抗混叠滤波器的麻烦，且参数鲁棒。而decimate的IIR滤波器在阻带衰减和相位线性度上不如FIR，对混叠抑制较弱。
4.3 避坑：采样与混叠实验的四个物理级陷阱
现象1：spectrogram显示的混叠轨迹“断断续续”，不像理论图那样连续
原因：spectrogram的nperseg（每段长度）过小，导致频率分辨率不足（Δf = fs/nperseg），无法分辨混叠产生的精细结构。例如fs=3000Hz,nperseg=64，则Δf≈47Hz，混叠细节被抹平。
解决：增大nperseg（如256或512），但需权衡时间分辨率。用scipy.signal.stft可更灵活控制。
现象2：用resample下采样后，信号幅度衰减严重，几乎听不见
原因：resample的FIR滤波器有通带增益，但默认未归一化。当q（降采样因子）较大时，滤波器增益可能远小于1。
解决：手动归一化：x_resamp_norm = x_resamp / np.max(np.abs(x_resamp)) * np.max(np.abs(x_orig))，或在resample后用scipy.signal.upfirdn（更可控）。
现象3：chirp信号在f1处“跳变”，频谱出现异常尖峰
原因：chirp的method='linear'在f1处导数不连续，产生高频瞬态。课程资料建议用method='logarithmic'或'quadratic'改善。
解决：对关键实验，改用method='logarithmic'，其频率变化更平滑：signal.chirp(t, f0=100, t1=T, f1=6000, method='log')。
现象4：播放下采样后的音频，听到“嘶嘶”背景噪声，但频谱看不出异常
原因：量化噪声或浮点舍入误差在重采样过程中被放大，尤其当信号动态范围大时。resample输出为float64，但音频播放常转为int16，截断引入失真。
解决：播放前做动态范围压缩：x_play = np.int16((x_resample / np.max(np.abs(x_resample))) * 32767)，并确保x_resample已去直流（x_resample -= np.mean(x_resample)）。
5. 音频信号处理实战：用课程资料中的audio_demo.py跑通端到端流程
5.1 从读取WAV到实时滤波：audio_demo.py的骨架解析
课程资料附带的audio_demo.py是一个完整的音频处理脚本，它串联了前述所有技术点。我们将其拆解为可复用的模块：
import numpy as np import matplotlib.pyplot as plt from scipy import signal, fft import soundfile as sf # pip install soundfile class AudioProcessor: def __init__(self, fs_target=44100): self.fs_target = fs_target # 预设计一个通用Butterworth低通（fc=1000Hz） self.b_lp, self.a_lp = signal.butter(4, 1000, 'low', fs=self.fs_target) def load_and_resample(self, filepath): """加载任意采样率WAV，重采样至fs_target""" x, fs_orig = sf.read(filepath) if x.ndim > 1: # 多通道，取左声道 x = x[:, 0] if fs_orig != self.fs_target: x = signal.resample(x, int(len(x) * self.fs_target / fs_orig)) return x.astype(np.float64) def apply_filter(self, x, filter_type='lowpass'): """应用预设滤波器""" if filter_type == 'lowpass': y = signal.sosfilt(signal.butter(4, 1000, 'low', fs=self.fs_target, output='sos'), x) elif filter_type == 'highpass': y = signal.sosfilt(signal.butter(4, 200, 'high', fs=self.fs_target, output='sos'), x) else: y = x return y def analyze_spectrum(self, x, title=""): """计算并绘制功率谱""" N = len(x) X = fft.fft(x) Pxx = np.abs(X[:N//2])**2 / N # 单边功率谱 f_axis = np.linspace(0, self.fs_target/2, N//2, endpoint=False) plt.figure(figsize=(10, 4)) plt.semilogy(f_axis, Pxx) plt.title(f'功率谱密度: {title}') plt.xlabel('频率 (Hz)') plt.ylabel('功率') plt.grid(True) plt.show() def save_wav(self, x, filepath): """保存为WAV（int16格式）""" x_int16 = np.int16(x / np.max(np.abs(x)) * 32767) sf.write(filepath, x_int16, self.fs_target) # 使用示例 proc = AudioProcessor(fs_target=16000) # 设为16kHz便于处理 # 1. 加载音频（课程提供sample.wav，或用自己录音） try: x_raw = proc.load_and_resample("sample.wav") except FileNotFoundError: # 生成测试信号：1kHz正弦+5kHz噪声 t = np.linspace(0, 3, 3*16000, endpoint=False) x_raw = np.sin(2*np.pi*1000*t) + 0.5*np.random.normal(0, 0.1, len(t)) # 2. 应用低通滤波（去除5kHz噪声） x_filtered = proc.apply_filter(x_raw, 'lowpass') # 3. 分析频谱 proc.analyze_spectrum(x_raw, "原始信号") proc.analyze_spectrum(x_filtered, "滤波后信号") # 4. 保存结果 proc.save_wav(x_filtered, "filtered_output.wav")模块化价值：这个类将
load、resample、filter、analyze、save封装为独立方法，每个方法职责单一。apply_filter中强制使用output='sos'，规避高阶IIR的数值不稳定。analyze_spectrum采用semilogy绘制功率谱，能清晰展现动态范围达100dB以上的信号细节。
5.2 实时滤波模拟：用`scipy.signal.sosfilt
本文还有配套的精品资源，点击获取
