简介：面向中南大学数字信号处理课程的课程设计任务书，适合该校相关专业学生完成DSP课程设计、复习理论与实践结合知识时参考。文档系统梳理了三大核心选题：连续信号采样与DFT谱分析及参数选择、周期方波信号滤波（要求滤除40Hz后分量并处理噪声）、音乐信号处理（设计单回声、多重回声与全通混响器并观察频谱变化），并给出GUI界面演示与报告撰写要求。资源包共1个doc文档，大小约318KB，内容覆盖设计目的、设计内容、设计要求、巴特沃斯/FIR/IIR滤波器原理、程序设计思路、测试输出结果、总结与参考文献等模块，便于对照完整流程撰写报告或准备验收。已有397人学习下载，适合需要按任务书逐项落实算法仿真、参数分析和滤波实验的本科阶段读者。
1. 从课程设计任务书看数字信号处理实验的完整闭环
中南大学这份数字信号处理课程设计任务书，表面上是三道习题，实际上覆盖了DSP实验课最核心的三条主线：频谱分析怎么做才不漏峰、滤波器设计怎么绕过工具箱函数、音频效果器怎么从差分方程落地成可听的Demo。和其他学校动辄要求“设计一个语音压缩系统”的题目不同，这份任务书把重点压在DFT参数选择和滤波器结构理解上，适合用来检验学生对采样定理、频率分辨率、FIR/IIR结构差异的掌握程度，也适合作为电子信息类学生的综合实训模板。任务书要求用MATLAB或其它高级语言实现，并明确要求尽量避免现成工具箱函数——这意味着核心算法要自己写，GUI只是最后一层壳。下面按理论、实现、验证的顺序拆开讲，代码基于MATLAB R2016a版本，全部可运行。
2. DFT谱分析：频率分辨率、栅栏效应与参数选择的MATLAB复现
2.1 分辨率公式背后的采样率陷阱
任务书第一大题给了一个经典信号：
x(n) = sin(2π·0.125·n) + cos(2π·(0.125+Δf)·n)
n = 0,1,...,N-1
这个信号的特点是：两个频率分量间隔恰好是Δf，而且都归一化到了数字频率。第一问要求N=16时分别取Δf=1/16和Δf=1/64观察频谱，第二问要求N=128时保持Δf不变再观察。这里有个初学者容易忽略的关键点：数字频率分辨率是2π/N，对应模拟频率分辨率是fs/N。当Δf=1/64而N=16时，频率间隔小于分辨率，频谱上两个峰必然合成一个宽包络，这不是窗函数选取的问题，而是DFT本身的观测能力极限。
代码实现上，要注意频谱显示必须用实际频率或归一化频率，不能用原始k值。我一般这样写：
N = 16; df = 1/64; n = 0:N-1; xn = sin(2*pi*0.125*n) + cos(2*pi*(0.125+df)*n); Xk = fft(xn, N); mag = abs(Xk(1:N/2+1)); % 只取单边谱 f_axis = (0:N/2)/N; % 归一化频率轴，单位是cycles/sample plot(f_axis, mag, '-o'); xlabel('归一化频率 f/fs'); ylabel('幅度'); grid on;这段代码里，fft(xn, N)的第二个参数明确指定了FFT点数，当长度不足时自动补零。实际做实验时你会发现，即使N=16、Δf=1/64时频谱已经无法分辨两个峰，但补零到1024点后曲线变得光滑了——注意这是插值，不是分辨率提升。这个区分在课程设计报告里要写清楚，否则答辩容易被问住。
2.2 补零操作与分辨率的关系验证
继续验证N=128时的表现。当N增大到128，频率分辨率变为1/128，Δf=1/64时两个谱线恰好间隔2个分辨单元，可以明显看到双峰。但要做对比实验，建议把四种情况放在同一张图上：
figure; cases = {16, 1/16; 16, 1/64; 128, 1/16; 128, 1/64}; for i = 1:4 N_cur = cases{i,1}; df_cur = cases{i,2}; xn = sin(2*pi*0.125*(0:N_cur-1)) + cos(2*pi*(0.125+df_cur)*(0:N_cur-1)); Xk = fft(xn, N_cur); mag = abs(Xk(1:N_cur/2+1)); f_axis = (0:N_cur/2)/N_cur; subplot(2,2,i); stem(f_axis, mag, 'filled', 'MarkerSize', 3); title(sprintf('N=%d, Δf=1/%d', N_cur, 1/df_cur)); xlabel('归一化频率'); ylabel('幅度'); grid on; end观察这四张子图，能直接验证任务书里“提高N可以降低频谱泄露但旁瓣相对幅度不减小”的结论。实际操作中，如果波形出现明显的不对称或毛刺，优先检查是否漏了fftshift。多数教材里的频谱图以零频为中心，而fft的结果是0到fs，需要fftshift才能正确显示负频率部分。
2.3 频谱泄露与窗函数选择
截断效应导致的频谱泄露是这份任务书第一个隐含考点。矩形窗主瓣窄但旁瓣高，汉宁窗主瓣宽但旁瓣低。课程设计如果时间充裕，可以加一个对比实验：
N = 64; n = 0:N-1; x = cos(2*pi*0.125*n); w_rect = rectwin(N)'; w_hann = hann(N)'; X1 = abs(fft(x.*w_rect, 1024)); X2 = abs(fft(x.*w_hann, 1024)); figure; plot((0:1023)/1024, 20*log10(X1/max(X1))); hold on; plot((0:1023)/1024, 20*log10(X2/max(X2))); legend('矩形窗', '汉宁窗'); ylim([-80, 10]);运行后能看到矩形窗第一旁瓣约-13dB，汉宁窗约-31dB，这就是为什么频谱分析时更推荐后者。注意补零后频谱变光滑，但主瓣宽度没有变窄，这正是区分“补零插值”和“增加N提高分辨率”的最好演示。
3. 周期方波滤波与巴特沃斯低通滤波器设计
3.1 方波频谱分析：傅里叶级数与DFT的对应关系
第二题需要生成10Hz基频的周期方波并对它做频谱分析。方波的傅里叶级数只含奇次谐波，即10Hz、30Hz、50Hz、70Hz……这意味着如果采样率取1000Hz，50Hz以远的谐波都落在阻带内，“滤除40Hz以后的频率分量”等价于保留基波和三次谐波，衰减五次以上谐波。
生成方波并分析频谱的代码：
Fs = 1000; T = 1/10; % 基频10Hz，周期0.1s t = 0:1/Fs:T*10; % 取10个周期 x = square(2*pi*10*t); % square函数生成周期方波，幅度±1 N = length(x); Xk = fft(x); mag = abs(fftshift(Xk))/N; f_axis = (-N/2:N/2-1)*Fs/N; subplot(211); plot(t, x); xlabel('t/s'); title('10Hz周期方波时域波形'); subplot(212); stem(f_axis, mag, 'MarkerSize', 2); xlim([0, 200]); xlabel('频率/Hz'); ylabel('幅度');这里有三点需要注意。第一，square(2*pi*10*t)产生的是-1到1的方波，占空比默认50%；第二，频谱用fftshift后零频在中间，显示时用xlim截取0-200Hz便于观察奇次谐波的衰减规律；第三，方波频谱的包络按1/f衰减，第n次谐波幅度约4/(nπ)，这在报告数据分析部分要体现出来。
3.2 巴特沃斯滤波器参数归一化流程
确定滤波器指标时，任务书给出的读取方式是：从频谱图上读出通带边界频率40Hz和阻带截止频率50Hz，通带最大衰减0.7dB，阻带最小衰减0.1dB。注意到MATLAB的buttord函数输入要求是归一化频率，所以40Hz和50Hz都要除以奈奎斯特频率500Hz（Fs/2），得到0.08和0.1。
滤波器设计代码：
Fs = 1000; Wp = 40/(Fs/2); % 通带边界归一化：40Hz/500Hz = 0.08 Ws = 50/(Fs/2); % 阻带截止归一化：50Hz/500Hz = 0.1 Rp = 0.7; % 通带最大衰减 dB Rs = 15; % 阻带最小衰减 dB [n, Wn] = buttord(Wp, Ws, Rp, Rs); [b, a] = butter(n, Wn); fprintf('滤波器阶数: %d, 截止频率: %.4f\\n', n, Wn); y = filter(b, a, x); Yk = abs(fftshift(fft(y)))/N; figure; subplot(211); plot(t, x); hold on; plot(t, y, 'r', 'LineWidth', 1.5); legend('原始方波', '滤波后'); xlabel('t/s'); subplot(212); plot(f_axis, mag); hold on; plot(f_axis, Yk, 'r'); xlim([0, 200]); xlabel('频率/Hz'); ylabel('幅度'); legend('滤波前频谱', '滤波后频谱');参数说明：buttord返回最小阶数n和3dB截止频率Wn，butter(n, Wn)设计低通滤波器并返回传递函数系数b和a。对M点数据用filter(b,a,x)做零状态滤波，得到输出y。滤波后50Hz以上的谐波明显被压制，时域波形从方波变成接近正弦的形态——这正是基波+三次谐波叠加的效果，五次以上谐波已经被衰减到可以忽略的程度。
3.3 设计过程中避不开的几个坑
这个题有个版本差异问题：老版MATLAB的buttord输入参数是模拟频率，用的是[n,Wn] = buttord(Wp,Ws,Rp,Rs,'s')这种带s参数的写法，而新版直接给归一化数字频率即可。课程设计报告里如果混用新旧版本，很容易出现“滤波器设计出来完全不对”的情况，现象是滤波后信号幅度异常或完全没有滤波效果。遇到这种情况，第一步先检查buttord返回的n是否合理——如果n返回1或2，多半是归一化频率写错了。
另一个常见问题是filter和filtfilt的选择。filter会产生相位延迟，滤波后的波形会向右偏移；filtfilt是零相位滤波，没有延迟但计算量翻倍且不能用于实时系统。课程设计用filter就好，但报告里要解释为什么波形有延时。
4. 回声与混响效果器：从差分方程到音频处理实现
4.1 单回声滤波器：FIR梳状滤波器的参数选择
第三大题要求设计单回声、多重回声和全通混响器。单回声的差分方程是：
y[n] = x[n] + α·x[n-R]
传输函数H(z) = 1 + α·z^(-R)
当α绝对值小于1时系统稳定。这是个梳状滤波器，频响有等间隔的峰谷，峰谷位置由延迟R决定，深度由α控制。在MATLAB里实现非常直接：
[x_audio, Fs] = audioread('music.wav'); if size(x_audio, 2) > 1 x_audio = mean(x_audio, 2); % 双声道转单声道 end R = 0.1 * Fs; % 延迟0.1秒 alpha = 0.6; % 回声衰减系数 N = length(x_audio); y = zeros(N + R, 1); y(1:N) = x_audio; y(1+R:N+R) = y(1+R:N+R) + alpha * x_audio; soundsc(y(1:N+R), Fs);这段实现的思路是：把原信号放在输出数组的前N个位置，然后把衰减后的原信号叠加到偏移R个采样点的位置。注意如果音频本身是双声道的，要先用mean转成单声道，否则filter操作会在声道维度上出错。α的取值是关键——大于0.7时回声过重，听起来像山洞里说话；小于0.3时几乎感觉不到回声效果，自己试听调整即可。
幅频特性可以这样绘制：
[h, w] = freqz([1, zeros(1, R-1), alpha], 1, 1024, Fs); figure; plot(w, 20*log10(abs(h))); xlabel('频率/Hz'); ylabel('幅度/dB'); title('单回声滤波器幅频特性');这里freqz(B, A, N, Fs)返回的频率单位是Hz而不是rad/sample，显示更直观。梳状特性的峰谷间距是Fs/R Hz，R越大峰谷越密，这个参数直接决定了回声的音色效果。
4.2 多重回声：IIR递归结构与衰减因子的平衡
多重回声用IIR滤波器实现，传递函数为：
H(z) = 1 / (1 - α·z^(-R))
对应的差分方程是 y[n] = x[n] + α·y[n-R]
注意这里是反馈结构，当前输出是当前输入加上R个采样点之前的输出乘以α。这个结构的记忆长度是无限长的，只要α绝对值小于1就稳定。实现代码：
R = 0.15 * Fs; % 延迟0.15秒 alpha = 0.5; % 反馈系数，必须小于1 b = 1; a = zeros(1, R+1); a(1) = 1; a(R+1) = -alpha; y_multi = filter(b, a, x_audio); soundsc(y_multi, Fs);这段代码里，a向量有R+1个元素，a(1)=1、a(R+1)=-α，中间全是零。filter(b,a,x)用这个稀疏系数向量实现反馈延迟结构。运行后能听到间隔约0.15秒、幅度按指数衰减的多重回声，类似于在山谷里的回声效果。
注意一个容易被忽略的问题：R值比较大时a向量很长（例如Fs=44100、R=0.15秒时R≈6615），filter的运行速度会明显变慢。优化做法是把大延迟拆成多个小延迟级联。
4.3 全通混响器的结构改进与听感验证
任务书最精彩的部分在3.3节——多重回声滤波器的幅频特性不是常数，会产生“染色”效应，听起来不自然。全通混响器的传递函数为：
H(z) = (z^(-R) - g) / (1 - g·z^(-R))
|g| < 1
这个结构的幅频特性恒为1，对所有频率的增益相同，因此不会改变音色，只会增加回声密度。实现：
R = 0.05 * Fs; % 混响器延迟通常较短 g = 0.7; % 全通系数 b_ap = [-g, zeros(1, R-1), 1]; a_ap = [1, zeros(1, R-1), -g]; y_ap = filter(b_ap, a_ap, x_audio); soundsc(y_ap, Fs);b和a的构造逻辑：b向量是z^(-R)系数减g；a向量是1减g·z^(-R)。全通滤波后每个输入脉冲会产生无限多个衰减回声，但因为每个频率的延迟时间不同，听感比单纯多重回声更自然。
课程设计要求给出一组级联结构：4个多重回声滤波器和2个全通混响器串联。级联顺序对听感影响很大，我通常的做法是多重回声在前（建立节奏感），全通在后（增加密度，模拟房间反射）。每个滤波器级联时都要重新归一化输出幅度，否则经过6级滤波器后信号会溢出削波。
5. GUI集成：把三个实验装进一个可演示的程序
5.1 菜单结构规划与回调函数分工
课程设计报告要求GUI集成打包，这本质上是把前面三段程序用菜单控件串联起来。用GUIDE新建空白界面，通过菜单编辑器创建三个一级菜单项“第一大题”、“第二大题”、“第三大题”，每个菜单项挂一个回调函数。整个GUI的组织逻辑很简单：菜单项的回调里调用对应的绘图脚本，绘图脚本的输出figures独立弹出，GUI本身只作为入口。
GUIDE自动生成的模板里，菜单回调函数名是Untitled_1_Callback这种形式，不方便管理。建议手动改名，例如：
function menu_dft_analysis_Callback(hObject, eventdata, handles) % hObject handle to menu_dft_analysis % 回调内容：调用DFT分析脚本 run_dft_analysis();注意要同步修改test6.m中gui_State结构里的gui_Callback字段设置，否则菜单点击没有响应。
5.2 音频播放与路径问题的处理
音频处理部分涉及文件读取和播放，有几个容易踩的坑。audioread的输入路径必须是当前工作目录下或绝对路径，建议用fileparts(mfilename('fullpath'))获取脚本所在目录，再拼接音频文件路径：
scriptPath = fileparts(mfilename('fullpath')); audioPath = fullfile(scriptPath, 'music.wav'); [x_audio, Fs] = audioread(audioPath);soundsc播放后要加pause，否则MATLAB会立即结束播放导致听不到声音。如果需要在GUI中播放，用audioplayer对象更稳定：
player = audioplayer(y, Fs); playblocking(player); % 阻塞直到播放完成5.3 数据复现与参数对比技巧
课程设计报告和答辩展示时，同一个参数下可能有多个图表需要对比。我建议写一个参数扫描函数，用inputdlg弹窗接收参数值而不是硬编码，这样演示时可以直接证明“参数改变、现象随之改变”的逻辑关系：
prompt = {'采样点数N:', '频率间隔Δf:'}; dlgtitle = 'DFT参数设置'; dims = [1 35]; definput = {'128', '1/64'}; answer = inputdlg(prompt, dlgtitle, dims, definput); N = str2double(answer{1}); df = str2double(answer{2});一个完整可演示的课程设计系统，最后应该包含：频谱分析模块（参数可调）、滤波模块（显示滤波前后对比）、音频效果器模块（含试听按钮）。只要这三块跑通，答辩时按顺序演示，基本能覆盖任务书的全部考核点。
本文还有配套的精品资源，点击获取
