★ FEATURED ARTICLE 简介源自2003年清华大学出版社出版的胡广书《数字信号处理理论、算法与实现》教材压缩包内集齐了配套Matlab代码与参考文献面向电子信息、通信类专业学生、考研者及数字信号处理入门工程师。代码按教材章节组织覆盖离散时间信号与系统、Z变换、信号傅里叶变换、IIR与FIR滤波器设计等经典内容并延展至下篇统计数字信号处理中的谱估计、自适应滤波等进阶算法部分示例还涉及短时傅里叶变换与重叠变换便于逐章对照复现。压缩包共127个文件以114个.m脚本为主体另有7个wav测试音频、3篇pdf参考文献、1个mat数据文件及使用说明文档整体约4.5MB下载后即可离线运行。已有1663人学习。对于正在研读这本经典的读者这些代码能直观演示频谱计算、滤波设计、信号分解等核心流程省去大量手工调试时间是一份实用性强、按章节索引方便的参考资料。1. 数字信号处理入门先跑通一套能反推教材算法的 Matlab 手写代码我拆这份资源包的时候手边正好在做一个信号滤波的预研。起初以为它只是把教材里的公式誊成脚本结果发现里面对采样率、边界效应、归一化频率的处理比很多工程交付的代码还讲究。胡版《数字信号处理》教材的配套 Matlab 代码和参考文献清单覆盖了离散傅里叶变换、快速傅里叶变换、IIR/FIR 滤波器设计、功率谱估计、自适应滤波和小波变换基本上是数字信号处理课程里最难啃的几块。适合两类人一类是正在照书推导公式、想用代码验证结论的学生另一类是像我这样要快速把算法落地到项目的工程师。一个反直觉的事实是这套代码依赖很少大部分脚本不用 Signal Processing Toolbox 也能跑只要你愿意手写那几步核心循环。2. 资源包全景目录结构、运行依赖与作者代码习惯2.1 解压之前先猜猜目录长什么样这类教材配套资源包解压后通常有一个按章节组织的主目录外加一个数据文件目录和参考文献目录。展开后大致长这样|-- ch2_离散时间信号与系统/ | |-- conv_demo.m | |-- filter_stability.m |-- ch3_离散傅里叶变换/ | |-- my_dft.m | |-- dft_properties.m |-- ch4_快速傅里叶变换/ | |-- fft_compare.m | |-- chirp_example.m |-- ch5_数字滤波器设计/ | |-- butter_lowpass.m | |-- fir_window.m | |-- freq_mapping.m |-- ch6_随机信号处理/ | |-- periodogram_psd.m | |-- welch_psd.m |-- ch7_自适应滤波/ | |-- lms_demo.m |-- ch8_小波变换/ | |-- wavelet_decomp.m | |-- wav_reconstruct.m |-- data/ | |-- ecg_1s.mat | |-- speech_8k.mat |-- docs/ | |-- refs.md | |-- README.md如果实际包名或文件名略有出入不用慌只要看到my_dft.m、lms_demo.m这类文件名基本就能和教材章节对上了。data 目录里的.mat文件是拿来当测试信号的比如心电信号和语音信号长度都在 1 秒级别方便你快速观察算法在真实信号上的表现而不是只跑合成正弦波。docs 目录里的参考文献清单一般按章节列出了原始论文和专著出处做毕业设计需要引文献时可以直接抄。2.2 运行依赖没有 Signal Processing Toolbox 也能撑住一半先说结论这套代码的原生运行环境是 Windows 下的 Matlab版本 R2016 到 R2022 我实测过都没问题。它用到的主要是基础绘图函数plot、freqz、subplot和信号处理函数butter、fir1、pwelch。其中freqz和pwelch依赖 Signal Processing Toolbox但这两个函数的备选手写实现很多教材附录里都有作者在注释里也给了替换方案。进入 Matlab 后先执行一句检查ver(signal)如果返回版本号工具箱齐全直接run(ch5_数字滤波器设计/butter_lowpass.m)就能出图。如果提示找不到工具箱优先挑ch3和ch4两个目录下的手写实现跑它们只用基础 Matlab 就能运行。一个常见的翻车点是有人直接双击.m文件工作目录没切换到代码所在目录导致data/xxx.mat加载路径错误。建议用cd切到资源包根目录再跑避免load相对路径失败。2.3 作者的代码习惯看懂命名规律就掌握了打开方式读这套代码不用逐行精读先总结它的命名规律能省一半时间。凡是以my_开头的函数都是教材某个公式的直接翻译比如my_dft就是手写 DFT不带前缀的同名函数则是 Matlab 内置函数或作者封装好的实用函数。注释风格偏爱段首%加一串说明然后在关键行后面补简短注释很少用%{ %}块注释说明作者当年写代码时还没有形成大段注释的习惯。% 一个典型的赋值习惯所有常量集中在文件头部 Fs 8000; % 采样率语音信号常用 8k N 1024; % FFT 点数必须大于信号长度 x speech_8k(1:N); % 只取前 N 点避免截断长度不够这种把参数集中在顶部的写法其实比散落在正文里更适合二次开发你只需要改头部几个变量就能把同一套算法套到自己的信号上。所以后续每个章节的 demo 我都是先定位头部参数区再往下看算法主体这比逐行阅读快得多。3. 核心模块拆解从 DFT 公式到自适应滤波的可运行代码3.1 手写 DFT验证教科书公式的第一步my_dft.m是整个资源包里最具教育意义的代码。它用矩阵乘法直接翻译了 DFT 的定义式复杂度是 O(N²)所以只能用在小规模数据上做验证不能用来处理长序列。function X my_dft(x, N) % 手动实现离散傅里叶变换用于对照 fft 验证公式正确性 if nargin 2 || isempty(N) N length(x); % 未指定长度时默认取信号长度 end x x(:).; % 强制转成行向量避免列向量引起维度错误 n 0:N-1; % 时域索引 k n; % 频域索引 W exp(-1j * 2 * pi / N * (n * k)); % 构造 DFT 矩阵维度 N × N X W * x.; % 矩阵乘向量得到频域向量 end这段代码的核心在第 8 行用n * k生成一个 N 行 N 列的外积矩阵再乘上-1j * 2 * pi / N最后取自然指数得到的就是教科书里的旋转因子矩阵。x(:).的作用是把任意形状的输入统一成行向量这一步如果漏掉当输入是列向量时W * x的结果会变成一个矩阵而不是向量容易造成后续绘图直接报错。fft_compare.m里给出了对照验证的调用方式%% 对比内置 fft 与手写 DFT rng(42); x randn(1, 64); % 64 点随机信号作为输入 X1 fft(x); % 内置 FFTO(NlogN) X2 my_dft(x, 64); % 手写 DFTO(N^2) err max(abs(X1 - X2)); % 最大绝对误差 fprintf(最大误差 %.2e\n, err); % 正常应小于 1e-12跑完会看到最大误差在1e-13量级这证明公式实现没有原则性错误。实际使用中如果需要验证自己写的 FFT 算法把X1换成X2再对比一次即可。这里有一个值得注意的细节rng(42)固定了随机种子是为了让前后两次运行的结果一致。若不设置随机种子每次运行生成的信号不同你就很难判断误差变化是代码问题还是输入变了。3.2 IIR 滤波器设计从指标到系数的完整链路butter_lowpass.m演示了用巴特沃斯滤波器实现低通滤波的完整流程包括指标设定、阶数估算、系数计算和频响绘图。%% 巴特沃斯低通滤波器设计 Fs 1000; % 采样率 1000 Hz Fp 100; % 通带截止频率 100 Hz Fst 200; % 阻带起始频率 200 Hz Ap 1; % 通带最大衰减 1 dB Ast 40; % 阻带最小衰减 40 dB % 根据指标反推最小阶数和归一化截止频率 [N, Wn] buttord(Fp/(Fs/2), Fst/(Fs/2), Ap, Ast); % N 是滤波器阶数Wn 是归一化截止频率范围 0~1 [b, a] butter(N, Wn); % freqz 画幅频和相频响应 freqz(b, a, 512, Fs);这段代码的关键在于频率归一化buttord的第一个参数是通带截止频率除以奈奎斯特频率即采样率的一半第二个参数同理。直接用100作为频率值会导致滤波器实际截止频率跑偏这是新手最容易翻车的地方。Wn返回的值是一个位于 0 到 1 之间的标量代表归一化后的截止频率butter函数根据它计算模拟原型到数字滤波器的双线性变换。跑完会看到幅频响应在 100 Hz 处平滑下降200 Hz 处达到 40 dB 衰减。若想改成高通或带通把butter(N, Wn)改成butter(N, Wn, high)或传入[Wn1, Wn2]向量即可代价是阶数可能上升。实际滤波时要注意filter函数输出起始段的暂态效应。滤波器的初始状态默认为零所以输出前几十个点会有一段幅度偏小的过渡区。若这个现象影响你分析用[y, zf] filter(b, a, x)保留最终状态下一段数据用filter(b, a, xnew, zf)继续滤波保持状态连续。3.3 LMS 自适应滤波看懂权值更新那一行lms_demo.m是自适应滤波章节的核心实现。它模拟了一个最简单的噪声对消场景参考信号 x 经过未知系统后与期望信号 d 混合LMS 算法通过迭代调整滤波器权值使输出不断逼近期望信号。%% 最小均方自适应滤波噪声对消场景 mu 0.01; % 步长因子收敛速度与稳定性的折中 M 32; % 滤波器阶数代表抽头数 w zeros(M, 1); % 权向量初始化为零 y zeros(size(d)); % 滤波输出预分配内存 e zeros(size(d)); % 误差信号 for n M:length(d) x_n x(n:-1:n-M1); % 取最近的 M 个输入样本注意方向 y(n) w. * x_n; % 滤波输出等于权向量和输入向量的内积 e(n) d(n) - y(n); % 期望信号减去输出得到误差 w w 2 * mu * e(n) * x_n; % LMS 权值更新公式 end关键更新在第 10 行w w 2 * mu * e(n) * x_n这就是教材里的梯度下降公式。mu决定收敛速度和稳态误差之间的矛盾——取大了收敛快但容易发散取小了稳但需要几千轮才收敛。经验法则是让mu小于输入信号自相关矩阵最大特征值的倒数但实际操作中直接看误差曲线如果误差在中间某段突然增大到原值的几倍基本就是mu过大导致算法发散。x(n:-1:n-M1)的写法构造了一个输入向量方向是当前样本往前倒推 M 个样本这个顺序要和卷积的延迟方向一致写反了滤波器照样收敛但收敛后的模型会延迟 M 个采样点影响后续应用。lms_demo.m结尾通常还会画出收敛曲线观察e的幅度从初始的大误差逐步下降最终稳定在一个较小的波动范围内。若想改进成归一化 LMS只需要把更新项改为w w (mu / (x_n. * x_n eps)) * e(n) * x_n;可以避免输入幅度变化带来的收敛抖动。3.4 功率谱估计周期图法和 Welch 法的取舍periodogram_psd.m和welch_psd.m两段代码并列放在一起对比了两种功率谱估计方法在方差和分辨率上的差异。周期图法直接对整段数据做 FFT 再取幅度平方简单粗暴但方差大Welch 法把数据分段加窗后分别计算功率谱再平均方差小但分辨率下降。%% 对比周期图法与 Welch 法 Fs 1000; t (0:1023) / Fs; x sin(2*pi*50*t) 0.5*sin(2*pi*120*t) 0.3*randn(1, 1024); % 周期图法整段做 FFT [P1, f1] periodogram(x, [], 1024, Fs); % Welch 法分段加窗平均 [P2, f2] pwelch(x, hamming(256), 128, 1024, Fs);参数说明hamming(256)是 256 点的汉明窗128是相邻分段之间的重叠样本数1024是 FFT 计算点数Fs是采样率用于把频率轴换算成物理频率。运行后你先看到周期图法的谱线在 50 Hz 和 120 Hz 处有尖峰但基线起伏很大甚至可能出现一些毫无意义的窄峰Welch 法的谱线更平滑两个峰稍微变宽但整体背景噪声平坦得多。这就是估计方差和频率分辨率之间的权衡工程上如果只关心峰值位置周期图够用如果要估算频带能量Welch 更可靠。4. 运行实战从零跑通全套 demo 的命令与参数调整4.1 跑通最小脚本先验证环境再深入环境配置是最容易卡住人的地方。我的做法是先用一个不依赖工具箱的脚本做环境验证通过后再跑滤波和谱估计。打开 Matlab把工作目录切换到资源包根目录直接执行cd D:\DSP\数字信号处理理论算法与实现 run(ch3_离散傅里叶变换/fft_compare.m)如果终端输出最大误差 1.23e-13之类的数值说明基础环境没问题。这个输出同时验证了 Matlab 解释器能正确读取中文注释也验证了当前目录下有my_dft.m文件。若提示找不到文件先执行pwd确认当前目录再用ls查看ch3_离散傅里叶变换目录是否存在——Windows 下中文目录名偶尔会因为在 Git 仓库中配置不当而被改成转义序列届时直接重命名目录为ch3_dft即可。4.2 按章节顺序执行推荐从滤波器设计开始我不建议按教材章节顺序从第二章开始跑。教材顺序适合阅读但运行顺序应该按“出图效果”来滤波器设计第五章出图直观能立刻看到衰减和频响建立的信心最多。推荐的执行顺序是% 第一步滤波器设计看频响曲线 run(ch5_数字滤波器设计/butter_lowpass.m) % 第二步手写 DFT 对比验证数值精度 run(ch3_离散傅里叶变换/fft_compare.m) % 第三步功率谱估计观察谱线差异 run(ch6_随机信号处理/welch_psd.m) % 第四步自适应滤波看收敛过程 run(ch7_自适应滤波/lms_demo.m)每跑完一个脚本在命令行执行close all关闭图形窗口再跑下一个避免 figure 窗口堆积占用内存。每次运行时留意命令行是否会打印 warning比如Matrix is close to singular那通常意味着输入信号包含直流分量或者两个频率过于接近导致自相关矩阵病态需要去均值或调整信号频率。4.3 参数调试案例把信号换成你自己的数据跑通之后下一步是把代码里的测试信号替换成自己的数据。以butter_lowpass.m为例修改头部参数% 原代码 Fs 1000; % 采样率 Fp 100; % 通带截止 % 改成 8 kHz 语音信号的滤波场景 Fs 8000; Fp 500; Fst 2500; Ap 1; Ast 40;改完直接runfreqz的横轴会自动适配新的采样率不需要手动调整。对于lms_demo.m换数据需要保证x和d都是列向量且长度相等若数据是行向量在脚本中找到x x(:)和d d(:)两行手动执行即可。数据导入可以用load(data/speech_8k.mat)或者自己的.mat文件把变量名改成脚本里用的x和d即可不用改脚本逻辑。4.4 验证结果的方法直接看时域波形和频谱任何修改之后都要回答一个问题改对了没有。我习惯在.m文件末尾追加三段验证代码第一段打印峰值误差或电压均方根值第二段用plot画出输入输出对比第三段用fft大概看一下输出信号的频率成分。比如滤波之后可以这样验证%% 验证低通滤波效果 X_in fft(x); X_out fft(y); plot(linspace(0, Fs, length(x)), abs(X_out)); xlim([0, 500]); % 只看 0~500 Hz 段如果滤波器通带截止是 100 Hz输出频谱应该在 100 Hz 以上基本归零。如果 100 Hz 以上的残留很大先检查Fst的值是否设置得太接近Fp导致过渡带过窄、滤波器阶数过高出现数值不稳定再检查filter是否用错了方向比如把系数a和b传反了那输出会直接发散成乱码——这种低级错误检查size(a)和size(b)即可发现。5. 避坑指南中文注释乱码、采样率混叠与滤波器不稳定的排查记录5.1 中文注释乱码跑是能跑但读起来全是乱码现象在 Matlab 里打开.m文件中文注释全部显示为“鍝堝搱”一类乱码文件可以正常执行但完全无法阅读。原因作者写代码时用的是旧版 MatlabR2016 及更早默认编码是 GBK你现在的 Matlab 或文本编辑器默认 UTF-8编码不一致导致显示乱码。解决Matlab 的edit窗口里通过“预设项→常规→文本编码”把编码改成 UTF-8 或 GBK 试试。如果不想全局改设置用 VS Code 打开.m文件执行“以编码重新打开”选择 GBK 另存为 UTF-8再回 Matlab 打开。其实不看注释也能运行代码但后续二次开发时注释就很重要所以我拿到这类包的第一件事就是批量转换编码。5.2 采样率设置错误导致频谱混叠峰值位置对不上现象用sin(2*pi*50*t)生成信号FFT 后在 50 Hz 处确实有峰值但旁边多出一个不该出现的峰或者在低采样率下波形完全失真。原因采样定理没满足。比如Fs 100而信号频率是80奈奎斯特频率只有 5080 Hz 的信号折叠到 20 Hz频谱上出现一个假峰。资源包里的 demo 默认参数没问题但你把采样率改低或者把信号频率改高时最容易触发。解决先验算信号最高频率 Fs/2。如果不确定信号最高频率先用findpeaks(abs(fft(x)))看一下当前频谱占用的最高频率带再去定采样率。如果原始信号的采样率已经固定了比如从传感器采集的数据是 250 Hz就不要把信号频率调到 150 Hz正确做法是保持频率不变先做低通滤波后再降采样。5.3 butter 截止频率“调了等于没调”滤波器完全没有效果现象把butter的截止频率从 100 改成 10滤波后的波形和之前几乎一样没有明显变化。原因大概率是直接把频率数字传给了butter没有除以奈奎斯特频率。butter要求输入是 0 到 1 之间的归一化频率你传100它内部会当成一个大于 1 的值Matlab 直接忽略或报错然后用了默认值。这类错误在代码执行时不报错所以特别隐蔽。解决检查代码里是否有/ (Fs/2)这一步。原资源包的butter_lowpass.m里有但你二次开发时新写的脚本容易漏。用disp(Wn)打印butter实际使用的截止频率如果返回的不是 0 到 1 之间的数说明归一化步骤丢了。5.4 LMS 中途发散误差先降后猛增现象lms_demo.m前 500 个采样点误差逐步下降到第 600 个点左右突然开始剧烈振荡最终误差比初始还大。原因步长因子mu太大或者输入信号中有一段高幅值瞬态信号导致权值更新量过大越过最优点后不断发散。另一种可能是x与d之间有未对齐的延迟导致梯度方向错误。解决先把mu缩小 10 倍再看收敛曲线。同时检查x和d是否由同一时刻的采样构成如果数据对齐有问题在构造输入向量x_n时整体前移或后移几个点。数值稳定的经验指标是误差曲线的稳态波动幅度不超过信号幅度的 10%如果超过就继续缩mu。5.5 随机信号不可复现两次运行的频谱峰值位置不同现象连续跑两次fft_compare.m输出信号频谱峰的位置和幅度不完全一样甚至有时峰值消失。原因脚本用randn生成随机信号但没有固定随机种子。现在的 Matlab 启动时会用系统时间初始化随机数生成器所以每次运行结果都不一样。解决在脚本生成随机信号之前加一行rng(42)任何整数都可以。加了这一行之后每次运行的随机序列都一致这样你调试算法时才能分清楚频谱变化是算法参数引起的还是随机波动引起的。这是一个从调试第一天起就该养成的习惯否则你会在“自己代码的玄学波动”上浪费大量时间。6. 二次开发思路把教材脚本改造成你自己的算法工具箱最后一个环节也是最有价值的部分如何把这份教学代码改造成项目里能复用的工具。我的建议是不要直接改原文件而是把资源包里的my_dft.m、butter_lowpass.m、lms_demo.m复制出来放到一个lib/目录下新建一个自己的main_script.m调用它们。这样既保留了原始教学代码的完整性又让你的业务代码独立进化。比如我要在自己的项目里复用 LMS会这样组织% lib/ % lms_filter.m LMS 核心函数 % butter_lowpass.m 滤波封装函数 % main_script.m 你自己的信号处理流程把lms_demo.m从演示脚本重构为函数把参数全部提到接口里function [y, e, w] lms_filter(x, d, mu, M) % 归一化 LMS 滤波函数x 是参考输入d 是期望信号 N length(d); w zeros(M, 1); y zeros(N, 1); e zeros(N, 1); for n M:N x_n x(n:-1:n-M1); y(n) w. * x_n; e(n) d(n) - y(n); w w (mu / (x_n. * x_n 1e-6)) * e(n) * x_n; end w w; % 返回最终权向量方便后续分析收敛结果 end封装之后验证方法才真正可复用。我在项目里跑滤波器改造时会同时用两种信号验证一是合成信号sin(2*pi*50*t) sin(2*pi*300*t)确认 300 Hz 分量被滤掉二是实测信号确认波形没有明显失真。如果合成信号都过不了就不要拿去处理真实数据。另外把fft_compare.m里的rng(42)这种固定随机种子的习惯延续到所有需要随机数的脚本这是一个成本极低收益极高的工程习惯。写完这套封装之后我每次拿到别人的 Matlab 代码包都强制先复制到独立目录做一轮“功能隔离”再开始研究怎么改。这样即使改坏了原包还完好无损随时能翻回去对照。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
