cleanup, comments, better handling of clipping in AEC.
[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       -25
81 #define ECHO_SUPPRESS_DEFAULT        -45
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    agc_enabled;
189    float  agc_level;
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 #ifndef FIXED_POINT
229    float *loudness_weight;   /**< Perceptual loudness curve */
230    float  loudness;          /**< loudness estimate */
231    float  loudness2;         /**< loudness estimate */
232    int    nb_loudness_adapt; /**< Number of frames used for loudness adaptation so far */
233 #endif
234    int    nb_adapt;          /**< Number of frames used for adaptation so far */
235    int    was_speech;
236    int    min_count;         /**< Number of frames processed so far */
237    void  *fft_lookup;        /**< Lookup table for the FFT */
238 #ifdef FIXED_POINT
239    int    frame_shift;
240 #endif
241 };
242
243
244 static void conj_window(spx_word16_t *w, int len)
245 {
246    int i;
247    for (i=0;i<len;i++)
248    {
249       spx_word16_t tmp;
250       spx_word16_t x = DIV32_16(MULT16_16(QCONST16(4.f,13),i),len);
251       int inv=0;
252       if (x<QCONST16(1.f,13))
253       {
254       } else if (x<QCONST16(2.f,13))
255       {
256          x=QCONST16(2.f,13)-x;
257          inv=1;
258       } else if (x<QCONST16(3.f,13))
259       {
260          x=x-QCONST16(2.f,13);
261          inv=1;
262       } else {
263          x=QCONST16(2.f,13)-x+QCONST16(2.f,13); /* 4 - x */
264       }
265       x = MULT16_16_Q14(QCONST16(1.271903f,14), x);
266       tmp = SQR16_Q15(QCONST16(.5f,15)-MULT16_16_P15(QCONST16(.5f,15),spx_cos_norm(QCONST32(x,2))));
267       if (inv)
268          tmp=SUB16(Q15_ONE,tmp);
269       w[i]=spx_sqrt(SHL32(EXTEND32(tmp),15));
270    }
271 }
272
273       
274 #ifdef FIXED_POINT
275 /* This function approximates the gain function 
276    y = gamma(1.25)^2 * M(-.25;1;-x) / sqrt(x)  
277    which multiplied by xi/(1+xi) is the optimal gain
278    in the loudness domain ( sqrt[amplitude] )
279    Input in Q11 format, output in Q15
280 */
281 static inline spx_word32_t hypergeom_gain(spx_word32_t xx)
282 {
283    int ind;
284    spx_word16_t frac;
285    /* Q13 table */
286    static const spx_word16_t table[21] = {
287        6730,  8357,  9868, 11267, 12563, 13770, 14898,
288       15959, 16961, 17911, 18816, 19682, 20512, 21311,
289       22082, 22827, 23549, 24250, 24931, 25594, 26241};
290       ind = SHR32(xx,10);
291       if (ind<0)
292          return Q15_ONE;
293       if (ind>19)
294          return ADD32(EXTEND32(Q15_ONE),EXTEND32(DIV32_16(QCONST32(.1296,23), SHR32(xx,EXPIN_SHIFT-SNR_SHIFT))));
295       frac = SHL32(xx-SHL32(ind,10),5);
296       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);
297 }
298
299 static inline spx_word16_t qcurve(spx_word16_t x)
300 {
301    x = MAX16(x, 1);
302    return DIV32_16(SHL32(EXTEND32(32767),9),ADD16(512,MULT16_16_Q15(QCONST16(.60f,15),DIV32_16(32767,x))));
303 }
304
305 /* Compute the gain floor based on different floors for the background noise and residual echo */
306 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)
307 {
308    int i;
309    
310    if (noise_suppress > effective_echo_suppress)
311    {
312       spx_word16_t noise_gain, gain_ratio;
313       noise_gain = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(0.11513,11),noise_suppress)),1)));
314       gain_ratio = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(.2302585f,11),effective_echo_suppress-noise_suppress)),1)));
315
316       /* gain_floor = sqrt [ (noise*noise_floor + echo*echo_floor) / (noise+echo) ] */
317       for (i=0;i<len;i++)
318          gain_floor[i] = MULT16_16_Q15(noise_gain,
319                                        spx_sqrt(SHL32(EXTEND32(DIV32_16_Q15(PSHR32(noise[i],NOISE_SHIFT) + MULT16_32_Q15(gain_ratio,echo[i]),
320                                              (1+PSHR32(noise[i],NOISE_SHIFT) + echo[i]) )),15)));
321    } else {
322       spx_word16_t echo_gain, gain_ratio;
323       echo_gain = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(0.11513,11),effective_echo_suppress)),1)));
324       gain_ratio = EXTRACT16(MIN32(Q15_ONE,SHR32(spx_exp(MULT16_16(QCONST16(.2302585f,11),noise_suppress-effective_echo_suppress)),1)));
325
326       /* gain_floor = sqrt [ (noise*noise_floor + echo*echo_floor) / (noise+echo) ] */
327       for (i=0;i<len;i++)
328          gain_floor[i] = MULT16_16_Q15(echo_gain,
329                                        spx_sqrt(SHL32(EXTEND32(DIV32_16_Q15(MULT16_32_Q15(gain_ratio,PSHR32(noise[i],NOISE_SHIFT)) + echo[i],
330                                              (1+PSHR32(noise[i],NOISE_SHIFT) + echo[i]) )),15)));
331    }
332 }
333
334 #else
335 /* This function approximates the gain function 
336    y = gamma(1.25)^2 * M(-.25;1;-x) / sqrt(x)  
337    which multiplied by xi/(1+xi) is the optimal gain
338    in the loudness domain ( sqrt[amplitude] )
339 */
340 static inline spx_word32_t hypergeom_gain(spx_word32_t xx)
341 {
342    int ind;
343    float integer, frac;
344    float x;
345    static const float table[21] = {
346       0.82157f, 1.02017f, 1.20461f, 1.37534f, 1.53363f, 1.68092f, 1.81865f,
347       1.94811f, 2.07038f, 2.18638f, 2.29688f, 2.40255f, 2.50391f, 2.60144f,
348       2.69551f, 2.78647f, 2.87458f, 2.96015f, 3.04333f, 3.12431f, 3.20326f};
349       x = EXPIN_SCALING_1*xx;
350       integer = floor(2*x);
351       ind = (int)integer;
352       if (ind<0)
353          return FRAC_SCALING;
354       if (ind>19)
355          return FRAC_SCALING*(1+.1296/x);
356       frac = 2*x-integer;
357       return FRAC_SCALING*((1-frac)*table[ind] + frac*table[ind+1])/sqrt(x+.0001f);
358 }
359
360 static inline spx_word16_t qcurve(spx_word16_t x)
361 {
362    return 1.f/(1.f+.15f/(SNR_SCALING_1*x));
363 }
364
365 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)
366 {
367    int i;
368    float echo_floor;
369    float noise_floor;
370
371    noise_floor = exp(.2302585f*noise_suppress);
372    echo_floor = exp(.2302585f*effective_echo_suppress);
373
374    /* Compute the gain floor based on different floors for the background noise and residual echo */
375    for (i=0;i<len;i++)
376       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]);
377 }
378
379 #endif
380 SpeexPreprocessState *speex_preprocess_state_init(int frame_size, int sampling_rate)
381 {
382    int i;
383    int N, N3, N4, M;
384
385    SpeexPreprocessState *st = (SpeexPreprocessState *)speex_alloc(sizeof(SpeexPreprocessState));
386    st->frame_size = frame_size;
387
388    /* Round ps_size down to the nearest power of two */
389 #if 0
390    i=1;
391    st->ps_size = st->frame_size;
392    while(1)
393    {
394       if (st->ps_size & ~i)
395       {
396          st->ps_size &= ~i;
397          i<<=1;
398       } else {
399          break;
400       }
401    }
402    
403    
404    if (st->ps_size < 3*st->frame_size/4)
405       st->ps_size = st->ps_size * 3 / 2;
406 #else
407    st->ps_size = st->frame_size;
408 #endif
409
410    N = st->ps_size;
411    N3 = 2*N - st->frame_size;
412    N4 = st->frame_size - N3;
413    
414    st->sampling_rate = sampling_rate;
415    st->denoise_enabled = 1;
416    st->agc_enabled = 0;
417    st->agc_level = 8000;
418    st->vad_enabled = 0;
419    st->dereverb_enabled = 0;
420    st->reverb_decay = 0;
421    st->reverb_level = 0;
422    st->noise_suppress = NOISE_SUPPRESS_DEFAULT;
423    st->echo_suppress = ECHO_SUPPRESS_DEFAULT;
424    st->echo_suppress_active = ECHO_SUPPRESS_ACTIVE_DEFAULT;
425
426    st->speech_prob_start = SPEECH_PROB_START_DEFAULT;
427    st->speech_prob_continue = SPEECH_PROB_CONTINUE_DEFAULT;
428
429    st->echo_state = NULL;
430    
431    st->nbands = NB_BANDS;
432    M = st->nbands;
433    st->bank = filterbank_new(M, sampling_rate, N, 1);
434    
435    st->frame = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));
436    st->window = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));
437    st->ft = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));
438    
439    st->ps = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
440    st->noise = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
441    st->echo_noise = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
442    st->residual_echo = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
443    st->reverb_estimate = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
444    st->old_ps = (spx_word32_t*)speex_alloc((N+M)*sizeof(spx_word32_t));
445    st->prior = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
446    st->post = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
447    st->gain = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
448    st->gain2 = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
449    st->gain_floor = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
450    st->zeta = (spx_word16_t*)speex_alloc((N+M)*sizeof(spx_word16_t));
451    
452    st->S = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t));
453    st->Smin = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t));
454    st->Stmp = (spx_word32_t*)speex_alloc(N*sizeof(spx_word32_t));
455    st->update_prob = (int*)speex_alloc(N*sizeof(int));
456    
457    st->inbuf = (spx_word16_t*)speex_alloc(N3*sizeof(spx_word16_t));
458    st->outbuf = (spx_word16_t*)speex_alloc(N3*sizeof(spx_word16_t));
459
460    conj_window(st->window, 2*N3);
461    for (i=2*N3;i<2*st->ps_size;i++)
462       st->window[i]=Q15_ONE;
463    
464    if (N4>0)
465    {
466       for (i=N3-1;i>=0;i--)
467       {
468          st->window[i+N3+N4]=st->window[i+N3];
469          st->window[i+N3]=1;
470       }
471    }
472    for (i=0;i<N+M;i++)
473    {
474       st->noise[i]=QCONST32(1.f,NOISE_SHIFT);
475       st->reverb_estimate[i]=0;
476       st->old_ps[i]=1;
477       st->gain[i]=Q15_ONE;
478       st->post[i]=SHL16(1, SNR_SHIFT);
479       st->prior[i]=SHL16(1, SNR_SHIFT);
480    }
481
482    for (i=0;i<N;i++)
483       st->update_prob[i] = 1;
484    for (i=0;i<N3;i++)
485    {
486       st->inbuf[i]=0;
487       st->outbuf[i]=0;
488    }
489 #ifndef FIXED_POINT
490    st->loudness_weight = (float*)speex_alloc(N*sizeof(float));
491    for (i=0;i<N;i++)
492    {
493       float ff=((float)i)*.5*sampling_rate/((float)N);
494       st->loudness_weight[i] = .35f-.35f*ff/16000.f+.73f*exp(-.5f*(ff-3800)*(ff-3800)/9e5f);
495       if (st->loudness_weight[i]<.01f)
496          st->loudness_weight[i]=.01f;
497       st->loudness_weight[i] *= st->loudness_weight[i];
498    }
499    st->loudness = pow(6000,LOUDNESS_EXP);
500    st->loudness2 = 6000;
501    st->nb_loudness_adapt = 0;
502 #endif
503    st->was_speech = 0;
504
505    st->fft_lookup = spx_fft_init(2*N);
506
507    st->nb_adapt=0;
508    st->min_count=0;
509    return st;
510 }
511
512 void speex_preprocess_state_destroy(SpeexPreprocessState *st)
513 {
514    speex_free(st->frame);
515    speex_free(st->ft);
516    speex_free(st->ps);
517    speex_free(st->gain2);
518    speex_free(st->gain_floor);
519    speex_free(st->window);
520    speex_free(st->noise);
521    speex_free(st->reverb_estimate);
522    speex_free(st->old_ps);
523    speex_free(st->gain);
524    speex_free(st->prior);
525    speex_free(st->post);
526 #ifndef FIXED_POINT
527    speex_free(st->loudness_weight);
528 #endif
529    speex_free(st->echo_noise);
530    speex_free(st->residual_echo);
531
532    speex_free(st->S);
533    speex_free(st->Smin);
534    speex_free(st->Stmp);
535    speex_free(st->update_prob);
536    speex_free(st->zeta);
537
538    speex_free(st->inbuf);
539    speex_free(st->outbuf);
540
541    spx_fft_destroy(st->fft_lookup);
542    filterbank_destroy(st->bank);
543    speex_free(st);
544 }
545
546 /* FIXME: The AGC doesn't work yet with fixed-point*/
547 #ifndef FIXED_POINT
548 static void speex_compute_agc(SpeexPreprocessState *st)
549 {
550    int i;
551    int N = st->ps_size;
552    float scale=.5f/N;
553    float agc_gain;
554    int freq_start, freq_end;
555    float active_bands = 0;
556
557    freq_start = (int)(300.0f*2*N/st->sampling_rate);
558    freq_end   = (int)(2000.0f*2*N/st->sampling_rate);
559    for (i=freq_start;i<freq_end;i++)
560    {
561       if (st->S[i] > 20.f*st->Smin[i]+1000.f)
562          active_bands+=1;
563    }
564    active_bands /= (freq_end-freq_start+1);
565
566    if (active_bands > .2f)
567    {
568       float loudness=0.f;
569       float rate, rate2=.2f;
570       st->nb_loudness_adapt++;
571       rate=2.0f/(1+st->nb_loudness_adapt);
572       if (rate < .05f)
573          rate = .05f;
574       if (rate < .1f && pow(loudness, LOUDNESS_EXP) > st->loudness)
575          rate = .1f;
576       if (rate < .2f && pow(loudness, LOUDNESS_EXP) > 3.f*st->loudness)
577          rate = .2f;
578       if (rate < .4f && pow(loudness, LOUDNESS_EXP) > 10.f*st->loudness)
579          rate = .4f;
580
581       for (i=2;i<N;i++)
582       {
583          loudness += scale*st->ps[i] * FRAC_SCALING_1*FRAC_SCALING_1*st->gain2[i] * st->gain2[i] * st->loudness_weight[i];
584       }
585       loudness=sqrt(loudness);
586       /*if (loudness < 2*pow(st->loudness, 1.0/LOUDNESS_EXP) &&
587         loudness*2 > pow(st->loudness, 1.0/LOUDNESS_EXP))*/
588       st->loudness = (1-rate)*st->loudness + (rate)*pow(loudness, LOUDNESS_EXP);
589       
590       st->loudness2 = (1-rate2)*st->loudness2 + rate2*pow(st->loudness, 1.0f/LOUDNESS_EXP);
591
592       loudness = pow(st->loudness, 1.0f/LOUDNESS_EXP);
593
594       /*fprintf (stderr, "%f %f %f\n", loudness, st->loudness2, rate);*/
595    }
596    
597    agc_gain = st->agc_level/st->loudness2;
598    /*fprintf (stderr, "%f %f %f %f\n", active_bands, st->loudness, st->loudness2, agc_gain);*/
599    if (agc_gain>200)
600       agc_gain = 200;
601
602    for (i=0;i<N;i++)
603       st->gain2[i] *= agc_gain;
604    
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       for (i=0;i<N;i++)
740          st->echo_noise[i] = MAX32(MULT16_32_Q15(QCONST16(.6f,15),st->echo_noise[i]), st->residual_echo[i]);
741       filterbank_compute_bank32(st->bank, st->echo_noise, st->echo_noise+N);
742    } else {
743       for (i=0;i<N+M;i++)
744          st->echo_noise[i] = 0;
745    }
746    preprocess_analysis(st, x);
747
748    update_noise_prob(st);
749
750    /* Noise estimation always updated for the 10 first frames */
751    /*if (st->nb_adapt<10)
752    {
753       for (i=1;i<N-1;i++)
754          st->update_prob[i] = 0;
755    }
756    */
757    
758    /* Update the noise estimate for the frequencies where it can be */
759    for (i=0;i<N;i++)
760    {
761       if (!st->update_prob[i] || st->ps[i] < PSHR32(st->noise[i], NOISE_SHIFT))
762          st->noise[i] = MAX32(EXTEND32(0),MULT16_32_Q15(beta_1,st->noise[i]) + MULT16_32_Q15(beta,SHL32(st->ps[i],NOISE_SHIFT)));
763    }
764    filterbank_compute_bank32(st->bank, st->noise, st->noise+N);
765
766    /* Special case for first frame */
767    if (st->nb_adapt==1)
768       for (i=0;i<N+M;i++)
769          st->old_ps[i] = ps[i];
770
771    /* Compute a posteriori SNR */
772    for (i=0;i<N+M;i++)
773    {
774       spx_word16_t gamma;
775       
776       /* Total noise estimate including residual echo and reverberation */
777       spx_word32_t tot_noise = ADD32(ADD32(ADD32(EXTEND32(1), PSHR32(st->noise[i],NOISE_SHIFT)) , st->echo_noise[i]) , st->reverb_estimate[i]);
778       
779       /* A posteriori SNR = ps/noise - 1*/
780       st->post[i] = SUB16(DIV32_16_Q8(ps[i],tot_noise), QCONST16(1.f,SNR_SHIFT));
781       st->post[i]=MIN16(st->post[i], QCONST16(100.f,SNR_SHIFT));
782       
783       /* Computing update gamma = .1 + .9*(old/(old+noise))^2 */
784       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))));
785       
786       /* A priori SNR update = gamma*max(0,post) + (1-gamma)*old/noise */
787       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));
788       st->prior[i]=MIN16(st->prior[i], QCONST16(100.f,SNR_SHIFT));
789    }
790
791    /*print_vec(st->post, N+M, "");*/
792
793    /* Recursive average of the a priori SNR. A bit smoothed for the psd components */
794    st->zeta[0] = PSHR32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[0]), MULT16_16(QCONST16(.3f,15),st->prior[0])),15);
795    for (i=1;i<N-1;i++)
796       st->zeta[i] = PSHR32(ADD32(ADD32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[i]), MULT16_16(QCONST16(.15f,15),st->prior[i])),
797                            MULT16_16(QCONST16(.075f,15),st->prior[i-1])), MULT16_16(QCONST16(.075f,15),st->prior[i+1])),15);
798    for (i=N-1;i<N+M;i++)
799       st->zeta[i] = PSHR32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[i]), MULT16_16(QCONST16(.3f,15),st->prior[i])),15);
800
801    /* Speech probability of presence for the entire frame is based on the average filterbank a priori SNR */
802    Zframe = 0;
803    for (i=N;i<N+M;i++)
804       Zframe = ADD32(Zframe, EXTEND32(st->zeta[i]));
805    Pframe = QCONST16(.1f,15)+MULT16_16_Q15(QCONST16(.899f,15),qcurve(DIV32_16(Zframe,st->nbands)));
806    
807    effective_echo_suppress = EXTRACT16(PSHR32(ADD32(MULT16_16(SUB16(Q15_ONE,Pframe), st->echo_suppress), MULT16_16(Pframe, st->echo_suppress_active)),15));
808    
809    compute_gain_floor(st->noise_suppress, effective_echo_suppress, st->noise+N, st->echo_noise+N, st->gain_floor+N, M);
810          
811    /* Compute Ephraim & Malah gain speech probability of presence for each critical band (Bark scale) 
812       Technically this is actually wrong because the EM gaim assumes a slightly different probability 
813       distribution */
814    for (i=N;i<N+M;i++)
815    {
816       /* See EM and Cohen papers*/
817       spx_word32_t theta;
818       /* Gain from hypergeometric function */
819       spx_word32_t MM;
820       /* Weiner filter gain */
821       spx_word16_t prior_ratio;
822       /* a priority probability of speech presence based on Bark sub-band alone */
823       spx_word16_t P1;
824       /* Speech absence a priori probability (considering sub-band and frame) */
825       spx_word16_t q;
826 #ifdef FIXED_POINT
827       spx_word16_t tmp;
828 #endif
829       
830       prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1,SNR_SHIFT)));
831       theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT));
832
833       MM = hypergeom_gain(theta);
834       /* Gain with bound */
835       st->gain[i] = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM)));
836       /* Save old Bark power spectrum */
837       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]);
838
839       P1 = QCONST16(.199f,15)+MULT16_16_Q15(QCONST16(.8f,15),qcurve (st->zeta[i]));
840       q = Q15_ONE-MULT16_16_Q15(Pframe,P1);
841 #ifdef FIXED_POINT
842       theta = MIN32(theta, EXTEND32(32767));
843 /*Q8*/tmp = MULT16_16_Q15((SHL32(1,SNR_SHIFT)+st->prior[i]),EXTRACT16(MIN32(Q15ONE,SHR32(spx_exp(-EXTRACT16(theta)),1))));
844       tmp = MIN16(QCONST16(3.,SNR_SHIFT), tmp); /* Prevent overflows in the next line*/
845 /*Q8*/tmp = PSHR(MULT16_16(PDIV32_16(SHL32(EXTEND32(q),8),(Q15_ONE-q)),tmp),8);
846       st->gain2[i]=DIV32_16(SHL(EXTEND32(32767),SNR_SHIFT), ADD16(256,tmp));
847 #else
848       st->gain2[i]=1/(1.f + (q/(1.f-q))*(1+st->prior[i])*exp(-theta));
849 #endif
850    }
851    /* Convert the EM gains and speech prob to linear frequency */
852    filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2);
853    filterbank_compute_psd16(st->bank,st->gain+N, st->gain);
854    
855    /* Use 1 for linear gain resolution (best) or 0 for Bark gain resolution (faster) */
856    if (1)
857    {
858       filterbank_compute_psd16(st->bank,st->gain_floor+N, st->gain_floor);
859    
860       /* Compute gain according to the Ephraim-Malah algorithm -- linear frequency */
861       for (i=0;i<N;i++)
862       {
863          spx_word32_t MM;
864          spx_word32_t theta;
865          spx_word16_t prior_ratio;
866          spx_word16_t tmp;
867          spx_word16_t p;
868          spx_word16_t g;
869          
870          /* Wiener filter gain */
871          prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1,SNR_SHIFT)));
872          theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT));
873
874          /* Optimal estimator for loudness domain */
875          MM = hypergeom_gain(theta);
876          /* EM gain with bound */
877          g = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM)));
878          /* Interpolated speech probability of presence */
879          p = st->gain2[i];
880                   
881          /* Constrain the gain to be close to the Bark scale gain */
882          if (MULT16_16_Q15(QCONST16(.333f,15),g) > st->gain[i])
883             g = MULT16_16(3,st->gain[i]);
884          st->gain[i] = g;
885          
886          /* Save old power spectrum */
887          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]);
888          
889          /* Apply gain floor */
890          if (st->gain[i] < st->gain_floor[i])
891             st->gain[i] = st->gain_floor[i];
892
893          /* Exponential decay model for reverberation (unused) */
894          /*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];*/
895          
896          /* Take into account speech probability of presence (loudness domain MMSE estimator) */
897          /* gain2 = [p*sqrt(gain)+(1-p)*sqrt(gain _floor) ]^2 */
898          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)));
899          st->gain2[i]=SQR16_Q15(tmp);
900
901          /* Use this if you want a log-domain MMSE estimator instead */
902          /*st->gain2[i] = pow(st->gain[i], p) * pow(st->gain_floor[i],1.f-p);*/
903       }
904    } else {
905       for (i=N;i<N+M;i++)
906       {
907          spx_word16_t tmp;
908          spx_word16_t p = st->gain2[i];
909          st->gain[i] = MAX16(st->gain[i], st->gain_floor[i]);         
910          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)));
911          st->gain2[i]=SQR16_Q15(tmp);
912       }
913       filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2);
914    }
915    
916    /* If noise suppression is off, don't apply the gain (but then why call this in the first place!) */
917    if (!st->denoise_enabled)
918    {
919       for (i=0;i<N+M;i++)
920          st->gain2[i]=Q15_ONE;
921    }
922    
923    /*FIXME: This *will* not work for fixed-point */
924 #ifndef FIXED_POINT
925    if (st->agc_enabled)
926       speex_compute_agc(st);
927 #endif
928    
929    /* Apply computed gain */
930    for (i=1;i<N;i++)
931    {
932       st->ft[2*i-1] = MULT16_16_P15(st->gain2[i],st->ft[2*i-1]);
933       st->ft[2*i] = MULT16_16_P15(st->gain2[i],st->ft[2*i]);
934    }
935    st->ft[0] = MULT16_16_P15(st->gain2[0],st->ft[0]);
936    st->ft[2*N-1] = MULT16_16_P15(st->gain2[N-1],st->ft[2*N-1]);
937    
938    /* Inverse FFT with 1/N scaling */
939    spx_ifft(st->fft_lookup, st->ft, st->frame);
940    /* Scale back to original (lower) amplitude */
941    for (i=0;i<2*N;i++)
942       st->frame[i] = PSHR16(st->frame[i], st->frame_shift);
943
944    /*FIXME: This *will* not work for fixed-point */
945 #ifndef FIXED_POINT
946    if (st->agc_enabled)
947    {
948       float max_sample=0;
949       for (i=0;i<2*N;i++)
950          if (fabs(st->frame[i])>max_sample)
951             max_sample = fabs(st->frame[i]);
952       if (max_sample>28000.f)
953       {
954          float damp = 28000.f/max_sample;
955          for (i=0;i<2*N;i++)
956             st->frame[i] *= damp;
957       }
958    }
959 #endif
960    
961    /* Synthesis window (for WOLA) */
962    for (i=0;i<2*N;i++)
963       st->frame[i] = MULT16_16_Q15(st->frame[i], st->window[i]);
964
965    /* Perform overlap and add */
966    for (i=0;i<N3;i++)
967       x[i] = st->outbuf[i] + st->frame[i];
968    for (i=0;i<N4;i++)
969       x[N3+i] = st->frame[N3+i];
970    
971    /* Update outbuf */
972    for (i=0;i<N3;i++)
973       st->outbuf[i] = st->frame[st->frame_size+i];
974
975    /* FIXME: This VAD is a kludge */
976    if (st->vad_enabled)
977    {
978       if (Pframe > st->speech_prob_start || (st->was_speech && Pframe > st->speech_prob_continue))
979       {
980          st->was_speech=1;
981          return 1;
982       } else
983       {
984          st->was_speech=0;
985          return 0;
986       }
987    } else {
988       return 1;
989    }
990 }
991
992 void speex_preprocess_estimate_update(SpeexPreprocessState *st, spx_int16_t *x, spx_int32_t *echo)
993 {
994    int i;
995    int N = st->ps_size;
996    int N3 = 2*N - st->frame_size;
997    int M;
998    spx_word32_t *ps=st->ps;
999
1000    M = st->nbands;
1001    st->min_count++;
1002    
1003    preprocess_analysis(st, x);
1004
1005    update_noise_prob(st);
1006    
1007    for (i=1;i<N-1;i++)
1008    {
1009       if (!st->update_prob[i] || st->ps[i] < PSHR32(st->noise[i],NOISE_SHIFT))
1010       {
1011          st->noise[i] = MULT16_32_Q15(QCONST16(.95f,15),st->noise[i]) + MULT16_32_Q15(QCONST16(.05f,15),SHL32(st->ps[i],NOISE_SHIFT));
1012       }
1013    }
1014
1015    for (i=0;i<N3;i++)
1016       st->outbuf[i] = MULT16_16_Q15(x[st->frame_size-N3+i],st->window[st->frame_size+i]);
1017
1018    /* Save old power spectrum */
1019    for (i=0;i<N+M;i++)
1020       st->old_ps[i] = ps[i];
1021
1022    for (i=0;i<N;i++)
1023       st->reverb_estimate[i] = MULT16_32_Q15(st->reverb_decay, st->reverb_estimate[i]);
1024 }
1025
1026
1027 int speex_preprocess_ctl(SpeexPreprocessState *state, int request, void *ptr)
1028 {
1029    int i;
1030    SpeexPreprocessState *st;
1031    st=(SpeexPreprocessState*)state;
1032    switch(request)
1033    {
1034    case SPEEX_PREPROCESS_SET_DENOISE:
1035       st->denoise_enabled = (*(spx_int32_t*)ptr);
1036       break;
1037    case SPEEX_PREPROCESS_GET_DENOISE:
1038       (*(spx_int32_t*)ptr) = st->denoise_enabled;
1039       break;
1040 #ifndef FIXED_POINT
1041    case SPEEX_PREPROCESS_SET_AGC:
1042       st->agc_enabled = (*(spx_int32_t*)ptr);
1043       break;
1044    case SPEEX_PREPROCESS_GET_AGC:
1045       (*(spx_int32_t*)ptr) = st->agc_enabled;
1046       break;
1047
1048    case SPEEX_PREPROCESS_SET_AGC_LEVEL:
1049       st->agc_level = (*(float*)ptr);
1050       if (st->agc_level<1)
1051          st->agc_level=1;
1052       if (st->agc_level>32768)
1053          st->agc_level=32768;
1054       break;
1055    case SPEEX_PREPROCESS_GET_AGC_LEVEL:
1056       (*(float*)ptr) = st->agc_level;
1057       break;
1058 #endif
1059    case SPEEX_PREPROCESS_SET_VAD:
1060       speex_warning("The VAD has been replaced by a hack pending a complete rewrite");
1061       st->vad_enabled = (*(spx_int32_t*)ptr);
1062       break;
1063    case SPEEX_PREPROCESS_GET_VAD:
1064       (*(spx_int32_t*)ptr) = st->vad_enabled;
1065       break;
1066    
1067    case SPEEX_PREPROCESS_SET_DEREVERB:
1068       st->dereverb_enabled = (*(spx_int32_t*)ptr);
1069       for (i=0;i<st->ps_size;i++)
1070          st->reverb_estimate[i]=0;
1071       break;
1072    case SPEEX_PREPROCESS_GET_DEREVERB:
1073       (*(spx_int32_t*)ptr) = st->dereverb_enabled;
1074       break;
1075
1076    case SPEEX_PREPROCESS_SET_DEREVERB_LEVEL:
1077       st->reverb_level = (*(float*)ptr);
1078       break;
1079    case SPEEX_PREPROCESS_GET_DEREVERB_LEVEL:
1080       (*(float*)ptr) = st->reverb_level;
1081       break;
1082    
1083    case SPEEX_PREPROCESS_SET_DEREVERB_DECAY:
1084       st->reverb_decay = (*(float*)ptr);
1085       break;
1086    case SPEEX_PREPROCESS_GET_DEREVERB_DECAY:
1087       (*(float*)ptr) = st->reverb_decay;
1088       break;
1089
1090    case SPEEX_PREPROCESS_SET_PROB_START:
1091       *(spx_int32_t*)ptr = MIN32(Q15_ONE,MAX32(0, *(spx_int32_t*)ptr));
1092       st->speech_prob_start = DIV32_16(MULT16_16(32767,*(spx_int32_t*)ptr), 100);
1093       break;
1094    case SPEEX_PREPROCESS_GET_PROB_START:
1095       (*(spx_int32_t*)ptr) = MULT16_16_Q15(st->speech_prob_start, 100);
1096       break;
1097
1098    case SPEEX_PREPROCESS_SET_PROB_CONTINUE:
1099       *(spx_int32_t*)ptr = MIN32(Q15_ONE,MAX32(0, *(spx_int32_t*)ptr));
1100       st->speech_prob_continue = DIV32_16(MULT16_16(32767,*(spx_int32_t*)ptr), 100);
1101       break;
1102    case SPEEX_PREPROCESS_GET_PROB_CONTINUE:
1103       (*(spx_int32_t*)ptr) = MULT16_16_Q15(st->speech_prob_continue, 100);
1104       break;
1105
1106    case SPEEX_PREPROCESS_SET_NOISE_SUPPRESS:
1107       st->noise_suppress = -ABS(*(spx_int32_t*)ptr);
1108       break;
1109    case SPEEX_PREPROCESS_GET_NOISE_SUPPRESS:
1110       (*(spx_int32_t*)ptr) = st->noise_suppress;
1111       break;
1112    case SPEEX_PREPROCESS_SET_ECHO_SUPPRESS:
1113       st->echo_suppress = -ABS(*(spx_int32_t*)ptr);
1114       break;
1115    case SPEEX_PREPROCESS_GET_ECHO_SUPPRESS:
1116       (*(spx_int32_t*)ptr) = st->echo_suppress;
1117       break;
1118    case SPEEX_PREPROCESS_SET_ECHO_SUPPRESS_ACTIVE:
1119       st->echo_suppress_active = -ABS(*(spx_int32_t*)ptr);
1120       break;
1121    case SPEEX_PREPROCESS_GET_ECHO_SUPPRESS_ACTIVE:
1122       (*(spx_int32_t*)ptr) = st->echo_suppress_active;
1123       break;
1124    case SPEEX_PREPROCESS_SET_ECHO_STATE:
1125       st->echo_state = (SpeexEchoState*)ptr;
1126       break;
1127    case SPEEX_PREPROCESS_GET_ECHO_STATE:
1128       ptr = (void*)st->echo_state;
1129       break;
1130
1131    default:
1132       speex_warning_int("Unknown speex_preprocess_ctl request: ", request);
1133       return -1;
1134    }
1135    return 0;
1136 }