Some AGC improvements, noise/echo suppression slightly less aggressive by
[speexdsp.git] / libspeex / preprocess.c
1 /* Copyright (C) 2003 Epic Games (written by Jean-Marc Valin)
2    Copyright (C) 2004-2006 Epic Games 
3    
4    File: preprocess.c
5    Preprocessor with denoising based on the algorithm by Ephraim and Malah
6
7    Redistribution and use in source and binary forms, with or without
8    modification, are permitted provided that the following conditions are
9    met:
10
11    1. Redistributions of source code must retain the above copyright notice,
12    this list of conditions and the following disclaimer.
13
14    2. Redistributions in binary form must reproduce the above copyright
15    notice, this list of conditions and the following disclaimer in the
16    documentation and/or other materials provided with the distribution.
17
18    3. The name of the author may not be used to endorse or promote products
19    derived from this software without specific prior written permission.
20
21    THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
22    IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
23    OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24    DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
25    INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
26    (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27    SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28    HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
29    STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31    POSSIBILITY OF SUCH DAMAGE.
32 */
33
34
35 /*
36    Recommended papers:
37    
38    Y. Ephraim and D. Malah, "Speech enhancement using minimum mean-square error
39    short-time spectral amplitude estimator". IEEE Transactions on Acoustics, 
40    Speech and Signal Processing, vol. ASSP-32, no. 6, pp. 1109-1121, 1984.
41    
42    Y. Ephraim and D. Malah, "Speech enhancement using minimum mean-square error
43    log-spectral amplitude estimator". IEEE Transactions on Acoustics, Speech and 
44    Signal Processing, vol. ASSP-33, no. 2, pp. 443-445, 1985.
45    
46    I. Cohen and B. Berdugo, "Speech enhancement for non-stationary noise environments".
47    Signal Processing, vol. 81, no. 2, pp. 2403-2418, 2001.
48
49    Stefan Gustafsson, Rainer Martin, Peter Jax, and Peter Vary. "A psychoacoustic 
50    approach to combined acoustic echo cancellation and noise reduction". IEEE 
51    Transactions on Speech and Audio Processing, 2002.
52    
53    J.-M. Valin, J. Rouat, and F. Michaud, "Microphone array post-filter for separation
54    of simultaneous non-stationary sources". In Proceedings IEEE International 
55    Conference on Acoustics, Speech, and Signal Processing, 2004.
56 */
57
58 #ifdef HAVE_CONFIG_H
59 #include "config.h"
60 #endif
61
62 #include <math.h>
63 #include "speex/speex_preprocess.h"
64 #include "speex/speex_echo.h"
65 #include "misc.h"
66 #include "fftwrap.h"
67 #include "filterbank.h"
68 #include "math_approx.h"
69
70 #ifndef M_PI
71 #define M_PI 3.14159263
72 #endif
73
74 #define LOUDNESS_EXP 2.5
75
76 #define NB_BANDS 24
77
78 #define SPEECH_PROB_START_DEFAULT       QCONST16(0.35f,15)
79 #define SPEECH_PROB_CONTINUE_DEFAULT    QCONST16(0.20f,15)
80 #define NOISE_SUPPRESS_DEFAULT       -15
81 #define ECHO_SUPPRESS_DEFAULT        -40
82 #define ECHO_SUPPRESS_ACTIVE_DEFAULT -15
83
84 #ifndef NULL
85 #define NULL 0
86 #endif
87
88 #define SQR(x) ((x)*(x))
89 #define SQR16(x) (MULT16_16((x),(x)))
90 #define SQR16_Q15(x) (MULT16_16_Q15((x),(x)))
91
92 #ifdef FIXED_POINT
93 static inline spx_word16_t DIV32_16_Q8(spx_word32_t a, spx_word32_t b)
94 {
95    if (SHR32(a,7) >= b)
96    {
97       return 32767;
98    } else {
99       if (b>=QCONST32(1,23))
100       {
101          a = SHR32(a,8);
102          b = SHR32(b,8);
103       }
104       if (b>=QCONST32(1,19))
105       {
106          a = SHR32(a,4);
107          b = SHR32(b,4);
108       }
109       if (b>=QCONST32(1,15))
110       {
111          a = SHR32(a,4);
112          b = SHR32(b,4);
113       }
114       a = SHL32(a,8);
115       return PDIV32_16(a,b);
116    }
117    
118 }
119 static inline spx_word16_t DIV32_16_Q15(spx_word32_t a, spx_word32_t b)
120 {
121    if (SHR32(a,15) >= b)
122    {
123       return 32767;
124    } else {
125       if (b>=QCONST32(1,23))
126       {
127          a = SHR32(a,8);
128          b = SHR32(b,8);
129       }
130       if (b>=QCONST32(1,19))
131       {
132          a = SHR32(a,4);
133          b = SHR32(b,4);
134       }
135       if (b>=QCONST32(1,15))
136       {
137          a = SHR32(a,4);
138          b = SHR32(b,4);
139       }
140       a = SHL32(a,15)-a;
141       return DIV32_16(a,b);
142    }
143 }
144 #define SNR_SCALING 256.f
145 #define SNR_SCALING_1 0.0039062f
146 #define SNR_SHIFT 8
147
148 #define FRAC_SCALING 32767.f
149 #define FRAC_SCALING_1 3.0518e-05
150 #define FRAC_SHIFT 1
151
152 #define EXPIN_SCALING 2048.f
153 #define EXPIN_SCALING_1 0.00048828f
154 #define EXPIN_SHIFT 11
155 #define EXPOUT_SCALING_1 1.5259e-05
156
157 #define NOISE_SHIFT 7
158
159 #else
160
161 #define DIV32_16_Q8(a,b) ((a)/(b))
162 #define DIV32_16_Q15(a,b) ((a)/(b))
163 #define SNR_SCALING 1.f
164 #define SNR_SCALING_1 1.f
165 #define SNR_SHIFT 0
166 #define FRAC_SCALING 1.f
167 #define FRAC_SCALING_1 1.f
168 #define FRAC_SHIFT 0
169 #define NOISE_SHIFT 0
170
171 #define EXPIN_SCALING 1.f
172 #define EXPIN_SCALING_1 1.f
173 #define EXPOUT_SCALING_1 1.f
174
175 #endif
176
177 /** Speex pre-processor state. */
178 struct SpeexPreprocessState_ {
179    /* Basic info */
180    int    frame_size;        /**< Number of samples processed each time */
181    int    ps_size;           /**< Number of points in the power spectrum */
182    int    sampling_rate;     /**< Sampling rate of the input/output */
183    int    nbands;
184    FilterBank *bank;
185    
186    /* Parameters */
187    int    denoise_enabled;
188    int    vad_enabled;
189    int    dereverb_enabled;
190    spx_word16_t  reverb_decay;
191    spx_word16_t  reverb_level;
192    spx_word16_t speech_prob_start;
193    spx_word16_t speech_prob_continue;
194    int    noise_suppress;
195    int    echo_suppress;
196    int    echo_suppress_active;
197    SpeexEchoState *echo_state;
198    
199    /* DSP-related arrays */
200    spx_word16_t *frame;      /**< Processing frame (2*ps_size) */
201    spx_word16_t *ft;         /**< Processing frame in freq domain (2*ps_size) */
202    spx_word32_t *ps;         /**< Current power spectrum */
203    spx_word16_t *gain2;      /**< Adjusted gains */
204    spx_word16_t *gain_floor; /**< Minimum gain allowed */
205    spx_word16_t *window;     /**< Analysis/Synthesis window */
206    spx_word32_t *noise;      /**< Noise estimate */
207    spx_word32_t *reverb_estimate; /**< Estimate of reverb energy */
208    spx_word32_t *old_ps;     /**< Power spectrum for last frame */
209    spx_word16_t *gain;       /**< Ephraim Malah gain */
210    spx_word16_t *prior;      /**< A-priori SNR */
211    spx_word16_t *post;       /**< A-posteriori SNR */
212
213    spx_word32_t *S;          /**< Smoothed power spectrum */
214    spx_word32_t *Smin;       /**< See Cohen paper */
215    spx_word32_t *Stmp;       /**< See Cohen paper */
216    int *update_prob;       /**< Propability of speech presence for noise update */
217
218    spx_word16_t *zeta;       /**< Smoothed a priori SNR */
219    spx_word32_t *echo_noise;
220    spx_word32_t *residual_echo;
221
222    /* Misc */
223    spx_word16_t *inbuf;      /**< Input buffer (overlapped analysis) */
224    spx_word16_t *outbuf;     /**< Output buffer (for overlap and add) */
225
226    /* AGC stuff, only for floating point for now */
227 #ifndef FIXED_POINT
228    int    agc_enabled;
229    float  agc_level;
230    float *loudness_weight;   /**< Perceptual loudness curve */
231    float  loudness;          /**< Loudness estimate */
232    float  agc_gain;          /**< Current AGC gain */
233    int    nb_loudness_adapt; /**< Number of frames used for loudness adaptation so far */
234    float  max_gain;          /**< Maximum gain allowed */
235    float  max_increase_step; /**< Maximum increase in gain from one frame to another */
236    float  max_decrease_step; /**< Maximum decrease in gain from one frame to another */
237    float  prev_loudness;     /**< Loudness of previous frame */
238    float  init_max;          /**< Current gain limit during initialisation */
239 #endif
240    int    nb_adapt;          /**< Number of frames used for adaptation so far */
241    int    was_speech;
242    int    min_count;         /**< Number of frames processed so far */
243    void  *fft_lookup;        /**< Lookup table for the FFT */
244 #ifdef FIXED_POINT
245    int    frame_shift;
246 #endif
247 };
248
249
250 static void conj_window(spx_word16_t *w, int len)
251 {
252    int i;
253    for (i=0;i<len;i++)
254    {
255       spx_word16_t tmp;
256       spx_word16_t x = DIV32_16(MULT16_16(QCONST16(4.f,13),i),len);
257       int inv=0;
258       if (x<QCONST16(1.f,13))
259       {
260       } else if (x<QCONST16(2.f,13))
261       {
262          x=QCONST16(2.f,13)-x;
263          inv=1;
264       } else if (x<QCONST16(3.f,13))
265       {
266          x=x-QCONST16(2.f,13);
267          inv=1;
268       } else {
269          x=QCONST16(2.f,13)-x+QCONST16(2.f,13); /* 4 - x */
270       }
271       x = MULT16_16_Q14(QCONST16(1.271903f,14), x);
272       tmp = SQR16_Q15(QCONST16(.5f,15)-MULT16_16_P15(QCONST16(.5f,15),spx_cos_norm(QCONST32(x,2))));
273       if (inv)
274          tmp=SUB16(Q15_ONE,tmp);
275       w[i]=spx_sqrt(SHL32(EXTEND32(tmp),15));
276    }
277 }
278
279       
280 #ifdef FIXED_POINT
281 /* This function approximates the gain function 
282    y = gamma(1.25)^2 * M(-.25;1;-x) / sqrt(x)  
283    which multiplied by xi/(1+xi) is the optimal gain
284    in the loudness domain ( sqrt[amplitude] )
285    Input in Q11 format, output in Q15
286 */
287 static inline spx_word32_t hypergeom_gain(spx_word32_t xx)
288 {
289    int ind;
290    spx_word16_t frac;
291    /* Q13 table */
292    static const spx_word16_t table[21] = {
293        6730,  8357,  9868, 11267, 12563, 13770, 14898,
294       15959, 16961, 17911, 18816, 19682, 20512, 21311,
295       22082, 22827, 23549, 24250, 24931, 25594, 26241};
296       ind = SHR32(xx,10);
297       if (ind<0)
298          return Q15_ONE;
299       if (ind>19)
300          return ADD32(EXTEND32(Q15_ONE),EXTEND32(DIV32_16(QCONST32(.1296,23), SHR32(xx,EXPIN_SHIFT-SNR_SHIFT))));
301       frac = SHL32(xx-SHL32(ind,10),5);
302       return SHL32(DIV32_16(PSHR32(MULT16_16(Q15_ONE-frac,table[ind]) + MULT16_16(frac,table[ind+1]),7),(spx_sqrt(SHL32(xx,15)+6711))),7);
303 }
304
305 static inline spx_word16_t qcurve(spx_word16_t x)
306 {
307    x = MAX16(x, 1);
308    return DIV32_16(SHL32(EXTEND32(32767),9),ADD16(512,MULT16_16_Q15(QCONST16(.60f,15),DIV32_16(32767,x))));
309 }
310
311 /* Compute the gain floor based on different floors for the background noise and residual echo */
312 static void compute_gain_floor(int noise_suppress, int effective_echo_suppress, spx_word32_t *noise, spx_word32_t *echo, spx_word16_t *gain_floor, int len)
313 {
314    int i;
315    
316    if (noise_suppress > effective_echo_suppress)
317    {
318       spx_word16_t noise_gain, gain_ratio;
319       noise_gain = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(0.11513,11),noise_suppress)),1)));
320       gain_ratio = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(.2302585f,11),effective_echo_suppress-noise_suppress)),1)));
321
322       /* gain_floor = sqrt [ (noise*noise_floor + echo*echo_floor) / (noise+echo) ] */
323       for (i=0;i<len;i++)
324          gain_floor[i] = MULT16_16_Q15(noise_gain,
325                                        spx_sqrt(SHL32(EXTEND32(DIV32_16_Q15(PSHR32(noise[i],NOISE_SHIFT) + MULT16_32_Q15(gain_ratio,echo[i]),
326                                              (1+PSHR32(noise[i],NOISE_SHIFT) + echo[i]) )),15)));
327    } else {
328       spx_word16_t echo_gain, gain_ratio;
329       echo_gain = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(0.11513,11),effective_echo_suppress)),1)));
330       gain_ratio = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(.2302585f,11),noise_suppress-effective_echo_suppress)),1)));
331
332       /* gain_floor = sqrt [ (noise*noise_floor + echo*echo_floor) / (noise+echo) ] */
333       for (i=0;i<len;i++)
334          gain_floor[i] = MULT16_16_Q15(echo_gain,
335                                        spx_sqrt(SHL32(EXTEND32(DIV32_16_Q15(MULT16_32_Q15(gain_ratio,PSHR32(noise[i],NOISE_SHIFT)) + echo[i],
336                                              (1+PSHR32(noise[i],NOISE_SHIFT) + echo[i]) )),15)));
337    }
338 }
339
340 #else
341 /* This function approximates the gain function 
342    y = gamma(1.25)^2 * M(-.25;1;-x) / sqrt(x)  
343    which multiplied by xi/(1+xi) is the optimal gain
344    in the loudness domain ( sqrt[amplitude] )
345 */
346 static inline spx_word32_t hypergeom_gain(spx_word32_t xx)
347 {
348    int ind;
349    float integer, frac;
350    float x;
351    static const float table[21] = {
352       0.82157f, 1.02017f, 1.20461f, 1.37534f, 1.53363f, 1.68092f, 1.81865f,
353       1.94811f, 2.07038f, 2.18638f, 2.29688f, 2.40255f, 2.50391f, 2.60144f,
354       2.69551f, 2.78647f, 2.87458f, 2.96015f, 3.04333f, 3.12431f, 3.20326f};
355       x = EXPIN_SCALING_1*xx;
356       integer = floor(2*x);
357       ind = (int)integer;
358       if (ind<0)
359          return FRAC_SCALING;
360       if (ind>19)
361          return FRAC_SCALING*(1+.1296/x);
362       frac = 2*x-integer;
363       return FRAC_SCALING*((1-frac)*table[ind] + frac*table[ind+1])/sqrt(x+.0001f);
364 }
365
366 static inline spx_word16_t qcurve(spx_word16_t x)
367 {
368    return 1.f/(1.f+.15f/(SNR_SCALING_1*x));
369 }
370
371 static void compute_gain_floor(int noise_suppress, int effective_echo_suppress, spx_word32_t *noise, spx_word32_t *echo, spx_word16_t *gain_floor, int len)
372 {
373    int i;
374    float echo_floor;
375    float noise_floor;
376
377    noise_floor = exp(.2302585f*noise_suppress);
378    echo_floor = exp(.2302585f*effective_echo_suppress);
379
380    /* Compute the gain floor based on different floors for the background noise and residual echo */
381    for (i=0;i<len;i++)
382       gain_floor[i] = FRAC_SCALING*sqrt(noise_floor*PSHR32(noise[i],NOISE_SHIFT) + echo_floor*echo[i])/sqrt(1+PSHR32(noise[i],NOISE_SHIFT) + echo[i]);
383 }
384
385 #endif
386 SpeexPreprocessState *speex_preprocess_state_init(int frame_size, int sampling_rate)
387 {
388    int i;
389    int N, N3, N4, M;
390
391    SpeexPreprocessState *st = (SpeexPreprocessState *)speex_alloc(sizeof(SpeexPreprocessState));
392    st->frame_size = frame_size;
393
394    /* Round ps_size down to the nearest power of two */
395 #if 0
396    i=1;
397    st->ps_size = st->frame_size;
398    while(1)
399    {
400       if (st->ps_size & ~i)
401       {
402          st->ps_size &= ~i;
403          i<<=1;
404       } else {
405          break;
406       }
407    }
408    
409    
410    if (st->ps_size < 3*st->frame_size/4)
411       st->ps_size = st->ps_size * 3 / 2;
412 #else
413    st->ps_size = st->frame_size;
414 #endif
415
416    N = st->ps_size;
417    N3 = 2*N - st->frame_size;
418    N4 = st->frame_size - N3;
419    
420    st->sampling_rate = sampling_rate;
421    st->denoise_enabled = 1;
422    st->vad_enabled = 0;
423    st->dereverb_enabled = 0;
424    st->reverb_decay = 0;
425    st->reverb_level = 0;
426    st->noise_suppress = NOISE_SUPPRESS_DEFAULT;
427    st->echo_suppress = ECHO_SUPPRESS_DEFAULT;
428    st->echo_suppress_active = ECHO_SUPPRESS_ACTIVE_DEFAULT;
429
430    st->speech_prob_start = SPEECH_PROB_START_DEFAULT;
431    st->speech_prob_continue = SPEECH_PROB_CONTINUE_DEFAULT;
432
433    st->echo_state = NULL;
434    
435    st->nbands = NB_BANDS;
436    M = st->nbands;
437    st->bank = filterbank_new(M, sampling_rate, N, 1);
438    
439    st->frame = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));
440    st->window = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));
441    st->ft = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));
442    
443    st->ps = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
444    st->noise = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
445    st->echo_noise = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
446    st->residual_echo = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
447    st->reverb_estimate = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
448    st->old_ps = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
449    st->prior = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
450    st->post = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
451    st->gain = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
452    st->gain2 = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
453    st->gain_floor = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
454    st->zeta = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
455    
456    st->S = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t));
457    st->Smin = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t));
458    st->Stmp = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t));
459    st->update_prob = (int*)speex_alloc(N*sizeof(int));
460    
461    st->inbuf = (spx_word16_t*)speex_alloc(N3*sizeof(spx_word16_t));
462    st->outbuf = (spx_word16_t*)speex_alloc(N3*sizeof(spx_word16_t));
463
464    conj_window(st->window, 2*N3);
465    for (i=2*N3;i<2*st->ps_size;i++)
466       st->window[i]=Q15_ONE;
467    
468    if (N4>0)
469    {
470       for (i=N3-1;i>=0;i--)
471       {
472          st->window[i+N3+N4]=st->window[i+N3];
473          st->window[i+N3]=1;
474       }
475    }
476    for (i=0;i<N+M;i++)
477    {
478       st->noise[i]=QCONST32(1.f,NOISE_SHIFT);
479       st->reverb_estimate[i]=0;
480       st->old_ps[i]=1;
481       st->gain[i]=Q15_ONE;
482       st->post[i]=SHL16(1, SNR_SHIFT);
483       st->prior[i]=SHL16(1, SNR_SHIFT);
484    }
485
486    for (i=0;i<N;i++)
487       st->update_prob[i] = 1;
488    for (i=0;i<N3;i++)
489    {
490       st->inbuf[i]=0;
491       st->outbuf[i]=0;
492    }
493 #ifndef FIXED_POINT
494    st->agc_enabled = 0;
495    st->agc_level = 8000;
496    st->loudness_weight = (float*)speex_alloc(N*sizeof(float));
497    for (i=0;i<N;i++)
498    {
499       float ff=((float)i)*.5*sampling_rate/((float)N);
500       /*st->loudness_weight[i] = .5f*(1.f/(1.f+ff/8000.f))+1.f*exp(-.5f*(ff-3800.f)*(ff-3800.f)/9e5f);*/
501       st->loudness_weight[i] = .35f-.35f*ff/16000.f+.73f*exp(-.5f*(ff-3800)*(ff-3800)/9e5f);
502       if (st->loudness_weight[i]<.01f)
503          st->loudness_weight[i]=.01f;
504       st->loudness_weight[i] *= st->loudness_weight[i];
505    }
506    st->loudness = pow(st->agc_level,LOUDNESS_EXP);
507    st->agc_gain = 1;
508    st->nb_loudness_adapt = 0;
509    st->max_gain = 10;
510    st->max_increase_step = exp(0.11513f * 6.*st->frame_size / st->sampling_rate);
511    st->max_decrease_step = exp(-0.11513f * 30.*st->frame_size / st->sampling_rate);
512    st->prev_loudness = 1;
513    st->init_max = 1;
514 #endif
515    st->was_speech = 0;
516
517    st->fft_lookup = spx_fft_init(2*N);
518
519    st->nb_adapt=0;
520    st->min_count=0;
521    return st;
522 }
523
524 void speex_preprocess_state_destroy(SpeexPreprocessState *st)
525 {
526    speex_free(st->frame);
527    speex_free(st->ft);
528    speex_free(st->ps);
529    speex_free(st->gain2);
530    speex_free(st->gain_floor);
531    speex_free(st->window);
532    speex_free(st->noise);
533    speex_free(st->reverb_estimate);
534    speex_free(st->old_ps);
535    speex_free(st->gain);
536    speex_free(st->prior);
537    speex_free(st->post);
538 #ifndef FIXED_POINT
539    speex_free(st->loudness_weight);
540 #endif
541    speex_free(st->echo_noise);
542    speex_free(st->residual_echo);
543
544    speex_free(st->S);
545    speex_free(st->Smin);
546    speex_free(st->Stmp);
547    speex_free(st->update_prob);
548    speex_free(st->zeta);
549
550    speex_free(st->inbuf);
551    speex_free(st->outbuf);
552
553    spx_fft_destroy(st->fft_lookup);
554    filterbank_destroy(st->bank);
555    speex_free(st);
556 }
557
558 /* FIXME: The AGC doesn't work yet with fixed-point*/
559 #ifndef FIXED_POINT
560 static void speex_compute_agc(SpeexPreprocessState *st, spx_word16_t Pframe, spx_word16_t *ft)
561 {
562    int i;
563    int N = st->ps_size;
564    float target_gain;
565    float loudness=1.f;
566    float rate;
567    
568    for (i=2;i<N;i++)
569    {
570       loudness += 2.f*N*st->ps[i]* st->loudness_weight[i];
571    }
572    loudness=sqrt(loudness);
573       /*if (loudness < 2*pow(st->loudness, 1.0/LOUDNESS_EXP) &&
574    loudness*2 > pow(st->loudness, 1.0/LOUDNESS_EXP))*/
575    if (Pframe>.5)
576    {
577       st->nb_loudness_adapt++;
578       rate=2.0f/(1+st->nb_loudness_adapt);
579       st->loudness = (1-rate)*st->loudness + (rate)*pow(loudness, LOUDNESS_EXP);
580       if (st->init_max < st->max_gain)
581          st->init_max *= 1.1;
582    }
583    /*printf ("%f %f %f %f\n", Pframe, loudness, pow(st->loudness, 1.0f/LOUDNESS_EXP), st->loudness2);*/
584    
585    target_gain = st->agc_level*pow(st->loudness, -1.0f/LOUDNESS_EXP);
586
587    if (Pframe>.5 || target_gain < st->agc_gain)
588    {
589       if (target_gain > st->max_increase_step*st->agc_gain)
590          target_gain = st->max_increase_step*st->agc_gain;
591       if (target_gain < st->max_decrease_step*st->agc_gain && loudness < 10*st->prev_loudness)
592          target_gain = st->max_decrease_step*st->agc_gain;
593       if (target_gain > st->max_gain)
594          target_gain = st->max_gain;
595       if (target_gain > st->init_max)
596          target_gain = st->init_max;
597    
598       st->agc_gain = target_gain;
599    }
600    /*fprintf (stderr, "%f %f %f\n", loudness, st->loudness2, rate);*/
601       
602    for (i=0;i<2*N;i++)
603       ft[i] *= st->agc_gain;
604    st->prev_loudness = loudness;
605 }
606 #endif
607
608 static void preprocess_analysis(SpeexPreprocessState *st, spx_int16_t *x)
609 {
610    int i;
611    int N = st->ps_size;
612    int N3 = 2*N - st->frame_size;
613    int N4 = st->frame_size - N3;
614    spx_word32_t *ps=st->ps;
615
616    /* 'Build' input frame */
617    for (i=0;i<N3;i++)
618       st->frame[i]=st->inbuf[i];
619    for (i=0;i<st->frame_size;i++)
620       st->frame[N3+i]=x[i];
621    
622    /* Update inbuf */
623    for (i=0;i<N3;i++)
624       st->inbuf[i]=x[N4+i];
625
626    /* Windowing */
627    for (i=0;i<2*N;i++)
628       st->frame[i] = MULT16_16_Q15(st->frame[i], st->window[i]);
629
630 #ifdef FIXED_POINT
631    {
632       spx_word16_t max_val=0;
633       for (i=0;i<2*N;i++)
634          max_val = MAX16(max_val, ABS16(st->frame[i]));
635       st->frame_shift = 14-spx_ilog2(EXTEND32(max_val));
636       for (i=0;i<2*N;i++)
637          st->frame[i] = SHL16(st->frame[i], st->frame_shift);
638    }
639 #endif
640    
641    /* Perform FFT */
642    spx_fft(st->fft_lookup, st->frame, st->ft);
643          
644    /* Power spectrum */
645    ps[0]=MULT16_16(st->ft[0],st->ft[0]);
646    for (i=1;i<N;i++)
647       ps[i]=MULT16_16(st->ft[2*i-1],st->ft[2*i-1]) + MULT16_16(st->ft[2*i],st->ft[2*i]);
648    for (i=0;i<N;i++)
649       st->ps[i] = PSHR32(st->ps[i], 2*st->frame_shift);
650
651    filterbank_compute_bank32(st->bank, ps, ps+N);
652 }
653
654 static void update_noise_prob(SpeexPreprocessState *st)
655 {
656    int i;
657    int min_range;
658    int N = st->ps_size;
659
660    for (i=1;i<N-1;i++)
661       st->S[i] =  MULT16_32_Q15(QCONST16(.8f,15),st->S[i]) + MULT16_32_Q15(QCONST16(.05f,15),st->ps[i-1]) 
662                       + MULT16_32_Q15(QCONST16(.1f,15),st->ps[i]) + MULT16_32_Q15(QCONST16(.05f,15),st->ps[i+1]);
663    st->S[0] =  MULT16_32_Q15(QCONST16(.8f,15),st->S[0]) + MULT16_32_Q15(QCONST16(.2f,15),st->ps[0]);
664    st->S[N-1] =  MULT16_32_Q15(QCONST16(.8f,15),st->S[N-1]) + MULT16_32_Q15(QCONST16(.2f,15),st->ps[N-1]);
665    
666    if (st->nb_adapt==1)
667    {
668       for (i=0;i<N;i++)
669          st->Smin[i] = st->Stmp[i] = 0;
670    }
671
672    if (st->nb_adapt < 100)
673       min_range = 15;
674    else if (st->nb_adapt < 1000)
675       min_range = 50;
676    else if (st->nb_adapt < 10000)
677       min_range = 150;
678    else
679       min_range = 300;
680    if (st->min_count > min_range)
681    {
682       st->min_count = 0;
683       for (i=0;i<N;i++)
684       {
685          st->Smin[i] = MIN32(st->Stmp[i], st->S[i]);
686          st->Stmp[i] = st->S[i];
687       }
688    } else {
689       for (i=0;i<N;i++)
690       {
691          st->Smin[i] = MIN32(st->Smin[i], st->S[i]);
692          st->Stmp[i] = MIN32(st->Stmp[i], st->S[i]);      
693       }
694    }
695    for (i=0;i<N;i++)
696    {
697       if (MULT16_32_Q15(QCONST16(.4f,15),st->S[i]) > ADD32(st->Smin[i],EXTEND32(20)))
698          st->update_prob[i] = 1;
699       else
700          st->update_prob[i] = 0;
701       /*fprintf (stderr, "%f ", st->S[i]/st->Smin[i]);*/
702       /*fprintf (stderr, "%f ", st->update_prob[i]);*/
703    }
704
705 }
706
707 #define NOISE_OVERCOMPENS 1.
708
709 void speex_echo_get_residual(SpeexEchoState *st, spx_word32_t *Yout, int len);
710
711 int speex_preprocess(SpeexPreprocessState *st, spx_int16_t *x, spx_int32_t *echo)
712 {
713    return speex_preprocess_run(st, x);
714 }
715
716 int speex_preprocess_run(SpeexPreprocessState *st, spx_int16_t *x)
717 {
718    int i;
719    int M;
720    int N = st->ps_size;
721    int N3 = 2*N - st->frame_size;
722    int N4 = st->frame_size - N3;
723    spx_word32_t *ps=st->ps;
724    spx_word32_t Zframe;
725    spx_word16_t Pframe;
726    spx_word16_t beta, beta_1;
727    spx_word16_t effective_echo_suppress;
728    
729    st->nb_adapt++;
730    st->min_count++;
731    
732    beta = MAX16(QCONST16(.03,15),DIV32_16(Q15_ONE,st->nb_adapt));
733    beta_1 = Q15_ONE-beta;
734    M = st->nbands;
735    /* Deal with residual echo if provided */
736    if (st->echo_state)
737    {
738       speex_echo_get_residual(st->echo_state, st->residual_echo, N);
739 #ifndef FIXED_POINT
740       /* If there are NaNs or ridiculous values, it'll show up in the DC and we just reset everything to zero */
741       if (!(st->residual_echo[0] >=0 && st->residual_echo[0]<N*1e9f))
742       {
743          for (i=0;i<N;i++)
744             st->residual_echo[i] = 0;
745       }
746 #endif
747       for (i=0;i<N;i++)
748          st->echo_noise[i] = MAX32(MULT16_32_Q15(QCONST16(.6f,15),st->echo_noise[i]), st->residual_echo[i]);
749       filterbank_compute_bank32(st->bank, st->echo_noise, st->echo_noise+N);
750    } else {
751       for (i=0;i<N+M;i++)
752          st->echo_noise[i] = 0;
753    }
754    preprocess_analysis(st, x);
755
756    update_noise_prob(st);
757
758    /* Noise estimation always updated for the 10 first frames */
759    /*if (st->nb_adapt<10)
760    {
761       for (i=1;i<N-1;i++)
762          st->update_prob[i] = 0;
763    }
764    */
765    
766    /* Update the noise estimate for the frequencies where it can be */
767    for (i=0;i<N;i++)
768    {
769       if (!st->update_prob[i] || st->ps[i] < PSHR32(st->noise[i], NOISE_SHIFT))
770          st->noise[i] = MAX32(EXTEND32(0),MULT16_32_Q15(beta_1,st->noise[i]) + MULT16_32_Q15(beta,SHL32(st->ps[i],NOISE_SHIFT)));
771    }
772    filterbank_compute_bank32(st->bank, st->noise, st->noise+N);
773
774    /* Special case for first frame */
775    if (st->nb_adapt==1)
776       for (i=0;i<N+M;i++)
777          st->old_ps[i] = ps[i];
778
779    /* Compute a posteriori SNR */
780    for (i=0;i<N+M;i++)
781    {
782       spx_word16_t gamma;
783       
784       /* Total noise estimate including residual echo and reverberation */
785       spx_word32_t tot_noise = ADD32(ADD32(ADD32(EXTEND32(1), PSHR32(st->noise[i],NOISE_SHIFT)) , st->echo_noise[i]) , st->reverb_estimate[i]);
786       
787       /* A posteriori SNR = ps/noise - 1*/
788       st->post[i] = SUB16(DIV32_16_Q8(ps[i],tot_noise), QCONST16(1.f,SNR_SHIFT));
789       st->post[i]=MIN16(st->post[i], QCONST16(100.f,SNR_SHIFT));
790       
791       /* Computing update gamma = .1 + .9*(old/(old+noise))^2 */
792       gamma = QCONST16(.1f,15)+MULT16_16_Q15(QCONST16(.89f,15),SQR16_Q15(DIV32_16_Q15(st->old_ps[i],ADD32(st->old_ps[i],tot_noise))));
793       
794       /* A priori SNR update = gamma*max(0,post) + (1-gamma)*old/noise */
795       st->prior[i] = EXTRACT16(PSHR32(ADD32(MULT16_16(gamma,MAX16(0,st->post[i])), MULT16_16(Q15_ONE-gamma,DIV32_16_Q8(st->old_ps[i],tot_noise))), 15));
796       st->prior[i]=MIN16(st->prior[i], QCONST16(100.f,SNR_SHIFT));
797    }
798
799    /*print_vec(st->post, N+M, "");*/
800
801    /* Recursive average of the a priori SNR. A bit smoothed for the psd components */
802    st->zeta[0] = PSHR32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[0]), MULT16_16(QCONST16(.3f,15),st->prior[0])),15);
803    for (i=1;i<N-1;i++)
804       st->zeta[i] = PSHR32(ADD32(ADD32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[i]), MULT16_16(QCONST16(.15f,15),st->prior[i])),
805                            MULT16_16(QCONST16(.075f,15),st->prior[i-1])), MULT16_16(QCONST16(.075f,15),st->prior[i+1])),15);
806    for (i=N-1;i<N+M;i++)
807       st->zeta[i] = PSHR32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[i]), MULT16_16(QCONST16(.3f,15),st->prior[i])),15);
808
809    /* Speech probability of presence for the entire frame is based on the average filterbank a priori SNR */
810    Zframe = 0;
811    for (i=N;i<N+M;i++)
812       Zframe = ADD32(Zframe, EXTEND32(st->zeta[i]));
813    Pframe = QCONST16(.1f,15)+MULT16_16_Q15(QCONST16(.899f,15),qcurve(DIV32_16(Zframe,st->nbands)));
814    
815    effective_echo_suppress = EXTRACT16(PSHR32(ADD32(MULT16_16(SUB16(Q15_ONE,Pframe), st->echo_suppress), MULT16_16(Pframe, st->echo_suppress_active)),15));
816    
817    compute_gain_floor(st->noise_suppress, effective_echo_suppress, st->noise+N, st->echo_noise+N, st->gain_floor+N, M);
818          
819    /* Compute Ephraim & Malah gain speech probability of presence for each critical band (Bark scale) 
820       Technically this is actually wrong because the EM gaim assumes a slightly different probability 
821       distribution */
822    for (i=N;i<N+M;i++)
823    {
824       /* See EM and Cohen papers*/
825       spx_word32_t theta;
826       /* Gain from hypergeometric function */
827       spx_word32_t MM;
828       /* Weiner filter gain */
829       spx_word16_t prior_ratio;
830       /* a priority probability of speech presence based on Bark sub-band alone */
831       spx_word16_t P1;
832       /* Speech absence a priori probability (considering sub-band and frame) */
833       spx_word16_t q;
834 #ifdef FIXED_POINT
835       spx_word16_t tmp;
836 #endif
837       
838       prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1,SNR_SHIFT)));
839       theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT));
840
841       MM = hypergeom_gain(theta);
842       /* Gain with bound */
843       st->gain[i] = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM)));
844       /* Save old Bark power spectrum */
845       st->old_ps[i] = MULT16_32_P15(QCONST16(.2f,15),st->old_ps[i]) + MULT16_32_P15(MULT16_16_P15(QCONST16(.8f,15),SQR16_Q15(st->gain[i])),ps[i]);
846
847       P1 = QCONST16(.199f,15)+MULT16_16_Q15(QCONST16(.8f,15),qcurve (st->zeta[i]));
848       q = Q15_ONE-MULT16_16_Q15(Pframe,P1);
849 #ifdef FIXED_POINT
850       theta = MIN32(theta, EXTEND32(32767));
851 /*Q8*/tmp = MULT16_16_Q15((SHL32(1,SNR_SHIFT)+st->prior[i]),EXTRACT16(MIN32(Q15ONE,SHR32(spx_exp(-EXTRACT16(theta)),1))));
852       tmp = MIN16(QCONST16(3.,SNR_SHIFT), tmp); /* Prevent overflows in the next line*/
853 /*Q8*/tmp = EXTRACT16(PSHR32(MULT16_16(PDIV32_16(SHL32(EXTEND32(q),8),(Q15_ONE-q)),tmp),8));
854       st->gain2[i]=DIV32_16(SHL32(EXTEND32(32767),SNR_SHIFT), ADD16(256,tmp));
855 #else
856       st->gain2[i]=1/(1.f + (q/(1.f-q))*(1+st->prior[i])*exp(-theta));
857 #endif
858    }
859    /* Convert the EM gains and speech prob to linear frequency */
860    filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2);
861    filterbank_compute_psd16(st->bank,st->gain+N, st->gain);
862    
863    /* Use 1 for linear gain resolution (best) or 0 for Bark gain resolution (faster) */
864    if (1)
865    {
866       filterbank_compute_psd16(st->bank,st->gain_floor+N, st->gain_floor);
867    
868       /* Compute gain according to the Ephraim-Malah algorithm -- linear frequency */
869       for (i=0;i<N;i++)
870       {
871          spx_word32_t MM;
872          spx_word32_t theta;
873          spx_word16_t prior_ratio;
874          spx_word16_t tmp;
875          spx_word16_t p;
876          spx_word16_t g;
877          
878          /* Wiener filter gain */
879          prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1,SNR_SHIFT)));
880          theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT));
881
882          /* Optimal estimator for loudness domain */
883          MM = hypergeom_gain(theta);
884          /* EM gain with bound */
885          g = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM)));
886          /* Interpolated speech probability of presence */
887          p = st->gain2[i];
888                   
889          /* Constrain the gain to be close to the Bark scale gain */
890          if (MULT16_16_Q15(QCONST16(.333f,15),g) > st->gain[i])
891             g = MULT16_16(3,st->gain[i]);
892          st->gain[i] = g;
893          
894          /* Save old power spectrum */
895          st->old_ps[i] = MULT16_32_P15(QCONST16(.2f,15),st->old_ps[i]) + MULT16_32_P15(MULT16_16_P15(QCONST16(.8f,15),SQR16_Q15(st->gain[i])),ps[i]);
896          
897          /* Apply gain floor */
898          if (st->gain[i] < st->gain_floor[i])
899             st->gain[i] = st->gain_floor[i];
900
901          /* Exponential decay model for reverberation (unused) */
902          /*st->reverb_estimate[i] = st->reverb_decay*st->reverb_estimate[i] + st->reverb_decay*st->reverb_level*st->gain[i]*st->gain[i]*st->ps[i];*/
903          
904          /* Take into account speech probability of presence (loudness domain MMSE estimator) */
905          /* gain2 = [p*sqrt(gain)+(1-p)*sqrt(gain _floor) ]^2 */
906          tmp = MULT16_16_P15(p,spx_sqrt(SHL32(EXTEND32(st->gain[i]),15))) + MULT16_16_P15(SUB16(Q15_ONE,p),spx_sqrt(SHL32(EXTEND32(st->gain_floor[i]),15)));
907          st->gain2[i]=SQR16_Q15(tmp);
908
909          /* Use this if you want a log-domain MMSE estimator instead */
910          /*st->gain2[i] = pow(st->gain[i], p) * pow(st->gain_floor[i],1.f-p);*/
911       }
912    } else {
913       for (i=N;i<N+M;i++)
914       {
915          spx_word16_t tmp;
916          spx_word16_t p = st->gain2[i];
917          st->gain[i] = MAX16(st->gain[i], st->gain_floor[i]);         
918          tmp = MULT16_16_P15(p,spx_sqrt(SHL32(EXTEND32(st->gain[i]),15))) + MULT16_16_P15(SUB16(Q15_ONE,p),spx_sqrt(SHL32(EXTEND32(st->gain_floor[i]),15)));
919          st->gain2[i]=SQR16_Q15(tmp);
920       }
921       filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2);
922    }
923    
924    /* If noise suppression is off, don't apply the gain (but then why call this in the first place!) */
925    if (!st->denoise_enabled)
926    {
927       for (i=0;i<N+M;i++)
928          st->gain2[i]=Q15_ONE;
929    }
930       
931    /* Apply computed gain */
932    for (i=1;i<N;i++)
933    {
934       st->ft[2*i-1] = MULT16_16_P15(st->gain2[i],st->ft[2*i-1]);
935       st->ft[2*i] = MULT16_16_P15(st->gain2[i],st->ft[2*i]);
936    }
937    st->ft[0] = MULT16_16_P15(st->gain2[0],st->ft[0]);
938    st->ft[2*N-1] = MULT16_16_P15(st->gain2[N-1],st->ft[2*N-1]);
939    
940    /*FIXME: This *will* not work for fixed-point */
941 #ifndef FIXED_POINT
942    if (st->agc_enabled)
943       speex_compute_agc(st, Pframe, st->ft);
944 #endif
945
946    /* Inverse FFT with 1/N scaling */
947    spx_ifft(st->fft_lookup, st->ft, st->frame);
948    /* Scale back to original (lower) amplitude */
949    for (i=0;i<2*N;i++)
950       st->frame[i] = PSHR16(st->frame[i], st->frame_shift);
951
952    /*FIXME: This *will* not work for fixed-point */
953 #ifndef FIXED_POINT
954    if (st->agc_enabled)
955    {
956       float max_sample=0;
957       for (i=0;i<2*N;i++)
958          if (fabs(st->frame[i])>max_sample)
959             max_sample = fabs(st->frame[i]);
960       if (max_sample>28000.f)
961       {
962          float damp = 28000.f/max_sample;
963          for (i=0;i<2*N;i++)
964             st->frame[i] *= damp;
965       }
966    }
967 #endif
968    
969    /* Synthesis window (for WOLA) */
970    for (i=0;i<2*N;i++)
971       st->frame[i] = MULT16_16_Q15(st->frame[i], st->window[i]);
972
973    /* Perform overlap and add */
974    for (i=0;i<N3;i++)
975       x[i] = st->outbuf[i] + st->frame[i];
976    for (i=0;i<N4;i++)
977       x[N3+i] = st->frame[N3+i];
978    
979    /* Update outbuf */
980    for (i=0;i<N3;i++)
981       st->outbuf[i] = st->frame[st->frame_size+i];
982
983    /* FIXME: This VAD is a kludge */
984    if (st->vad_enabled)
985    {
986       if (Pframe > st->speech_prob_start || (st->was_speech && Pframe > st->speech_prob_continue))
987       {
988          st->was_speech=1;
989          return 1;
990       } else
991       {
992          st->was_speech=0;
993          return 0;
994       }
995    } else {
996       return 1;
997    }
998 }
999
1000 void speex_preprocess_estimate_update(SpeexPreprocessState *st, spx_int16_t *x)
1001 {
1002    int i;
1003    int N = st->ps_size;
1004    int N3 = 2*N - st->frame_size;
1005    int M;
1006    spx_word32_t *ps=st->ps;
1007
1008    M = st->nbands;
1009    st->min_count++;
1010    
1011    preprocess_analysis(st, x);
1012
1013    update_noise_prob(st);
1014    
1015    for (i=1;i<N-1;i++)
1016    {
1017       if (!st->update_prob[i] || st->ps[i] < PSHR32(st->noise[i],NOISE_SHIFT))
1018       {
1019          st->noise[i] = MULT16_32_Q15(QCONST16(.95f,15),st->noise[i]) + MULT16_32_Q15(QCONST16(.05f,15),SHL32(st->ps[i],NOISE_SHIFT));
1020       }
1021    }
1022
1023    for (i=0;i<N3;i++)
1024       st->outbuf[i] = MULT16_16_Q15(x[st->frame_size-N3+i],st->window[st->frame_size+i]);
1025
1026    /* Save old power spectrum */
1027    for (i=0;i<N+M;i++)
1028       st->old_ps[i] = ps[i];
1029
1030    for (i=0;i<N;i++)
1031       st->reverb_estimate[i] = MULT16_32_Q15(st->reverb_decay, st->reverb_estimate[i]);
1032 }
1033
1034
1035 int speex_preprocess_ctl(SpeexPreprocessState *state, int request, void *ptr)
1036 {
1037    int i;
1038    SpeexPreprocessState *st;
1039    st=(SpeexPreprocessState*)state;
1040    switch(request)
1041    {
1042    case SPEEX_PREPROCESS_SET_DENOISE:
1043       st->denoise_enabled = (*(spx_int32_t*)ptr);
1044       break;
1045    case SPEEX_PREPROCESS_GET_DENOISE:
1046       (*(spx_int32_t*)ptr) = st->denoise_enabled;
1047       break;
1048 #ifndef FIXED_POINT
1049    case SPEEX_PREPROCESS_SET_AGC:
1050       st->agc_enabled = (*(spx_int32_t*)ptr);
1051       break;
1052    case SPEEX_PREPROCESS_GET_AGC:
1053       (*(spx_int32_t*)ptr) = st->agc_enabled;
1054       break;
1055
1056    case SPEEX_PREPROCESS_SET_AGC_LEVEL:
1057       st->agc_level = (*(float*)ptr);
1058       if (st->agc_level<1)
1059          st->agc_level=1;
1060       if (st->agc_level>32768)
1061          st->agc_level=32768;
1062       break;
1063    case SPEEX_PREPROCESS_GET_AGC_LEVEL:
1064       (*(float*)ptr) = st->agc_level;
1065       break;
1066 #endif
1067    case SPEEX_PREPROCESS_SET_VAD:
1068       speex_warning("The VAD has been replaced by a hack pending a complete rewrite");
1069       st->vad_enabled = (*(spx_int32_t*)ptr);
1070       break;
1071    case SPEEX_PREPROCESS_GET_VAD:
1072       (*(spx_int32_t*)ptr) = st->vad_enabled;
1073       break;
1074    
1075    case SPEEX_PREPROCESS_SET_DEREVERB:
1076       st->dereverb_enabled = (*(spx_int32_t*)ptr);
1077       for (i=0;i<st->ps_size;i++)
1078          st->reverb_estimate[i]=0;
1079       break;
1080    case SPEEX_PREPROCESS_GET_DEREVERB:
1081       (*(spx_int32_t*)ptr) = st->dereverb_enabled;
1082       break;
1083
1084    case SPEEX_PREPROCESS_SET_DEREVERB_LEVEL:
1085       st->reverb_level = (*(float*)ptr);
1086       break;
1087    case SPEEX_PREPROCESS_GET_DEREVERB_LEVEL:
1088       (*(float*)ptr) = st->reverb_level;
1089       break;
1090    
1091    case SPEEX_PREPROCESS_SET_DEREVERB_DECAY:
1092       st->reverb_decay = (*(float*)ptr);
1093       break;
1094    case SPEEX_PREPROCESS_GET_DEREVERB_DECAY:
1095       (*(float*)ptr) = st->reverb_decay;
1096       break;
1097
1098    case SPEEX_PREPROCESS_SET_PROB_START:
1099       *(spx_int32_t*)ptr = MIN32(Q15_ONE,MAX32(0, *(spx_int32_t*)ptr));
1100       st->speech_prob_start = DIV32_16(MULT16_16(32767,*(spx_int32_t*)ptr), 100);
1101       break;
1102    case SPEEX_PREPROCESS_GET_PROB_START:
1103       (*(spx_int32_t*)ptr) = MULT16_16_Q15(st->speech_prob_start, 100);
1104       break;
1105
1106    case SPEEX_PREPROCESS_SET_PROB_CONTINUE:
1107       *(spx_int32_t*)ptr = MIN32(Q15_ONE,MAX32(0, *(spx_int32_t*)ptr));
1108       st->speech_prob_continue = DIV32_16(MULT16_16(32767,*(spx_int32_t*)ptr), 100);
1109       break;
1110    case SPEEX_PREPROCESS_GET_PROB_CONTINUE:
1111       (*(spx_int32_t*)ptr) = MULT16_16_Q15(st->speech_prob_continue, 100);
1112       break;
1113
1114    case SPEEX_PREPROCESS_SET_NOISE_SUPPRESS:
1115       st->noise_suppress = -ABS(*(spx_int32_t*)ptr);
1116       break;
1117    case SPEEX_PREPROCESS_GET_NOISE_SUPPRESS:
1118       (*(spx_int32_t*)ptr) = st->noise_suppress;
1119       break;
1120    case SPEEX_PREPROCESS_SET_ECHO_SUPPRESS:
1121       st->echo_suppress = -ABS(*(spx_int32_t*)ptr);
1122       break;
1123    case SPEEX_PREPROCESS_GET_ECHO_SUPPRESS:
1124       (*(spx_int32_t*)ptr) = st->echo_suppress;
1125       break;
1126    case SPEEX_PREPROCESS_SET_ECHO_SUPPRESS_ACTIVE:
1127       st->echo_suppress_active = -ABS(*(spx_int32_t*)ptr);
1128       break;
1129    case SPEEX_PREPROCESS_GET_ECHO_SUPPRESS_ACTIVE:
1130       (*(spx_int32_t*)ptr) = st->echo_suppress_active;
1131       break;
1132    case SPEEX_PREPROCESS_SET_ECHO_STATE:
1133       st->echo_state = (SpeexEchoState*)ptr;
1134       break;
1135    case SPEEX_PREPROCESS_GET_ECHO_STATE:
1136       ptr = (void*)st->echo_state;
1137       break;
1138
1139    default:
1140       speex_warning_int("Unknown speex_preprocess_ctl request: ", request);
1141       return -1;
1142    }
1143    return 0;
1144 }