尧图建网站 尧图建网站 YAOTU WEB BUILD 免费咨询
ARTICLE DETAIL

资讯详情

深耕网站建设与建站编程的一线实战洞察。

MATLAB实现MFCC特征提取:从原理到实践的完整指南

MATLAB实现MFCC特征提取:从原理到实践的完整指南 1. 项目概述从声音到数字的桥梁做语音处理的朋友对MFCC这个词一定不陌生。它就像是声音世界的“指纹”能把一段连续、复杂的语音波形转换成一串能够代表其本质特性的数字序列。无论是语音识别、说话人识别还是语音情感分析MFCC特征都是最基础、最核心的“原材料”。这个项目就是带你亲手用MATLAB搭建一个完整的MFCC特征提取流程从一段.wav文件开始一步步拆解直到得到那组神秘的数字矩阵。我自己在早期做语音项目时第一个坎儿就是特征提取网上代码片段很多但要么缺斤少两要么原理讲得云里雾里跑通了也不知道为什么。所以这次我们不只给代码更要把每一步背后的“为什么”讲透让你不仅能仿真出来更能理解每一个参数变动对结果的影响真正掌握这门手艺。简单来说这个教程适合三类人一是刚接触语音信号处理的在校学生需要一个能跑通、能理解的入门案例二是从事相关领域开发但想夯实基础、理解细节的工程师三是任何对“机器如何听懂声音”感到好奇的技术爱好者。通过这个案例你将获得一个可复现的MATLAB仿真程序并深刻理解MFCC特征提取的完整链条包括预加重、分帧加窗、短时傅里叶变换、Mel滤波器组、对数能量、离散余弦变换等核心环节。我们不止步于得到结果更要探究结果背后的物理意义和数学逻辑。2. MFCC特征提取的核心原理与设计思路为什么是MFCC这得从人耳的听觉特性说起。人耳对声音频率的感知并不是线性的在低频部分比如100Hz到1000Hz非常敏感分辨率很高而在高频部分比如3000Hz以上敏感度会下降分辨率变低。这种非线性的频率感知尺度就是Mel尺度。MFCCMel-Frequency Cepstral Coefficients梅尔频率倒谱系数的设计初衷就是模仿人耳的听觉特性对语音信号进行一种“仿生”处理从而提取出更符合人类听觉感知、也更有利于机器区分的特征。整个MFCC提取流程是一个精心设计的管道每一步都承担着特定的信号处理任务。其核心设计思路可以概括为将时域的非平稳信号语音转化为相对稳定的短时帧再映射到模仿人耳听觉的Mel频域最后通过倒谱分析剥离出声源激励和声道响应的信息取其中相对稳定的声道响应部分作为特征。这个思路决定了我们技术方案的选择必须采用短时分析方法必须引入Mel滤波器组进行频率弯折必须使用离散余弦变换DCT来近似倒谱分析。在方案选型上我们坚持使用经典且被广泛验证的流程。有些简化方案会跳过预加重或倒谱提升但对于教学和夯实基础而言完整的流程更能体现精髓。我们选择用三角滤波器组来模拟Mel尺度因为它的计算简单效果接近人耳的带通滤波特性。选择DCT而不是真正的逆傅里叶变换来求倒谱是因为DCT具有优秀的能量压缩特性前几个系数就包含了绝大部分信息非常有利于后续的建模和识别。这个方案的优势在于流程清晰、每步可解释性强、结果稳定可靠是工业界和学术界的通用标准。它避免了直接使用FFT频谱的缺点维度高、对频率变化过于敏感也为我们后续可能的操作比如差分系数计算、归一化留下了清晰的接口。2.1 完整流程拆解与环节作用让我们把MFCC提取这个“黑箱”彻底打开看看里面究竟有哪些关键环节以及它们各自扮演什么角色。理解这个你就能在调试时有的放矢。预加重语音信号尤其是浊音其频谱在高频部分的能量通常比低频部分衰减得更快。预加重就是一个一阶高通滤波器目的是提升高频分量使信号的频谱变得平坦便于后续频谱分析。从公式上看它很简单s s(n) - a * s(n-1)其中a通常取0.97左右。这个操作相当于给频谱乘以一个随频率升高而增大的斜坡补偿了高频衰减。分帧语音信号是短时平稳的即在10-30毫秒的时间内其特性如共振峰可以认为是基本不变的。因此我们需要将整个语音流切割成许多这样的小片段称为“帧”。帧长通常取20-30ms如256点8kHz采样或512点16kHz采样。为了保持帧与帧之间的连续性相邻帧之间会有重叠帧移通常为帧长的一半如10-15ms。这个操作是将一个全局非平稳问题转化为一系列局部平稳问题来处理的基础。加窗分帧本质上是给信号加了一个矩形窗这会在频谱中引入很高的旁瓣造成“频谱泄漏”即一个频率的能量会泄漏到其他频率上。为了减少泄漏我们给每一帧乘以一个窗函数如汉明窗。汉明窗在两端平滑地衰减到零可以有效地抑制旁瓣。加窗的代价是主瓣稍微变宽但用频谱分辨率的轻微损失换来泄漏的大幅减少在语音处理中是值得的。快速傅里叶变换对每一帧加窗后的信号进行FFT将其从时域变换到频域得到短时幅度谱。这一步让我们能观察到该帧信号在不同频率上的能量分布。我们通常取FFT点数大于等于帧长并通过补零来增加频率分辨率即频谱的平滑度。Mel滤波器组这是模仿人耳的关键一步。我们在频率轴上0Hz到奈奎斯特频率之间设置一组三角带通滤波器。这些滤波器的中心频率在Mel尺度上是均匀分布的但在线性Hz尺度上则是低频密集、高频稀疏。每个滤波器的作用是计算该频带内的能量总和。通过这组滤波器我们将高维的线性频谱可能512个点映射到低维的Mel频谱通常20-40个维度。这个过程模拟了人耳对不同频带敏感度的差异。取对数能量计算每个Mel滤波器输出的能量并对其取以10为底的对数。这么做有两个主要原因一是人耳对声音强度的感知近似对数关系分贝就是对数单位二是取对数可以将乘性噪声如信道效应转化为加性噪声便于后续处理。这一步得到的是Log-Mel频谱。离散余弦变换对上面得到的Log-Mel频谱向量进行DCT得到MFCC系数。DCT在这里起到了“倒谱分析”的作用。在理想的同态信号处理中我们希望将语音声源激励与声道滤波的卷积结果解卷积开。取对数后卷积变加法再做逆傅里叶变换即倒谱分析可以将信号分离。我们使用DCT来近似这个逆傅里叶变换因为它计算高效且能产生高度不相关的特征这对很多统计模型如高斯混合模型非常有利。通常我们只保留前12-13个DCT系数因为它们包含了声道形状的主要信息而更高的系数往往代表细微的细节或噪声。动态特征提取静态的MFCC系数只描述了一帧的静态特性。为了捕捉语音的动态变化如音素的过渡我们还会计算一阶差分系数Delta和二阶差分系数Delta-Delta。它们分别近似于MFCC系数对时间的一阶和二阶导数能够显著提升识别性能。注意整个流程中参数的选择如帧长、帧移、滤波器个数、保留的DCT系数个数需要根据具体任务和语音特性进行调整没有一成不变的“最优值”。例如对于音高较高的语音或音乐可能需要调整Mel滤波器的上限频率。2.2 关键参数选择背后的考量参数不是随便填的数字每一个背后都有其声学或工程上的理由。这里我们深入探讨几个最核心的参数理解它们如何影响最终的特征。采样频率与帧长/帧移这是所有参数的起点。假设采样率Fs 16000 Hz。根据短时平稳性假设帧长通常取20-30ms。我们取frame_len 0.025 * Fs 400个采样点。为了帧间平滑过渡帧移取帧长的一半即frame_shift 0.0125 * Fs 200点。这里为什么是0.025和0.012520-30ms是大量实验证实的语音准平稳段长度。帧移取一半是权衡重叠太少会导致帧间变化剧烈丢失信息重叠太多则计算冗余。一半是一个经验上的甜点。FFT点数为了计算FFT我们需要将每帧信号400点通过补零扩展到2的整数次幂通常取NFFT 512。这保证了FFT计算效率并且NFFT frame_len的补零操作相当于对频谱进行了插值让频谱图看起来更平滑便于后续Mel滤波器组的能量求和计算。如果取NFFT 400得到的频谱是“粗糙”的每个频率点间隔大用三角滤波器去求和时误差会更大。Mel滤波器个数通常取20-40个。滤波器个数决定了Mel频谱的维度即特征的粗细粒度。个数太少如10个则频率分辨率太低无法区分细致的频谱形状个数太多如60个则特征维度高且相邻滤波器高度相关增加了计算量和模型过拟合的风险。对于大多数语音识别任务26个滤波器是一个稳健的起点。它能在保留足够频谱细节和控制维度之间取得良好平衡。保留的MFCC系数个数通常取12-13个再加上第0个系数对数能量构成13维的静态MFCC。为什么是12-13DCT具有能量压缩特性前几个系数包含了Log-Mel频谱中的大部分信息尤其是表征声道形状的平滑、低频分量。更高的系数往往对应频谱的快速变化和细节这些部分更容易受到噪声和个体差异的影响对识别贡献小且不稳定。因此截断高阶系数实际上起到了一种降维和去噪的作用。预加重系数通常取0.97。这个值接近1实现了一个轻微的高通滤波。其传递函数为H(z) 1 - a*z^{-1}。系数a决定了高频提升的程度。0.97是一个广泛使用的经验值它能有效平衡高频提升和数值稳定性。如果a太接近1如0.99滤波效果过强可能放大高频噪声如果太小如0.90则提升效果不足。3. 基于MATLAB的MFCC特征提取实现详解理论说得再多不如动手实现一遍。下面我将结合代码详细讲解如何在MATLAB中一步步实现上述流程。我们会从读取一个语音文件开始完整地走完整个流程并可视化中间结果让你对每一步的输入输出有直观的认识。首先我们需要准备环境。确保你的MATLAB安装了信号处理工具箱因为我们需要用到hamming,fft等函数。当然这些基础函数也可以自己写但为了效率和可靠性直接使用工具箱是更好的选择。3.1 语音读取与预处理实操我们从一个干净的语音文件开始。这里假设你有一个名为speech.wav的16kHz单声道语音文件。% 1. 读取语音文件 [audio, Fs] audioread(speech.wav); % Fs为采样频率audio为信号向量 % 确保是单声道 if size(audio, 2) 1 audio mean(audio, 2); % 如果是立体声取平均转为单声道 end % 2. 预加重 pre_emphasis_coeff 0.97; emphasized_audio filter([1, -pre_emphasis_coeff], 1, audio); % 可视化原始信号和预加重后信号可选 t (0:length(audio)-1) / Fs; figure; subplot(2,1,1); plot(t, audio); title(原始语音信号); xlabel(时间 (s)); ylabel(幅度); subplot(2,1,2); plot(t, emphasized_audio); title(预加重后语音信号); xlabel(时间 (s)); ylabel(幅度);提示filter([1, -a], 1, x)是实现y(n) x(n) - a*x(n-1)的标准方法。[1, -a]是分子系数1是分母系数代表一个FIR滤波器。分帧与加窗这是将一维信号转化为二维矩阵帧数 x 帧长的关键步骤。% 3. 分帧参数设置 frame_length round(0.025 * Fs); % 25ms的帧长以采样点计 frame_step round(0.010 * Fs); % 10ms的帧移通常小于帧长以实现重叠 signal_length length(emphasized_audio); num_frames floor((signal_length - frame_length) / frame_step) 1; % 初始化帧矩阵 frames zeros(frame_length, num_frames); for i 0:num_frames-1 start_index i * frame_step 1; end_index start_index frame_length - 1; % 确保索引不超出范围最后一帧可能不够长需要补零 if end_index signal_length frame emphasized_audio(start_index:signal_length); frames(:, i1) [frame; zeros(end_index - signal_length, 1)]; else frames(:, i1) emphasized_audio(start_index:end_index); end end % 4. 加窗汉明窗 hamming_window hamming(frame_length); % 对每一帧点乘窗函数 windowed_frames frames .* hamming_window;这里有一个实操心得在循环中直接索引拼接矩阵在MATLAB中对于长语音可能效率不高。更高效也更MATLAB风格的向量化方法是使用buffer函数来自信号处理工具箱或通过构造索引矩阵来实现。但为了代码清晰易懂我们使用了直观的循环。在实际处理大批量数据时可以考虑优化。3.2 频谱计算与Mel滤波器组应用接下来我们对每一帧加窗信号进行FFT并计算其功率谱。% 5. 计算功率谱 NFFT 512; % FFT点数通常取2的幂且大于帧长 mag_frames abs(fft(windowed_frames, NFFT)); % 计算幅度谱 pow_frames ((1/NFFT) * (mag_frames .^ 2)); % 计算功率谱1/NFFT是归一化因子 % 由于频谱对称我们只取前半部分0~Nyquist频率 num_fft_bins NFFT/2 1; pow_frames pow_frames(1:num_fft_bins, :);现在到了核心环节设计并应用Mel滤波器组。我们需要在Hz频率和Mel频率之间进行转换。% 6. 定义Mel滤波器组 num_filters 26; % 滤波器个数 low_freq_mel 0; high_freq_mel 2595 * log10(1 (Fs/2) / 700); % 将最高频率Fs/2转换为Mel mel_points linspace(low_freq_mel, high_freq_mel, num_filters 2); % 在Mel尺度上均匀取点 hz_points 700 * (10.^(mel_points / 2595) - 1); % 将Mel点转换回Hz频率 % 将Hz频率点映射到FFT的bin索引上 fft_bin_indices floor((NFFT 1) * hz_points / Fs); % 创建滤波器组矩阵 filter_bank zeros(num_filters, num_fft_bins); for m 1:num_filters f_left fft_bin_indices(m); f_center fft_bin_indices(m 1); f_right fft_bin_indices(m 2); for k f_left:f_center filter_bank(m, k1) (k - f_left) / (f_center - f_left); end for k f_center:f_right filter_bank(m, k1) (f_right - k) / (f_right - f_center); end end % 可视化Mel滤波器组 figure; plot(0:(num_fft_bins-1), filter_bank.); xlabel(FFT Bin索引); ylabel(幅度); title(Mel三角滤波器组Hz频率轴); % 为了更直观可以尝试将x轴转换为Hz频率x_Hz (0:(num_fft_bins-1)) * Fs / NFFT;注意公式mel 2595 * log10(1 f/700)是将线性频率f(Hz) 转换为Mel频率的常用近似公式。其反变换为f 700 * (10^(mel/2595) - 1)。不同的文献可能使用不同的系数如1127但2595和700是最常见的组合效果差异不大。应用滤波器组并取对数% 7. 应用滤波器组并取对数能量 filter_bank_energies filter_bank * pow_frames; % 矩阵乘法每个滤波器对每一帧的频谱能量求和 % 防止出现log(0)的情况加一个很小的数 filter_bank_energies max(filter_bank_energies, 1e-10); log_filter_bank_energies log10(filter_bank_energies); % 得到Log-Mel频谱3.3 离散余弦变换与动态特征计算最后一步对Log-Mel频谱进行DCT得到静态的MFCC系数。% 8. 离散余弦变换 (DCT) 获取MFCC系数 num_ceps 12; % 需要保留的MFCC系数个数不包括第0阶 mfcc dct(log_filter_bank_energies); mfcc mfcc(1:(num_ceps1), :); % 保留0~12阶系数 % 通常我们会去掉第0阶系数对数能量因为它容易受录音条件影响波动较大。 % 但有时也会保留作为一个单独的特征。这里我们先保留后续可根据需求选择。 % mfcc mfcc(2:end, :); % 如果想去掉第0阶 % 转置一下使得每一行代表一帧每一列代表一个MFCC系数更常见的表示 mfcc mfcc;至此我们得到了最基本的静态MFCC特征矩阵其大小为[num_frames, num_ceps1]。为了捕捉动态信息我们计算一阶和二阶差分Delta和Delta-Delta。% 9. 计算一阶差分Delta系数 % 使用一个简单的回归器近似导数 delta_window 2; % 通常取2 delta_coeffs [-delta_window:delta_window]; denom sum(delta_coeffs.^2); delta zeros(size(mfcc)); for i 1:size(mfcc, 2) % 对每个系数维度 padded [repmat(mfcc(1, i), delta_window, 1); mfcc(:, i); repmat(mfcc(end, i), delta_window, 1)]; for t 1:size(mfcc, 1) delta(t, i) sum(delta_coeffs .* padded(t:t2*delta_window)) / denom; end end % 10. 计算二阶差分Delta-Delta系数方法同上 delta_delta zeros(size(delta)); for i 1:size(delta, 2) padded_d [repmat(delta(1, i), delta_window, 1); delta(:, i); repmat(delta(end, i), delta_window, 1)]; for t 1:size(delta, 1) delta_delta(t, i) sum(delta_coeffs .* padded_d(t:t2*delta_window)) / denom; end end % 11. 特征拼接静态MFCC Delta Delta-Delta % 假设我们决定使用静态MFCC不含能量加上动态特征 static_features mfcc(:, 2:end); % 去掉第0阶能量 final_features [static_features, delta(:, 2:end), delta_delta(:, 2:end)];现在final_features就是一个完整的MFCC特征矩阵通常维度是[num_frames, (num_ceps * 3)]例如12维静态MFCC 12维Delta 12维Delta-Delta 36维。这个矩阵就可以作为后续语音识别或说话人识别模型的输入了。4. 仿真结果分析与可视化解读代码跑通了输出了一堆数字怎么判断我们提取的特征是“好”还是“坏”呢可视化是关键。通过将中间结果和最终特征画出来我们可以直观地检查流程是否合理特征是否具有区分度。首先我们可以绘制语音的波形图、频谱图以及MFCC特征图进行对比观察。% 绘制原始语音波形、频谱图、MFCC特征图 figure(Position, [100, 100, 1200, 800]); % 子图1: 原始语音波形 subplot(3, 1, 1); plot((0:length(audio)-1)/Fs, audio); xlabel(时间 (s)); ylabel(幅度); title(原始语音信号波形); xlim([0, length(audio)/Fs]); grid on; % 子图2: 语谱图 (Spectrogram) subplot(3, 1, 2); spectrogram(audio, hamming(256), 250, 512, Fs, yaxis); title(语音信号的语谱图); colorbar; % 子图3: MFCC特征热图 subplot(3, 1, 3); imagesc(1:size(static_features,2), (0:size(static_features,1)-1)*frame_step/Fs, static_features); xlabel(MFCC系数索引); ylabel(时间 (s)); title(MFCC特征序列 (静态部分)); colorbar; colormap(jet); % 使用jet颜色图便于观察对比 axis xy; % 确保y轴方向正确时间向下增长通过对比这三个图你可以看到波形图显示了信号的时域振幅变化语谱图频谱随时间的变化显示了信号频率成分的时变特性你能看到共振峰能量集中的频带形成的条纹而MFCC特征图则是一种高度压缩和抽象后的表示它抹去了精细的谐波结构突出了由声道形状决定的平滑频谱包络。不同音素如元音/a/、/i/对应的MFCC图案会有明显的不同。如何解读MFCC特征图纵轴时间代表语音的进程从左到右是一帧一帧的。横轴MFCC系数索引通常0号系数是能量1号系数代表频谱包络的“整体倾斜度”低频vs高频能量2号、3号系数代表更细致的频谱形状比如第一、第二共振峰的位置信息。系数索引越大代表频谱包络中越快速变化的细节。颜色代表该MFCC系数值的大小。你可以看到在发不同音的时候MFCC系数的模式会发生明显变化。为了更定量地分析我们可以计算同一说话人发相同元音、不同元音以及不同说话人发相同元音的MFCC特征并观察它们的距离。% 假设我们有两段语音的特征矩阵 feat1 和 feat2已按帧对齐或计算了统计量如均值 % 常用动态时间规整 (DTW) 计算序列距离或简单计算均值向量的欧氏距离 % 这里演示计算均值向量的距离 mean_feat1 mean(static_features_from_audio1, 1); mean_feat2 mean(static_features_from_audio2, 1); % 计算欧氏距离 dist_euclidean sqrt(sum((mean_feat1 - mean_feat2).^2)); fprintf(两段语音MFCC均值向量的欧氏距离为: %.4f\n, dist_euclidean); % 也可以计算余弦相似度 cos_sim dot(mean_feat1, mean_feat2) / (norm(mean_feat1) * norm(mean_feat2)); fprintf(两段语音MFCC均值向量的余弦相似度为: %.4f\n, cos_sim);在理想情况下同一个说话人说同一个词的MFCC距离应该很小相似度高说不同词的MFCC距离应该较大而不同说话人说同一个词的MFCC距离介于两者之间。这验证了MFCC特征的有效性。5. 常见问题、调试技巧与性能优化在实际实现和调试MFCC提取代码时你肯定会遇到各种问题。下面我整理了一些常见坑点和解决技巧很多都是我在项目里踩过坑才总结出来的。5.1 数值不稳定与边界处理问题1对数运算遇到零或负值。在计算log10(filter_bank_energies)时如果某个滤波器的输出能量为零或极小由于数值计算误差可能为负取对数会得到-Inf或复数导致后续计算崩溃。解决技巧在取对数前一定要加一个极小值的下限。filter_bank_energies max(filter_bank_energies, 1e-10); % 确保所有值大于一个很小的正数问题2分帧时最后一帧长度不足。如果语音总长度不是帧移的整数倍最后一帧的采样点会不够frame_length。解决技巧采用补零Zero-Padding的方式。就像我们代码中做的判断end_index是否超出信号长度如果超出则用零补齐。这能保证所有帧长度一致方便矩阵运算。虽然补零会引入一点频谱畸变但对整体影响通常可接受。问题3Mel滤波器组在低频或高频的bin索引越界。在计算fft_bin_indices floor((NFFT 1) * hz_points / Fs)时由于数值计算第一个点可能为0最后一个点可能超过num_fft_bins。解决技巧在创建滤波器循环中对索引k进行边界检查确保其在[1, num_fft_bins]范围内。更稳健的做法是在计算hz_points时就将最低频率设为大于0如50Hz以避免直流偏移最高频率设为略小于Fs/2。5.2 特征维度与参数调优问题4MFCC特征看起来“太平滑”缺乏区分度。这可能是因为Mel滤波器个数 (num_filters) 太少或者保留的DCT系数 (num_ceps) 太多。滤波器太少导致频谱信息压缩过度保留的高阶DCT系数太多引入了噪声而非有效信息。调试建议首先可视化你的Mel滤波器组和Log-Mel频谱。确保滤波器覆盖了有效的频率范围通常是0到Fs/2并且形状正确。尝试增加num_filters例如从26增加到40观察MFCC特征图是否显示出更丰富的纵向条纹。尝试减少num_ceps例如从12减少到8看看是否核心的区分信息仍然保留。通常前6-8个系数已经包含了绝大部分声道信息。问题5不同录音条件下的特征差异巨大。MFCC特征尤其是第0阶能量对录音音量、麦克风距离、背景噪声非常敏感。这会导致模型难以学习到本质的语音内容。解决技巧进行特征归一化。常见的有两种倒谱均值减在整个语句层面上计算每一维MFCC系数的均值然后从每一帧中减去这个均值。这可以消除信道麦克风、传输线路的固定影响。mfcc_normalized mfcc - mean(mfcc, 1); % 按列求均值滑动窗归一化对于实时流式处理可以使用一个滑动窗口计算窗口内特征的均值和方差进行标准化。 归一化操作应在计算动态特征Delta之前进行因为动态特征对均值偏移也很敏感。5.3 MATLAB代码性能优化我们之前的示例代码为了清晰大量使用了循环。在处理长语音或大批量数据时效率可能成为瓶颈。优化1向量化分帧操作。可以使用buffer函数信号处理工具箱一次性完成分帧和补零。% 使用buffer函数分帧 (更高效) windowed_frames buffer(emphasized_audio, frame_length, frame_length-frame_step, nodelay); % 注意buffer函数输出的矩阵每一列是一帧但帧的顺序可能需要调整最后一列可能是补零的帧 % 需要根据实际情况进行裁剪和转置优化2矩阵运算替代循环应用滤波器组。我们已经用矩阵乘法filter_bank * pow_frames实现了滤波器组的应用这本身就是向量化的。确保你的filter_bank是[num_filters, num_fft_bins]pow_frames是[num_fft_bins, num_frames]这样一次乘法就得到所有帧所有滤波器的能量。优化3使用内置DCT函数。MATLAB的dct函数默认对矩阵的每一列进行操作。我们的log_filter_bank_energies是[num_filters, num_frames]每一列是一帧的Log-Mel频谱。直接调用mfcc dct(log_filter_bank_energies)会对每一帧每一列分别做DCT这正是我们需要的无需循环。一个综合性的、经过部分优化的MFCC提取函数框架可能长这样function [mfcc_features, delta, delta_delta] extract_mfcc(audio, Fs, varargin) % 解析输入参数设置默认值 p inputParser; addParameter(p, FrameLength, 0.025, isnumeric); % 秒 addParameter(p, FrameShift, 0.010, isnumeric); % 秒 addParameter(p, NumFilters, 26, isnumeric); addParameter(p, NumCeps, 12, isnumeric); addParameter(p, PreEmph, 0.97, isnumeric); addParameter(p, NFFT, 512, isnumeric); parse(p, varargin{:}); params p.Results; % 1. 预加重 emphasized_audio filter([1, -params.PreEmph], 1, audio(:)); % 2. 分帧 (使用buffer高效) frame_len_samples round(params.FrameLength * Fs); frame_step_samples round(params.FrameShift * Fs); frames buffer(emphasized_audio, frame_len_samples, frame_len_samples-frame_step_samples, nodelay); num_frames size(frames, 2); % 3. 加窗 window hamming(frame_len_samples); windowed_frames frames .* window; % 4. 计算功率谱 (向量化) mag_spec abs(fft(windowed_frames, params.NFFT)); pow_spec (1/params.NFFT) * (mag_spec(1:params.NFFT/21, :) .^ 2); % 5. 创建并应用Mel滤波器组 [filter_bank, ~] mel_filter_bank(params.NumFilters, params.NFFT, Fs); filter_bank_energies max(filter_bank * pow_spec, 1e-10); log_energies log10(filter_bank_energies); % 6. DCT得到静态MFCC mfcc_static dct(log_energies); mfcc_static mfcc_static(1:params.NumCeps1, :); % 保留0~NumCeps阶 % 7. (可选) 倒谱均值减 % mfcc_static mfcc_static - mean(mfcc_static, 2); % 8. 计算动态特征 [delta, delta_delta] compute_delta_features(mfcc_static(2:end, :)); % 去掉能量再计算 % 9. 拼接特征 mfcc_features [mfcc_static(2:end, :); delta; delta_delta]; % 转置帧为行 end function [filter_bank, f] mel_filter_bank(num_filters, nfft, fs) % 创建Mel滤波器组的子函数 % ... 实现代码同上文 ... end function [delta, delta_delta] compute_delta_features(features) % 计算Delta和Delta-Delta的子函数 % ... 实现代码同上文可向量化优化 ... end把这个函数封装好以后提取特征就是一行代码的事features extract_mfcc(audio, Fs, NumFilters, 40, NumCeps, 13);。参数可以灵活调整适应不同的任务需求。最后再分享一个调试心得当你怀疑MFCC特征提取的某个环节有问题时最好的办法是找一个已知的、简单的信号比如一个纯净的正弦波或方波作为输入然后逐步检查每个环节的输出。例如对于一个440Hz的正弦波它的功率谱应该只在440Hz处有一个尖峰。经过Mel滤波器组后只有中心频率覆盖440Hz的那个滤波器输出能量很大其他滤波器输出接近零。这样一步步验证能快速定位问题所在。语音信号处理就是这样理论、代码和调试经验三者缺一不可。
返回列表