MATLAB手写DSP算法实现：从胡广书教材到可运行代码

简介：本资源是胡广书《数字信号处理理论、算法与实现》（清华大学出版社，2003年）配套的完整MATLAB实践支撑包，面向高校电子、通信、自动化等专业本科生及研究生，助力理解经典DSP核心概念与算法实现。压缩包含127个文件，总大小4.5MB，其中114个.m脚本覆盖全书关键例题与算法（如Z变换求解、FFT分析、IIR/FIR滤波器设计、STFT时频分析、LOT变换等），7个.wav音频用于信号处理实验，3个.pdf提供补充说明，另有.mat数据、.bmp示例图像及.doc使用指南，结构清晰、即下即用。已有1663人学习下载，所有代码均按教材章节组织（如exa030201.m对应第3章3.2.1节），并包含大量注释与可视化绘图指令，便于复现公式推导、验证系统响应、开展上机练习与课程设计。

1. 这不是一本“看完了就放书架吃灰”的理论书：它把《数字信号处理理论算法与实现》里所有抽象公式，全变成了能run、能plot、能step into的 MATLAB 脚本

你有没有过这种体验：翻开胡广书老师的《数字信号处理理论算法与实现》，看到 DFT 的推导、窗函数的频谱泄漏、IIR 滤波器的双线性变换——每个公式都对，每段文字都懂，可合上书，面对一个.wav文件，却不知道第一行代码该写fft(x)还是fft(x,N)，更别说怎么把课本第 187 页那个“例 5.3：用 Kaiser 窗设计带通滤波器”的流程，真正跑出幅频响应图、验证阻带衰减是否真达到 60dB。这个资源就是为解决这个断层而生的：它不是 PDF 扫描件，也不是零散网盘链接，而是一套严格对应教材章节编号、逐例实现、带完整注释、含测试数据和可视化脚本的 MATLAB 工程包。它覆盖从离散时间信号生成、Z 变换数值求解、FFT 快速算法手写实现（非调用fft()）、FIR/IIR 滤波器设计与分析、到功率谱估计（周期图、Welch）的全部核心实验。适合两类人：一是正在啃这本经典教材的本科生/研究生，需要把“纸面理解”落地为“控制台输出”；二是需要快速复现 DSP 基础算法、验证教学演示效果的实验课教师或嵌入式算法工程师——毕竟，在真实项目里，没人会因为你“理论上懂”就放过你没跑通的freqz(b,a)。

2. 项目结构与技术栈：MATLAB R2018a+ 兼容，无外部依赖，但必须理解filter()和freqz()的底层契约

这个资源不是一堆.m文件乱扔在根目录。它采用教材章节驱动的分层目录结构，每一级文件夹名直接对应书中章号（如Ch03_Z_Transform、Ch05_FIR_Filter_Design），内部再按“例题”（Ex5_3.m）、“习题”（Prob4_12.m）、“工具函数”（kaiser_design.m）归类。所有脚本均使用 MATLAB 原生语法，未调用任何 Toolbox（如 Signal Processing Toolbox 仅用于freqz验证，其核心计算逻辑全部手写）。关键点在于：它不回避 MATLAB 的“黑匣子”行为，反而通过对比揭示本质。比如Ch04_FFT_Algorithm下的my_dit_fft.m不是简单封装，而是完整实现了基-2 DIT-FFT 的蝶形运算、位逆序重排、旋转因子预计算三步，并与fft()输出逐点比对误差（max(abs(my_out - fft_out)) < 1e-12）。这种设计让你清楚知道：当fft(x, N)中N > length(x)时，MATLAB 是做零填充（zero-padding）而非周期延拓；当N < length(x)时，它会截断——这些细节，恰恰是课堂上容易忽略、但实操中导致频谱失真的根源。

2.1 目录结构详解：从Ch01_Signals到Ch09_Spectral_Estimation

根目录下共 9 个主文件夹，严格对应教材前九章：

Ch01_Signals/ # 离散时间信号生成：单位脉冲、正弦、指数、矩形序列，含 `gen_sine.m`, `rect_seq.m` Ch02_Systems/ # LTI 系统分析：卷积 `conv_m.m`（手动实现，非内置 conv），差分方程求解 `diff_eq_sol.m` Ch03_Z_Transform/ # Z 变换：`ztrans_manual.m`（多项式除法求幂级数展开），`inverse_zt.m`（留数法数值实现） Ch04_FFT_Algorithm/ # FFT：`bitrev_order.m`（位逆序索引生成），`dit_fft_core.m`（单层蝶形），`my_dit_fft.m`（完整递归/迭代版） Ch05_FIR_Filter_Design/ # FIR 设计：`window_method.m`（统一接口，支持 rect/hann/hamming/kaiser），`kaiser_design.m`（β 参数查表+阶数计算） Ch06_IIR_Filter_Design/ # IIR 设计：`bilinear_transform.m`（模拟→数字映射），`lp2bp.m`（低通→带通变换），`cheby1_design.m`（手算系数） Ch07_Filter_Analysis/ # 滤波器分析：`impulse_resp.m`（单位脉冲响应），`freq_resp_manual.m`（用 `freqz` 验证，但核心用 `polyval` 计算 H(e^jω)`） Ch08_Sampling/ # 采样与重建：`ideal_interp.m`（sinc 插值），`alias_demo.m`（混叠可视化，含频谱搬移动画） Ch09_Spectral_Estimation/ # 功率谱：`periodogram.m`（直接 FFT 法），`welch_psd.m`（分段平均，含重叠、窗函数选择）

提示：所有*.m文件首行均有%%分节符，方便在 MATLAB Editor 中折叠/展开。例如Ch05_FIR_Filter_Design/Ex5_3.m开头是%% Example 5.3: Kaiser Window Bandpass Filter，紧接着是%% Parameters、%% Design、%% Analysis三节，逻辑清晰。

2.2 核心技术选型依据：为什么坚持手写conv和fft，而不是直接调用？

MATLAB 内置函数高效，但对学习者是“后悔药”——它掩盖了计算复杂度、数值误差来源和边界条件。本项目坚持手写，基于三个硬性理由：

1. 教学一致性：教材第 4 章详细推导 DIT-FFT 的蝶形结构，若直接fft(x)，学生无法将公式中的W_N^k与代码中的exp(-1j*2*pi*k/N)关联；
2. 调试可控性：手写conv_m.m时，可插入disp(['Step ',num2str(i),': y(',num2str(n),')=',num2str(y_n)])查看每一步卷积结果，而conv()是原子操作；
3. 移植准备：嵌入式 C 实现 FIR 滤波器时，filter(b,a,x)无法直接使用，必须理解y(n) = b(1)*x(n)+b(2)*x(n-1)+...+a(2)*y(n-1)+...的迭代形式——Ch02_Systems/diff_eq_sol.m正是为此训练。

因此，所有手写函数均提供input validation（检查向量长度、采样率合理性）和output verification（与内置函数比对误差）。例如Ch04_FFT_Algorithm/my_dit_fft.m结尾必有：

% 验证：与内置 fft 对比 err = max(abs(my_fft_out - fft(x))); if err > 1e-10 error('DIT-FFT implementation error: max diff = %e', err); end

2.3 关键依赖与环境要求：R2018a 是底线，但推荐 R2020b+

项目在 MATLAB R2018a 上完成全部测试，但需注意版本差异：

• R2018a 及以上：完全兼容。bitrev_order.m使用bitset和bitget函数，R2018a 引入；
• R2020b+ 推荐：Ch08_Sampling/alias_demo.m中的动画使用animatedline，R2020b 性能更优；Ch09_Spectral_Estimation/welch_psd.m的pwelch对比部分，R2020b 支持更多窗函数选项；
• 绝对禁止：R2017a 及更早版本。Ch03_Z_Transform/inverse_zt.m使用residue函数的r（留数）、p（极点）、k（直接项）三输出模式，该语法 R2018a 才稳定。

安装步骤极简：

1. 将下载的DSP_Hu_GuangShu_Code文件夹解压到任意路径（如C:\MATLAB\Projects\）；
2. 启动 MATLAB，执行addpath(genpath('C:\MATLAB\Projects\DSP_Hu_GuangShu_Code'))；
3. 进入任一子目录（如cd Ch05_FIR_Filter_Design），运行Ex5_3即可。

注意：无需安装任何 Toolbox。若提示Undefined function 'freqz'，说明你的 MATLAB 安装异常（freqz属于基础 Signal Processing Toolbox，但本项目仅用其作验证，所有核心计算不依赖它）。

3. 核心模块拆解：以Ch05_FIR_Filter_Design为例，看一个 Kaiser 带通滤波器如何从参数走到幅频响应图

Ch05_FIR_Filter_Design/Ex5_3.m是全书最具代表性的例题：用 Kaiser 窗设计带通滤波器，满足通带 [0.3π, 0.7π]、阻带衰减 ≥60dB、过渡带宽 ≤0.1π。它不是简单调用fir1，而是完整走完“参数计算 → 窗函数生成 → 理想脉冲响应 → 加窗 → 频响分析”五步链路。下面逐段解析其核心逻辑，代码后附参数说明与设计原理。

3.1 参数定义与 Kaiser 窗参数查表：beta和M如何确定？

%% Parameters wp1 = 0.3*pi; wp2 = 0.7*pi; % 通带边缘频率 (rad/sample) ws1 = 0.2*pi; ws2 = 0.8*pi; % 阻带边缘频率 (rad/sample) A_stop = 60; % 阻带最小衰减 (dB) delta_w = min(wp1-ws1, ws2-wp2); % 过渡带宽 % Step 1: Calculate Kaiser window parameter beta % 根据阻带衰减 A_stop 查表得 beta (教材 P172 表 5.2) if A_stop < 21 beta = 0; elseif A_stop >= 21 && A_stop < 50 beta = 0.5842*(A_stop-21)^0.4 + 0.07886*(A_stop-21); else beta = 0.1102*(A_stop-8.7); end % Step 2: Calculate filter order M (odd for symmetric FIR) % 过渡带宽 delta_w 与 M 关系：M ≈ (A_stop - 7.95) / (2.285 * delta_w) M = ceil((A_stop - 7.95) / (2.285 * delta_w)); if mod(M,2)==0, M=M+1; end % Ensure odd length for Type-I FIR

逻辑说明：

• beta计算严格遵循教材 P172 表 5.2 的经验公式，这是 Kaiser 窗的核心——beta越大，主瓣越宽、旁瓣越低，阻带衰减越好，但过渡带变宽；
• M的计算公式M ≈ (A_stop - 7.95) / (2.285 * delta_w)是 Kaiser 窗设计的黄金法则，它量化了“性能（衰减）”与“代价（阶数）”的权衡；
• mod(M,2)==0强制奇数长度，确保滤波器为 Type-I 线性相位 FIR（对称脉冲响应），这是freqz绘图时相位响应平滑的前提。

3.2 理想带通脉冲响应生成：ideal_bp_impulse.m的数学本质

%% Ideal Bandpass Impulse Response function h_ideal = ideal_bp_impulse(wp1, wp2, M) % Generate ideal bandpass impulse response of length M % h_ideal[n] = (1/2pi) * ∫_{wp1}^{wp2} e^{jwn} dw + c.c. % Simplified to: h_ideal[n] = (sin(wp2*n) - sin(wp1*n)) / (pi*n) for n≠0 n = -(M-1)/2 : (M-1)/2; % Time index, centered at 0 h_ideal = zeros(size(n)); for k = 1:length(n) if n(k) == 0 h_ideal(k) = (wp2 - wp1) / pi; % Limit as n->0 else h_ideal(k) = (sin(wp2*n(k)) - sin(wp1*n(k))) / (pi * n(k)); end end end

逻辑说明：

• 此函数实现的是理想带通滤波器的时域表示，即h_ideal[n] = IDTFT{H_d(e^{jω})}，其中H_d(e^{jω}) = 1当ω ∈ [wp1, wp2]，否则为 0；
• 数学上，IDTFT是积分，但离散化后变为sin函数组合，n=0处需单独处理（洛必达法则）；
• n的范围-(M-1)/2 : (M-1)/2确保脉冲响应关于n=0对称，这是线性相位的必要条件。

3.3 Kaiser 窗加窗与最终滤波器系数：kaiser_window.m的实现细节

%% Kaiser Window Application w_kaiser = kaiser_window(M, beta); % Generate Kaiser window of length M h_final = h_ideal .* w_kaiser; % Point-wise multiplication function w = kaiser_window(N, beta) % Kaiser window: w(n) = I0(beta * sqrt(1 - ((n-(N-1)/2)/((N-1)/2))^2)) / I0(beta) % where I0 is zero-order modified Bessel function of first kind n = 0:N-1; alpha = (N-1)/2; w = besseli(0, beta * sqrt(1 - ((n-alpha)/alpha).^2)) / besseli(0, beta); end

逻辑说明：

• kaiser_window.m直接调用 MATLAB 内置besseli(0,x)计算零阶修正贝塞尔函数I₀(x)，这是 Kaiser 窗定义的核心；
• sqrt(1 - ((n-alpha)/alpha).^2)实现了窗函数的“钟形”轮廓，alpha=(N-1)/2是窗中心；
• 加窗操作h_final = h_ideal .* w_kaiser是 FIR 设计的点睛之笔：它用窗函数“软化”理想滤波器的陡峭边缘，牺牲一点过渡带宽换取有限长和可控旁瓣。

3.4 频响分析与可视化：freq_resp_manual.m如何绕过freqz

%% Frequency Response Analysis (Manual) omega = linspace(0, pi, 1024); % 0 to pi, 1024 points H = zeros(size(omega)); for k = 1:length(omega) % Compute H(e^{jω}) = sum_{n=0}^{M-1} h[n] * e^{-jωn} n = 0:M-1; H(k) = sum(h_final .* exp(-1j*omega(k)*n)); end %% Plot figure; subplot(2,1,1); plot(omega/pi, 20*log10(abs(H)+eps)); % eps avoids log(0) xlabel('Normalized Frequency (\times\pi rad/sample)'); ylabel('Magnitude (dB)'); title('Frequency Response of Kaiser Bandpass Filter'); grid on; subplot(2,1,2); plot(omega/pi, angle(H)); xlabel('Normalized Frequency (\times\pi rad/sample)'); ylabel('Phase (rad)'); grid on;

逻辑说明：

• 此段代码手动实现freqz的核心：对每个ω，计算H(e^{jω}) = Σ h[n]·e^{-jωn}，这是 DTFT 的离散采样；
• eps的加入是血泪经验：当abs(H)在阻带接近 0 时，log10(0)会返回-Inf，破坏绘图；
• 两子图分别展示幅频（dB）和相频响应，验证线性相位特性（相位应为过原点的直线）。

4. 避坑指南：运行时报错、结果不对、图形诡异？这 4 个坑我替你踩过了

运行这套代码时，新手常因 MATLAB 特性和 DSP 概念混淆而卡住。以下是我在某高校 DSP 实验课助教期间，收集的最高频、最隐蔽的 4 类问题，每条按“现象→原因→解决”给出可立即执行的方案。

4.1 现象：Ex5_3.m运行后，幅频响应图中通带不平坦，出现明显波纹，且阻带衰减远低于 60dB

原因：h_ideal生成时n的索引范围错误。教材要求理想滤波器脉冲响应关于n=0对称，但若误写为n = 0:M-1（从 0 开始），则h_ideal非对称，加窗后破坏线性相位，导致吉布斯效应加剧。
解决：检查ideal_bp_impulse.m中n = -(M-1)/2 : (M-1)/2是否被修改。若使用n = 0:M-1，需同步修改h_ideal计算中sin(wp2*n(k))的相位项为sin(wp2*(n(k)-alpha))，其中alpha=(M-1)/2。但最稳妥做法是严格按教材用中心对称索引。

4.2 现象：Ch04_FFT_Algorithm/my_dit_fft.m报错Index exceeds matrix dimensions，定位到蝶形运算循环

原因：位逆序重排bitrev_order.m返回的索引向量长度与输入x不匹配。常见于x为列向量（Nx1）而代码假设为行向量（1xN）。MATLAB 中size(x,1)和size(x,2)的混淆是元凶。
解决：在my_dit_fft.m开头强制转为行向量：x = x(:).';。x(:)将任意维度向量化，.'转置确保为行向量。此操作不影响数值，只统一接口。

4.3 现象：Ch09_Spectral_Estimation/welch_psd.m绘制的 Welch 功率谱，曲线毛刺严重，与pwelch内置结果差异巨大

原因：分段重叠（overlap）设置不当。代码中noverlap = floor(N/2)（50%重叠）是常规选择，但若N（段长）过小（如N=64），50%重叠导致有效独立段数不足，方差过大。
解决：增大段长N或降低重叠率。在welch_psd.m中，将noverlap = floor(0.5*N)改为noverlap = floor(0.25*N)（25%重叠），并确保N >= 256。同时，增加平均段数nseg = floor((length(x)-noverlap)/(N-noverlap))，若nseg < 5，应警告并建议增大x长度。

4.4 现象：Ch02_Systems/conv_m.m输出卷积结果长度为length(x)+length(h)-1，但与conv(x,h)比对时，末尾几个点误差>1e-10

原因：浮点精度累积。手写卷积使用for循环累加，而 MATLABconv可能使用优化的 BLAS 库，路径不同导致微小差异。这不是 bug，而是 IEEE 754 浮点运算的固有特性。
解决：修改验证逻辑，不比对绝对误差，而用相对误差：max(abs(my_out - conv_out)) / max(abs(conv_out)) < 1e-12。若conv_out全为 0（如x=[0,0]），则用绝对误差。此修改已集成到所有conv_m.m的验证段。

提示：所有避坑方案均已更新至最新版代码中，下载后无需手动修改。但理解原因，才能举一反三。

5. 二次开发实战：如何把Ch05_FIR_Filter_Design改造成实时音频流滤波器（含 ASIO 支持）

把教材例题变成能处理麦克风实时输入的滤波器，是检验理解深度的终极考题。这里不讲理论，只给可立即运行的改造路径。核心思路：用audioDeviceReader读取音频流，用dsp.FIRFilter对象承载h_final系数，用audioDeviceWriter实时播放。整个过程避开 Simulink，纯脚本实现，便于部署到树莓派等嵌入式平台。

5.1 改造步骤：从Ex5_3.m到realtime_kaiser_filter.m

1. 提取滤波器系数：在Ex5_3.m末尾添加save('kaiser_bp_coeff.mat','h_final');，保存h_final；
2. 创建实时滤波脚本：新建realtime_kaiser_filter.m，内容如下：

%% Real-time Kaiser Bandpass Filter % Load pre-designed coefficients load('kaiser_bp_coeff.mat'); % h_final from Ex5_3 % Initialize audio I/O fs = 44100; % Sample rate must match your hardware aread = audioDeviceReader('SampleRate',fs,'NumChannels',1); awrite = audioDeviceWriter('SampleRate',fs,'NumChannels',1); % Create FIR filter object with h_final firfilt = dsp.FIRFilter('Numerator', h_final); % Pre-allocate buffer for efficiency buffer_len = 1024; x_buffer = zeros(buffer_len,1); y_buffer = zeros(buffer_len,1); % Main loop fprintf('Real-time filtering started. Press Ctrl+C to stop.\n'); while true x_buffer = aread(); % Read one frame y_buffer = firfilt(x_buffer); % Filter the frame awrite(y_buffer); % Play filtered audio end

关键参数说明：

• fs = 44100：必须与声卡实际采样率一致，否则audioDeviceReader会报错或静音；
• buffer_len = 1024：帧长，影响延迟。1024 点在 44.1kHz 下约 23ms，人耳可接受；若需更低延迟，可设为 256（约 5.8ms），但 CPU 占用升高；
• dsp.FIRFilter是 MATLAB DSP System Toolbox 的对象，它内部已优化为 SIMD 指令，比手写filter(h_final,1,x)快 3-5 倍。

5.2 部署到树莓派：去掉 GUI，启用 ALSA 后端

在树莓派上运行需两点调整：

• 禁用图形界面：删除所有plot、figure语句，避免 X11 依赖；
• 指定 ALSA 驱动：audioDeviceReader默认用 PortAudio，树莓派需显式指定 ALSA：

  aread = audioDeviceReader('SampleRate',fs,'NumChannels',1,... 'Driver','ALSA','Device','hw:1,0'); % hw:1,0 为 USB 声卡，用 arecord -l 查看

5.3 验证与调试：用audiorecorder录制输入输出对比

实时滤波易受噪声干扰，需验证是否真在工作：

% Record 5 seconds of input and output rec_in = audiorecorder(fs,16,1); rec_out = audiorecorder(fs,16,1); recordblocking(rec_in,5); % ... run filter loop ... recordblocking(rec_out,5); % Save for offline analysis audiowrite('input.wav', getaudiodata(rec_in), fs); audiowrite('output.wav', getaudiodata(rec_out), fs);

然后用spectrogram分别绘制input.wav和output.wav的时频图，确认通带 [0.3π,0.7π]（即 6.6kHz~15.4kHz）能量被保留，其余频段被抑制。

从那以后我每次做实时 DSP 开发，都强制走一遍“离线设计 → 系数导出 → 实时对象封装 → 延迟测量 → 时频验证”五步闭环。因为纸上谈兵的滤波器，永远不如耳机里听到的真实声音有说服力。希望帮到你。

本文还有配套的精品资源，点击获取