fe_sigproc.c

CMU大名鼎鼎的SPHINX－3大词汇量连续语音识别系统

/* ====================================================================
 * Copyright (c) 1996-2004 Carnegie Mellon University.  All rights 
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer. 
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 * This work was supported in part by funding from the Defense Advanced 
 * Research Projects Agency and the National Science Foundation of the 
 * United States of America, and the CMU Sphinx Speech Consortium.
 *
 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * ====================================================================
 *
 */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "fe.h"
#include "fe_internal.h"


/*
  31 Jan 00 mseltzer - changed rounding of filter edges to -not- use 
                        rint() function. 
   3 Dec 99 mseltzer - corrected inverse DCT-2 
                        period is 1/NumFilts not 1/(2*NumFilts)
			added "beta" factor in summation
	  	     - changed mel filter bank construction so that 
                        left,right,center freqs are rounded to DFT 
                        points before filter is constructed
  
*/
                        

int32 fe_build_melfilters(melfb_t *MEL_FB)
{    
    int32 i, whichfilt, start_pt;
    float32 leftfr, centerfr, rightfr, fwidth, height, *filt_edge;
    float32 melmax, melmin, dmelbw, freq, dfreq, leftslope,rightslope;

    /*estimate filter coefficients*/
    MEL_FB->filter_coeffs = (float32 **)fe_create_2d(MEL_FB->num_filters, MEL_FB->fft_size, sizeof(float32));
    MEL_FB->left_apex = (float32 *) calloc(MEL_FB->num_filters,sizeof(float32));
    MEL_FB->width = (int32 *) calloc(MEL_FB->num_filters,sizeof(int32));
    
    if (MEL_FB->doublewide==ON)
	filt_edge = (float32 *) calloc(MEL_FB->num_filters+4,sizeof(float32));
    else	
	filt_edge = (float32 *) calloc(MEL_FB->num_filters+2,sizeof(float32));

    if (MEL_FB->filter_coeffs==NULL || MEL_FB->left_apex==NULL || MEL_FB->width==NULL || filt_edge==NULL){
	fprintf(stderr,"memory alloc failed in fe_build_mel_filters()\n...exiting\n");
	exit(0);
    }
    
    dfreq = MEL_FB->sampling_rate/(float32)MEL_FB->fft_size;
    
    melmax = fe_mel(MEL_FB->upper_filt_freq);
    melmin = fe_mel(MEL_FB->lower_filt_freq);
    dmelbw = (melmax-melmin)/(MEL_FB->num_filters+1);

    if (MEL_FB->doublewide==ON){
	melmin = melmin-dmelbw;
	melmax = melmax+dmelbw;
	if ((fe_melinv(melmin)<0) ||
	    (fe_melinv(melmax)>MEL_FB->sampling_rate/2)){
	    fprintf(stderr,"Out of Range: low  filter edge = %f (%f)\n",fe_melinv(melmin),0.0);
	    fprintf(stderr,"              high filter edge = %f (%f)\n",fe_melinv(melmax),MEL_FB->sampling_rate/2);
	    fprintf(stderr,"exiting...\n");
	    exit(0);
	}
    }    
    
    if (MEL_FB->doublewide==ON){
        for (i=0;i<=MEL_FB->num_filters+3; ++i){
	    filt_edge[i] = fe_melinv(i*dmelbw + melmin);
        }
    }
    else {
	for (i=0;i<=MEL_FB->num_filters+1; ++i){
	    filt_edge[i] = fe_melinv(i*dmelbw + melmin);   
	}
    }
    
    for (whichfilt=0;whichfilt<MEL_FB->num_filters; ++whichfilt) {
      /*line triangle edges up with nearest dft points... */
      if (MEL_FB->doublewide==ON){
	leftfr   = (float32)((int32)((filt_edge[whichfilt]/dfreq)+0.5))*dfreq;
	centerfr = (float32)((int32)((filt_edge[whichfilt+2]/dfreq)+0.5))*dfreq;
	rightfr  = (float32)((int32)((filt_edge[whichfilt+4]/dfreq)+0.5))*dfreq;
      }else{
	leftfr   = (float32)((int32)((filt_edge[whichfilt]/dfreq)+0.5))*dfreq;
	centerfr = (float32)((int32)((filt_edge[whichfilt+1]/dfreq)+0.5))*dfreq;
	rightfr  = (float32)((int32)((filt_edge[whichfilt+2]/dfreq)+0.5))*dfreq;
      }
      MEL_FB->left_apex[whichfilt] = leftfr;
      fwidth = rightfr - leftfr;
      
      /* 2/fwidth for triangles of area 1 */
      height = 2/(float32)fwidth;
      leftslope = height/(centerfr-leftfr);
      rightslope = height/(centerfr-rightfr);
      
      start_pt = 1 + (int32)(leftfr/dfreq);
      freq = (float32)start_pt*dfreq;
      i=0;
      
      while (freq<=centerfr){
	MEL_FB->filter_coeffs[whichfilt][i] = (freq-leftfr)*leftslope;	    
	freq += dfreq;
	i++;
      }
      while (freq<rightfr){
	MEL_FB->filter_coeffs[whichfilt][i] = (freq-rightfr)*rightslope;
	freq += dfreq;
	i++;
      }
      
      MEL_FB->width[whichfilt] = i;
    }
    
    free(filt_edge);
    return(0);
}

int32 fe_compute_melcosine(melfb_t *MEL_FB)
{

    float32 period, freq;
    int32 i,j;
    
    period = (float32)2*MEL_FB->num_filters;

    if ((MEL_FB->mel_cosine = (float32 **) fe_create_2d(MEL_FB->num_cepstra,MEL_FB->num_filters,
					      sizeof(float32)))==NULL){
	fprintf(stderr,"memory alloc failed in fe_compute_melcosine()\n...exiting\n");
	exit(0);
    }
    
    
    for (i=0; i<MEL_FB->num_cepstra; i++) {
	freq = 2*(float32)M_PI*(float32)i/period;
	for (j=0;j< MEL_FB->num_filters;j++)
	    MEL_FB->mel_cosine[i][j] = (float32)cos((float64)(freq*(j+0.5)));	
    }    

    return(0);
	
}


float32 fe_mel(float32 x)
{
    return (float32)(2595.0*log10(1.0+x/700.0));
}

float32 fe_melinv(float32 x)
{
    return (float32)(700.0*(pow(10.0,x/2595.0) - 1.0));
}


void fe_pre_emphasis(int16 const *in, float64 *out, int32 len, float32
		     factor, int16 prior)
{
    int32 i;
  
    out[0] = (float64)in[0]-factor*(float64)prior;
    for (i=1; i<len;i++) {
	out[i] = (float64)in[i] - factor*(float64)in[i-1];
    }
 
}

void fe_short_to_double(int16 const *in, float64 *out, int32 len)
{
    int32 i;
    
    for (i=0;i<len;i++)
	out[i] = (float64)in[i];
}

    
void fe_create_hamming(float64 *in, int32 in_len)
{
    int i;
     
    if (in_len>1){
	for (i=0; i<in_len; i++)
	    in[i] = 0.54 - 0.46*cos(2*M_PI*i/((float64)in_len-1.0));
    }
    return;
    
}


void fe_hamming_window(float64 *in, float64 *window, int32 in_len)
{
    int i;
    
    if (in_len>1){
	for (i=0; i<in_len; i++)
	    in[i] *= window[i];
    }
    return;
    
}


int32 fe_frame_to_fea(fe_t *FE, float64 *in, float64 *fea)
{
    float64 *spec, *mfspec;
    int32 returnValue = FE_SUCCESS;
    
    if (FE->FB_TYPE == MEL_SCALE){
	spec = (float64 *)calloc(FE->FFT_SIZE, sizeof(float64));
	mfspec = (float64 *)calloc(FE->MEL_FB->num_filters, sizeof(float64));

	if (spec==NULL || mfspec==NULL){
	    fprintf(stderr,"memory alloc failed in fe_frame_to_fea()\n...exiting\n");
	    exit(0);
	}
	
 	fe_spec_magnitude(in, FE->FRAME_SIZE, spec, FE->FFT_SIZE);
	fe_mel_spec(FE, spec, mfspec);
	returnValue = fe_mel_cep(FE, mfspec, fea);

	free(spec);
	free(mfspec);	
    }
    else {
	fprintf(stderr,"MEL SCALE IS CURRENTLY THE ONLY IMPLEMENTATION!\n");
	exit(0);
    }
    return returnValue;    
}



void fe_spec_magnitude(float64 const *data, int32 data_len, float64 *spec, int32 fftsize)
{
    int32  j,wrap;
    complex  *FFT, *IN;
    
    /*fftsize defined at top of file*/
    FFT = (complex *) calloc(fftsize,sizeof(complex));
    IN = (complex *) calloc(fftsize,sizeof(complex));
    
    if (FFT==NULL || IN==NULL){
	fprintf(stderr,"memory alloc failed in fe_spec_magnitude()\n...exiting\n");
	exit(0);
    }
	
    if (data_len > fftsize)  /*aliasing */
    {
	
	for (j=0; j<fftsize;j++) {
	    IN[j].r = data[j];
	    IN[j].i = 0.0;
	}
	for (wrap=0; j<data_len; wrap++,j++) {
	    IN[wrap].r += data[j];
	    IN[wrap].i += 0.0;
	}
    }
    else
    {
	for (j=0; j < data_len; j++){
		IN[j].r = data[j];
		IN[j].i = 0.0;
	}
        for ( ;j<fftsize;j++) {  /*pad zeros if necessary */
		IN[j].r = 0.0;
		IN[j].i = 0.0;
	}
    }
    

    fe_fft(IN,FFT,fftsize,FORWARD_FFT);
    
    for (j=0; j <= fftsize/2; j++)
    {	
	spec[j] = FFT[j].r*FFT[j].r + FFT[j].i*FFT[j].i;
    }

    free(FFT);
    free(IN);
    return;
}

void fe_mel_spec(fe_t *FE, float64 const *spec, float64 *mfspec)
{
    int32 whichfilt, start, i;
    float32 dfreq;
    
    dfreq = FE->SAMPLING_RATE/(float32)FE->FFT_SIZE;