简介围绕胡广书《数字信号处理理论、算法与实现》经典教材编写的MATLAB代码与参考文献合集面向电子信息、通信工程等专业师生及数字信号处理开发者覆盖从入门到进阶的应用场景。书籍内容涵盖离散时间信号与系统、Z变换、傅里叶变换、经典数字信号处理与统计数字信号处理等核心模块配套M文件与各章节例题、习题对应兼具理论推导与算法实现双重参考价值。压缩包共127个文件包含114个MATLAB脚本、7个WAV语音样例、3个PDF参考文献、1个MAT数据文件及使用说明文档整体约4.5MB体量精简且目录结构清晰便于按章节索引与快速运行验证。内容预览可见短时傅里叶变换、重叠变换等专题实现脚本能辅助课程学习、考研复习以及科研算法验证。目前已有1663人浏览学习适合系统研读教材并需动手实践算法的读者深入使用。1. 数字信号处理为什么公式全会一写代码就翻车很多人在学《数字信号处理理论、算法与实现》时都会撞上同一个怪圈课本公式推导能看懂一到上机写数字信号处理相关的 MATLAB 代码结果就完全不听指挥。频谱峰值不在预期位置滤波器输出开头乱成一团时移后的信号像幽灵一样绕回原点。这不是你数学没学明白而是从连续域公式到离散域代码之间隔着一层没人点破的“翻译约定”。频率怎么归一化、下标从几开始、序列怎么补零这些细节写在代码里时每个都有理论依据。这套教材配套的 MATLAB 代码价值恰恰在于把每一个公式都翻译成可运行的实现把那些课本里一带而过的转换层摊开给你看。搞懂这层翻译机制才算真正理解数字信号处理的“可计算形态”。2. 从公式到代码连续世界与离散计算之间的三层翻译2.1 归一化频率数字信号里的“π”为什么反复出现打开配套代码几乎每个频谱分析脚本都会出现这样几行fs 8000; % 采样率单位 Hz N 1024; % FFT 点数 f (0:N-1) / N * fs; % 频率轴从 0 到 fs - fs/N单位 Hz这一行里最关键的是“除以 N 再乘以 fs”。FFT 输出的是 N 个等间隔频率采样点相邻两个点的频率间隔是 fs/N这个值就是频率分辨率。(0:N-1)/N把整数 bin 序号映射成 0 到 1 之间的比例再乘 fs 变成物理频率。很多人用findpeaks找到峰值索引后直接用索引去对比频率或者忘了除 N结果频率轴整体错位。这里有个更底层的设计动机离散系统里的频率天然以采样率为基准做归一化单位是 rad/sample范围只有 [0, π]π 对应奈奎斯特频率 fs/2。为什么非要这么设计因为算法要和采样率解耦。同一组滤波器系数在 fs8000 和 fs44100 下都能工作只是归一化频率对应的物理频率不同。配套代码里fir1这类滤波器的截止频率参数 Wn 全部是归一化值Wn 截止频率 / (fs/2)。比如设计 4000 Hz 低通采样率 8000 HzWn 就是 0.5不是 4000。很多翻车现场就是把模拟频率直接塞进数字滤波器设计函数结果滤波器形状完全不对。还有一点容易被忽略真正可区分的最小频率间隔由时域窗长度 T 决定等于 1/T。补零只让频谱曲线更平滑不会提高分辨率。N决定的是频域采样密度不是分辨率的物理上限。源码里那些“先补零到 4096 点再做 FFT”的操作本质是插值显示很多人误以为补零能分辨更近的两个频率这是对频率轴含义的误读。2.2 DFT 的隐含周期性与频谱泄漏DFT 公式看起来只是对 N 个点求加权和但它隐含了一个非常强的假设这 N 点数据是某个周期信号的一个周期。也就是说把这段有限长序列周期延拓之后边界处的采样点是连续衔接的。问题在于现实里你截取的数据段往往不满足这个条件。配套代码里几乎都会有一个对照实验Fs 1000; t (0:999)/Fs; x1 sin(2*pi*50*t); % 50Hz 正弦整数周期 x2 sin(2*pi*47*t); % 47Hz 正弦非整数周期 X1 fft(x1); X2 fft(x2);x1 在一个采样秒内恰好完成 50 个完整周期首尾相切周期延拓后边界连续频谱是干净的单根谱线。x2 在 1000 个点里只完成 47 个周期数据段首尾值不相等周期延拓后在边界处形成一个阶跃。这个跳变在频域里表现为能量从真实频率处向所有 bin 扩散就是频谱泄漏。换一种更本质的说法对数据段做 DFT 等价于先对无限长信号乘一个矩形窗再做周期延拓。时域相乘对应频域卷积矩形窗的频谱是 sinc 函数真实谱线被它卷积后主瓣展宽、旁瓣出现。理解这个机制再去看代码里的加窗操作就很清晰了。Hamming、Blackman 这些窗函数的本质是让数据段边缘平滑衰减到零压低周期延拓边界处的阶跃幅度代价是主瓣变宽。主瓣宽度决定过渡带旁瓣高度决定阻带衰减这就是为什么滤波器设计里永远存在“过渡带陡峭”和“阻带衰减大”这对矛盾。你把这段推理串起来就会发现频谱泄漏不是数值 bug而是有限数据长度的数学必然。2.3 索引约定从 1 开始的数组下标如何让公式“移位”MATLAB 的数组下标从 1 开始但 DFT 的数学定义里 n 和 k 都从 0 数到 N-1。这个错位是配套代码里最常见的隐性 bug 来源。直接写exp(-2i*pi*k*(1:N)/N)看起来和公式很像但序列起点从 n0 变成了 n1整体多乘了一个 e^{-j2πk/N} 的相位因子。对幅度谱没影响因为取模把相位信息丢掉了但一旦你做相位分析或者用 ifft 重建信号重建出来的序列会比原来延迟一个点表现为整体时移。配套代码里通常这样生成时间轴n 0:N-1; % 时间索引从 0 开始符合 DFT 定义 k 0:N-1; % 频率索引从 0 开始 x cos(2*pi*5*n/N); % 正好 5 个周期的余弦 X x * exp(-2i*pi/N * (n*k)); % 用矩阵乘法直接算 DFTn*k生成的是 N×N 的外积矩阵每一行对应一个时间点 n每一列对应一个频率 bin k整个矩阵就是 DFT 的旋转因子矩阵。这种写法把 DFT 从“循环求和”变成“矩阵乘”和书上公式的对应关系一目了然X[k] 等于 x 与第 k 列旋转因子做内积。用内置fft函数验证两者结果完全一致只是内置实现利用了旋转因子的周期性和对称性把复杂度从 O(N²) 降到 O(N log N)。从这段代码你能学到一个习惯写数字信号处理的迭代代码前先把 n、k 的定义域写清楚下标从 0 开始生成序列再考虑索引的偏移。这个习惯能避免大量相位类的“玄学”问题。3. 源码级拆解卷积、相关和 FFT 在代码里到底做了什么3.1 线性卷积为什么要补零到 LM-1卷积是线性时不变系统的核心运算但 DFT 的卷积定理给的是循环卷积不是线性卷积。这两个东西在什么条件下等价是配套代码里 FIR 滤波实现的关键前提。一个长度为 L 的信号 x 和长度为 M 的滤波器 h线性卷积结果的长度是 LM-1。如果直接把 x 和 h 做 N 点 FFT、频域相乘、再 IFFT得到的是 N 点循环卷积。当 N 小于 LM-1 时线性卷积结果中本该分开的部分会在 N 点圆环上缠绕重叠输出看起来像信号回声混杂。验证一下就非常直观x [1, 2, 3, 4]; h [1, 1, 1]; N length(x) length(h) - 1; % 结果长度 6 yc ifft(fft(x, N) .* fft(h, N)); % 补零到 N 点再做频域相乘 yl conv(x, h); % 直接线性卷积yc和yl完全一致。机制上可以这样理解DFT 隐含周期延拓补零相当于把原本会缠绕的圆环“撑大”让线性卷积的非零部分互不重叠。N 只要大于等于 LM-1 即可再大也只是频域采样更密结果不会变。这个补零条件不是工程习惯而是循环卷积与线性卷积等价的充要条件。配套代码里凡是出现fft(x, N) .* fft(h, N)的地方N 的取值必然满足这个约束。工程上用 FFT 做卷积图的是 O(N log N) 比直接卷积的 O(N²) 快得多。但省掉的复杂度也带来了边界条件你必须记住 N 至少是 LM-1否则结果就是错的。很多人在长信号滤波里拿默认长度直接相乘得到一段看似平滑但首尾异常的输出多半就是栽在这个补零长度上。3.2 相关与卷积共轭倒序的时间和代码相关运算和卷积在教科书里是两个符号在代码里却是同一个计算路径。互相关函数 r[m] Σ x[nm]·conj(y[n])用卷积实现时需要把 y 做时间翻转再加共轭。这个关系在配套代码中频繁出现尤其是时延估计和匹配滤波的例程。x [1, 2, 3, 4, 5]; y [0, 0, 1, 2, 3]; % y 看起来是 x 延迟了 2 个点 [c, lags] xcorr(x, y); [~, idx] max(c); delay lags(idx); % 理想情况下等于 2xcorr内部返回的 lags 从 -(M-1) 到 M-1峰值位置对应两个信号最相似的相对位移。为什么相关要翻转因为卷积的定义是 Σ x[k]·h[n-k]把 h 翻转后再滑动相关是 Σ x[nm]·conj(y[n])不翻转直接滑动。要让它们共用一个卷积引擎就必须把 y 提前翻转一次翻转之后再按卷积滑动效果就和相关一致。频域实现里这个翻转表现为 FFT 后取共轭时域翻转加共轭对应频域共轭。匹配滤波器本质上就是让接收信号与发送波形做相关所以雷达、声呐里的时延估计全部归结为找相关峰。配套代码里演示相关时你会发现用卷积实现和直接算相关结果一致但边界模式有差异——这是实现细节不是数学差异。3.3 overlap-add长信号滤波时为什么要分段FFT 卷积虽然快但一次处理几十万点意味着整段时延不可控内存开销也大。实时系统里信号是源源不断进来的不可能等整段收齐再做 FFT。overlap-add 方法把长信号切成块逐块做 FFT 滤波再把输出按重叠区相加解决了“无限长输入 有限长滤波器”的矛盾。M 64; % 滤波器长度 B 256; % 分块长度 L M B - 1; % FFT 长度 % H 是滤波器 h 的 L 点 FFT预先算好 y zeros(1, total_len M - 1); for idx 1:B:total_len block x(idx:min(idxB-1, total_len)); block [block, zeros(1, L - length(block))]; % 补零到 L 点 yseg ifft(fft(block, L) .* H); y(idx:idxL-1) y(idx:idxL-1) yseg; % 重叠相加 end每个分块的长度是 B滤波器长度是 M。对每个块做 L 点 FFT 卷积得到的线性卷积长度是 BM-1比本块多出的 M-1 个点属于滤波器的“拖尾”必须叠加到下一块输出上。这就是为什么叫 overlap-add——相邻输出块之间重叠了 M-1 个点。常见错误是忘了把重合区加回去或者块长和 FFT 长度没满足 L ≥ BM-1结果块边界出现明显的接缝噪声。你还会看到 overlap-save 的变体它用重叠输入块加丢弃尾部的方式规避相加操作两者取舍在实现复杂度上各有优劣。配套代码选用哪一种不重要重要的是理解分块没有改变线性卷积的数学本质只是把一个大卷积拆成了互不干扰的小卷积再拼回去。4. 滤波器设计源码从理想频响到可实现系统的多次妥协4.1 窗函数法为什么“截断”能等同于“调制”理想低通滤波器的冲激响应是无限长的 sinc 函数非因果、不可实现。窗函数法用一个有限长窗 w[n] 去截断它h[n] h_d[n] · w[n]。时域相乘对应频域卷积理想频响与窗谱卷积的结果就是实际频响——矩形窗的陡峭边缘被抹平成过渡带旁瓣变成阻带波纹。配套代码里fir1就是窗函数法的封装。一个典型调用h fir1(50, 0.25); % 50 阶低通归一化截止频率 0.25 [H, w] freqz(h, 1, 1024); % H 是复频响w 是数字角频率 0~pifreqz返回的 w 从 0 到 π对应物理频率 0 到 fs/2。画出幅度谱你会发现过渡带宽度和阻带衰减这两个指标被窗函数直接绑定Hamming 窗主瓣宽度约 8π/N阻带衰减约 53 dBBlackman 窗主瓣更宽换来的阻带衰减可以到 58 dB 以上。这不是玄学而是窗谱卷积的定量结果。窗谱的主瓣宽度决定过渡带因为理想频响的阶跃边缘被窗谱主瓣“涂抹”开旁瓣高度决定阻带波纹幅度因为泄漏过去的能量压不下来。配套代码里常用一个参数表来对比不同窗的设计结果。窗函数主瓣宽度(近似)阻带衰减(近似)过渡带代价矩形窗4π/N21 dB最窄但波纹大Hamming8π/N53 dB过渡带加宽一倍Blackman12π/N58 dB过渡带最宽选择哪个窗本质上是在过渡带陡峭和阻带干净之间做交易。理解这一点你再看代码里fir1的默认 Hamming 窗就不会觉得它是个随便选的经验值。4.2 阶数估算公式为什么要“宁可多一拍”设计滤波器时先估算阶数再验证指标这是配套例程里的标准流程。凯泽窗法的阶数估算公式大致形式是 N ≈ (A - 8) / (2.285 · Δω)其中 A 是期望阻带衰减 dB 值Δω 是归一化过渡带宽度。这个公式给出的是满足指标的近似最小阶数。实际代码里有个微妙的处理估算结果几乎总是要向上取整并且可能再加一两拍。原因有三。第一公式推导时假设的是理想化连续频率变量而实际 FFT、量化都会引入边界偏差。第二线性相位 FIR 滤波器要求冲激响应严格对称阶数 N 必须是偶数才能构成第一类线性相位滤波器估算出来是奇数就要加一。第三公式给的是“临界值”在临界点上设计出的滤波器刚好压线通过指标任何数值误差都可能让它掉出指标。与其反复试错不如直接多给几阶代价只是一点点计算量。配套代码里你常看到这样的处理N ceil(N_est); % 向上取整 if rem(N, 2) 1 N N 1; % 保证线性相位 FIR 的对称性 end这个“保守一拍”的习惯看起来不优雅却是工程里非常实用的设计取舍。滤波器阶数多一阶过渡带只会略微变窄但指标余量大了很多。4.3 零相位滤波器filtfilt 为什么好看却不能用于实时配套代码的离线分析例程里经常用filtfilt代替filtery filtfilt(b, a, x); % 先正向后反向等效相位为零filtfilt的原理是先让信号正向通过滤波器然后把输出序列在时间上翻转再通过同一个滤波器最后再翻转回来。两次滤波的相位延迟一正一负完全抵消总相位响应恒为零幅度响应变成 |H(e^jω)|²阻带衰减翻倍。优点是波形不失真看起来非常漂亮。但这里有一个根本性的代价非因果。第二次滤波时信号的“未来”部分已经作为输入参与计算了输出依赖整段信号的未来信息。实时系统里未来数据还没到达根本没法执行这种操作。所以filtfilt只能用于离线批处理。配套代码里如果要做在线滤波演示用的必然是filter而且起始段会有明显的瞬态响应。很多初学者看完filtfilt的漂亮输出就把在线系统也换成它等到低延迟场景一测发现算法根本无法实时运行这就是没分清因果性和非因果性的差别。filter的起始瞬态可以通过filtic设置初始状态来压低但这是另一个话题了。4.4 参考文献在代码里的价值参数从哪来边界条件在哪定义这套资源里除了 MATLAB 代码还有一份参考文献列表。很多人忽略它觉得文献是“论文里装点门面用的”其实不然。代码里很多看似凭空冒出来的参数都能回溯到文献中的推导。凯泽窗的 β 参数与阻带衰减的近似关系式、各类窗函数的频谱解析式、滤波器阶数估算公式的系数都来自原始论文中的推导或图表。配套代码注释里标注“见参考文献[x]”的地方就是在告诉你这个常量不是调参调出来的而是某个定理或近似公式的产物。明文规定的边界条件也一样。线性相位 FIR 的对称性分类、IIR 滤波器稳定性条件、窗函数设计法中通带/阻带纹波的定义方式这些在教材正文里分散在各处代码的注释往往只写结论。当你的代码行为和书本推导冲突时最优的排查路径不是盲调参数而是翻开文献目录里对应条目回到定义源头核对。参考文献在代码世界里扮演的是“设计意图注释”的角色——它解释这个公式从哪里来、在什么条件下成立比代码里的注释更接近算法的原始定义。5. 避坑排查信号处理代码里最容易翻车的五个细节5.1 高频错误的“现象、原因、解决”清单频谱幅值对不上现象对正弦信号做 FFT 后峰值幅度和信号真实幅度差很多有时候是一半有时候还要乘个系数才对。原因FFT 结果包含正负频率两个对称分量实信号的能量被分成两半。传统画法用单边谱要把正频部分取出来乘 2再除以 N 才得到真实幅度。P2 abs(X / N); P1 P2(1:N/21); P1(2:end-1) 2*P1(2:end-1); % 直流分量不用乘 2 f (0:N/2) / N * fs;解决先打印 P2 在峰值处的原始值验证是不是幅值一半再检查除没除 N。直流分量是一个例外它没有负频对应项不能乘 2。这条做对了幅度再去查相位。滤波器输出开头有一段混乱现象filter的输出前几十个点不跟随输入规律过一会儿才正常。原因滤波器内部状态初始为零输出要经历一个暂态过程才进入稳态。这不是滤波器坏了而是零初始条件下的自然响应。IIR 滤波器反馈越深暂态持续越长。解决在线场景里预热一组数据再取有效输出离线场景可以用filtic估算稳态初值让起始状态直接跳到工作点。千万别在没搞清楚状态变量的情况下改滤波器系数那只会让稳态也跟着坏掉。频率轴差一格或多出镜像现象画出abs(fft(x))后横轴从 0 到 fs图像有两个对称峰或者你只想要 0 到 fs/2 的部分峰值却落在最后一格。原因FFT 输出 k0 对应直流kN/2 对应奈奎斯特频率 fs/2再往后是负频率部分。没做单边截断时频谱自然是对称的双边谱。解决明确画图目的。分析真实信号的频谱用单边截断f 轴取 (0:N/2)/N*fs想展示周期延拓特性用fftshift把零频移到中心画 -fs/2 到 fs/2。时移操作不生效现象把信号频谱乘上 e^{-j2πk·m/N} 之后做 ifft波形确实移动了但移出去的部分从另一端冒了出来像是环绕一周。原因DFT 的时移定理给出的是循环时移序列在 N 点圆环上移动越过边界就绕回起点。这是 DFT 的周期延拓假设的直接后果不是代码写错。解决理解“DFT 下的时移天然是循环的”这件事。想要线性延迟要么用circshift配合索引取模运算要么在时域直接补零后整体移动。验证定理本身时接受环绕现象就好。用 fir1 设计完滤波器增益不对现象设计低通滤波器期望通带增益 1实际freqz看到的幅度不是 1或者带内整体变小。原因fir1默认会自动缩放滤波器系数让通带中心频率处的增益精确为 1。当设计带通或者其他需要任意增益时默认缩放不符合预期。解决查看代码里是否有scale参数或者显式的增益归一化步骤。需要自定义增益就自行对 h 乘以目标增益值再用freqz验证。这条规则很隐蔽因为 FIR 线性相位滤波器的系数和增益之间的比例关系不影响形状但影响幅度。5.2 排查顺序先验证坐标轴再验证算法路径遇到频谱异常很多人的第一反应是把 FFT 换成更高阶的算法或者怀疑浮点数精度。血泪经验告诉我先别碰算法先验证坐标轴和增益约定。可以用一个已知信号把整条链路打通fs 1000; N 1024; t (0:N-1) / fs; x 1.5 * sin(2*pi*50*t); % 幅度 1.5频率 50Hz X fft(x); f (0:N-1) / N * fs; [val, idx] max(abs(X)); fprintf(峰值频率 %.2f Hz, 幅度 %.3f\n, f(idx), 2*val/N);如果这个测试输出接近“峰值频率 50 Hz幅度 1.500”说明 FFT 这条路没问题问题出在输入信号的构造或滤波环节。排查顺序我一般固定为先确认采样率和频率刻度对不对再检查补零长度是否满足卷积约束然后用fftshift确认频谱排列方式最后才怀疑算法实现本身。数字信号处理的代码里整体波形“看起来差不多”很容易掩盖坐标轴的错误最有效的防呆手段是拿已知信号做一次全链路对照。这也解释了为什么配套代码每个 demo 前都要造一个简单的正弦波——它不是教学仪式而是调试基准。6. 设计取舍教学代码的价值在“可读的映射”不在性能这套配套代码里大量使用 for 循环和逐点操作而不是向量化矩阵运算。从性能角度看这些代码跑起来比工程优化版慢不少但这个取舍是故意的。循环结构的每一步都对应课本公式里的一步一个 for k 循环对应一个求和号一个移位操作对应一个延迟代码和公式之间的映射几乎是逐字对应的。教学代码的首要目标不是运行速度而是让人能在公式和实现之间来回对照。如果一开始就用矩阵化的奇异技巧阅读者反而看不出它在算哪个公式。工程上把这些循环改成矩阵乘法或 FFT 分块通常能获得一个数量级的加速但代价是映射关系变模糊。我有过一次印象深刻的教训把某个相关性计算的循环展开成矩阵乘法后速度确实提升了二十倍但共轭取反那一步写错输出全成了噪声更麻烦的是矩阵版本的结果看起来太规整第一眼根本没发现异常。从那以后我每次改完计算路径都强制拿一段已知信号跑一遍和参考文献中的原始公式核对差值而不是盯着最终输出看“对不对”。把可读性版本和优化版本做差分测试比在结果图上看形状可靠得多。这套资源真正稀缺的地方不是代码本身能跑多快而是它把公式、代码、参考文献三者串成了同一条线索。当你碰到“代码行为和书不一致”的时候文献目录就是地图带你回定义源头重新核对边界条件。学会用这套映射去阅读他人代码、核对自己的实现比记住任何一个具体函数都值钱。希望帮到你。本文还有配套的精品资源点击获取