MATLAB实现声纹识别特征提取

声纹识别特征概述

声纹识别依赖于从语音信号中提取能够唯一表征说话人身份的声学特征。这些特征需要具有良好的区分性、稳定性和抗噪性。

核心特征类型

• 频谱特征：MFCC、LPCC、PLP
• 基音特征：基频、共振峰
• ** prosodic特征**：能量、语速、节奏
• 高级特征：GFCC、IMFCC

特征提取流程框架

语音信号 → 预处理 → 分帧加窗 → 特征提取 → 特征后处理 → 特征向量

MATLAB实现

1. 主函数框架

function [features, featureNames] = voiceprintFeatureExtraction(audioData, fs, varargin)
% 声纹识别特征提取主函数
% 输入：
%   audioData - 音频数据
%   fs - 采样率
% 输出：
%   features - 提取的特征向量
%   featureNames - 特征名称

    % 参数解析
    p = inputParser;
    addParameter(p, 'frameLength', 0.025, @isnumeric); % 帧长(秒)
    addParameter(p, 'frameOverlap', 0.010, @isnumeric); % 帧重叠(秒)
    addParameter(p, 'mfccCoeffs', 13, @isnumeric); % MFCC系数个数
    addParameter(p, 'lpOrder', 12, @isnumeric); % LPC阶数
    addParameter(p, 'includeDelta', true, @islogical); % 是否包含差分特征
    addParameter(p, 'includeEnergy', true, @islogical); % 是否包含能量特征
    parse(p, varargin{:});
    params = p.Results;
    
    fprintf('开始声纹特征提取...\n');
    fprintf('采样率: %d Hz, 音频长度: %.2f 秒\n', fs, length(audioData)/fs);
    
    % 1. 语音预处理
    processedAudio = preprocessAudio(audioData, fs);
    
    % 2. 分帧加窗
    frames = frameBlocking(processedAudio, fs, params.frameLength, params.frameOverlap);
    
    % 3. 多维度特征提取
    featureSet = extractAllFeatures(frames, fs, params);
    
    % 4. 特征后处理和统计聚合
    features = postProcessFeatures(featureSet);
    featureNames = generateFeatureNames(params);
    
    fprintf('特征提取完成，特征维度: %d\n', length(features));
end

2. 语音预处理函数

function processedAudio = preprocessAudio(audioData, fs)
% 语音信号预处理
    
    % 转换为单声道（如果是立体声）
    if size(audioData, 2) > 1
        audioData = mean(audioData, 2);
        fprintf('转换为单声道信号\n');
    end
    
    % 预加重：提升高频分量，H(z) = 1 - 0.97z^(-1)
    preemph = [1, -0.97];
    processedAudio = filter(preemph, 1, audioData);
    
    % 静音检测和端点检测（可选）
    voicedFrames = voiceActivityDetection(processedAudio, fs);
    processedAudio = processedAudio(voicedFrames);
    
    % 幅度归一化
    processedAudio = processedAudio / max(abs(processedAudio));
    
    fprintf('预处理完成，有效语音长度: %.2f 秒\n', length(processedAudio)/fs);
end

function voicedSegments = voiceActivityDetection(audio, fs)
% 基于能量的端点检测
    
    frameLen = round(0.025 * fs); % 25ms帧
    overlap = round(0.015 * fs);  % 15ms重叠
    step = frameLen - overlap;
    
    numFrames = floor((length(audio) - frameLen) / step) + 1;
    energy = zeros(numFrames, 1);
    
    % 计算每帧能量
    for i = 1:numFrames
        startIdx = (i-1)*step + 1;
        endIdx = startIdx + frameLen - 1;
        frame = audio(startIdx:endIdx);
        energy(i) = sum(frame.^2);
    end
    
    % 基于能量阈值进行静音检测
    energyThresh = 0.1 * max(energy);
    voicedFrames = energy > energyThresh;
    
    % 找到有声段的起始和结束
    voicedSegments = findVoicedSegments(voicedFrames, step, frameLen, length(audio));
end

3. 分帧加窗处理

function frames = frameBlocking(audio, fs, frameLength, frameOverlap)
% 分帧和加窗处理
    
    frameSamples = round(frameLength * fs);
    overlapSamples = round(frameOverlap * fs);
    stepSamples = frameSamples - overlapSamples;
    
    numFrames = floor((length(audio) - frameSamples) / stepSamples) + 1;
    frames = zeros(frameSamples, numFrames);
    
    % 汉明窗
    window = hamming(frameSamples);
    
    for i = 1:numFrames
        startIdx = (i-1)*stepSamples + 1;
        endIdx = startIdx + frameSamples - 1;
        
        if endIdx <= length(audio)
            frame = audio(startIdx:endIdx) .* window;
            frames(:, i) = frame;
        end
    end
    
    fprintf('分帧完成: %d帧, 帧长%d点, 重叠%d点\n', ...
        numFrames, frameSamples, overlapSamples);
end

4. MFCC特征提取

function mfccFeatures = extractMFCC(frames, fs, numCoeffs)
% 提取MFCC特征
    
    [frameLength, numFrames] = size(frames);
    mfccFeatures = zeros(numCoeffs, numFrames);
    
    % 梅尔滤波器组参数
    lowFreq = 0;
    highFreq = fs/2;
    numFilters = 26; % 梅尔滤波器数量
    
    % 创建梅尔滤波器组
    melFilterBank = createMelFilterBank(frameLength, fs, numFilters, lowFreq, highFreq);
    
    for i = 1:numFrames
        frame = frames(:, i);
        
        % 1. 预加重（已在预处理完成）
        
        % 2. 计算功率谱
        fftFrame = fft(frame);
        powerSpectrum = abs(fftFrame(1:frameLength/2+1)).^2;
        
        % 3. 梅尔滤波器组处理
        melSpectrum = melFilterBank * powerSpectrum;
        logMelSpectrum = log(melSpectrum + eps); % 加eps防止log(0)
        
        % 4. DCT变换得到MFCC
        mfcc = dct(logMelSpectrum);
        mfccFeatures(:, i) = mfcc(1:numCoeffs);
    end
    
    % 提升MFCC（可选）
    mfccFeatures = lifterMFCC(mfccFeatures, 22);
end

function melFilterBank = createMelFilterBank(fftSize, fs, numFilters, lowFreq, highFreq)
% 创建梅尔滤波器组
    
    % 频率转换为梅尔尺度
    lowMel = 2595 * log10(1 + lowFreq/700);
    highMel = 2595 * log10(1 + highFreq/700);
    melPoints = linspace(lowMel, highMel, numFilters+2);
    
    % 梅尔尺度转回频率
    freqPoints = 700 * (10.^(melPoints/2595) - 1);
    
    % 转换为FFT bin索引
    fftBins = floor((fftSize/2 + 1) * freqPoints / (fs/2));
    
    % 创建三角形滤波器
    melFilterBank = zeros(numFilters, fftSize/2+1);
    
    for i = 1:numFilters
        leftBin = fftBins(i);
        centerBin = fftBins(i+1);
        rightBin = fftBins(i+2);
        
        % 上升部分
        if leftBin ~= centerBin
            melFilterBank(i, leftBin:centerBin) = ...
                linspace(0, 1, centerBin - leftBin + 1);
        end
        
        % 下降部分
        if centerBin ~= rightBin
            melFilterBank(i, centerBin:rightBin) = ...
                linspace(1, 0, rightBin - centerBin + 1);
        end
    end
end

function liftedMFCC = lifterMFCC(mfcc, liftCoeff)
% MFCC提升，增强高阶系数
    [numCoeffs, numFrames] = size(mfcc);
    lift = 1 + (liftCoeff/2) * sin(pi * (0:numCoeffs-1)' / liftCoeff);
    liftedMFCC = mfcc .* lift;
end

5. LPCC特征提取

function lpccFeatures = extractLPCC(frames, fs, lpOrder)
% 提取LPCC特征（线性预测倒谱系数）
    
    [frameLength, numFrames] = size(frames);
    lpccFeatures = zeros(lpOrder, numFrames);
    
    for i = 1:numFrames
        frame = frames(:, i);
        
        % 计算LPC系数
        lpcCoeffs = lpc(frame, lpOrder);
        
        % LPC转LPCC
        lpcc = lpc2lpcc(lpcCoeffs, lpOrder);
        lpccFeatures(:, i) = lpcc;
    end
end

function lpcc = lpc2lpcc(lpcCoeffs, numCoeffs)
% LPC系数转换为LPCC系数
    
    p = length(lpcCoeffs) - 1; % LPC阶数
    lpcc = zeros(numCoeffs, 1);
    
    lpcc(1) = log(lpcCoeffs(1)); % 增益项
    
    for n = 2:numCoeffs
        if n <= p
            lpcc(n) = -lpcCoeffs(n);
            for k = 1:n-1
                lpcc(n) = lpcc(n) - (1 - k/n) * lpcCoeffs(k) * lpcc(n-k);
            end
        else
            lpcc(n) = 0;
            for k = 1:p
                lpcc(n) = lpcc(n) - (1 - k/n) * lpcCoeffs(k) * lpcc(n-k);
            end
        end
    end
end

6. 基音和共振峰特征

function [pitch, formants] = extractProsodicFeatures(frames, fs)
% 提取基音和共振峰特征
    
    [frameLength, numFrames] = size(frames);
    pitch = zeros(1, numFrames);
    formants = zeros(3, numFrames); % 提取前三个共振峰
    
    for i = 1:numFrames
        frame = frames(:, i);
        
        % 基音提取（自相关法）
        pitch(i) = computePitch(frame, fs);
        
        % 共振峰提取（LPC谱峰值）
        formants(:, i) = computeFormants(frame, fs);
    end
end

function pitch = computePitch(frame, fs)
% 基于自相关的基音检测
    
    % 预加重
    preemphFrame = filter([1 -0.97], 1, frame);
    
    % 自相关计算
    correlation = xcorr(preemphFrame, 'coeff');
    correlation = correlation(length(frame):end); % 取正延迟部分
    
    % 寻找基音峰值（在合理基音范围内）
    minPitch = 50; % Hz
    maxPitch = 400; % Hz
    
    minLag = round(fs / maxPitch);
    maxLag = round(fs / minPitch);
    
    [peaks, locations] = findpeaks(correlation(minLag:maxLag));
    
    if ~isempty(peaks)
        [~, maxPeakIdx] = max(peaks);
        pitchLag = locations(maxPeakIdx) + minLag - 1;
        pitch = fs / pitchLag;
    else
        pitch = 0; % 无声段
    end
end

function formants = computeFormants(frame, fs)
% 基于LPC的共振峰提取
    
    lpcOrder = 12; % LPC阶数
    lpcCoeffs = lpc(frame, lpcOrder);
    
    % 计算LPC谱
    [h, w] = freqz(1, lpcCoeffs, 1024, fs);
    lpcSpectrum = 20*log10(abs(h));
    
    % 寻找频谱峰值（共振峰）
    [peaks, locations] = findpeaks(lpcSpectrum, 'SortStr', 'descend');
    
    numFormants = min(3, length(peaks));
    formants = zeros(3, 1);
    
    if numFormants > 0
        formantFreqs = w(locations(1:numFormants));
        formants(1:numFormants) = formantFreqs;
    end
end

7. 多维度特征整合

function featureSet = extractAllFeatures(frames, fs, params)
% 提取所有特征并整合
    
    [frameLength, numFrames] = size(frames);
    featureSet = struct();
    
    fprintf('提取MFCC特征...\n');
    featureSet.mfcc = extractMFCC(frames, fs, params.mfccCoeffs);
    
    fprintf('提取LPCC特征...\n');
    featureSet.lpcc = extractLPCC(frames, fs, params.lpOrder);
    
    fprintf('提取基音和共振峰特征...\n');
    [featureSet.pitch, featureSet.formants] = extractProsodicFeatures(frames, fs);
    
    fprintf('提取能量特征...\n');
    featureSet.energy = extractEnergy(frames);
    
    % 差分特征（一阶和二阶差分）
    if params.includeDelta
        featureSet.mfccDelta = computeDelta(featureSet.mfcc);
        featureSet.mfccDeltaDelta = computeDelta(featureSet.mfccDelta);
        featureSet.lpccDelta = computeDelta(featureSet.lpcc);
        featureSet.lpccDeltaDelta = computeDelta(featureSet.lpccDelta);
    end
    
    fprintf('特征提取完成: MFCC(%dx%d), LPCC(%dx%d), 基音(%dx%d)\n', ...
        size(featureSet.mfcc,1), size(featureSet.mfcc,2), ...
        size(featureSet.lpcc,1), size(featureSet.lpcc,2), ...
        1, size(featureSet.pitch,2));
end

function energy = extractEnergy(frames)
% 计算帧能量
    energy = sum(frames.^2, 1);
end

function delta = computeDelta(staticFeatures)
% 计算差分特征（一阶导数）
    [numCoeffs, numFrames] = size(staticFeatures);
    delta = zeros(numCoeffs, numFrames);
    
    for t = 2:numFrames-1
        delta(:, t) = (staticFeatures(:, t+1) - staticFeatures(:, t-1)) / 2;
    end
    
    % 边界处理
    delta(:, 1) = delta(:, 2);
    delta(:, end) = delta(:, end-1);
end

8. 特征后处理

function finalFeatures = postProcessFeatures(featureSet)
% 特征后处理和统计聚合
    
    featureMatrix = [];
    
    % MFCC特征统计量（均值、标准差）
    mfccMean = mean(featureSet.mfcc, 2);
    mfccStd = std(featureSet.mfcc, 0, 2);
    featureMatrix = [featureMatrix; mfccMean; mfccStd];
    
    % LPCC特征统计量
    lpccMean = mean(featureSet.lpcc, 2);
    lpccStd = std(featureSet.lpcc, 0, 2);
    featureMatrix = [featureMatrix; lpccMean; lpccStd];
    
    % 基音特征统计量
    pitchMean = mean(featureSet.pitch(featureSet.pitch > 0)); % 只统计有声段
    pitchStd = std(featureSet.pitch(featureSet.pitch > 0));
    featureMatrix = [featureMatrix; pitchMean; pitchStd];
    
    % 共振峰特征统计量
    formantMeans = mean(featureSet.formants, 2);
    formantStds = std(featureSet.formants, 0, 2);
    featureMatrix = [featureMatrix; formantMeans; formantStds];
    
    % 能量特征统计量
    energyMean = mean(featureSet.energy);
    energyStd = std(featureSet.energy);
    featureMatrix = [featureMatrix; energyMean; energyStd];
    
    % 差分特征统计量
    if isfield(featureSet, 'mfccDelta')
        mfccDeltaMean = mean(featureSet.mfccDelta, 2);
        mfccDeltaStd = std(featureSet.mfccDelta, 0, 2);
        featureMatrix = [featureMatrix; mfccDeltaMean; mfccDeltaStd];
    end
    
    % 去除NaN值
    finalFeatures = featureMatrix(~isnan(featureMatrix));
end

function featureNames = generateFeatureNames(params)
% 生成特征名称（用于可解释性）
    featureNames = {};
    
    % MFCC特征名
    for i = 1:params.mfccCoeffs
        featureNames{end+1} = sprintf('MFCC_%d_Mean', i);
    end
    for i = 1:params.mfccCoeffs
        featureNames{end+1} = sprintf('MFCC_%d_Std', i);
    end
    
    % LPCC特征名
    for i = 1:params.lpOrder
        featureNames{end+1} = sprintf('LPCC_%d_Mean', i);
    end
    for i = 1:params.lpOrder
        featureNames{end+1} = sprintf('LPCC_%d_Std', i);
    end
    
    % 其他特征名
    featureNames = [featureNames, {'Pitch_Mean', 'Pitch_Std'}];
    featureNames = [featureNames, {'Formant1_Mean', 'Formant1_Std'}];
    featureNames = [featureNames, {'Formant2_Mean', 'Formant2_Std'}];
    featureNames = [featureNames, {'Formant3_Mean', 'Formant3_Std'}];
    featureNames = [featureNames, {'Energy_Mean', 'Energy_Std'}];
    
    if params.includeDelta
        for i = 1:params.mfccCoeffs
            featureNames{end+1} = sprintf('MFCC_Delta_%d_Mean', i);
        end
        for i = 1:params.mfccCoeffs
            featureNames{end+1} = sprintf('MFCC_Delta_%d_Std', i);
        end
    end
end

9. 可视化分析工具

function visualizeFeatures(audioData, fs, features, featureNames)
% 特征可视化分析
    
    figure('Position', [100, 100, 1400, 800]);
    
    % 原始语音波形
    subplot(3, 3, 1);
    t = (0:length(audioData)-1) / fs;
    plot(t, audioData);
    title('原始语音信号');
    xlabel('时间(s)'); ylabel('幅度');
    grid on;
    
    % 语谱图
    subplot(3, 3, 2);
    spectrogram(audioData, 256, 128, 256, fs, 'yaxis');
    title('语谱图');
    
    % MFCC特征示例（第一帧）
    subplot(3, 3, 3);
    mfccExample = features(1:13); % 假设前13个是MFCC均值
    bar(mfccExample);
    title('MFCC特征示例');
    xlabel('系数索引'); ylabel('幅度');
    grid on;
    
    % 特征分布直方图
    subplot(3, 3, 4);
    histogram(features, 30);
    title('特征值分布');
    xlabel('特征值'); ylabel('频数');
    grid on;
    
    % 特征重要性排序
    subplot(3, 3, 5);
    [~, idx] = sort(abs(features), 'descend');
    topFeatures = features(idx(1:min(10, length(features))));
    bar(topFeatures);
    title('Top 10特征值');
    xlabel('特征索引'); ylabel('幅度');
    grid on;
    
    % 特征相关性矩阵（部分特征）
    subplot(3, 3, 6);
    numShow = min(15, length(features));
    corrMatrix = corrcoef([features(1:numShow)']);
    imagesc(corrMatrix);
    colorbar;
    title('特征相关性矩阵');
    xlabel('特征索引'); ylabel('特征索引');
    
    % 特征维度约减可视化（PCA）
    subplot(3, 3, 7);
    if length(features) > 2
        [coeff, score] = pca(features');
        if size(score, 2) >= 2
            scatter(score(:,1), score(:,2));
            title('PCA可视化');
            xlabel('主成分1'); ylabel('主成分2');
            grid on;
        end
    end
    
    % 特征统计摘要
    subplot(3, 3, 8);
    featureStats = [mean(features), std(features), min(features), max(features)];
    bar(featureStats);
    title('特征统计摘要');
    set(gca, 'XTickLabel', {'均值', '标准差', '最小值', '最大值'});
    grid on;
    
    % 特征提取流水线示意图
    subplot(3, 3, 9);
    text(0.1, 0.5, sprintf('特征维度: %d\n特征类型: %d种\n音频长度: %.2fs', ...
        length(features), length(unique(regexp(strjoin(featureNames), '_', 'split'))), ...
        length(audioData)/fs), 'FontSize', 12);
    title('特征提取摘要');
    axis off;
end

批量处理和数据库构建

function featureDatabase = buildVoiceprintDatabase(dataFolder, outputFile)
% 构建声纹特征数据库
    
    % 查找所有音频文件
    audioFiles = dir(fullfile(dataFolder, '*.wav'));
    numSpeakers = length(audioFiles);
    
    featureDatabase = struct();
    
    fprintf('开始构建声纹数据库，共%d个说话人...\n', numSpeakers);
    
    for i = 1:numSpeakers
        audioFile = fullfile(dataFolder, audioFiles(i).name);
        
        % 读取音频文件
        [audioData, fs] = audioread(audioFile);
        
        % 提取特征
        [features, featureNames] = voiceprintFeatureExtraction(audioData, fs);
        
        % 存储到数据库
        [~, speakerID, ~] = fileparts(audioFiles(i).name);
        featureDatabase(i).speakerID = speakerID;
        featureDatabase(i).features = features;
        featureDatabase(i).featureNames = featureNames;
        featureDatabase(i).audioFile = audioFiles(i).name;
        featureDatabase(i).fs = fs;
        
        fprintf('完成: %s (%d/%d)\n', speakerID, i, numSpeakers);
    end
    
    % 保存数据库
    if nargin > 1 && ~isempty(outputFile)
        save(outputFile, 'featureDatabase');
        fprintf('数据库已保存到: %s\n', outputFile);
    end
    
    fprintf('声纹数据库构建完成，共%d个说话人\n', numSpeakers);
end

参考代码 声纹识别特征提取程序 www.3dddown.com/zha/79557.html

实际应用示例

示例1：单个音频文件特征提取

% 读取音频文件
[audio, fs] = audioread('speaker1.wav');

% 提取特征
[features, featureNames] = voiceprintFeatureExtraction(audio, fs, ...
    'mfccCoeffs', 13, ...
    'lpOrder', 12, ...
    'includeDelta', true);

% 可视化分析
visualizeFeatures(audio, fs, features, featureNames);

% 显示特征统计
fprintf('提取特征数: %d\n', length(features));
fprintf('特征范围: [%.4f, %.4f]\n', min(features), max(features));

示例2：批量处理构建数据库

% 构建声纹数据库
dataFolder = 'audio_database/';
featureDatabase = buildVoiceprintDatabase(dataFolder, 'voiceprint_database.mat');

% 分析数据库统计信息
numSpeakers = length(featureDatabase);
featureDims = arrayfun(@(x) length(x.features), featureDatabase);

fprintf('数据库统计:\n');
fprintf('  说话人数: %d\n', numSpeakers);
fprintf('  特征维度: %d ± %d\n', mean(featureDims), std(featureDims));

注意事项和优化建议

1. 参数调优建议

% 针对不同语音类型的优化参数
if strcmp(speechType, 'telephone')
    params.mfccCoeffs = 12; % 电话语音频带较窄
    params.lowFreq = 300;   % 调整低频截止
    params.highFreq = 3400; % 调整高频截止
elseif strcmp(speechType, 'wideband')
    params.mfccCoeffs = 16; % 宽带语音可用更多系数
    params.highFreq = 8000; % 更高频率范围
end

2. 鲁棒性增强

% 添加噪声鲁棒性处理
function enhancedFeatures = addRobustness(features, method)
    switch method
        case 'CMN'
            % Cepstral Mean Normalization
            enhancedFeatures = features - mean(features, 2);
        case 'MVN'
            % Mean-Variance Normalization
            enhancedFeatures = (features - mean(features, 2)) ./ std(features, 0, 2);
        case 'RASTA'
            % RASTA滤波
            enhancedFeatures = rastaplp(features);
    end
end

3. 实时处理优化

% 流式处理版本（用于实时应用）
function streamingFeatureExtraction()
    % 初始化缓冲区
    bufferSize = 2048;
    audioBuffer = zeros(bufferSize, 1);
    
    while hasAudioData()
        % 获取新音频数据
        newData = getAudioFrame();
        
        % 更新缓冲区
        audioBuffer = [audioBuffer(bufferSize/2+1:end); newData];
        
        % 实时特征提取
        features = voiceprintFeatureExtraction(audioBuffer, fs, ...
            'frameLength', 0.020, ... % 更短的帧长
            'frameOverlap', 0.010);
        
        % 实时更新声纹模型
        updateVoiceprintModel(features);
    end
end

高级特征扩展

1. GFCC特征（Gammatone频率倒谱系数）

function gfccFeatures = extractGFCC(frames, fs, numCoeffs)
% 提取GFCC特征 - 更适合噪声环境
    % 实现Gammatone滤波器组
    % 类似于MFCC但使用Gammatone滤波器
end

2. 深度特征提取

function deepFeatures = extractDeepFeatures(audioData, fs)
% 使用预训练的深度学习模型提取高级特征
    % 加载预训练模型
    net = load('pretrained_voice_model.mat');
    
    % 提取深度特征
    deepFeatures = predict(net, audioData);
end