MFCC音频特征提取流程demo

MFCC的核心思想是模仿人耳对声音的非线性感知特性（人耳对低频声音比高频声音更敏感），从而更好地表示语音特征。其流程可以概括为以下步骤：

分帧

• 目的：语音信号是非平稳信号，其特性是随时间变化的。但在一个非常短的时间段内，可以近似认为是平稳的。分帧就是将长信号切分成许多短片段来分析。
• 操作：使用一个固定长度的窗口以一定的帧移沿着信号滑动并截取数据。

加窗

• 目的：减少每一帧信号在其两端处的频谱泄漏，使帧两端平滑地衰减到零，从而降低后续傅里叶变换后旁瓣的强度。
• 操作：将每一帧信号乘上一个窗函数。

傅里叶变换和功率谱计算

• 目的：将信号从时域转换到频域。
• 操作：
  • 对每一帧加窗后的信号进行FFT，得到复数频谱。
  • 然后计算其功率谱（取模的平方）。P = |FFT(frame)|² / N （N是FFT点数）

梅尔滤波器组滤波

• 目的：模拟人耳的听觉特性。人耳对于不同频率的感知能力是不同的，在低频区域区分度高，在高频区域区分度低。梅尔刻度是一种将实际频率转换为更符合人耳感知的频率刻度。
• 操作：
  • 定义一组三角带通滤波器（梅尔滤波器组），这些滤波器在梅尔刻度上是等宽的，但在线性频率刻度上是不等宽的（低频区域滤波器窄且密集，高频区域宽且稀疏）。
  • 将上一步得到的功率谱与每一个梅尔滤波器进行积分（实际上是相乘并求和），得到每个滤波器输出的能量。这样就将频谱压缩为n个更能反映人耳听觉特性的滤波器组能量。

离散余弦变换

• 目的：对梅尔滤波器组能量取对数后，进行离散余弦变换 (DCT)。这一步也称为倒谱分析。

  • 取对数：因为人耳对声音强度的感知也是对数的。同时，对数操作可以将卷积信号（声源激励和声道滤波）转化为可加的。

  • DCT：对取对数后的滤波器组能量进行DCT，相当于对频谱进行“压缩”。DCT的结果中，低维（前几个）系数代表了包络信息（反映了声道形状，是语音内容的关键），而高维系数代表了细节信息（如激励源）。

最终的输出通常只保留DCT后的前12-13个系数，再加上一个能量项（第0项或直接从时域帧中计算），共同构成最终的MFCC特征向量。

这个流程为音频信号的每一帧都计算出了一个固定维度的MFCC特征向量，这些向量序列就可以用于后续的语音识别、说话人识别等机器学习任务。

 

下面是在一个项目中使用MFCC方法进行特征提取的预处理部分，这里用的不是标准MFCC，省去了最后的DCT部分，更直接简单。这个处理流程是在python仿真后迁移的，python中运用了librosa库，不过该库里面有一些隐式处理，因此在转换成C版本代码时遇到了较多不匹配的情况，不过调试中也是逐步改进了。

分帧函数如下，实际使用时调整传参使其能够保持60%的重合率，让每一帧平滑过渡，不过这个函数只计算帧数，真正的分帧操作在下面的流式处理里面循环执行。

// 分帧函数，使用局部缓存处理每一帧数据
int frame_signal(const int16_t* y, int y_len, int frame_length, int hop_length, int* n_frames) {
    *n_frames = 1 + (y_len - frame_length) / hop_length;
    return 0; // 假设不会出错
}
//n_frames 经过处理后为实际帧数
//实际运行时，hoplength和framlength定义如下
// #define WIN_LENGTH 400 窗长400
// #define HOP_LENGTH 160 帧移160

 

手动FFT实现，经典的Cooley-Tukey迭代FFT算法，采用了位反转置换和蝶形运算：

// Complex结构体定义
typedef struct {
    float real;
    float imag;
} Complex;


void fft(Complex* x) {
    // 预先计算旋转因子
    static Complex W[N_FFT/2];
    static int initialized = 0;
    if (!initialized) {
        for (int k = 0; k < N_FFT/2; k++) {
            float angle = -2 * PI * k / N_FFT;
            W[k].real = cosf(angle);
            W[k].imag = sinf(angle);
        }
        initialized = 1;
    }

    // 位反转置换
    int j = 0;
    for (int i = 0; i < N_FFT - 1; i++) {
        if (i < j) {
            Complex temp = x[i];
            x[i] = x[j];
            x[j] = temp;
        }
        int k = N_FFT >> 1;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }

    // 迭代FFT
    for (int stride = 2; stride <= N_FFT; stride <<= 1) {
        int half_stride = stride >> 1;
        int W_step = N_FFT / stride;
        
        for (int k = 0; k < N_FFT; k += stride) {
            for (int m = 0; m < half_stride; m++) {
                int idx_e = k + m;
                int idx_o = idx_e + half_stride;
                int W_idx = m * W_step;
                
                Complex t = {
                    W[W_idx].real * x[idx_o].real - W[W_idx].imag * x[idx_o].imag,
                    W[W_idx].real * x[idx_o].imag + W[W_idx].imag * x[idx_o].real
                };
                
                x[idx_o].real = x[idx_e].real - t.real;
                x[idx_o].imag = x[idx_e].imag - t.imag;
                x[idx_e].real += t.real;
                x[idx_e].imag += t.imag;
            }
        }
    }
}

 

Mel功率谱计算，函数中的汉明窗与Mel滤波器组都是事先生成的。

// 计算Mel功率谱
void manual_melspectrogram(const int16_t* y, int y_len, int sr,
                           float* mel_spec, int* n_mels_out,
                           int* n_frames_out) {
    int n_frames;

    // 计算最大绝对值（用于归一化）
    float max_val = 0.0f;
    for (int i = 0; i < y_len; i++) {
        float abs_val = fabsf((float)y[i]);
        if (abs_val > max_val) {
            max_val = abs_val;
        }
    }

    // 分配并归一化为 float 缓冲区
    static float y_norm[8000];  // 假设最大长度8000
    if (y_len > 8000) {
        *n_frames_out = 0;
        return;  // 防止越界
    }
    if (max_val > 0.0f) {
        for (int i = 0; i < y_len; i++) {
            y_norm[i] = (float)y[i] / max_val;
        }
    } else {
        for (int i = 0; i < y_len; i++) {
            y_norm[i] = 0.0f;
        }
    }

    // 计算帧数
    frame_signal(y, y_len, WIN_LENGTH, HOP_LENGTH, &n_frames);  
    // 注意这里使用原 y，因为只用于 n_frames 计算
    *n_mels_out = N_MELS;
    *n_frames_out = n_frames;

    // 建立FFT缓冲区
    static Complex fft_input[N_FFT];
    static float power_spectrum[MEL_BINS];

    for (int t = 0; t < n_frames; t++) {
        // 加窗处理（改为使用 y_norm）
        for (int n = 0; n < WIN_LENGTH; n++) {
            fft_input[n].real = y_norm[t * HOP_LENGTH + n] * HAMMING_WINDOW[n];
            fft_input[n].imag = 0.0f;
        }
        // 清零剩余部分
        for (int n = WIN_LENGTH; n < N_FFT; n++) {
            fft_input[n].real = 0.0f;
            fft_input[n].imag = 0.0f;
        }

        // FFT计算
        fft(fft_input);

        // 计算功率谱 (直接存储到power_spectrum)
        for (int k = 0; k < MEL_BINS; k++) {
            float re = fft_input[k].real / N_FFT;
            float im = fft_input[k].imag / N_FFT;
            power_spectrum[k] = re * re + im * im;
        }

        // 梅尔滤波 (直接计算到输出mel_spec)
        for (int m = 0; m < N_MELS; m++) {
            float mel_energy = 0.0f;
            for (int k = 0; k < MEL_BINS; k++) {
                mel_energy += power_spectrum[k] * mel_filter[m][k];
            }
            mel_spec[m * n_frames + t] = mel_energy;
        }
    }
}

 

最终处理函数，这里是将上面流程输出的值每帧取均值后放入后续神经网络的input。

#include "algo.h"
#include "MFCC.h"
#include <stdint.h>
#define FRAME_LENGTH 8000
#define N_MELS 40

// 静态缓冲区，用于存储帧、梅尔频谱和对数梅尔特征
//static float frame_global[FRAME_LENGTH];
float mel_spec_global[N_MELS * MAX_FRAMES];
//static float mean_log_mel_global[N_MELS];
int speaker_recognition_run(int16_t* pcm_data, int len) {
    // 确保只有一帧数据
    if (len != FRAME_LENGTH) {
        printf("lengtherror: need %d ，current %d\n", FRAME_LENGTH, len);
        return -1;  // 如果数据长度不匹配，返回错误
    }
    // 计算梅尔频谱
    int sr = 8000;
    int n_mels, n_frames;
    // 直接在全局缓冲区上操作，避免分配新内存
    manual_melspectrogram(pcm_data, len, sr, mel_spec_global, &n_mels, &n_frames);
    // 计算对数梅尔特征
    float* log_mel_features = compute_log_mel_mean(mel_spec_global, n_mels, n_frames);
    if (!log_mel_features) {
        // 处理错误
        printf("log error\n");
        return 1;
    }

    // Prepare input (1 sample, 40 features)
    float input[BATCH_SIZE * SEQ_LENGTH * INPUT_DIM];
    for (int i = 0; i < INPUT_DIM; i++) {
        input[i] = log_mel_features[i];
    }

    free(log_mel_features);

    // Output array for the classification results
    float output[BATCH_SIZE * NUM_CLASSES];

    // Forward pass through the model
    forward(input, output, BATCH_SIZE, SEQ_LENGTH);

    // FindMax
    float max = output[0];
    int out = 0;
    for (int j = 1; j < NUM_CLASSES; j++) {
        if (max < output[j]) {
            max = output[j];  // 更新最大值
            out = j;           // 更新输出的索引
        }
    }

    return out;
}

其中compute_log_mel_mean函数如下定义：

// 计算对数梅尔特征并取均值
float* compute_log_mel_mean(float* mel_spec, int n_mels, int n_frames) {
    if (n_mels < 40) {
        fprintf(stderr, "The audio is too short to extract Mel features.\n");
        return NULL;
    }

    float* mean_log_mel = (float*)malloc(40 * sizeof(float));  // 保留40个梅尔频率带

    // 计算每个梅尔频率带的对数均值
    for (int i = 0; i < n_mels; ++i) {
        float sum = 0.0f;
        for (int j = 0; j < n_frames; ++j) {
            float value = mel_spec[i * n_frames + j]; // 按列优先排列
            value = log10f(value + 1e-6);  // 加上一个小常数以避免对数的零值
            sum += value;
        }
        mean_log_mel[i] = sum / n_frames * SCALE_FACTOR;  // 按帧数取均值并进行缩放
    }

    return mean_log_mel;
}

其中SCALE_FACTOR 10.0 是比例因子，用于补偿与 Python 输出的差异。

当时仿真效果还可以，就没有再考虑添加DCT转换为标准MFCC，主打能跑就行。