noise_est.c

ears-0.32, linux下有用的语音信号处理工具包

/****************************************************************************
*                                                                           *
*    ROUTINES IN THIS FILE:                                                 *
*                                                                           *
*        comp_noisepower() : computes the noise power of a speech signal    *
*        creat_filt() : find the coefficients of the filter for high        *
*                       energy speech detection                             *
*        det() :      for detecting high energy speech                      * 
*                     
****************************************************************************/

#include <stdio.h>
#include "others/mymath.h"
#include "rasta.h"
#include "functions.h"

/*---------------------------------------------------------------------------
  comp_noisepower
  
  This routine computes the noise power of a speech signal. In the first 100ms,
  it is assumed that the first imcoming speech samples of a recording are 
  related to the noise only. Thus the average signal energy is used as the first
  estimation of the noise power. After 200ms, the noise level in a certain   
  subband is estimated by a statistical analysis of a segment of the
  magnitude spectral envelope("Estimation of Noise Spectrum and its Application
  to SNR-Estimation and Speech Enhancement" by H. Guenter Hirsch) 
  Given a spectral envelope and the corresponding distribution density function
  in a certain subband, the most frequently occurring spectral magnitude value
  is taken as an estimation for the noise level inside this band. These noise
  levels for different subbands are squared and then the average of these   
  squared values gives the noise power estimate. The noise power is computed
  using this histogram method every 100ms. In Hirsch's paper, FFT subbands 
  were used for the analysis. In this routine, we have chosen to use critical
  subbands for the analysis. Some modifications are required. Sometimes the
  most frequently occurring value has an unrealistically  high value. A simple
  detector is used to exclude some of the unrealistically high values 
  due to speech from the histograms. 
  --------------------------------------------------------------------------*/       


/* local function prototypes */
void create_filt(float,  int , float , struct fvec *);
int det(struct fmat *, int ,  struct fvec *, struct fvec *);   


struct fvec *comp_noisepower( struct fhistory *hptr, const struct param *pptr, struct fvec *aspectrum)
{
	char *funcname;
	static int i_call = 0;  /*  Frame number */
	static struct fvec *noiseptr = NULL;  /* the noise power (i.e. the return value) */
	static struct fvec *max = NULL;  /* Most frequently occurring value in a subband */  
	static struct svec *new_max = NULL; /*  A flag. Flag = 1 Current spectral magnitude  
											in a subband is a new maximum 
											Flag = 0 otherwise                       */
	static struct svec *new_max_index = NULL;  /* index of the current maximum */  
	static struct fvec *noise_level_lt = NULL; /* estimated noise level in a subband */
	static struct fvec *waste_trajec = NULL;   /* oldest value in the histogram are
												  thrown away  */ 
	static struct fmat *all_trajec = NULL;     /* matrix containing all trajectories */
	static struct fmat *histo_values = NULL;   /* matrix containing the entries for
												  the histograms for each band       */
	static struct svec *histo_index = NULL;  /* index for the entries of the histograms */ 
	static int   stat_lt[MAXFILTS][HISTO_RESOLUTION];  /* the histograms, the first index is 
														  the subband number, the second index
														  is the histogram bin number */ 
	static struct fvec *filter = NULL;   /* filter coefficients for high energy speech detector  */
	static struct fmat *buffer = NULL;   /* buffer for the high energy speech detector */  
	
	int i,j,k;
	int t;
	static int m1;
	int lastfilt;
	float average, sum = 0.0;
	static int forget_frame = 1; /* forget first frame (building-up of highpass filter) */
	static float powerOLD;     /* noise power of history file */
	static float noisepower;   /* the estimated noise power */
	static int histo_length;   /* number of entries in each histograms */
	static int histo_step_framesize; /* The histograms are re-evaluated
										every histo_step_framesize */  
	int mmax_lp;     
	int mindex = 0;
	int length;    
	float mmax;     
	int mmax_index;  
	float average1;  
	
	/*---------------------------------------------------------------*/ 
	static FILE *outdebug;
	/*---------------------------------------------------------------*/
	
	funcname = "comp_noisepower";
	if (noiseptr == (struct fvec *)NULL)
	{
        noiseptr = alloc_fvec(1);
	}
	
	if(forget_frame == 1)
	{
		forget_frame = 0;
		if(pptr->HPfilter == 1) /* forget first frame */
		{
			noiseptr->values[0] = 1/(C_CONSTANT * 1.0e-7);
			return(noiseptr);
		}
	} 
	
	i_call++;
	lastfilt = pptr->nfilts - pptr->first_good; 
	
	/*----------------------------------------------------------------*/ 
	if ((i_call == 1) && (pptr->debug==TRUE))
	{
        if ((outdebug = fopen("rasta.debug","wt")) == NULL)
        {
			fprintf(stderr,"Cannot open output file rasta.debug\n");
			exit(-1);
        }
	}
	/*------------------------------------------------------------------*/ 
	
	/* Initialization before the first frame */
	if (i_call  == 1)
	{
        /* There are histo_length entries for each histogram */
        histo_length = HISTO_TIME/(pptr->stepsize);  
        /* The histograms are evaluated every histo_step_framesize frames */                               
        histo_step_framesize = HISTO_STEPTIME/(pptr->stepsize);  
        max = alloc_fvec(pptr->nfilts);
        new_max = alloc_svec(pptr->nfilts);
        new_max_index = alloc_svec(pptr->nfilts);
        noise_level_lt = alloc_fvec(pptr->nfilts);
        waste_trajec = alloc_fvec(pptr->nfilts);
        all_trajec = alloc_fmat(pptr->nfilts, histo_length);
        histo_values =  alloc_fmat(pptr->nfilts, histo_length);
        histo_index = alloc_svec(pptr->nfilts);
        filter = alloc_fvec(FILT_LENGTH);
        buffer = alloc_fmat(pptr->nfilts, FILT_LENGTH);
        /* Find coefficients of the filter for the high energy speech detector */
        create_filt(OVER_EST_FACTOR, FILT_LENGTH, TIME_CONSTANT, filter); 
        for (i = 0; i < pptr->nfilts; i++)
        {
			max->values[i] = 0;
			new_max->values[i] = RESET;
			new_max_index->values[i] = 0;
			noise_level_lt->values[i] = 0;
			histo_index->values[i] = -1;
			for (j=0; j<histo_length; j++)
			{
				all_trajec->values[i][j] = 0;
				histo_values->values[i][j] = 0;
			}
	    }
		
		if(hptr->normInit == FALSE)
		{
			for(i=pptr->first_good, powerOLD = 0.0; i<lastfilt; i++)
				powerOLD += (hptr->noiseOLD[i] * hptr->noiseOLD[i]);
			powerOLD /= (float)(lastfilt - (pptr->first_good));
		}
	}   /*   end of if (i_call == 1)  */
	
	t = ((i_call-1)% histo_length);  
	
	/* copy the spectral magnitudes in the matrix all_trajec */
	if((i_call <= histo_step_framesize) && (hptr->normInit == FALSE))
	{
		for (i=pptr->first_good, sum=0.0; i<lastfilt; i++)
			sum += aspectrum->values[i];
		sum /= (float)(lastfilt - (pptr->first_good));
		if(sum <= powerOLD)
		{
			for (i=pptr->first_good; i<lastfilt; i++)
				all_trajec->values[i][t] = (float)sqrt((double)aspectrum->values[i]);
		}
		else
		{
			for (i=pptr->first_good; i<lastfilt; i++)
				all_trajec->values[i][t] =	hptr->noiseOLD[i];
		}
	}
	else
	{
		for (i=pptr->first_good; i<lastfilt; i++)
			all_trajec->values[i][t] = (float)sqrt((double)aspectrum->values[i]);
	}
	
	/* A simple detector is used in each subband. When the value returned by the subroutine
	   det is "1" for a certain subband, the spectral magnitude in that subband is not entered
	   into the histogram. Otherwise the spectral magnitude is entered into the matrix 
	   histo_values and histo_index is incremented
	   For high energy speech detection, there is a FIR filter for each subband. The filter for each
	   subband is identical 
	   The buffers of the filters for high energy speech detection are initialized  with spectral magnitude
	   values. The first spectral magnitude value of each subband is assigned to the "histo_step_framesize"
	   th element (f th) of the corresponding buffer. The 2rd to f th values are assigned to the (f-1)th
	   to 1st elements of the buffer. The rest of the elements are initialized with the average spectral  
	   magnitude values */
	
	if (i_call <= histo_step_framesize)
	{
        for(i = pptr->first_good; i<lastfilt; i++)
        {
            buffer->values[i][histo_step_framesize - i_call] = all_trajec->values[i][t];
            histo_values->values[i][t] = all_trajec->values[i][t];
            histo_index->values[i] ++;
        }
	} /* end of if (i_call <= histo_step_framesize) */
	
	if (i_call == histo_step_framesize)
	{
        for (i= pptr->first_good; i<lastfilt; i++)
        {
			mmax = 0;
			mmax_index = 0;
			for (j= 0; j< histo_step_framesize; j++)
			{
				if (buffer->values[i][j] > mmax)
				{
					mmax = buffer->values[i][j];
					mmax_index = j;
				}
			}
			average1 = 0.;
			for (j = 0; j < histo_step_framesize; j++)
			{