enhance.c

remove background noise from speech signal

	
	float   *p;    //define a float pointer which will used to rotate M1[] to M4[]


	/* work out fraction of available CPU time used by algorithm */    
	cpufrac = ((float) (io_ptr & (FRAMEINC - 1)))/FRAMEINC;  
		
	/* wait until io_ptr is at the start of the current frame */ 	
	while((io_ptr/FRAMEINC) != frame_ptr); 
	
	/* then increment the framecount (wrapping if required) */ 
	if (++frame_ptr >= (CIRCBUF/FRAMEINC)) frame_ptr=0;
 	
 	/* save a pointer to the position in the I/O buffers (inbuffer/outbuffer) where the 
 	data should be read (inbuffer) and saved (outbuffer) for the purpose of processing */
 	io_ptr0=frame_ptr * FRAMEINC;
	
	/* copy input data from inbuffer into inframe (starting from the pointer position) */ 
	 
	m=io_ptr0;
    for (k=0;k<FFTLEN;k++)
	{                           
		inframe[k] = inbuffer[m] * inwin[k]; 
		if (++m >= CIRCBUF) m=0; /* wrap if required */
	} 
	
	/************************* DO PROCESSING OF FRAME  HERE **************************/
	
	/******************************** FFT processing *********************************/
	/* copy data frome inframe[] to inframe_copy[] */	
	for (k=0;k<FFTLEN;k++)
	{
		inframe_copy[k]=cmplx(inframe[k],0);
	}

    /* do FFT on inframe_copy[] */
	fft(FFTLEN,inframe_copy);

	/* calculate magnitude of the spectrum and store them in inframe_mag[] */
	for (k=0;k<FFTLEN;k++)
	{
		inframe_mag[k]=cabs(inframe_copy[k]);
	}
	     
	
	/**************************** noise estimation *************************************/
    /* different versions of |X(w)|(3 versions totally) */
	/* all the versions of |X(w)| will be eventually set to Px[], 
	   use "X_opt" to select which version will be selected */
	/* (enhancement 1 ; 2) */
	switch(X_opt)
		{
			case 0: //Px(w)=|X(w)|
				for (k=0;k<FFTLEN;k++)
					Px[k]=inframe_mag[k];    //Px[] is the version of |X(w)|
				break;
			case 1: //Px(w) low pass version of |X(w)|
				for (k=0;k<FFTLEN;k++)
					Px[k]=inframe_mag[k]*(1-K_lowpass)+Px[k]*K_lowpass;
				break;
			case 2: //Px(w) low pass version of |X(w)|^2
				for (k=0;k<FFTLEN;k++)
					Px[k]=sqrt((1-K_lowpass)*(inframe_mag[k]*inframe_mag[k])+K_lowpass*(Px[k]*Px[k]));
				break;
		}
		

	/* get M1[] */
	/* there are 3 different versions of M1[] as well. The reason is that Px[] might be
	copied to M1[] and Px[] has different version */
	if (num==0)
		/* at the beginning of each 2.5sec(i.e. num=0), set the initial values of M1[] 
	    to be inframe_mag[] */
	{
		for (k=0;k<FFTLEN;k++)
			M1[k]=Px[k];	
	}
	else
	{
    	for (k=0;k<FFTLEN;k++)
		{
			/* overwrite M1[k] if inframe_mag[k]<M1[k] */
			M1[k]=min(Px[k],M1[k]);		
		} 
	}
	/* increase the number of frames having been processed within 2.5sec by 1 */
	num++;   
	
	
	/* noise spectrum estimation*/
	/* when M1[] to M4[] have been filled with their corrsponding minimum noise 
	magnitude, compare M1[k] to M4[k] to find the minimum valued during last 10 sec*/
	for (k=0;k<FFTLEN;k++)
	{
		/* find out minimum magnitude of each frequency bin from M1[] to M4[] */
		noise[k]=min(M1[k],M_min[k]);  // M_min[] stores the minimum value from M2[] to M4[]
		/* different versions of noise (use "N_opt" to select different version of noise);
		   use N_opt to select whether to use the low pass version of noise or not */
	    if (N_opt==3)  //low-pass version of noise
			Pn[k]=(1-K_lowpass) * noise[k]+K_lowpass * Pn[k];
		else
			Pn[k]=noise[k];
		
		/* different value of a corresponding to different SNR (enhancement 6);
		   use "a_opt" to select whether to do enhancement 6 or not */
		if (a_opt==1)	 
		{   
			SNR[k]=20.0 * log10(Px[k]-Pn[k])/Pn[k];
			/* calculate a_SNR[] according to SNR[k] */
			if (SNR[k]>20.0)
				a_SNR[k]=1.0;
			else if (SNR[k]<-5.0)
				a_SNR[k]=4.75;
			else if (SNR[k]>=-5.0 && SNR[k]<=20.0)
				a_SNR[k]=(-0.15) * SNR[k]+4.0;
			a= a_SNR[k];
		}
		/* correct underestimation of noise */
		Pn[k]=a*Pn[k];  
	}


	/* if NFRAME frames have been processed, then rotate M1[] to M4[] */
	if (num==NFRAME)
	{
		p  = M1;
	    M1 = M4;
	    M4 = M3;
        M3 = M2;
	    M2 = p;
		/*set num (number of frames having been processed in 2.5 sec) to be zero */
		num= 0; 
		/* find the minimum value of noise stored in M2 to M4 */
		for (k=0;k<FFTLEN;k++)
			M_min[k]=min(min(M2[k],M3[k]),M4[k]);
	}


   /********************************* spectral substraction ***************************/
   /* use "G_opt" to select different expressions of G[k]*/
	for (k=0;k<FFTLEN;k++)
	{
		/* different version of G[k] (enhancement 4;5) */
		switch(G_opt)
		{
		    case 0: 
				G[k]=max(LAMDA,1-Pn[k]/Px[k]);
				break;
			case 1: 
				G[k]=max(LAMDA*Pn[k]/Px[k],(1-Pn[k]/Px[k]));
				break;
			case 2:
				LPx[k]=inframe_mag[k]*(1-K_lowpass)+LPx[k]*K_lowpass;
				G[k]=max(LAMDA*LPx[k]/Px[k],(1-Pn[k]/Px[k]));
				break;
			case 3:
				LPx[k]=inframe_mag[k]*(1-K_lowpass)+LPx[k]*K_lowpass;
				G[k]=max(LAMDA*Pn[k]/LPx[k],(1-Pn[k]/Px[k]));
				break;
			case 4:
				LPx[k]=inframe_mag[k]*(1-K_lowpass)+LPx[k]*K_lowpass;
				G[k]=max(LAMDA,(1-Pn[k]/LPx[k]));
				break;
			case 5:
				G[k]=max(LAMDA,sqrt(1-Pn[k]*Pn[k]/Px[k]*Px[k]));
				break;
		}
		
		/* substract noise spectrum */
		/* enhancement 8 */
		/* use "Y_opt" to select whether to do "enhancement 8" or not */
		switch (Y_opt)
		{
			case 0:  // not do enhancement 8
				inframe_copy[k].r=inframe_copy[k].r*G[k];
				inframe_copy[k].i=inframe_copy[k].i*G[k];
				break;
			case 1:  // enhancement 8
				/* calculate the current Y(w) of the current frame*/
				Y3[k].r=inframe_copy[k].r*G[k];
				Y3[k].i=inframe_copy[k].i*G[k];
				/* if Pn1[k]/Px1[k]>threshold in previous frame, do process below */
				if ((Pn1[k]/Px1[k])>threshold)
				{
					/* find the minimum value from Y1[] to Y3[], 
					and copy the minimum value to Y2[] */
					if ((Y1[k].r*Y1[k].r+Y1[k].i*Y1[k].i)<(Y2[k].r*Y2[k].r+Y2[k].i*Y2[k].i))
					{
						Y2[k].r=Y1[k].r;
						Y2[k].i=Y1[k].i;
					}
					if ((Y3[k].r*Y3[k].r+Y3[k].i*Y3[k].i)<(Y2[k].r*Y2[k].r+Y2[k].i*Y2[k].i))
					{
						Y2[k].r=Y3[k].r;
						Y2[k].i=Y3[k].i;
					}
				}
				/* copy processed data to inframe_copy[] */
				inframe_copy[k].r=Y2[k].r;
				inframe_copy[k].i=Y2[k].i;
				/* rotate Y1[] to Y3[] */
				Y1=Y2;
				Y2=Y3;
				/* set Pn1 to be Pn; Px1 to be Px*/
				Pn1[k]=Pn[k];
				Px1[k]=Px[k];
				break;
		}				
	}

	
	/********************************* IFFT processing *******************************/
	/* IFFT on inframe_copy[] */
	ifft(FFTLEN, inframe_copy);

	/* copy the real part of inframe_copy[] to inframe[]*/
	for (k=0;k<FFTLEN;k++)
	{
		inframe[k]=inframe_copy[k].r;
	}
						      	
										
    for (k=0;k<FFTLEN;k++)
	{                           
		outframe[k] = inframe[k];      /* copy input straight into output */ 
	} 


	
	
	/********************************************************************************/
	
    /* multiply outframe by output window and overlap-add into output buffer */  
                           
	m=io_ptr0;
    
    for (k=0;k<(FFTLEN-FRAMEINC);k++) 
	{    										/* this loop adds into outbuffer */                       
	  	outbuffer[m] = outbuffer[m]+outframe[k]*outwin[k];   
		if (++m >= CIRCBUF) m=0; /* wrap if required */
	}         
    for (;k<FFTLEN;k++) 
	{                           
		outbuffer[m] = outframe[k]*outwin[k];   /* this loop over-writes outbuffer */        
	    m++;
	}	                                   
}        
/*************************** INTERRUPT SERVICE ROUTINE  *****************************/

// Map this to the appropriate interrupt in the CDB file
   
void ISR_AIC(void)
{       
	short sample;
	/* Read and write the ADC and DAC using inbuffer and outbuffer */
	
	sample = mono_read_16Bit();
	inbuffer[io_ptr] = ((float)sample)*ingain;
		/* write new output data */
	mono_write_16Bit((int)(outbuffer[io_ptr]*outgain)); 
	
	/* update io_ptr and check for buffer wraparound */    
	
	if (++io_ptr >= CIRCBUF) io_ptr=0;
}

/************************************************************************************/