min_gamma = .6;
int speex_preprocess(SpeexPreprocessState *st, float *x, float *echo){   int i;   int is_speech=1;   float mean_post=0;   float mean_prior=0;   float energy;   int N = st->ps_size;   int N3 = 2*N - st->frame_size;   int N4 = st->frame_size - N3;   float scale=.5/N;   float *ps=st->ps;   /* 'Build' input frame */   for (i=0;i<N3;i++)      st->frame[i]=st->inbuf[i];   for (i=0;i<st->frame_size;i++)      st->frame[N3+i]=x[i];      /* Update inbuf */   for (i=0;i<N3;i++)      st->inbuf[i]=x[N4+i];   /* Windowing */   for (i=0;i<2*N;i++)      st->frame[i] *= st->window[i];   /* Perform FFT */   drft_forward(st->fft_lookup, st->frame);   /**************************************************************     *  Denoise in spectral domain using Ephraim-Malah algorithm  *    **************************************************************/   /* Power spectrum */   ps[0]=1;   for (i=1;i<N;i++)      ps[i]=1+st->frame[2*i-1]*st->frame[2*i-1] + st->frame[2*i]*st->frame[2*i];   energy=0;   for (i=1;i<N;i++)      energy += log(100+ps[i]);   energy /= 160;   st->last_energy[st->nb_preprocess%STABILITY_TIME]=energy;   if (st->nb_preprocess>=STABILITY_TIME)   {      float E=0, E2=0;      float std;      for (i=0;i<STABILITY_TIME;i++)      {         E+=st->last_energy[i];         E2+=st->last_energy[i]*st->last_energy[i];      }      E2=E2/STABILITY_TIME;      E=E/STABILITY_TIME;      std = sqrt(E2-E*E);      if (std<.15 && st->last_update>20)      {         update_noise(st, &st->last_ps[st->last_id*N], echo);      }      /*fprintf (stderr, "%f\n", std);*/   }   st->nb_preprocess++;   /* Noise estimation always updated for the 20 first times */   if (st->nb_adapt<20)      /*if (st->nb_adapt<25 && st->nb_adapt>15)*/   {      update_noise(st, ps, echo);      st->last_update=0;   }   /* Deal with residual echo if provided */   if (echo)      for (i=1;i<N;i++)         st->echo_noise[i] = (.7*st->echo_noise[i] + .3* echo[i]);   /* Compute a posteriori SNR */   for (i=1;i<N;i++)   {      st->post[i] = ps[i]/(1+st->noise[i]+st->echo_noise[i]) - 1;      if (st->post[i]>100)         st->post[i]=100;      /*if (st->post[i]<0)        st->post[i]=0;*/      mean_post+=st->post[i];   }   mean_post /= N;   if (mean_post<0)      mean_post=0;   /* Special case for first frame */   if (st->nb_adapt==1)      for (i=1;i<N;i++)         st->old_ps[i] = ps[i];   /* Compute a priori SNR */   {      /* A priori update rate */      float gamma;      float min_gamma=0.12;      gamma = 1.0/st->nb_preprocess;      /*Make update rate smaller when there's no speech*/#if 0      if (mean_post<3.5 && mean_prior < 1)         min_gamma *= (mean_post+.5);      else         min_gamma *= 4.;#else      min_gamma = .2*fabs(mean_prior - mean_post)*fabs(mean_prior - mean_post);      if (min_gamma>.6)         min_gamma = .6;      if (min_gamma<.01)         min_gamma = .01;#endif      min_gamma = .6;      if (gamma<min_gamma)         gamma=min_gamma;            for (i=1;i<N;i++)      {                  /* A priori SNR update */         st->prior[i] = gamma*max(0.0,st->post[i]) +         (1-gamma)*st->gain[i]*st->gain[i]*st->old_ps[i]/(1+st->noise[i]+st->echo_noise[i]);                  if (st->prior[i]>100)            st->prior[i]=100;                  mean_prior+=st->prior[i];      }   }   mean_prior /= N;#if 0   for (i=0;i<N;i++)   {      fprintf (stderr, "%f ", st->prior[i]);   }   fprintf (stderr, "\n");#endif   /*fprintf (stderr, "%f %f\n", mean_prior,mean_post);*/   if (st->nb_preprocess>=20)   {      int do_update = 0;      float noise_ener=0, sig_ener=0;      /* If SNR is low (both a priori and a posteriori), update the noise estimate*/      /*if (mean_prior<.23 && mean_post < .5)*/      if (mean_prior<.23 && mean_post < .5)         do_update = 1;      for (i=1;i<N;i++)      {         noise_ener += st->noise[i];         sig_ener += ps[i];      }      if (noise_ener > 3*sig_ener)         do_update = 1;      /*do_update = 0;*/      if (do_update)      {         st->consec_noise++;      } else {         st->consec_noise=0;      }   }   if (st->vad_enabled)      is_speech = speex_compute_vad(st, ps, mean_prior, mean_post);   if (st->consec_noise>=3)   {      update_noise(st, st->old_ps, echo);      st->last_update=0;   } else {      st->last_update++;   }   /* Compute gain according to the Ephraim-Malah algorithm */   for (i=1;i<N;i++)   {      float MM;      float theta;      float prior_ratio;      prior_ratio = st->prior[i]/(1.0001+st->prior[i]);      theta = (1+st->post[i])*prior_ratio;#if 0      /* Spectral magnitude estimator */      /* Approximation of:         exp(-theta/2)*((1+theta)*I0(theta/2) + theta.*I1(theta/2))         because I don't feel like computing Bessel functions      */      /*MM = -.22+1.155*sqrt(theta+1.1);*/      MM=-.22+1.163*sqrt(theta+1.1)-.0015*theta;      st->gain[i] = SQRT_M_PI_2*sqrt(prior_ratio/(1.0001+st->post[i]))*MM;      if (st->gain[i]>1)      {         st->gain[i]=1;      }#else      /* log-spectral magnitude estimator */      if (theta<6)         MM = 0.74082*pow(theta+1,.61)/sqrt(.0001+theta);      else         MM=1;      st->gain[i] = prior_ratio * MM;      /*Put some (very arbitraty) limit on the gain*/      if (st->gain[i]>2)      {         st->gain[i]=2;      }#endif      /*st->gain[i] = prior_ratio;*/   }   st->gain[0]=0;   st->gain[N-1]=0;   if (st->denoise_enabled)   {      for (i=1;i<N-1;i++)      {         st->gain2[i]=st->gain[i];         /* Limits noise reduction to -26 dB, put prevents some musical noise */         if (st->gain2[i]<.05)            st->gain2[i]=.05;      }      st->gain2[N-1]=0;   } else {      for (i=0;i<N;i++)         st->gain2[i] = 1;   }   if (st->agc_enabled)      speex_compute_agc(st, mean_prior);#if 0   if (!is_speech)   {      for (i=0;i<N;i++)         st->gain2[i] = 0;   }#if 0 else {      for (i=0;i<N;i++)         st->gain2[i] = 1;   }#endif#endif   /* Apply computed gain */   for (i=1;i<N;i++)   {      st->frame[2*i-1] *= st->gain2[i];      st->frame[2*i] *= st->gain2[i];   }   /* Get rid of the DC and very low frequencies */   st->frame[0]=0;   st->frame[1]=0;   st->frame[2]=0;   /* Nyquist frequency is mostly useless too */   st->frame[2*N-1]=0;   /* Inverse FFT with 1/N scaling */   drft_backward(st->fft_lookup, st->frame);   for (i=0;i<2*N;i++)      st->frame[i] *= scale;   {      float max_sample=0;      for (i=0;i<2*N;i++)         if (fabs(st->frame[i])>max_sample)            max_sample = fabs(st->frame[i]);      if (max_sample>28000)      {         float damp = 28000./max_sample;         for (i=0;i<2*N;i++)            st->frame[i] *= damp;      }   }   for (i=0;i<2*N;i++)      st->frame[i] *= st->window[i];   /* Perform overlap and add */   for (i=0;i<N3;i++)      x[i] = st->outbuf[i] + st->frame[i];   for (i=0;i<N4;i++)      x[N3+i] = st->frame[N3+i];      /* Update outbuf */   for (i=0;i<N3;i++)      st->outbuf[i] = st->frame[st->frame_size+i];   /* Save old power spectrum */   for (i=1;i<N;i++)      st->old_ps[i] = ps[i];   for (i=1;i<N;i++)      st->last_ps[st->last_id*N+i] = ps[i];   st->last_id++;   if (st->last_id>=NB_LAST_PS)      st->last_id=0;   return is_speech;}