2026/9/11 19:10:48

GMM说话人识别原理与MATLAB实现:从MFCC特征到EM训练

GMM说话人识别原理与MATLAB实现:从MFCC特征到EM训练 简介一套基于高斯混合模型GMM的说话人身份识别Matlab仿真资源面向语音信号处理、声纹识别方向的学生、研究人员以及算法初学者也适合作为课程设计或毕业设计的基础代码。资源覆盖从语音特征提取、GMM模型训练到说话人身份匹配的完整流程利用多种脚本完成类似MFCC的特征参数计算通过GMM初始化与期望最大化迭代完成模型参数估计再结合后验概率计算进行说话人判别整体算法链条清晰便于边读代码边理解识别原理。用户只需在MATLAB 2021a或更高版本中运行训练和识别主程序即可复现结果也可以修改高斯分量个数、更换实验语音或替换数据文件做进一步对比实验深入观察GMM参数对识别率的影响。压缩包共16个文件包含12个.m功能脚本、3个.mat数据文件训练数据、说话人模型和测试数据和1个avi操作演示视频整体大小仅2.85MB结构清晰、便于按需查阅。附带的录屏详细演示了运行环境配置和注意事项能有效降低上手门槛即使是初学者也能快速完成实验。目前已有330人学习使用内容紧凑实用可在其基础上二次开发。1. GMM说话人识别一段30秒录音如何锁定“你是谁”一段30秒的录音换个麦克风、隔两天重录频谱形态会有明显漂移但人耳依然能判断“还是这个人”。说话人识别要复现的正是这种能力高斯混合模型GMM的做法是给每个注册说话人构建一张概率密度地图用若干个高斯分量加权叠加来覆盖这个人语音特征的散布区域。测试时把一段语音的特征序列分别送进所有模型得分最高者即为识别结果。这个思路和模板匹配最大的区别在于GMM不追求“最像某一次录音”而是统计“这个人通常怎么发声”所以对麦克风差异、情绪波动、语速起伏都有天然容忍度。它适合几十人以内、30秒级注册语音的身份确认场景也是研究生做语音方向仿真时最先接触的基线系统。下面这套MATLAB流程覆盖特征、训练、识别、调参四段标题里提到的操作演示视频按同流程录制对照视频逐步跑通即可。2. GMM为什么是说话人识别的经典模型从混合分布到EM迭代2.1 声学特征天生“多峰”单个高斯分布必然欠拟合一段语音里包含元音、辅音、鼻音和停顿不同音素的共振峰位置相差很大。如果把MFCC特征投影到二维平面观察会看到多个密集区域而不是一个椭圆点云。单个高斯分布去拟合这些区域时均值会被推到几个峰之间的低概率区结果在真实特征频繁出现的位置模型给出的概率反而偏低。GMM用M个加权高斯分量覆盖这些区域每个分量大致对应一类发音状态或共振峰配置权重表示这类状态在一段话里出现的比例。M是唯一关键模型容量参数。M太小多个声学状态被迫挤进同一个分量模型像“近视眼”M太大某些分量只学到一个说话人的两三个孤立帧测试时分数方差激增。常见做法是在8到64之间做扫描扫描手法在第5章展开。相比深度模型GMM至今仍是说话人识别的可靠基线计算可控、样本需求低、结果可解释。每个分量的均值、方差、权重都有明确物理含义出问题时可以直接检查某个分量是否塌缩这是黑盒模型不容易做到的。2.2 MFCC特征转换把语音切片变成GMM的输入向量GMM本身不关心时域波形它消费的是特征向量序列。最常见的特征是MFCC。提取流程是预加重、分帧、加窗、FFT、Mel三角滤波器组求能量、取对数、DCT得倒谱系数。预加重系数0.97用于削弱低频能量、让高频共振峰更容易被建模帧长25ms、帧移10ms是语音处理的常用起点兼顾时域分辨率和频域分辨率。function mfcc computeMFCC(x, fs, numCep) % x: 单声道语音列向量, fs: 采样率, numCep: 倒谱系数个数 frameLen round(0.025 * fs); % 25ms帧长 frameShift round(0.010 * fs); % 10ms帧移 x filter([1 -0.97], 1, x); % 预加重, 提升高频 numFrames floor((length(x) - frameLen) / frameShift) 1; frames zeros(frameLen, numFrames); for n 1:numFrames st (n - 1) * frameShift 1; frames(:, n) x(st:st frameLen - 1); end frames frames .* hamming(frameLen, periodic); % 加窗抑制频谱泄漏 nfft 512; % FFT点数 nfilt 24; % Mel滤波器个数 melFilt melfb(nfilt, nfft, fs); % nfilt x nfft 三角滤波器组 spec abs(fft(frames, nfft, 1)); % 各帧幅度谱 logMel log(spec * melFilt eps); % 帧数 x nfilt mfcc dct(logMel); % 每行对应一帧 mfcc mfcc(:, 1:numCep); % 截取前numCep维 end代码里melfb是Mel滤波器组构造函数本质是在Mel刻度上等间距放置nfilt个三角窗自己写只需十几行。hamming和dct属于Signal Processing Toolbox完整版MATLAB默认包含精简版没有时hamming可用0.54-0.46*cos(2*pi*(0:n-1)/(n-1))替代dct也有十几行的手写实现。返回的mfcc每一行对应一帧语音的特征向量行数等于帧数。需要说明的是这里没有加能量维和一阶二阶差分。纯静态MFCC已经足够反映说话人音色差分信息对噪声更敏感。演示视频里建议先跑通13维静态版本再视识别率决定要不要扩到39维——39维特征维度翻三倍EM迭代和打分耗时也翻三倍对验证场景未必划算。2.3 EM迭代的核心分工E步算后验、M步更新参数GMM没有解析解要用期望最大化EM算法迭代逼近。E步根据当前参数计算每个特征帧属于每个分量的后验概率M步用这些后验作为权重重新估计均值、方差和混合系数重复直到对数似然增量小于阈值。为防止概率下溢所有计算都在对数域进行。function logPdf computeLogPdf(X, mu, sig2, w) % X: NxD特征矩阵, mu: MxD, sig2: MxD(对角协方差), w: 1xM N size(X, 1); D size(X, 2); M size(mu, 1); logPdf zeros(N, M); for i 1:M diff X - mu(i, :); % NxD logPdf(:, i) log(w(i)) - 0.5 * D * log(2 * pi) ... - 0.5 * sum(log(sig2(i, :)), 2) ... - 0.5 * sum(diff.^2 ./ sig2(i, :), 2); end end function y logsumexp(x, dim) % 数值稳定的 log(sum(exp(x))), 避免小数下溢为0 xmax max(x, [], dim); y xmax log(sum(exp(x - xmax), dim)); endcomputeLogPdf里三个减号分别对应高斯分布的对数归一化系数、各维度方差对数和、马氏距离。logsumexp是EM迭代和识别打分都要复用的工具函数先取最大值再移位避免exp(负几百)直接归零。对角协方差的设计是有意为之——每个维度独立计算方差的倒数M步不需要做矩阵求逆训练帧数只有几万时也能在几十轮内稳定收敛换全矩阵协方差参数个数随维度平方上涨小样本下很容易奇异。E步的职责是“评价现状”在参数固定的前提下算每个点的归属。M步的职责是“按评价结果改进参数”把每个分量收进来的点重新求均值和方差。重复这两步模型就会逐步靠近局部最优。这里提前说一个常见误区EM只保证收敛到局部最优所以初始化依赖kmeans否则16个分量的模型可能收敛出一堆权重趋近于0的空分量。下一章的训练函数直接采用kmeans初始化。3. 训练阶段MATLAB从wav到GMM模型库的完整流程3.1 训练数据目录与预处理约定训练前先把语音整理成固定目录结构data/spk01/、data/spk02/……每个子目录代表一个说话人目录内放该说话人多条wav单人总时长建议不少于30秒。audioread是MATLAB基础函数不需要额外工具箱。读取后对多声道做平均再统一重采样到16kHz这是语音识别领域最常用的采样率能覆盖电话语音和大部分麦克风录音的频谱范围。训练脚本里有一个常被忽略的步骤去掉特征序列的首尾各几帧。wav开头经常是静音段或点击声静音帧特征集中在坐标原点附近会让某个高斯分量专门建模“没说话”的状态挤占本应描述发音的分量。处理办法见下方trainSpeakerGMM第9行的切片。3.2 训练核心kmeans初始化加EM迭代function model trainSpeakerGMM(wavDir, spkId, M) files dir(fullfile(wavDir, *.wav)); feats []; for k 1:numel(files) [x, fs] audioread(fullfile(wavDir, files(k).name)); x mean(x, 2); % 多声道转单声道 x resample(x, 16000, fs); % 统一到16kHz f computeMFCC(x, 16000, 13); % 每行一帧, 13维 if size(f, 1) 40 % 语音太短则整段丢弃 feats [feats; f(10:end-10, :)];% 两端各删9帧避静音 end end [mu, sig2, w] trainGMM(feats, M, 50); model.mu mu; model.sig2 sig2; model.w w; model.spkId spkId; model.M M; endsize(f,1)40的判断保证删去首尾后还有足够帧数。feats累积的是该说话人所有语音的帧集合行数通常上万、列数13后续EM直接在这张大矩阵上迭代。训练的主循环在trainGMM函数里它先做kmeans初始化再进入EMfunction [mu, sig2, w] trainGMM(X, M, maxIter) % X: NxD, M: 分量数, maxIter: 最多迭代轮数 [idx, ~] kmeans(X, M, MaxIter, 50, ... EmptyAction, singleton); % 各簇至少1帧 N size(X, 1); D size(X, 2); mu zeros(M, D); sig2 zeros(M, D); w zeros(1, M); for i 1:M pts X(idx i, :); mu(i, :) mean(pts, 1); sig2(i, :) var(pts, 0, 1) 1e-4; % 方差加下限防奇异 w(i) size(pts, 1) / N; end for iter 1:maxIter logPdf computeLogPdf(X, mu, sig2, w); % NxM llkOld mean(logsumexp(logPdf, 2)); % 上一轮对数似然 gamma exp(logPdf - llkOld); % E步, 每行和为1 Nk sum(gamma, 1); % 每个分量的有效帧数 w Nk / N; % M步: 权重 mu (gamma * X) ./ Nk; % M步: 均值 for i 1:M diff X - mu(i, :); sig2(i, :) (gamma(:, i) * (diff.^2)) / Nk(i) 1e-4; end if iter 1 logPdfNew computeLogPdf(X, mu, sig2, w); llkNew mean(logsumexp(logPdfNew, 2)); if abs(llkNew - llkOld) 1e-3 break; end end end endE步直接取exp(logPdf - llkOld)得到后验概率矩阵因为每行的llkOld是一个标量均值所有分量共享同一个减数这样计算省掉一次logsumexp。M步里mu那行是矩阵运算gamma是MxN乘以X后除以Nk按后验加权求均值。sig2的更新用gamma第i列对diff.^2做加权平均后加上1e-4下限这个下限能避免某维方差退化成0——仿真中经常出现的“训练发散”有一半根因就出在这个数值边界上。kmeans来自Statistics and Machine Learning Toolbox若没装该工具箱可以手写20行K-means替代装有完整版时也可以后续用fitgmdist(X,M)和这段自写EM对比收敛代数。训练完成后把5个model存入cell数组作为模型库注意测试文件绝不能出现在训练脚本扫描的路径里spkNames {spk01, spk02, spk03, spk04, spk05}; models cell(numel(spkNames), 1); for i 1:numel(spkNames) models{i} trainSpeakerGMM(fullfile(baseDir, spkNames{i}), i, 16); end save(gmmModels.mat, models);16个分量的选择依据在第5章的扫描曲线上。这里先记住经验训练语音总量越多M可以适当加大单人只有10秒语音时M8更稳M32大概率出现空分量。提示训练集和测试集一旦混用识别率会虚高到没有参考价值。按说话人目录整理文件时训练目录与测试目录各放一套wav测试文件绝不能出现在trainSpeakerGMM扫描的路径里。3.3 训练阶段的三个必查点第一个必查点是特征矩阵里有没有NaN或Inf。wav尾部常有一段电平极低的采集噪声方差若除到0会直接把整列方差变成Inf。排查方式是在trainSpeakerGMM里临时加一行assert(~any(isnan(feats(:))))触发即停。第二个必查点是空分量。训练结束后打印w如果某个分量权重小于0.005它基本没有贡献。可以降M也可以给kmeans加Replicates参数多试几次随机初始中心代价是训练时间翻倍。第三个必查点是训练轮数。50轮对16分量、几万帧数据通常足够如果50轮后收敛判据还没触发先别急着加迭代上限回看llkNew是否在小幅波动——EM在接近峰值时对数似然曲线会变平1e-3的相对门限可能需要更密的帧数才能触发这时候改成看连续5轮差值平均更稳定。4. 识别阶段极大似然判定与识别率统计4.1 对每个模型算平均对数似然识别流程比训练简单把测试语音打散成特征帧每帧对每个模型算computeLogPdf用logsumexp得到该帧的log似然再对帧数做平均。平均是为了消除录音时长差异带来的偏差——2秒的测试语音帧数少、5秒的帧数多不平均时长的录音天然占优平均后分数才具有跨录音比较的意义。% 测试单条语音: testFile为wav完整路径 [xt, fst] audioread(testFile); xt mean(xt, 2); % 单声道 xt resample(xt, 16000, fst); ft computeMFCC(xt, 16000, 13); % 与训练完全同参数 score zeros(numel(models), 1); for i 1:numel(models) lp computeLogPdf(ft, models{i}.mu, models{i}.sig2, models{i}.w); score(i) mean(logsumexp(lp, 2)); % 平均对数似然 end [bestScore, pred] max(score);score的数值单位是“每帧的平均对数似然”通常为负。遍历模型库后pred就是识别出的说话人编号。最后两行的max隐含了各说话人先验相同的假设如果事先知道某人说话更频繁可以在score上加上log(先验概率)这是贝叶斯最小错误率准则的直接应用。4.2 阈值判定不认识的人如何拒绝严格说上述流程只解决“从注册者中选一个”不解决“这个人根本没注册过”。真实门禁系统需要加一道阈值完成拒识if bestScore threshold pred 0; % 0表示拒绝/未知说话人 endthreshold怎么定常见做法是在独立验证集上跑一遍画出“正确接受注册者”与“错误接受冒认者”的DET曲线取等错误率EER对应的分数作为阈值。阈值和M、训练语音时长都有关调参扫描时阈值要和M一起记录在同一张实验表里改动M后必须重新测阈值。对仿真验证来说能画出曲线并解释形状即可不必追求EER特别低。4.3 识别率统计与混淆矩阵把测试集完整遍历一遍统计正确率和混淆矩阵total 0; correct 0; confMat zeros(numel(models), numel(models)); for spk 1:numel(models) testDir fullfile(testBaseDir, spkNames{spk}); files dir(fullfile(testDir, *.wav)); for k 1:numel(files) total total 1; pred recognizeFile(fullfile(testDir, files(k).name), models); correct correct (pred spk); if pred 0 confMat(spk, pred) confMat(spk, pred) 1; end end end fprintf(识别率 %.2f%%\n, 100 * correct / total);这里recognizeFile是把4.1和4.2合并封装的函数输入wav路径和模型库返回pred实现细节就是第4.1节那段代码加上阈值判断。装有Statistics and Machine Learning Toolbox时可以用confusionchart(confMat)直接出图没有也不影响把confMat打印出来看对角线占比即可。对角线比例越高说明区分度越好某一行整体偏向另一列说明两个说话人音色接近需要增加该人训练语音或调大M。4.4 演示视频的对照检查顺序操作演示视频一般按下面顺序录制按此顺序对照自己机器上的运行结果可快速定位问题先在命令窗口调用computeMFCC验证特征维度和帧数是否符合预期再画出前两个维度的散点图观察有多少个明显的簇接着运行trainGMM观察每轮llk是否单调上升然后执行识别脚本打印测试语音的得分向量最后统计5个说话人的识别率。如果散点图看不出簇状结构多半是MFCC的滤波器组或预加重写错不要急着调GMM参数。5. GMM参数边界与仿真排错分量数、收敛性和验证技巧5.1 用分量数扫描曲线代替拍脑袋M从4、8、16、32、64扫描每个M训练5个说话人模型用同一段测试集统计识别率脚本只需在4.3的统计外面包一层for循环。扫描结果通常是抛物线M过小声学模式被压缩进少数几个分量区分度不足M过大训练数据被切得太碎测试特征落入分量边界区域的概率升高。中低档数据量下M16或32往往就是拐点。扫描时把训练耗时一并记下M翻倍训练时间大约也翻倍性能提升却只有零点几个百分点时就该退回到小一档。5.2 训练发散与NaN的定位清单“训练发散”在GMM上下文里不是梯度爆炸而是三种可识别的现象对数似然在一轮内骤降、部分方差变成Inf、权重大量趋零。第一类多半是computeLogPdf里的维度或sum作用方向写错打印logPdf矩阵看某一列是否全为极小的负值即可第二类出现在特征含NaN或方差下限加错位置时第三类则是kmeans初始化散不开先加Replicates再考虑降M。调试阶段在脚本开头加rng(42)固定随机种子让症状可复现修完问题再移除。5.3 可复现验证的实验记录模板做参数实验时按以下模板记录避免调了半天说不清哪组参数有效表内数字仅为示意结构可直接复用实验编号M每人训练语音时长阈值识别率备注1830s-45.284.3%欠拟合明显21630s-47.892.1%推荐参数33230s-49.691.6%与16差别不大41615s-46.085.5%训练时长敏感这张表同时回答“M选多少”和“加训练语音值不值”两个问题。阈值一栏的变化提示我们改动M后必须重测阈值否则识别率数字不可比。做视频演示时也建议按这张表逐行跑最后指着表格中的某一行作为仿真结论比笼统说“识别率90多”更有说服力。本文还有配套的精品资源点击获取