杂记

#pragma once  // 头文件保护指令，作用是防止该头文件被同一个编译单元重复包含。

#include <vector> // 引入动态数组容器


class MelFeatureExtractor 
{

public:

    MelFeatureExtractor(
        int sample_rate,
        int frame_length,
        int frame_shift,
        int n_fft,
        int n_mels
    );


    // 输入PCM，输出Mel谱
    std::vector<std::vector<float>>
    Compute(
        const std::vector<float>& waveform
    ); // 返回值为std::vector<std::vector<float>> 二维数组


private:

    int sample_rate_;

    int frame_length_;

    int frame_shift_;

    int n_fft_;

    int n_mels_;


    std::vector<float> window_;  // 私有成员，构造时传入


    // Mel滤波器组
    std::vector<
        std::vector<float>
    > mel_filterbank_;



    void InitWindow();  // 生成窗函数 存入windows_


    void InitMelFilterbank(); // 生成梅尔三角滤波器组并存入mel_filterbank_


};

#include "MelFeatureExtractor.h"
#include <cmath>
#include <algorithm>

using namespace std;

// 构造函数：初始化参数 + 预计算窗函数 + 预计算梅尔滤波器
MelFeatureExtractor::MelFeatureExtractor( // 类里面的成员函数实现
    int sample_rate,
    int frame_length,
    int frame_shift,
    int n_fft,
    int n_mels
)
    : sample_rate_(sample_rate),
      frame_length_(frame_length),
      frame_shift_(frame_shift),
      n_fft_(n_fft),
      n_mels_(n_mels)
{
    InitWindow();
    InitMelFilterbank(); //生成传函数和滤波器组
}

// 初始化汉明窗 Hamming Window
void MelFeatureExtractor::InitWindow()
{
    window_.resize(frame_length_);
    for (int i = 0; i < frame_length_; ++i) {
        window_[i] = 0.54f - 0.46f * cosf(2.0f * (float)M_PI * i / (frame_length_ - 1));
    } //窗函数代码
}

// 线性频率 → 梅尔频率
static float linear_to_mel(float f)
{
    return 2595.0f * log10f(1.0f + f / 700.0f);
} // 人耳频率映射

// 梅尔频率 → 线性频率
static float mel_to_linear(float mel)
{
    return 700.0f * (powf(10.0f, mel / 2595.0f) - 1.0f);
}

// 初始化梅尔三角滤波器组
void MelFeatureExtractor::InitMelFilterbank()
{
    int num_bins = n_fft_ / 2 + 1;
    float f_max = sample_rate_ / 2.0f;

    // 梅尔轴上均匀分点
    float mel_min = linear_to_mel(0.0f);
    float mel_max = linear_to_mel(f_max);

    // 生成 n_mels + 2 个中心点
    vector<float> mel_points(n_mels_ + 2);
    for (int i = 0; i < n_mels_ + 2; ++i) {
        mel_points[i] = mel_min + (mel_max - mel_min) * i / (n_mels_ + 1);
    }

    // 转回频率
    vector<float> freq_points(n_mels_ + 2);
    for (int i = 0; i < n_mels_ + 2; ++i) {
        freq_points[i] = mel_to_linear(mel_points[i]);
    }

    // 转换成 FFT bin 索引
    vector<int> bin_points(n_mels_ + 2);
    for (int i = 0; i < n_mels_ + 2; ++i) {
        bin_points[i] = roundf(freq_points[i] * n_fft_ / sample_rate_);
    }

    // 构造每个三角滤波器
    mel_filterbank_.resize(n_mels_, vector<float>(num_bins, 0.0f));
    for (int m = 0; m < n_mels_; ++m) {
        int left = bin_points[m];
        int center = bin_points[m + 1];
        int right = bin_points[m + 2];

        // 上升沿
        for (int k = left; k < center; ++k) {
            mel_filterbank_[m][k] = (float)(k - left) / (center - left);
        }
        // 下降沿
        for (int k = center; k < right; ++k) {
            mel_filterbank_[m][k] = (float)(right - k) / (right - center);
        }
    }
}

// 最简单的复数结构用于 FFT（仅实现本代码所需）
struct Complex {
    float re, im;
    Complex(float r = 0, float i = 0) : re(r), im(i) {}
    Complex operator+(const Complex& o) const { return Complex(re+o.re, im+o.im); }
    Complex operator-(const Complex& o) const { return Complex(re-o.re, im-o.im); }
    Complex operator*(const Complex& o) const { return Complex(re*o.re - im*o.im, re*o.im + im*o.re); }
};

// 极精简的 FFT 实现（仅用于演示，工程建议用 kissfft / FFTW ）
static void fft(vector<Complex>& data, bool invert) {
    int n = data.size();
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1)
            j ^= bit;
        j ^= bit;
        if (i < j)
            swap(data[i], data[j]);
    }

    for (int len = 2; len <= n; len <<= 1) {
        float ang = 2 * M_PI / len * (invert ? -1 : 1);
        Complex wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len) {
            Complex w(1);
            for (int j = 0; j < len/2; j++) {
                Complex u = data[i+j], v = data[i+j+len/2] * w;
                data[i+j] = u + v;
                data[i+j+len/2] = u - v;
                w = w * wlen;
            }
        }
    }
}

// 计算梅尔谱
vector<vector<float>> MelFeatureExtractor::Compute(const vector<float>& waveform)
{
    int num_samples = waveform.size();
    int num_frames = (num_samples - frame_length_) / frame_shift_ + 1;

    if (num_frames <= 0)
        return {};

    vector<vector<float>> mel_frames(num_frames, vector<float>(n_mels_, 0.0f));
    int num_bins = n_fft_ / 2 + 1;

    for (int t = 0; t < num_frames; ++t) {
        // 1. 取一帧
        vector<float> frame(frame_length_);
        int start = t * frame_shift_;
        for (int i = 0; i < frame_length_; ++i)
            frame[i] = waveform[start + i];

        // 2. 加窗
        for (int i = 0; i < frame_length_; ++i)
            frame[i] *= window_[i];

        // 3. 补零到 n_fft
        vector<Complex> fft_in(n_fft_);
        for (int i = 0; i < frame_length_; ++i)
            fft_in[i] = Complex(frame[i]);

        // 4. FFT
        fft(fft_in, false);

        // 5. 计算功率谱
        vector<float> power(num_bins);
        for (int i = 0; i < num_bins; ++i) {
            float re = fft_in[i].re;
            float im = fft_in[i].im;
            power[i] = (re*re + im*im) / n_fft_;
        }

        // 6. 梅尔滤波 + 取 log
        for (int m = 0; m < n_mels_; ++m) {
            float sum = 0.0f;
            for (int i = 0; i < num_bins; ++i)
                sum += power[i] * mel_filterbank_[m][i];

            // 对数梅尔谱（加小epsilon防 log(0)）
            mel_frames[t][m] = logf(sum + 1e-6f);
        }
    }

    return mel_frames;
}