C++实现MFCC语音特征提取：从原理到工程实践

1次阅读

共计 4790 个字符，预计需要花费 12 分钟才能阅读完成。

MFCC（Mel 频率倒谱系数）是语音识别中最常用的特征表示方法之一。它模拟人耳对不同频率声音的感知特性，将语音信号转换为一组能够有效表征语音内容的特征向量。MFCC 特征提取过程通过一系列信号处理步骤，将时域语音信号转换为更紧凑且易于分类的频域表示，为后续的语音识别模型提供高质量的输入特征。

Python 实现
优点：开发快速，有现成的库（如 librosa）可直接调用；适合原型验证和算法研究
缺点：运行效率较低，不适合实时处理；内存消耗较大
C++ 实现
优点：运行效率高，适合实时语音处理；内存管理更精细
缺点：开发周期较长，需要手动实现各处理步骤

数学原理 ：
预加重通过一个一阶高通滤波器来增强语音信号的高频部分，补偿语音信号在传输过程中高频分量的衰减。

// 使用 Eigen 实现预加重
Eigen::VectorXf preEmphasis(const Eigen::VectorXf& signal, float alpha = 0.97) {Eigen::VectorXf result(signal.size());
    result[0] = signal[0];
    for (int i = 1; i < signal.size(); ++i) {result[i] = signal[i] - alpha * signal[i-1];
    }
    return result;
}

数学原理 ：
语音信号是准稳态信号，需要分帧处理（通常 20-40ms 每帧），帧移一般为帧长的一半。

// 分帧实现
std::vector<Eigen::VectorXf> framing(const Eigen::VectorXf& signal, 
                                    int frame_size, int frame_stride) {
    std::vector<Eigen::VectorXf> frames;
    int num_frames = 1 + (signal.size() - frame_size) / frame_stride;
    frames.reserve(num_frames);

    for (int i = 0; i < num_frames; ++i) {
        int start = i * frame_stride;
        frames.emplace_back(signal.segment(start, frame_size));
    }
    return frames;
}

数学原理 ：
使用汉明窗减少频谱泄漏，公式：$w(n) = 0.54 – 0.46\cos\left(\frac{2πn}{N-1}\right)$

// 加窗实现
void applyHammingWindow(Eigen::VectorXf& frame) {
    static Eigen::VectorXf window;
    if (window.size() != frame.size()) {window.resize(frame.size());
        for (int i = 0; i < frame.size(); ++i) {window[i] = 0.54f - 0.46f * std::cos(2 * M_PI * i / (frame.size() - 1));
        }
    }
    frame = frame.cwiseProduct(window);
}

数学原理 ：
快速傅里叶变换将时域信号转换为频域表示，获取频谱信息。

// 使用 FFTW 库实现 FFT
Eigen::VectorXf computeFFT(const Eigen::VectorXf& frame) {fftwf_complex* out = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex) * (frame.size()/2 + 1));
    fftwf_plan plan = fftwf_plan_dft_r2c_1d(frame.size(), 
                                          const_cast<float*>(frame.data()), 
                                          out, FFTW_ESTIMATE);
    fftwf_execute(plan);

    Eigen::VectorXf power_spectrum(frame.size()/2 + 1);
    for (int i = 0; i <= frame.size()/2; ++i) {power_spectrum[i] = out[i][0]*out[i][0] + out[i][1]*out[i][1];
    }

    fftwf_destroy_plan(plan);
    fftwf_free(out);
    return power_spectrum;
}

数学原理 ：
将线性频率转换为梅尔频率，并在梅尔尺度上均匀分布三角滤波器。

// 梅尔滤波器组实现
Eigen::MatrixXf createMelFilterBank(int num_filters, int fft_size, 
                                   int sample_rate, float low_freq, float high_freq) {
    // 将频率转换为梅尔尺度
    auto hz2mel = [](float hz) {return 2595 * std::log10(1 + hz / 700.0); };
    auto mel2hz = [](float mel) {return 700 * (std::pow(10, mel / 2595.0) - 1); };

    // 计算滤波器边界
    float low_mel = hz2mel(low_freq);
    float high_mel = hz2mel(high_freq);
    Eigen::VectorXf mel_points = Eigen::VectorXf::LinSpaced(num_filters + 2, low_mel, high_mel);

    // 转换回 Hz 并映射到 FFT bins
    Eigen::VectorXf bin_points = (mel_points.unaryExpr(mel2hz) / sample_rate * (fft_size - 1)).array().round();

    // 构建滤波器组
    Eigen::MatrixXf filter_bank = Eigen::MatrixXf::Zero(num_filters, fft_size/2 + 1);
    for (int i = 0; i < num_filters; ++i) {int left = bin_points[i];
        int center = bin_points[i+1];
        int right = bin_points[i+2];

        // 上升斜坡
        for (int j = left; j <= center; ++j) {filter_bank(i, j) = (j - left) / float(center - left);
        }

        // 下降斜坡
        for (int j = center; j <= right; ++j) {filter_bank(i, j) = (right - j) / float(right - center);
        }
    }
    return filter_bank;
}

数学原理 ：
对梅尔滤波器组输出的对数能量进行 DCT 变换，保留前 N 个系数作为 MFCC 特征。

// DCT 实现
Eigen::VectorXf applyDCT(const Eigen::VectorXf& mel_energies, int num_ceps) {int num_filters = mel_energies.size();
    Eigen::MatrixXf dct_matrix(num_ceps, num_filters);

    for (int i = 0; i < num_ceps; ++i) {for (int j = 0; j < num_filters; ++j) {dct_matrix(i, j) = std::cos(M_PI * i * (j + 0.5) / num_filters);
        }
    }

    // 归一化处理
    dct_matrix.row(0) *= 1.0 / std::sqrt(num_filters);
    dct_matrix.block(1, 0, num_ceps-1, num_filters) *= std::sqrt(2.0 / num_filters);

    return dct_matrix * mel_energies;
}

使用 Eigen 的向量化运算或显式 SIMD 指令加速关键计算：

// 使用 Eigen 向量化运算
Eigen::VectorXf computeMelEnergies(const Eigen::MatrixXf& filter_bank, 
                                  const Eigen::VectorXf& power_spectrum) {return (filter_bank * power_spectrum).array().max(1e-10).log();}

对于实时处理，使用环形缓冲区减少内存分配和拷贝：

class CircularBuffer {
public:
    CircularBuffer(size_t size) : buffer(size), head(0), tail(0), full(false) {}

    void push(float value) {buffer[head] = value;
        if (full) tail = (tail + 1) % buffer.size();
        head = (head + 1) % buffer.size();
        full = head == tail;
    }

    bool getFrame(Eigen::VectorXf& frame, int frame_size) {if (size() < frame_size) return false;

        frame.resize(frame_size);
        for (int i = 0; i < frame_size; ++i) {frame[i] = buffer[(tail + i) % buffer.size()];
        }
        return true;
    }

private:
    std::vector<float> buffer;
    size_t head, tail;
    bool full;
};

预先计算常用范围内的对数结果，减少实时计算量：

class LogTable {
public:
    LogTable(float min_val = 1e-10, float max_val = 1e6, int size = 100000) 
        : min_val(min_val), max_val(max_val), step((max_val - min_val) / size) {table.resize(size);
        for (int i = 0; i < size; ++i) {table[i] = std::log(min_val + i * step);
        }
    }

    float operator()(float x) const {if (x <= min_val) return table[0];
        if (x >= max_val) return table.back();
        int idx = static_cast<int>((x - min_val) / step);
        return table[idx];
    }

private:
    std::vector<float> table;
    float min_val, max_val, step;
};