float loudness=0;
int denoise(DenoiseState *st, float *x){   int i;   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_denoise%STABILITY_TIME]=energy;   if (st->nb_denoise>=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]);      }      /*fprintf (stderr, "%f\n", std);*/   }   st->nb_denoise++;#if 0   if (st->nb_min_estimate<50)   {      float ener=0;      for (i=1;i<N;i++)         ener += ps[i];      /*fprintf (stderr, "%f\n", ener);*/      if (ener < st->min_ener || st->nb_min_estimate==0)      {         st->min_ener = ener;         for (i=1;i<N;i++)            st->min_ps[i] = ps[i];      }      st->nb_min_estimate++;   } else {      float noise_ener=0;      st->nb_min_estimate=0;      for (i=1;i<N;i++)         noise_ener += st->noise[i];      /*fprintf (stderr, "%f %f\n", noise_ener, st->min_ener);*/      if (0&&(st->last_update>50 && st->min_ener > 3*noise_ener) || st->last_update>50)      {         for (i=1;i<N;i++)         {            if (st->noise[i] < st->min_ps[i])               st->noise[i] = st->min_ps[i];         }         /*fprintf (stderr, "tata %d\n",st->last_update);*/         st->last_update=0;      } else {         /*fprintf (stderr, "+");*/      }   }#endif   /* Noise estimation always updated for the 20 first times */   if (st->nb_adapt<15)   {      update_noise(st, ps);      st->last_update=0;   }   /* Compute a posteriori SNR */   for (i=1;i<N;i++)   {      st->post[i] = ps[i]/(1+st->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.05;      gamma = 1.0/st->nb_denoise;      /*Make update rate smaller when there's no speech*/      if (mean_post<3)         min_gamma *= (mean_post+.1);      else         min_gamma *= 3.1;      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]/st->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 SNR is low (both a priori and a posteriori), update the noise estimate*/   if (mean_prior<.23 && mean_post < .5 && st->nb_adapt>=20)   {      st->consec_noise++;      /*fprintf (stderr, "noise\n");*/   } else {      st->consec_noise=0;   }   if (st->consec_noise>=3)   {      update_noise(st, st->old_ps);      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;      /* 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;      }      /*st->gain[i] = prior_ratio;*/   }   st->gain[0]=0;   st->gain[N-1]=0;   for (i=1;i<N-1;i++)   {      st->gain2[i]=st->gain[i];      if (st->gain2[i]<.1)         st->gain2[i]=.1;   }   st->gain2[N-1]=0;   if ((mean_prior>3&&mean_prior>3))   {      st->nb_loudness_adapt++;      float rate=2.0/(1+st->nb_loudness_adapt);      if (rate < .01)         rate = .01;      float loudness=0;      for (i=2;i<N;i++)      {         loudness += scale*st->ps[i] * st->gain2[i] * st->gain2[i] * st->loudness_weight[i];      }      loudness=sqrt(loudness);      /*if (loudness < 2*pow(st->loudness, 1.0/LOUDNESS_EXP) &&        loudness*2 > pow(st->loudness, 1.0/LOUDNESS_EXP))*/      st->loudness = (1-rate)*st->loudness + (rate)*pow(loudness, LOUDNESS_EXP);            st->loudness2 = (1-rate)*st->loudness2 + rate*pow(st->loudness, 1.0/LOUDNESS_EXP);      loudness = pow(st->loudness, 1.0/LOUDNESS_EXP);      /*fprintf (stderr, "%f %f %f\n", loudness, st->loudness2, rate);*/   }   for (i=0;i<N;i++)      st->gain2[i] *= 6000.0/st->loudness2;   /* 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*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 1;}