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