fe_epoch.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"


static const int SkipF = 16;
static const int WinSize = 320;
static const int MinPitchInterval = 40; // org=24.
static const int MaxPitchInterval = 140;
static const float VoiceThreshold = (float)0.05; // org=0.1
static const int SamplesBetweenPeaks = 40; // For male=40. For female=20.


int Fe::calc_epoch(float *sample, int sampleN, vector<short>& epoch, short* PitchPeriod)
{
	int i,j,k;
	valarray<float> dataA(sampleN);
	valarray<int> len(3);
	valarray<float> dwt[3];
	valarray<float> a[3];

	float mean = 0;
	for(i=0;i<sampleN;i++) mean += sample[i];
	mean /= sampleN;
	for(i=0;i<sampleN;i++) dataA[i] = sample[i]-mean;	

	/* init */
	for(k=0;k<3;k++){
		len[k] = SamplesBetweenPeaks*ipower(2,k);
		dwt[k].resize(sampleN+SamplesBetweenPeaks, 0);
		a[k].resize(WinSize+1, 0);
		for(i=0;i<=WinSize;i++){
			float x = M_PI*(i-len[k]/2)/(float)len[k];
			a[k][i] = ((i <= len[k]) ? (float)((0.48829*sin(2*x)+0.28256*sin(4*x)+0.03504*sin(6*x))*cos(x)*cos(x)) : (float)0);
		}
	}

	/* dwt */
	valarray<float> maxdwt((float)0, 3);
	for(k=0; k<3; k++){
		// first pass
		maxdwt[k] = (float)0;		
		for(i=len[k];i<sampleN-len[k];i+=SkipF){
			dwt[k][i] = calc_DWT(&dataA[i-len[k]/2], &a[k][0], len[k]);
			if(dwt[k][i] > maxdwt[k]) maxdwt[k] = dwt[k][i];
		}

		// second pass
		float th = (float)0.8*maxdwt[k];
		int delta = (SkipF/4 > 0) ? SkipF/4 : 1;
		for(i=len[k];i<sampleN-len[k];i+=SkipF){
			if(dwt[k][i] < th) continue;
			for(j=i+delta;j<i+SkipF;j+=delta){
				dwt[k][j] = calc_DWT(&dataA[j-len[k]/2], &a[k][0], len[k]);
				if(dwt[k][j] > maxdwt[k]) maxdwt[k] = dwt[k][j];
			}
		}
	}

	/* peak */
	int bCanceled = 0;
	float pc0,pc1,pc2;
	vector<int> pc;
	pc.push_back(0);
	int epochN = 1;
	for(i=0;i<MinPitchInterval/2+3;i++) epoch[i] = 0;
	for(i=sampleN-3; i<sampleN; i++) epoch[i] = 0;
	epoch[0] = 10; // To show epoch clearly!
	pc0 = VoiceThreshold*maxdwt[0];
	pc1 = VoiceThreshold*maxdwt[1];
	pc2 = VoiceThreshold*maxdwt[2];
	for(i=MinPitchInterval/2+3;i<sampleN-3;i++){
		if(!CheckWinMessage()) {
			bCanceled = 1;
			break;
		}
		ShowProgress((int)((i*100)/sampleN));
		if(dwt[0][i] == 0) dwt[0][i] = calc_DWT(&dataA[i-len[0]/2], &a[0][0], len[0]);
		if(dwt[0][i-1] == 0) dwt[0][i-1] = calc_DWT(&dataA[i-1-len[0]/2], &a[0][0], len[0]);
		if(dwt[0][i+1] == 0) dwt[0][i+1] = calc_DWT(&dataA[i+1-len[0]/2], &a[0][0], len[0]);
	
		if((dwt[0][i]>=pc0) && (dwt[0][i]>=dwt[0][i-1]) && (dwt[0][i]>=dwt[0][i+1])){
			for(j=i-MinPitchInterval/2;j<=i+MinPitchInterval/2;j++){
				if(dwt[1][j] == 0) dwt[1][j] = calc_DWT(&dataA[j-len[1]/2], &a[1][0], len[1]);
				if(dwt[1][j-1] == 0) dwt[1][j-1] = calc_DWT(&dataA[j-1-len[1]/2], &a[1][0], len[1]);
				if(dwt[1][j+1] == 0) dwt[1][j+1] = calc_DWT(&dataA[j+1-len[1]/2], &a[1][0], len[1]);

				if((dwt[1][j]>=pc1) && (dwt[1][j]>=dwt[1][j-1]) && (dwt[1][j]>=dwt[1][j+1])){
					if(j-MinPitchInterval/2>=0){
						for(k=j-MinPitchInterval/2;k<=j+MinPitchInterval/2;k++){
							if(dwt[2][k] == 0) dwt[2][k] = calc_DWT(&dataA[k-len[2]/2], &a[2][0], len[2]);
							if(dwt[2][k-1] == 0) dwt[2][k-1] = calc_DWT(&dataA[k-1-len[2]/2], &a[2][0], len[2]);
							if(dwt[2][k+1] == 0) dwt[2][k+1] = calc_DWT(&dataA[k+1-len[2]/2], &a[2][0], len[2]);

							if((dwt[2][k]>=pc2) && (dwt[2][k]>=dwt[2][k-1]) && (dwt[2][k]>=dwt[2][k+1])){
								if(i-pc[epochN-1]<=MinPitchInterval){
									//printf("Too narrow\n");
									if(dwt[0][i]>dwt[0][pc[epochN-1]]){
										epoch[i]=7;	
										epoch[pc[epochN-1]]=0;
										//printf("reset t=%f,j=%d,k=%d\n",pc[epochN-1]/16000.0,j,k);
										pc[epochN-1]=i;
										//printf("t=%f,j=%d,k=%d\n",i/16000.0,j,k);
										
									}
									else{
										epoch[i]=0;
										//printf("set zero t=%f,j=%d,k=%d\n",i/16000.0,j,k);
									}
								}
								else{
									epoch[i]=7;
									pc.push_back(i);
									epochN++;
									//printf("t=%f,j=%d,k=%d\n",i/16000.0,j,k);
								}
							}
						}
					}
				}
			}
		}
		else{
			epoch[i]=0;
		}
	}

	/* pitch */
	if(PitchPeriod){
		for(i=0;i<sampleN;i++) PitchPeriod[i] = 0;
		for(i=1;i<=epochN;i++){
			if((pc[i]-pc[i-1]>=MinPitchInterval) && (pc[i]-pc[i-1]<=MaxPitchInterval))
				for(k=pc[i-1]; k<pc[i]; k++)
					PitchPeriod[k] = pc[i]-pc[i-1];
		}
	}
	return sampleN;
}


float Fe::calc_DWT(float *sample, float *a, int len)
{
	int k;
	float sum = 0;
	for(k=0;k<len/2;k++)
		sum += sample[k]*a[k];
	if(sum < 0) return (float)(-len);
	for(k=len/2;k<=len;k++)
		sum += sample[k]*a[k];
	return sum/100;
}