fe_mfcc.cpp

这是一个语音特征提取的程序源码

///////////////////////////////////////////////////////////////////////////////
// This is a part of the Feature program.
// Version: 1.0
// Date: February 22, 2003
// Programmer: Oh-Wook Kwon
// Copyright(c) 2003 Oh-Wook Kwon. All rights reserved. [e-mail]
///////////////////////////////////////////////////////////////////////////////

#include "StdAfx.h"
#include "FE_feature.h"


int	Fe::mel_cepstrum_basic(short *sample, int frameSize, float *mel_cep, int ceporder, int fft_size)
{
	vector<float> frameA(frameSize);
	preprocessing(sample, frameSize, &frameA[0]);
	m_window.Windowing(&frameA[0], frameSize, WIN_HAMMING);
	_mel_cepstrum_basic(&frameA[0], frameSize, mel_cep, m_fbOrder, ceporder, fft_size);
	return 1;
}


int Fe::_mel_cepstrum_basic(float *sample, int frameSize, float *mel_cep, int fborder, int ceporder, int fft_size)
{
	int i;
	if(m_dctMatrix.size() != ceporder*fborder || ceporder != m_cepOrder){
		m_dctMatrix.resize((ceporder+1)*fborder);
		MfccInitDCTMatrix (&m_dctMatrix[0], ceporder, fborder);
		m_idctMatrix.resize((ceporder+1)*fborder);
		MfccInitIDCTMatrix (&m_idctMatrix[0], ceporder, fborder);
	}
	vector<float> filter_bank(m_fbOrder);
	_filterbank_basic(sample, frameSize, &filter_bank[0], fborder, fft_size, 0, 0);
	MfccDCT (&filter_bank[0], &m_dctMatrix[0], ceporder, fborder, mel_cep);

	/* scale down to be consistent with other kinds of cepstrum coefficients */
	float f=fborder/(float)fft_size;
	for(i=0;i<=ceporder;i++) mel_cep[i]*=f;

	return 1;
}


/* Supporting routine for MFCC */
#define MfccRound(x) ((int)((x)+0.5))
void Fe::MfccInitMelFilterBanks (float startingFrequency, float samplingRate, int fftLength, int numChannels)
{
    int i, k;
    vector<float> freq(numChannels+2);
	vector<int> bin(numChannels+2);
	float start_mel, fs_per_2_mel;
	
    /* Constants for calculation */
	freq[0] = startingFrequency;
    start_mel = (float)(2595.0 * log10 (1.0 + startingFrequency / 700.0));
	bin[0] = MfccRound(fftLength * freq[0] / samplingRate);
	freq[numChannels+1] = (float)(samplingRate / 2.0);
    fs_per_2_mel = (float)(2595.0 * log10 (1.0 + (samplingRate / 2.0) / 700.0));
	bin[numChannels+1] = MfccRound(fftLength * freq[numChannels+1] / samplingRate);
	
	/* Calculating mel-scaled frequency and the corresponding FFT-bin */
    /* number for the lower edge of the band                          */
	for (k = 1; k <= numChannels; k++) {
        freq[k] = (float)(700 * (pow (10, (start_mel + (float) k / (numChannels + 1) * (fs_per_2_mel - start_mel)) / 2595.0) - 1.0));
		bin[k] = MfccRound(fftLength * freq[k] / samplingRate);
	}
	
	/* This part is never used to compute MFCC coefficients */
	/* but initialized for completeness                     */
	m_MelWeight.resize(fftLength/2+1);
	for(i = 0; i<bin[0]; i++){
		m_MelWeight[i]=0;
	}
	m_MelWeight[fftLength/2]=1;
	
	/* Initialize low, center, high indices to FFT-bin */
	m_MelFB.resize(numChannels);
	for (k = 0; k <= numChannels; k++) {
		if(k<numChannels){
			m_MelFB[k].m_lowX=bin[k];
			m_MelFB[k].m_centerX=bin[k+1];
			m_MelFB[k].m_highX=bin[k+2];
		}
		for(i = bin[k]; i<bin[k+1]; i++){
			m_MelWeight[i]=(i-bin[k]+1)/(float)(bin[k+1]-bin[k]+1);
		}
    }

	m_MelCenterFreq.resize(numChannels+2);
	m_MelCenterIdx.resize(numChannels+2);
	for(k=0;k<=numChannels+1;k++){
		m_MelCenterFreq[k]=freq[k];
		m_MelCenterIdx[k]=bin[k];
	}
	
    return;
}


/* Be careful that the assumed ordering is [c0 c1 c2 ... c12] */
int Fe::MfccInitDCTMatrix (float *dctMatrix, int ceporder, int numChannels)
{
    int i, j;
    for (i = 0; i <= ceporder; i++){
        for (j = 0; j < numChannels; j++){
            dctMatrix[i * numChannels + j] = (float) cos (M_PI * (float) i / (float) numChannels * ((float) j + 0.5));
			if(i==0) dctMatrix[i * numChannels + j]*=(float)sqrt(1/(float)numChannels);
			else     dctMatrix[i * numChannels + j]*=(float)sqrt(2/(float)numChannels);
		}
	}
	return 1;
}


/* Be careful that the assumed ordering is [c0 c1 c2 ... c12] */
int Fe::MfccInitIDCTMatrix (float *idctMatrix, int ceporder, int numChannels)
{
    int i, j;
	for (j = 0; j < numChannels; j++){
		for (i = 0; i <= ceporder; i++){
            idctMatrix[j * (ceporder+1) + i] = (float) cos (M_PI * (float) i / (float) numChannels * ((float) j + 0.5));
			if(i==0) idctMatrix[j * (ceporder+1) + i]*=(float)sqrt(1/(float)numChannels);
			else     idctMatrix[j * (ceporder+1) + i]*=(float)sqrt(2/(float)numChannels);
		}
	}
	return 1;
}


void Fe::MfccDCT (float *x, float *dctMatrix, int ceporder, int numChannels, float *mel_cep)
{
    int i, j;
    for (i = 0; i <= ceporder; i++) {
        mel_cep[i] = 0.0;
        for (j = 0; j < numChannels; j++){
            mel_cep[i] += x[j] * dctMatrix[i * numChannels + j];
		}
    }	
    return;
}


void Fe::MfccIDCT (float *mel_cep, float *idctMatrix, int ceporder, int numChannels, float *x)
{
    int i, j;
	for (j = 0; j < numChannels; j++) {
		x[j] = 0;
		for (i = 0; i <= ceporder; i++) {
            x[j] += mel_cep[i] * idctMatrix[j * (ceporder+1) + i];
		}
    }	
    return;
}


void Fe::MfccMelFilterBank (float *sigFFT, int numChannels, float* output, int normalize)
{
    float sum, wsum;
    int i, k;
	
    for (k=0;k<numChannels;k++){
		MfccMelFB *melFB=&m_MelFB[k];
        sum = sigFFT[melFB->m_centerX];
		wsum=1;
        for (i = melFB->m_lowX; i < melFB->m_centerX; i++){
            sum += m_MelWeight[i] * sigFFT[i];
			wsum += m_MelWeight[i];
		}
		for (i = melFB->m_centerX+1; i <= melFB->m_highX; i++){
            sum += (1 - m_MelWeight[i-1]) * sigFFT[i];
			wsum += (1 - m_MelWeight[i-1]);
		}
        output[k] = sum;
		if(normalize) {
			assert(wsum>0);
			output[k] /= wsum;
		}
    }
    return;
}