简介：这是一份面向电子信息、通信工程等专业毕业设计的MATLAB语音信号处理系统设计资料，系统讲解语音信号从采集、加噪到滤波还原的完整流程，并融入图形用户界面（GUI）开发，适合课程设计或毕业设计参考。压缩包内为单个PDF文档，共1.45MB，正文包含关键程序代码和仿真波形图，便于对照实现与验证。文档从声卡采集语音信号出发，利用FFT完成时频域转换与频谱分析，模拟添加干扰噪声后，重点讲解巴特沃斯、切比雪夫滤波器的设计原理，以及双线性变换法、脉冲响应不变法在数字滤波器实现中的应用；滤波后通过重采样、回放对比和参数提取评估处理效果，并给出GUI界面加载波形、调节参数和显示频谱的构建思路。整体内容既覆盖理论推导，又给出可运行的MATLAB实现细节，可作为毕业设计说明书、答辩PPT素材和实验复现依据。目前已有2000人学习，是上手语音信号处理与滤波器设计的实用资料。
1. 语音信号处理课设资源：从 wavread 到 GUI 的一条完整链路
做课设最怕的不是不会写代码，而是拿到一套源码却不知道每段在干什么，改一个参数就全线崩盘。这份基于 MATLAB 的语音信号处理系统设计，把语音采集、FFT 频谱分析、加噪、数字滤波、重采样、短时参数提取和 GUI 图形界面整条串成了一条可运行的链路。它解决的核心问题是：如何用一段 wav 文件在 MATLAB 里走完从时域到频域再到滤波回放的全流程，同时给出能直接交差的仿真图。适合正在做相关课设或毕业设计的学生，也适合想快速上手 MATLAB 语音分析基础的从业者。程序加仿真图的结构很完整，拿到手能直接跑，关键是跑完之后要能看懂参数怎么改——这正是下面要拆的部分。
2. FFT 频谱分析：先看懂原始语音在频域长什么样
2.1 读入文件与 FFT 变换：这几行代码就是整个系统的基础
语音信号处理的第一步是拿到数据。课程设计里最常见的方式是准备一段 wav 文件放在当前目录，然后直接读入。原始资源里的核心代码是这样的：
[x, fs, nbits] = wavread('mei.wav'); % 读取语音文件，x为采样数据，fs为采样率，nbits为量化位数 sound(x, fs, nbits); % 播放原始语音 N = length(x); % 获取采样点总数 n = 0:N-1; % 样本序号向量 w = 2*n*pi/N; % 数字角频率，范围 0 ~ 2*pi y1 = fft(x); % 对原始信号做 FFT 变换 subplot(2,1,1); plot(n, x); % 时域波形 title('原始语音信号时域图'); xlabel('时间 t'); ylabel('幅值'); subplot(2,1,2); plot(w/pi, abs(y1)); % 频谱幅度，横轴归一化到 0~2 title('原始语音信号频谱'); xlabel('频率 Hz'); ylabel('幅度');逻辑说明：这段代码把语音文件读成向量x，然后直接对该向量做 FFT。FFT 的结果是复数，取模abs(y1)得到幅度谱，横轴用w/pi做归一化，画出来的频谱范围是 0 到 2。其中w = 2*n*pi/N是 MATLAB 课设里很经典的写法，它把 0 到 N-1 的样本序号映射到 0 到 2π 的数字角频率范围。
参数说明：fs是采样率，决定了频谱图的真实频率范围。比如fs = 16000表示每秒采样 16000 个点，奈奎斯特频率是 8000 Hz。nbits是量化位数，常见为 16。N是采样点数，FFT 的分辨率等于fs/N，也就是说 N 越大，频谱上能分辨的频率间隔越小。如果后面要分析特定频段的噪声，这个分辨率值要先算清楚。
2.2 频谱图横轴换算：w/pi 与真实频率 Hz 的对应关系
上面代码里横轴是w/pi，取值范围 0~2，但这个刻度不是 Hz。实际画频谱图时，很多初学者直接把abs(fft(x))的数组画出来，横轴从 0 到 N，这就完全错了。
常见的换算方式有两种。第一种是用w/pi * fs/2把归一化频率转成真实频率，第二种是直接用f = (0:N-1) * fs / N生成频率轴。这里给出一段更实用的画频谱的写法：
f = (0:N-1) * fs / N; % 真实频率轴，单位 Hz mag = abs(fft(x)) / N * 2; % 幅度归一化，单边谱 plot(f, mag); xlabel('频率 (Hz)'); ylabel('幅度');逻辑说明：f = (0:N-1) * fs / N是频率轴的经典构造方式，每个点的频率间隔是fs/N。幅度除以 N 再乘以 2，是为了把双边谱折叠成单边谱，让纵轴数值更符合实际信号幅度。如果你只是想看噪声集中在哪里，不做幅度归一化问题不大；如果想对比滤波前后的幅度变化，归一化就很有必要。
参数说明：当 N 是 2 的幂次时，FFT 计算最快。如果N不是 2 的幂，MATLAB 会自动补零到合适的长度，但这时频谱分辨率会变。所以如果你发现某次 FFT 之后频率轴对不上，先检查是不是 N 的奇偶性出了问题。
3. 加噪与数字滤波器设计：巴特沃斯链路与参数表
3.1 加干扰噪声：固定频率余弦叠加的写法与频率换算
原始语音读进来之后，下一步是人为加入干扰噪声。课程设计里加的是单一频率的余弦噪声，代码是这样的：
z = 0.02 * cos(10 * n); % 构造噪声，幅度0.02，角频率10 y = x + z; % 带噪信号 sound(y, fs, nbits); % 播放带噪语音逻辑说明：z = 0.02*cos(10*n)生成的是一个幅度 0.02、数字角频率 10 的余弦序列。注意这里的n是样本序号（0,1,2,...），不是时间 t，所以 10 的单位是弧度/样本，不是 Hz。这个噪声加到语音上之后，FFT 频谱会在对应频点出现一根明显的谱线，滤波器滤掉的就是这根线。
参数说明：这里最容易翻车的地方是 10 这个数字。如果采样率 fs = 44100，那么实际噪声频率是10*fs/(2*pi)，约 70 Hz，这是个非常低的频率，人耳几乎听不出区别，滤波效果也不明显。课设里要演示“能听见的降噪效果”，建议把噪声改成可听频段，比如 3000 Hz 左右的单频干扰。常见做法是把代码改成：
f0 = 3000; % 噪声频率，单位 Hz z = 0.02 * cos(2*pi*f0*n/fs); % 换算到数字域 y = x + z;这样改了之后，滤波前后的回放差异会非常明显，答辩时演示效果也好得多。
3.2 滤波器设计链路：从模拟原型到数字滤波器的四步走
资源里滤波器设计走的是完整链路：先设计模拟巴特沃斯低通，再用脉冲响应不变法转成数字滤波器。滤波器参数是：通带截止频率 0.6π，阻带起始频率 0.75π，通带最大衰减 1 dB，阻带最小衰减 30 dB。核心代码分三个函数：
function [b, a] = afd_butt(Wp, Ws, Rp, As) % 计算巴特沃斯滤波器阶数 N = ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1))) / (2*log10(Wp/Ws))); fprintf('\n Butterworth Filter Order=%2.0f\n', N); OmegaC = Wp / ((10^(Rp/10)-1)^(1/(2*N))); % 3dB截止频率 [b, a] = u_buttap(N, OmegaC); % 求模拟滤波器系数 endfunction [b, a] = imp_invr(c, d, T) % 脉冲响应不变法：模拟滤波器转数字滤波器 [R, p, k] = residue(c, d); % 部分分式展开 p = exp(p * T); % 极点映射 z = e^(sT) [b, a] = residuez(R, p, k); % 转成多项式之比 b = real(b'); % 取实部得到数字滤波器系数 a = real(a'); endfunction [db, mag, pha, w] = freqz_m(b, a) [H, w] = freqz(b, a, 1000, 'whole'); % 计算0~2*pi频率响应 H = (H(1:501))'; % 取前半段 w = (w(1:501))'; mag = abs(H); % 幅度响应 db = 20*log10((mag + eps) / max(mag)); % 增益，单位dB pha = angle(H); % 相位响应 end主程序把这些函数串起来：
wp = 0.6*pi; ws = 0.75*pi; Rp = 1; As = 30; T = 1; OmegaP = wp/T; OmegaS = ws/T; % 模拟域频率 [cs, ds] = afd_butt(OmegaP, OmegaS, Rp, As); [b, a] = imp_invr(cs, ds, T); % 转数字滤波器 [db, mag, pha, w] = freqz_m(b, a); figure(1); subplot(2,1,1); plot(w/pi, mag); % 幅度响应 title('digital filter Magnitude Response'); axis([0 1 0 1.1]); subplot(2,1,2); plot(w/pi, db); % 增益响应 title('digital filter Magnitude in DB'); axis([0 1 -40 5]);逻辑说明：afd_butt里的阶数公式是巴特沃斯滤波器设计的标准公式，N = ceil(...)向上取整保证满足指标。imp_invr用部分分式展开做模拟到数字的极点映射，residuez再把极点形式转回系数形式。这条链路的价值在于它完整呈现了教材里“模拟原型→数字化”的每一步，适合写进设计报告的论证部分。
参数说明：wp = 0.6*pi和ws = 0.75*pi是数字截止频率，单位是弧度/样本。换算成真实频率是wp*fs/(2*pi)。比如 fs = 16000 时，通带截止约 4800 Hz，阻带起始约 6000 Hz。Rp = 1表示通带内最大衰减 1 dB，As = 30表示阻带最小衰减 30 dB。T 是采样间隔，通常取 1 配合归一化频率。
这里有个更快的替代方案：如果只是要一个能用的低通滤波器，直接调用 MATLAB 内置函数就行：
Wn = 0.4; % 归一化截止频率，0~1，对应 fs/2 [b, a] = butter(8, Wn); % 8阶巴特沃斯低通但课设报告里如果只写这一行，设计依据很难展开。资源里那套完整链路更适合需要写推导过程、画设计步骤图的场景。
3.3 filter 函数滤波与降噪效果的前后对比
滤波器系数设计好之后，用filter函数实现对带噪信号的滤波：
y1 = filter(b, a, y); % 用设计好的滤波器滤除噪声 y2 = fft(y1); % 对滤波后信号做 FFT sound(y1, fs, nbits); % 回放滤波后语音 subplot(2,1,1); plot(n, y1); % 滤波后时域波形 title('滤波后信号时域图'); xlabel('时间 t'); ylabel('幅值'); subplot(2,1,2); plot(w/pi, abs(y2)); % 滤波后频谱 title('滤波后信号频谱'); xlabel('频率 Hz'); ylabel('幅度');逻辑说明：filter(b, a, y)是 IIR 滤波器的标准调用方式，b是前馈系数，a是反馈系数。滤波后你会看到时域波形变平滑，频谱上那根突出的噪声谱线被压下去。对比滤波前的频谱图，能直观看出设计参数是否达标。
参数说明：filter函数默认初始状态为零。如果信号很长，建议关注滤波后的起始段是否有瞬态效应。另外filter不改变数据长度，和卷积conv不同，这在对比频谱时不用对齐长度，比较方便。
4. 重采样与短时参数提取：把语音信号“切片”来看
4.1 resample 重采样：8kHz 与 4kHz 两档降采样
语音信号处理里重采样是常见操作。资源里的实现是：
[x, fs, nbits] = wavread('mei.wav'); x1 = resample(x, 1, 2); % 采样率降一半，比如 16kHz -> 8kHz sound(x1, fs/2, nbits); % 播放 8kHz 采样语音 x2 = resample(x, 1, 4); % 再降一半，8kHz -> 4kHz sound(x2, fs/4, nbits); % 播放 4kHz 采样语音 y1 = fft(x1); y2 = fft(x2); figure(1); subplot(2,1,1); plot(x1); title('8KHZ重采样语音信号时域图'); subplot(2,1,2); plot(abs(y1)); title('8KHZ重采样语音信号频谱'); figure(2); subplot(2,1,1); plot(x2); title('4KHZ重采样语音信号时域图'); subplot(2,1,2); plot(abs(y2)); title('4KHZ重采样语音信号频谱');逻辑说明：resample(x, p, q)的作用是把采样率变成原来的 p/q 倍。resample(x,1,2)是抽取操作，每隔一个点取一个，采样率减半。注意抽取前滤波器会自动做抗混叠处理，所以比直接写x(1:2:end)更稳妥。
参数说明：采样率降低后，信号最高可用频率也跟着降低。8kHz 采样时奈奎斯特频率是 4kHz，4kHz 采样时只有 2kHz。回放时能明显听出声音变闷——高频信息丢了。这是演示采样率对音质影响的好素材。
这里有个容易被疏忽的细节：resample之后新信号长度不是简单的N/2。MATLAB 会根据滤波器的群延迟调整输出长度，所以后面画图时最好用length(x1)动态获取长度，别直接用N/2硬编码。
4.2 短时能量、平均幅度与过零率：区分清音浊音的三个抓手
语音信号是非平稳的，但在一段短时间（通常 20~30 ms）内可以看成平稳的。所以要做分帧处理，再在每一帧上计算短时能量、短时平均幅度和短时过零率。资源里给了三个核心公式的实现基础：短时能量En = sum(xn.^2)，短时平均幅度Mn = sum(abs(xn))，短时过零率Zn = 1/2 * sum(abs(sign(xn) - sign(xn-1)))。
一个完整的短时参数计算代码可以这样写：
function [energy, amag, zcr] = stparam(x, win_len, step) % 短时能量、平均幅度、过零率计算 % x: 输入语音信号, win_len: 帧长, step: 帧移 N = length(x); nf = fix((N - win_len) / step) + 1; % 帧数 energy = zeros(1, nf); amag = zeros(1, nf); zcr = zeros(1, nf); for i = 1:nf start = (i-1)*step + 1; seg = x(start : start+win_len-1); % 取一帧 energy(i) = sum(seg .^ 2); % 短时能量 amag(i) = sum(abs(seg)); % 短时平均幅度 zcr(i) = sum(abs(sign(seg(2:end)) - sign(seg(1:end-1)))) / 2; end end逻辑说明：分帧的起点是(i-1)*step + 1，帧长win_len通常取 20 ms 对应的采样点数，帧移step取 10 ms。比如 fs = 16000 时，win_len = 320，step = 160。过零率的实现用了sign函数：相邻采样点符号不同就计 1，总计除以 2 得到过零次数。注意sign函数在值为 0 时返回 0，实际语音信号很少正好为零，影响不大。
参数说明：短时能量的作用很直接——浊音段幅度大，能量高；清音段幅度小，能量低。短时过零率正好反过来：浊音低频成分多，过零率低；清音高频成分多，过零率高。这两个参数配合，可以在信噪比不差的情况下做有声/无声判决，也可以粗略区分清浊音。如果你改了帧长，阈值也要跟着重新调，没有一劳永逸的参数。
4.3 短时自相关函数：估计基音周期的基本操作
自相关函数的作用是发现信号的周期性。语音里浊音有周期成分，自相关会在基音周期处出现峰值，所以可以用来估计基音频率。资源里的自相关定义是Rn(k) = sum(xn(m) * xn(m+k))，也就是把一帧信号和自己错位 k 个点做内积。
对应代码：
function [r] = autocorr_seg(seg, maxlag) % 计算一帧语音的短时自相关 % seg: 一帧语音, maxlag: 最大延迟点数 N = length(seg); r = zeros(1, maxlag+1); for k = 0:maxlag r(k+1) = sum(seg(1:N-k) .* seg(1+k:N)); % 错位相乘累加 end r = r / r(1); % 归一化，r(0)=1 end逻辑说明：核心语句seg(1:N-k) .* seg(1+k:N)实现的是错位 k 个点的內积。当 k 等于基音周期时，波形对齐程度高，乘积和出现峰值。归一化用r(1)做分母，让 r(0)=1，后面找峰时只看相对值。
参数说明：maxlag的取值要覆盖预期的基音周期范围。正常人说话基频约 80~400 Hz，fs = 16000 时周期对应 40~200 个采样点，所以maxlag取 200 左右比较合适。如果取太小，周期长的低音男声找不到峰；取太大，计算量变大且容易出现假峰。
5. 复现避坑：从 wavread 失效到滤波不出效果的五个典型问题
5.1 老函数报错：新版 MATLAB 里 wavread 已经被替代
现象：直接跑资源里的代码，第一行wavread就报错“未定义函数或变量”，程序完全跑不起来。
原因：wavread是旧版本（R2016a 之前）的音频读取函数，新版 MATLAB 用audioread取代了它。同理，sound(x, fs, nbits)在新版本里nbits参数也被忽略。
解决：做两处替换。wavread换成audioread，返回值变成[x, fs]，没有 nbits 了，播放直接写sound(x, fs)。改一行就能跑通：
[x, fs] = audioread('mei.wav'); sound(x, fs);5.2 重采样后长度取整：N/2 与 N/4 索引溢出的边界
现象：跑重采样那段代码时，画图报错“索引超出数组边界”，但单独执行resample不报错。
原因：resample(x,1,2)的输出长度不是严格等于N/2。MATLAB 内部重采样滤波器有群延迟，输出长度通常是ceil(N*p/q)附近的值。原代码里直接写n1 = 0:N/2-1是假设整除，一旦 N 是奇数或者长度对不上，n1 就越界了。
解决：用length(x1)动态生成索引，别用 N 做算术推导：
N1 = length(x1); n1 = 0:N1-1;5.3 加噪频率没有换算：10 不是频率，是样本序号
现象：照着资源加了z = 0.02*cos(10*n)，滤波前后听着没什么区别，频谱上看那根噪声谱线也压不掉。
原因：cos(10*n)里的 10 是数字角频率，单位是弧度/样本，不是 Hz。当采样率较高时，这个 10 对应的实际频率非常低，可能只有几十 Hz，加到语音上几乎感觉不到。更关键的是，如果你的滤波器截止频率远高于这个噪声频率，滤波器根本不需要出力，噪声本来就在通带里。
解决：按 3.1 里的写法，把噪声频率换算到可听频段：
f0 = 3000; z = 0.02 * cos(2*pi*f0*n/fs);改了之后滤波效果才会在音频回放里真实可听，频谱对比图也更有说服力。
5.4 脉冲响应不变法的混叠：高通带阻场景不要用这条链路
现象：抄了资源里的滤波器设计链路，改成高通滤波器之后，滤波结果出来一团糟，甚至整个信号都变形了。
原因：imp_invr（脉冲响应不变法）存在频谱混叠问题。它只适合低通和带通滤波器设计，因为这两种滤波器天然没有高频分量，混叠影响小。高通和带阻滤波器有大量高频通过，混叠后数字滤波器频率响应严重失真。
解决：高通、带阻场景改用双线性变换法。用bilinear函数替换imp_invr，或者直接用 MATLAB 的高层设计函数：
% 双线性变换法的高通设计示例 [b, a] = butter(4, 0.4, 'high');顺便说一句，双线性变换本身也有代价——频率轴非线性压缩，截止频率需要做预畸变。资源里没有展开讲这块，但如果你要改高通，这个坑一定会踩到。
5.5 GUI 绘图被覆盖：subplot 与 axes 句柄管理
现象：GUI 界面上点了“开语音”按钮，图确实画出来了，但再点一次，图上会叠加新曲线，或者旧图不消失，界面越来越乱。
原因：GUI 回调里直接调用subplot和plot，每次点击都会在当前坐标区叠画，没有先清除上一次的内容。
解决：每次画图前显式清空坐标轴。常见做法是拿到坐标轴句柄后调用cla，或者直接用plot(handles.axes1, ...)这种指定坐标区的写法：
axes(handles.axes1); % 选中GUI上的axes1 cla; % 清空原内容 plot(n, x); % 重新绘图资源里开语音的代码只写了 pure MATLAB 的subplot用法，这个在脚本里没问题，放进 GUI 回调里就必须补上句柄管理，否则交互体验很差。
6. GUI 整合与回放验证：把滤波前后听出区别的最后一公里
6.1 按钮回调与全局变量组织
资源里的 GUI 部分结构是：界面上放“打开语音”“播放语音”“加噪”“滤波”“重采样”“参数提取”等按钮，每个按钮对应一个回调函数。数据在回调之间传递用的是全局变量，这是课设 GUI 最常见的组织方式：
global x fs N n w; % 在多个回调里共享数据 % “打开语音”按钮回调 function open_push_Callback(hObject, eventdata, handles) global x fs N n w; [x, fs] = audioread('mei.wav'); N = length(x); n = 0:N-1; w = 2*n*pi/N; axes(handles.axes1); cla; plot(n, x); title('原始语音信号时域图'); end % “加噪”按钮回调 function noisy_push_Callback(hObject, eventdata, handles) global x fs N n w; f0 = 3000; z = 0.02 * cos(2*pi*f0*n/fs); y = x + z; axes(handles.axes2); cla; plot(n, y); title('加噪信号时域图'); end % “滤波”按钮回调 function filter_push_Callback(hObject, eventdata, handles) global x fs N n w; [b, a] = butter(8, 0.4); y1 = filter(b, a, x); axes(handles.axes3); cla; plot(n, y1); title('滤波后信号时域图'); sound(y1, fs); end这里的回调函数名用了_Callback后缀，实际 GUI 里按钮 Tag 要对应好。handles.axes1是在 GUI 编辑器里放置坐标轴时自动生成的句柄。用axes(handles.axesN); cla;的组合，就能保证每次点击都刷新而不是叠画。global写法够用，调试方便；如果你想更规范，可以用guidata或handles存数据，但课设阶段不必过度设计。
6.2 验证方法：用声音和波形双通道确认滤波效果
很多人做完滤波，看一眼频谱图觉得差不多就交了。我的习惯是强制走一遍双通道验证：通道一是听觉，通道二是波形与频谱对比。具体做法是先把滤波前后的音频都播一遍，再用一段小脚本把降噪前后的频谱峰值打出来：
% 对比滤波前后噪声频点的幅度 f0 = 3000; f_index = round(f0 * N / fs) + 1; before = abs(fft(y)); after = abs(fft(y1)); fprintf('噪声频点滤波前幅度: %.1f\n', before(f_index)); fprintf('噪声频点滤波后幅度: %.1f\n', after(f_index));逻辑说明：f_index = round(f0*N/fs) + 1是把 3000 Hz 换算到 FFT 数组下标。打印出来滤波前后的幅度，能给出一个数字化的降噪指标，不只是“看着像滤掉了”。如果前后两个数差了一个数量级以上，说明滤波器确实压住了这个频点。
参数说明：这段验证脚本依赖前面把y和y1都留在工作区。如果中途清了工作区，变量丢了脚本就报错。实战里我会把加噪、滤波、验证写进一个完整脚本，一次跑完，避免手动一步步点按钮导致中间变量丢失。从那以后我每次做完滤波器，都强制走一遍“播音频 + 打印频点幅度”这个双通道验证，不再只盯着波形图判断效果。希望帮到你。
本文还有配套的精品资源，点击获取
