梅爾頻率倒譜系數（MFCC）的提取過程與C++程式碼實現
阿新 • • 發佈：2019-02-11
MFCC引數提取步驟
——>預加重
——>分幀
——>對每一幀加窗
——>對每一幀補零
——>各幀訊號的FFT變換及其功率譜
——>梅爾濾波（通過40個濾波器）
——>取對數
——>DCT變換
——>歸一化
1.預加重
如果資料在低頻的強度大於高頻，就會不利於處理，因此需要通過一個傳遞函式為s[n]-a*s[n]的高通濾波器。濾去資料中的低頻成分，使高頻特性更加突現。
2.分幀
分幀就是將N個取樣點集合成一個觀測單位。我們設定每幀涵蓋的時間是25ms，因為取樣率是16000，所以得到每幀的樣本點個數是400。
另外，為了避免相鄰兩幀的變化過大，因此會讓兩相鄰幀之間有一段重疊區域。我們設定的重疊區域是15ms，所以就是每隔10ms取一幀。
3.對每一幀加窗
分幀後馬上進行FFT，由於轉換時會將幀內訊號當作週期訊號處理，所以在幀的兩個端點處會發生突變，轉換出來的頻譜與原訊號頻譜差別很大。所以要對每一幀加窗，使幀內訊號作FFT時的兩個端點處不會發生突變。
我們採用的窗是漢明窗：（M為幀長，即400）
4.對每一幀補零
我們要對每一幀訊號進行FFT，而FFT要求輸入資料長度一定是2^K，現在一幀為400個取樣點，所以補零至最接近的512位。
5.各幀訊號的FFT變換及其功率譜
對分幀加窗後的各幀訊號進行512點的FFT變換得到各幀的頻譜。並對語音訊號的頻譜取模平方得到語音訊號的功率譜。
6.梅爾濾波（通過40個濾波器）
40個三角濾波器在MEL譜上均勻分佈，每兩個濾波器間有50%的重疊部分。
所以要先把實際頻率轉換成梅爾頻率，實際頻率最小為0Hz，最大為16000 / 2 = 8000Hz
轉換成梅爾頻率後，我們要實現的是40個濾波器，所以計算這40個濾波器的梅爾頻率分佈，然後把梅爾頻率轉換成實際頻率
然後根據實際頻率的陣列，計算出f陣列，公式如下（h陣列就是40個濾波器的實際頻率分佈陣列）：
最後根據以下公式，計算濾波器的輸出（m為濾波器的個數）
以下為10個濾波器的圖：
7.取對數
三角窗濾波器組的輸出求取對數，可以得到近似於同態變換的結果。
8.DCT變換
對對數能量梅爾譜進行DCT變換，取前13維輸出，得到梅爾倒譜。
公式為：
9.歸一化
對所有的梅爾倒譜歸一化。
先求出所有倒譜向量的均值向量，公式為
再用每一個倒譜向量減去均值向量，即
C++程式碼實現
#include<iostream>
#include<fstream>
#include<complex>
using namespace std;
int filterNum = 40;
int sampleRate = 16000;
#define Win_Time 0.025//把25ms裡的所有點作為一個點分析
#define Hop_Time 0.01//每隔10ms分一次幀
#define Pi 3.1415927
//1.預加重
double* pre_emphasizing(double *sample, int len, double factor)
{
double *Sample = new double[len];
Sample[0] = sample[0];
for(int i = 1; i < len; i++)
{
//預加重過程
Sample[i] = sample[i] - factor * sample[i - 1];
}
return Sample;
}
void Hamming( double *hamWin, int hamWinSize )
{
for (int i = 0; i < hamWinSize; i++)
{
hamWin[i] = (double)(0.54 - 0.46 * cos(2 * Pi * (double)i / ((double)hamWinSize - 1) ));
}
}
//計算每一幀的功率譜
void mfccFFT(double *frameSample, double *FFTSample, int frameSize, int pos)
{
//對分幀加窗後的各幀訊號進行FFT變換得到各幀的頻譜
//並對語音訊號的頻譜取模平方得到語音訊號的功率譜
double dataR[frameSize];
double dataI[frameSize];
for(int i = 0; i < frameSize; i++)
{
dataR[i] = frameSample[i + pos];
dataI[i] = 0.0f;
}
int x0, x1, x2, x3, x4, x5, x6, xx, x7, x8;
int i, j, k, b, p, L;
float TR, TI, temp;
/********** following code invert sequence ************/
for(i = 0; i < frameSize; i++)
{
x0 = x1 = x2 = x3 = x4 = x5 = x6 = x7 = x8 = 0;
x0 = i & 0x01; x1 = (i / 2) & 0x01; x2 = (i / 4) & 0x01; x3 = (i / 8) & 0x01; x4 = (i / 16) & 0x01;
x5 = (i / 32) & 0x01; x6 = (i / 64) & 0x01; x7 = (i / 128) & 0x01; x8 = (i / 256) & 0x01;
xx = x0 * 256 + x1 * 128 + x2 * 64 + x3 * 32 + x4 * 16 + x5 * 8 + x6 * 4 + x7 * 2 + x8;
dataI[xx] = dataR[i];
}
for(i = 0; i < frameSize; i++)
{
dataR[i] = dataI[i]; dataI[i] = 0;
}
/************** following code FFT *******************/
for(L = 1; L <= 9; L++)
{ /* for(1) */
b = 1; i = L - 1;
while(i > 0)
{
b = b * 2; i--;
} /* b= 2^(L-1) */
for(j = 0; j <= b-1; j++) /* for (2) */
{
p = 1; i = 9 - L;
while(i > 0) /* p=pow(2,7-L)*j; */
{
p = p * 2; i--;
}
p = p * j;
for(k = j; k < 512; k = k + 2*b) /* for (3) */
{
TR = dataR[k]; TI = dataI[k]; temp = dataR[k + b];
dataR[k] = dataR[k] + dataR[k + b] * cos(2 * Pi * p / frameSize) + dataI[k + b] * sin(2 * Pi * p / frameSize);
dataI[k] = dataI[k] - dataR[k + b] * sin(2 * Pi * p / frameSize) + dataI[k + b] * cos(2 * Pi * p / frameSize);
dataR[k + b] = TR - dataR[k + b] * cos(2 * Pi * p / frameSize) - dataI[k + b] * sin(2 * Pi * p / frameSize);
dataI[k + b] = TI + temp * sin(2 * Pi * p / frameSize) - dataI[k + b] * cos(2 * Pi * p / frameSize);
} /* END for (3) */
} /* END for (2) */
} /* END for (1) */
for(i = 0; i < frameSize / 2; i++)
{
FFTSample[i + pos] = (dataR[i] * dataR[i] + dataI[i] * dataI[i]);
}
}
//引數說明：frameSample為處理之後的陣列，Sample為對樣本進行預加重之後的陣列
// len為Sample的長度，frameSize為每幀的樣本點個數，frameSampleLen為處理之後的長度
double* mfccFrame(double *frameSample, double *Sample, int *len, int frameSize, int &frameSampleLen)
{
double *hamWin;
int hamWinSize = sampleRate * Win_Time;
hamWin = new double[hamWinSize];
Hamming(hamWin, hamWinSize);//計算hamming窗
int hopStep = Hop_Time * sampleRate;
int frameNum = ceil(double(*len) / hopStep);//計算一共會有多少幀
frameSampleLen = frameNum * frameSize;//經過處理之後的長度
frameSample = new double[frameSampleLen];
for(int i = 0; i < frameSampleLen; i++)
frameSample[i] = 0;
double *FFTSample = new double[frameSampleLen];
for(int i = 0; i < frameSampleLen; i++)
FFTSample[i] = 0;
for(int i = 0; i * hopStep < *len; i++)//分幀
{
for(int j = 0; j < frameSize; j++)
{
if(j < hamWinSize && i * hopStep + j < *len)
frameSample[i * frameSize + j] = Sample[i * hopStep + j] * hamWin[j];
else
frameSample[i * frameSize + j] = 0;//補0
}
mfccFFT(frameSample, FFTSample, frameSize, i * frameSize);
}
ofstream fileFrame("C:\\Users\\john\\Desktop\\txt\\Frame.txt");
for(int i = 0; i < frameSize; i++)
fileFrame << frameSample[100 * frameSize + i] << endl;
delete []hamWin;
return FFTSample;
}
void DCT(double mel[400][40], double c[400][40], int frameNum)
{
for(int k = 0; k < frameNum; k++)
{
for(int i = 0; i < 13; i++)
{
for(int j = 0; j < filterNum; j++)
{
c[k][i] += mel[k][j] * cos(Pi * i / (2 * filterNum) * (2 * j + 1));
//if(k == 0 && i ==0)
//cout << c[0][0] << endl;
}
}
}
//cout << "c[0][0] = " << c[0][0] << endl;
}
void computeMel(double mel[400][40], int sampleRate, double *FFTSample, int frameNum, int frameSize)
{
double freMax = sampleRate / 2;//實際最大頻率
double freMin = 0;//實際最小頻率
double melFremax = 1125 * log(1 + freMax / 700);//將實際頻率轉換成梅爾頻率
double melFremin = 1125 * log(1 + freMin / 700);
double melFilters[filterNum][3];
double k = (melFremax - melFremin) / (filterNum + 1);
double *m = new double[filterNum + 2];
double *h = new double[filterNum + 2];
double *f = new double[filterNum + 2];
for(int i = 0; i < filterNum + 2; i++)
{
m[i] = melFremin + k * i;
h[i] = 700 * (exp(m[i] / 1125) - 1);//將梅爾頻率轉換成實際頻率
f[i] = floor((frameSize + 1) * h[i] / sampleRate);
}
delete[] m;
delete[] h;
//delete[] f;
for(int k = 0; k < frameNum; k++)
{
for(int j = 0; j < filterNum; j++)
mel[k][j] = 0;
}
//計算出每個三角濾波器的輸出: 對每一幀進行處理
for(int i = 0; i < frameNum; i++)
{
for(int j = 1; j <= filterNum; j++)
{
double temp = 0;
for(int z = 0; z < frameSize; z++)
{
if(z < f[j - 1])
temp = 0;
else if(z >= f[j - 1] && z <= f[j])
temp = (z - f[j - 1]) / (f[j] - f[j - 1]);
else if(z >= f[j] && z <= f[j + 1])
temp = (f[j + 1] - z) / (f[j + 1] - f[j]);
else if(z > f[j + 1])
temp = 0;
mel[i][j - 1] += FFTSample[i * frameSize + z] * temp;
}
}
}
ofstream fileMel("C:\\Users\\john\\Desktop\\txt\\Mel.txt");
for(int i = 1; i <= filterNum; i++)
fileMel << mel[100][i] << endl;
//取對數
for(int i = 0; i < frameNum; i++)
{
for(int j = 0; j < filterNum; j++)
{
if(mel[i][j] <= 0.00000000001 || mel[i][j] >= 0.00000000001)
mel[i][j] = log(mel[i][j]);
}
}
}
void writeToFile(int frameNum, int frameSize, int hopStep, double DCT[400][40], double *sample, double *Sample, double *frameSample, double *FFTSample, double mel[400][40])
{
ofstream fileDCT("C:\\Users\\john\\Desktop\\txt\\DCT.txt");
ofstream filesample("C:\\Users\\john\\Desktop\\txt\\sample.txt");
ofstream filePreemphasized("C:\\Users\\john\\Desktop\\txt\\Preemphasized.txt");
ofstream fileFft("C:\\Users\\john\\Desktop\\txt\\Fft.txt");
ofstream fileLogMel("C:\\Users\\john\\Desktop\\txt\\LogMel.txt");
for(int i = 0; i < frameSize; i++)
{
filesample << sample[100 * hopStep + i] << endl;
filePreemphasized << Sample[100 * hopStep + i] << endl;
fileFft << FFTSample[100 * frameSize + i] << endl;
}
for(int i = 0; i < filterNum; i++)
fileLogMel << mel[100][i] << endl;
//歸一化
/* for(int i = 0; i < 13; i++)
{
double sum = 0.0f;
double Mrecording = 0.0f;
for(int j = 0; j < frameNum; j++)
{
sum = sum + DCT[j][i];
}
Mrecording = sum / frameNum;
cout << Mrecording << endl;
for(int j = 0; j < frameNum; j++)
{
DCT[j][i] = abs(DCT[j][i] - Mrecording);
}
} */
for (int i = 0; i < 13; i++)//write DCT
{
for (int j = 0; j < frameNum; j++)
fileDCT << DCT[j][i] << " ";
fileDCT << endl;
}
}
//MFCC
void MFCC(double *sample, int len)
{
double factor = 0.95;//預加重引數
double *Sample;
//1.預加重
Sample = pre_emphasizing(sample, len, factor);
//計算出每幀有多少個點，然後算出最接近點的個數的2^k，使得每幀的點的個數為2^k，以便進行補0
int frameSize = (int)pow(2, ceil( log(Win_Time * sampleRate) / log(2.0) ));
double *frameSample = NULL, *FFTSample = NULL;
int frameSampleLen;
//分幀、加窗、功率譜
FFTSample = mfccFrame(frameSample, Sample, &len, frameSize, frameSampleLen);
int hopStep = Hop_Time * sampleRate;//隔hopStep個樣本點分一次幀
int frameNum = ceil(double(len) / hopStep);//計算所有樣本點一共有多少幀
double mel[400][40];
computeMel(mel, sampleRate, FFTSample, frameNum, frameSize);
double c[400][40];
for(int i = 0; i < 400; i++)
{
for(int j = 0; j < 40; j++)
c[i][j] = 0;
}
DCT(mel, c, frameNum);
writeToFile(frameNum, frameSize, hopStep, c, sample, Sample, frameSample, FFTSample, mel);
}
int main()
{
ifstream filedata("C:\\Users\\john\\Desktop\\txt\\Data.txt");
int len = 0;
//讀取wav檔案的數值
double *sample = new double[160000];
while(!filedata.eof())
{
filedata >> sample[len];
sample[len] = sample[len] * 1000;
len++;
}
cout << len << endl;
MFCC(sample, len);
delete [] sample;
return 0;
} Matlab畫出結果圖
1.feature矩陣的色圖
clc;
clear all;
original = importdata('C:\\Users\\john\\Desktop\\txt\\Data.txt');
figure;
plot(original * 1000);
DCT = importdata('C:\\Users\\john\\Desktop\\txt\\DCT.txt');
figure;
imagesc(DCT);
colorbar;
xlabel('cepstrum after DCT');
2.單幀資料變化圖
sample = importdata('C:\\Users\\john\\Desktop\\txt\\sample.txt');
subplot(2, 3, 1),
plot(sample);
axis([0, 400, -60, 60]);
xlabel('400 sample segment from 16kHz signal');
yujiazhong = importdata('C:\\Users\\john\\Desktop\\txt\\Preemphasized.txt');
subplot(2, 3, 2),
plot(yujiazhong);
axis([0, 400, -20, 20]);
xlabel('preemphasized');
frame = importdata('C:\\Users\\john\\Desktop\\txt\\Frame.txt');
subplot(2, 3, 3),
plot(frame);
axis([0, 400, -25, 25]);
xlabel('windowed')
fft = importdata('C:\\Users\\john\\Desktop\\txt\\Fft.txt');
subplot(2, 3, 4),
max(fft)
plot(fft);
axis([0, 256, 0, 1000000]);
xlabel('Power spectrum');
mel = importdata('C:\\Users\\john\\Desktop\\txt\\Mel.txt');
subplot(2, 3, 5),
plot(mel);
xlabel('40 point Mel spectrum');
duishu = importdata('C:\\Users\\john\\Desktop\\txt\\LogMel.txt');
subplot(2, 3, 6),
plot(duishu);
axis([0, 40, 0, 15]);
xlabel('Log Mel spectrum');
參考資料：
最後，引用龍應臺的一句話作為總結：
”有些事，只能一個人做；有些關，只能一個人過；有些路啊，只能一個人走。“