Added a higher precision cosine approximation and used it for the
[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    /*FIXME: Set ps[0] properly */
646    ps[0]=MULT16_16(st->ft[0],st->ft[0]);
647    for (i=1;i<N;i++)
648       ps[i]=MULT16_16(st->ft[2*i-1],st->ft[2*i-1]) + MULT16_16(st->ft[2*i],st->ft[2*i]);
649    for (i=0;i<N;i++)
650       st->ps[i] = PSHR32(st->ps[i], 2*st->frame_shift);
651
652    filterbank_compute_bank32(st->bank, ps, ps+N);
653 }
654
655 static void update_noise_prob(SpeexPreprocessState *st)
656 {
657    int i;
658    int min_range;
659    int N = st->ps_size;
660
661    for (i=1;i<N-1;i++)
662       st->S[i] =  MULT16_32_Q15(QCONST16(.8f,15),st->S[i]) + MULT16_32_Q15(QCONST16(.05f,15),st->ps[i-1]) 
663                       + MULT16_32_Q15(QCONST16(.1f,15),st->ps[i]) + MULT16_32_Q15(QCONST16(.05f,15),st->ps[i+1]);
664    st->S[0] =  MULT16_32_Q15(QCONST16(.8f,15),st->S[0]) + MULT16_32_Q15(QCONST16(.2f,15),st->ps[0]);
665    st->S[N-1] =  MULT16_32_Q15(QCONST16(.8f,15),st->S[N-1]) + MULT16_32_Q15(QCONST16(.2f,15),st->ps[N-1]);
666    
667    if (st->nb_adapt==1)
668    {
669       for (i=0;i<N;i++)
670          st->Smin[i] = st->Stmp[i] = 0;
671    }
672
673    if (st->nb_adapt < 100)
674       min_range = 15;
675    else if (st->nb_adapt < 1000)
676       min_range = 50;
677    else if (st->nb_adapt < 10000)
678       min_range = 150;
679    else
680       min_range = 300;
681    if (st->min_count > min_range)
682    {
683       st->min_count = 0;
684       for (i=0;i<N;i++)
685       {
686          st->Smin[i] = MIN32(st->Stmp[i], st->S[i]);
687          st->Stmp[i] = st->S[i];
688       }
689    } else {
690       for (i=0;i<N;i++)
691       {
692          st->Smin[i] = MIN32(st->Smin[i], st->S[i]);
693          st->Stmp[i] = MIN32(st->Stmp[i], st->S[i]);      
694       }
695    }
696    for (i=0;i<N;i++)
697    {
698       if (MULT16_32_Q15(QCONST16(.4f,15),st->S[i]) > ADD32(st->Smin[i],EXTEND32(20)))
699          st->update_prob[i] = 1;
700       else
701          st->update_prob[i] = 0;
702       /*fprintf (stderr, "%f ", st->S[i]/st->Smin[i]);*/
703       /*fprintf (stderr, "%f ", st->update_prob[i]);*/
704    }
705
706 }
707
708 #define NOISE_OVERCOMPENS 1.
709
710 void speex_echo_get_residual(SpeexEchoState *st, spx_word32_t *Yout, int len);
711
712 int speex_preprocess(SpeexPreprocessState *st, spx_int16_t *x, spx_int32_t *echo)
713 {
714    return speex_preprocess_run(st, x);
715 }
716
717 int speex_preprocess_run(SpeexPreprocessState *st, spx_int16_t *x)
718 {
719    int i;
720    int M;
721    int N = st->ps_size;
722    int N3 = 2*N - st->frame_size;
723    int N4 = st->frame_size - N3;
724    spx_word32_t *ps=st->ps;
725    spx_word32_t Zframe;
726    spx_word16_t Pframe;
727    spx_word16_t beta, beta_1;
728    spx_word16_t effective_echo_suppress;
729    
730    st->nb_adapt++;
731    st->min_count++;
732    
733    beta = MAX16(QCONST16(.03,15),DIV32_16(Q15_ONE,st->nb_adapt));
734    beta_1 = Q15_ONE-beta;
735    M = st->nbands;
736    /* Deal with residual echo if provided */
737    if (st->echo_state)
738    {
739       speex_echo_get_residual(st->echo_state, st->residual_echo, N);
740       for (i=0;i<N;i++)
741          st->echo_noise[i] = MAX32(MULT16_32_Q15(QCONST16(.6f,15),st->echo_noise[i]), st->residual_echo[i]);
742       filterbank_compute_bank32(st->bank, st->echo_noise, st->echo_noise+N);
743    } else {
744       for (i=0;i<N+M;i++)
745          st->echo_noise[i] = 0;
746    }
747    preprocess_analysis(st, x);
748
749    update_noise_prob(st);
750
751    /* Noise estimation always updated for the 10 first frames */
752    /*if (st->nb_adapt<10)
753    {
754       for (i=1;i<N-1;i++)
755          st->update_prob[i] = 0;
756    }
757    */
758    
759    /* Update the noise estimate for the frequencies where it can be */
760    for (i=0;i<N;i++)
761    {
762       if (!st->update_prob[i] || st->ps[i] < PSHR32(st->noise[i], NOISE_SHIFT))
763          st->noise[i] = MULT16_32_Q15(beta_1,st->noise[i]) + MULT16_32_Q15(beta,SHL32(st->ps[i],NOISE_SHIFT));
764    }
765    filterbank_compute_bank32(st->bank, st->noise, st->noise+N);
766
767    /* Special case for first frame */
768    if (st->nb_adapt==1)
769       for (i=0;i<N+M;i++)
770          st->old_ps[i] = ps[i];
771
772    /* Compute a posteriori SNR */
773    for (i=0;i<N+M;i++)
774    {
775       spx_word16_t gamma;
776       
777       /* Total noise estimate including residual echo and reverberation */
778       spx_word32_t tot_noise = ADD32(ADD32(ADD32(EXTEND32(1), PSHR32(st->noise[i],NOISE_SHIFT)) , st->echo_noise[i]) , st->reverb_estimate[i]);
779       
780       /* A posteriori SNR = ps/noise - 1*/
781       st->post[i] = SUB16(DIV32_16_Q8(ps[i],tot_noise), QCONST16(1.f,SNR_SHIFT));
782       st->post[i]=MIN16(st->post[i], QCONST16(100.f,SNR_SHIFT));
783       
784       /* Computing update gamma = .1 + .9*(old/(old+noise))^2 */
785       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))));
786       
787       /* A priori SNR update = gamma*max(0,post) + (1-gamma)*old/noise */
788       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));
789       st->prior[i]=MIN16(st->prior[i], QCONST16(100.f,SNR_SHIFT));
790    }
791
792    /*print_vec(st->post, N+M, "");*/
793
794    /* Recursive average of the a priori SNR. A bit smoothed for the psd components */
795    st->zeta[0] = PSHR32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[0]), MULT16_16(QCONST16(.3f,15),st->prior[0])),15);
796    for (i=1;i<N-1;i++)
797       st->zeta[i] = PSHR32(ADD32(ADD32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[i]), MULT16_16(QCONST16(.15f,15),st->prior[i])),
798                            MULT16_16(QCONST16(.075f,15),st->prior[i-1])), MULT16_16(QCONST16(.075f,15),st->prior[i+1])),15);
799    for (i=N-1;i<N+M;i++)
800       st->zeta[i] = PSHR32(ADD32(MULT16_16(QCONST16(.7f,15),st->zeta[i]), MULT16_16(QCONST16(.3f,15),st->prior[i])),15);
801
802    /* Speech probability of presence for the entire frame is based on the average filterbank a priori SNR */
803    Zframe = 0;
804    for (i=N;i<N+M;i++)
805       Zframe = ADD32(Zframe, EXTEND32(st->zeta[i]));
806    Pframe = QCONST16(.1f,15)+MULT16_16_Q15(QCONST16(.899f,15),qcurve(DIV32_16(Zframe,st->nbands)));
807    
808    effective_echo_suppress = EXTRACT16(PSHR32(ADD32(MULT16_16(SUB16(Q15_ONE,Pframe), st->echo_suppress), MULT16_16(Pframe, st->echo_suppress_active)),15));
809    
810    compute_gain_floor(st->noise_suppress, effective_echo_suppress, st->noise+N, st->echo_noise+N, st->gain_floor+N, M);
811          
812    /* Compute Ephraim & Malah gain speech probability of presence for each critical band (Bark scale) 
813       Technically this is actually wrong because the EM gaim assumes a slightly different probability 
814       distribution */
815    for (i=N;i<N+M;i++)
816    {
817       /* See EM and Cohen papers*/
818       spx_word32_t theta;
819       /* Gain from hypergeometric function */
820       spx_word32_t MM;
821       /* Weiner filter gain */
822       spx_word16_t prior_ratio;
823       /* a priority probability of speech presence based on Bark sub-band alone */
824       spx_word16_t P1;
825       /* Speech absence a priori probability (considering sub-band and frame) */
826       spx_word16_t q;
827 #ifdef FIXED_POINT
828       spx_word16_t tmp;
829 #endif
830       
831       prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1,SNR_SHIFT)));
832       theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT));
833
834       MM = hypergeom_gain(theta);
835       /* Gain with bound */
836       st->gain[i] = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM)));
837       /* Save old Bark power spectrum */
838       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]);
839
840       P1 = QCONST16(.199f,15)+MULT16_16_Q15(QCONST16(.8f,15),qcurve (st->zeta[i]));
841       q = Q15_ONE-MULT16_16_Q15(Pframe,P1);
842 #ifdef FIXED_POINT
843       theta = MIN32(theta, EXTEND32(32767));
844 /*Q8*/tmp = MULT16_16_Q15((SHL32(1,SNR_SHIFT)+st->prior[i]),EXTRACT16(MIN32(Q15ONE,SHR32(spx_exp(-EXTRACT16(theta)),1))));
845       tmp = MIN16(QCONST16(3.,SNR_SHIFT), tmp); /* Prevent overflows in the next line*/
846 /*Q8*/tmp = PSHR(MULT16_16(PDIV32_16(SHL32(EXTEND32(q),8),(Q15_ONE-q)),tmp),8);
847       st->gain2[i]=DIV32_16(SHL(EXTEND32(32767),SNR_SHIFT), ADD16(256,tmp));
848 #else
849       st->gain2[i]=1/(1.f + (q/(1.f-q))*(1+st->prior[i])*exp(-theta));
850 #endif
851    }
852    /* Convert the EM gains and speech prob to linear frequency */
853    filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2);
854    filterbank_compute_psd16(st->bank,st->gain+N, st->gain);
855    
856    /* Use 1 for linear gain resolution (best) or 0 for Bark gain resolution (faster) */
857    if (1)
858    {
859       filterbank_compute_psd16(st->bank,st->gain_floor+N, st->gain_floor);
860    
861       /* Compute gain according to the Ephraim-Malah algorithm -- linear frequency */
862       for (i=0;i<N;i++)
863       {
864          spx_word32_t MM;
865          spx_word32_t theta;
866          spx_word16_t prior_ratio;
867          spx_word16_t tmp;
868          spx_word16_t p;
869          spx_word16_t g;
870          
871          /* Wiener filter gain */
872          prior_ratio = PDIV32_16(SHL32(EXTEND32(st->prior[i]), 15), ADD16(st->prior[i], SHL32(1,SNR_SHIFT)));
873          theta = MULT16_32_P15(prior_ratio, QCONST32(1.f,EXPIN_SHIFT)+SHL32(EXTEND32(st->post[i]),EXPIN_SHIFT-SNR_SHIFT));
874
875          /* Optimal estimator for loudness domain */
876          MM = hypergeom_gain(theta);
877          /* EM gain with bound */
878          g = EXTRACT16(MIN32(Q15_ONE, MULT16_32_Q15(prior_ratio, MM)));
879          /* Interpolated speech probability of presence */
880          p = st->gain2[i];
881                   
882          /* Constrain the gain to be close to the Bark scale gain */
883          if (MULT16_16_Q15(QCONST16(.333f,15),g) > st->gain[i])
884             g = MULT16_16(3,st->gain[i]);
885          st->gain[i] = g;
886          
887          /* Save old power spectrum */
888          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]);
889          
890          /* Apply gain floor */
891          if (st->gain[i] < st->gain_floor[i])
892             st->gain[i] = st->gain_floor[i];
893
894          /* Exponential decay model for reverberation (unused) */
895          /*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];*/
896          
897          /* Take into account speech probability of presence (loudness domain MMSE estimator) */
898          /* gain2 = [p*sqrt(gain)+(1-p)*sqrt(gain _floor) ]^2 */
899          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)));
900          st->gain2[i]=SQR16_Q15(tmp);
901
902          /* Use this if you want a log-domain MMSE estimator instead */
903          /*st->gain2[i] = pow(st->gain[i], p) * pow(st->gain_floor[i],1.f-p);*/
904       }
905    } else {
906       for (i=N;i<N+M;i++)
907       {
908          spx_word16_t tmp;
909          spx_word16_t p = st->gain2[i];
910          st->gain[i] = MAX16(st->gain[i], st->gain_floor[i]);         
911          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)));
912          st->gain2[i]=SQR16_Q15(tmp);
913       }
914       filterbank_compute_psd16(st->bank,st->gain2+N, st->gain2);
915    }
916    
917    /* If noise suppression is off, don't apply the gain (but then why call this in the first place!) */
918    if (!st->denoise_enabled)
919    {
920       for (i=0;i<N+M;i++)
921          st->gain2[i]=Q15_ONE;
922    }
923    
924    /*FIXME: This *will* not work for fixed-point */
925 #ifndef FIXED_POINT
926    if (st->agc_enabled)
927       speex_compute_agc(st);
928 #endif
929    
930    /* Apply computed gain */
931    for (i=1;i<N;i++)
932    {
933       st->ft[2*i-1] = MULT16_16_P15(st->gain2[i],st->ft[2*i-1]);
934       st->ft[2*i] = MULT16_16_P15(st->gain2[i],st->ft[2*i]);
935    }
936    st->ft[0] = MULT16_16_P15(st->gain2[0],st->ft[0]);
937    st->ft[2*N-1] = MULT16_16_P15(st->gain2[N-1],st->ft[2*N-1]);
938    
939    /* Inverse FFT with 1/N scaling */
940    spx_ifft(st->fft_lookup, st->ft, st->frame);
941    /* Scale back to original (lower) amplitude */
942    for (i=0;i<2*N;i++)
943       st->frame[i] = PSHR16(st->frame[i], st->frame_shift);
944
945    /*FIXME: This *will* not work for fixed-point */
946 #ifndef FIXED_POINT
947    if (st->agc_enabled)
948    {
949       float max_sample=0;
950       for (i=0;i<2*N;i++)
951          if (fabs(st->frame[i])>max_sample)
952             max_sample = fabs(st->frame[i]);
953       if (max_sample>28000.f)
954       {
955          float damp = 28000.f/max_sample;
956          for (i=0;i<2*N;i++)
957             st->frame[i] *= damp;
958       }
959    }
960 #endif
961    
962    /* Synthesis window (for WOLA) */
963    for (i=0;i<2*N;i++)
964       st->frame[i] = MULT16_16_Q15(st->frame[i], st->window[i]);
965
966    /* Perform overlap and add */
967    for (i=0;i<N3;i++)
968       x[i] = st->outbuf[i] + st->frame[i];
969    for (i=0;i<N4;i++)
970       x[N3+i] = st->frame[N3+i];
971    
972    /* Update outbuf */
973    for (i=0;i<N3;i++)
974       st->outbuf[i] = st->frame[st->frame_size+i];
975
976    /* FIXME: This VAD is a kludge */
977    if (st->vad_enabled)
978    {
979       if (Pframe > st->speech_prob_start || (st->was_speech && Pframe > st->speech_prob_continue))
980       {
981          st->was_speech=1;
982          return 1;
983       } else
984       {
985          st->was_speech=0;
986          return 0;
987       }
988    } else {
989       return 1;
990    }
991 }
992
993 void speex_preprocess_estimate_update(SpeexPreprocessState *st, spx_int16_t *x, spx_int32_t *echo)
994 {
995    int i;
996    int N = st->ps_size;
997    int N3 = 2*N - st->frame_size;
998    int M;
999    spx_word32_t *ps=st->ps;
1000
1001    M = st->nbands;
1002    st->min_count++;
1003    
1004    preprocess_analysis(st, x);
1005
1006    update_noise_prob(st);
1007    
1008    for (i=1;i<N-1;i++)
1009    {
1010       if (!st->update_prob[i] || st->ps[i] < PSHR32(st->noise[i],NOISE_SHIFT))
1011       {
1012          st->noise[i] = MULT16_32_Q15(QCONST16(.95f,15),st->noise[i]) + MULT16_32_Q15(QCONST16(.05f,15),SHL32(st->ps[i],NOISE_SHIFT));
1013       }
1014    }
1015
1016    for (i=0;i<N3;i++)
1017       st->outbuf[i] = MULT16_16_Q15(x[st->frame_size-N3+i],st->window[st->frame_size+i]);
1018
1019    /* Save old power spectrum */
1020    for (i=0;i<N+M;i++)
1021       st->old_ps[i] = ps[i];
1022
1023    for (i=0;i<N;i++)
1024       st->reverb_estimate[i] = MULT16_32_Q15(st->reverb_decay, st->reverb_estimate[i]);
1025 }
1026
1027
1028 int speex_preprocess_ctl(SpeexPreprocessState *state, int request, void *ptr)
1029 {
1030    int i;
1031    SpeexPreprocessState *st;
1032    st=(SpeexPreprocessState*)state;
1033    switch(request)
1034    {
1035    case SPEEX_PREPROCESS_SET_DENOISE:
1036       st->denoise_enabled = (*(spx_int32_t*)ptr);
1037       break;
1038    case SPEEX_PREPROCESS_GET_DENOISE:
1039       (*(spx_int32_t*)ptr) = st->denoise_enabled;
1040       break;
1041 #ifndef FIXED_POINT
1042    case SPEEX_PREPROCESS_SET_AGC:
1043       st->agc_enabled = (*(spx_int32_t*)ptr);
1044       break;
1045    case SPEEX_PREPROCESS_GET_AGC:
1046       (*(spx_int32_t*)ptr) = st->agc_enabled;
1047       break;
1048
1049    case SPEEX_PREPROCESS_SET_AGC_LEVEL:
1050       st->agc_level = (*(float*)ptr);
1051       if (st->agc_level<1)
1052          st->agc_level=1;
1053       if (st->agc_level>32768)
1054          st->agc_level=32768;
1055       break;
1056    case SPEEX_PREPROCESS_GET_AGC_LEVEL:
1057       (*(float*)ptr) = st->agc_level;
1058       break;
1059 #endif
1060    case SPEEX_PREPROCESS_SET_VAD:
1061       speex_warning("The VAD has been replaced by a hack pending a complete rewrite");
1062       st->vad_enabled = (*(spx_int32_t*)ptr);
1063       break;
1064    case SPEEX_PREPROCESS_GET_VAD:
1065       (*(spx_int32_t*)ptr) = st->vad_enabled;
1066       break;
1067    
1068    case SPEEX_PREPROCESS_SET_DEREVERB:
1069       st->dereverb_enabled = (*(spx_int32_t*)ptr);
1070       for (i=0;i<st->ps_size;i++)
1071          st->reverb_estimate[i]=0;
1072       break;
1073    case SPEEX_PREPROCESS_GET_DEREVERB:
1074       (*(spx_int32_t*)ptr) = st->dereverb_enabled;
1075       break;
1076
1077    case SPEEX_PREPROCESS_SET_DEREVERB_LEVEL:
1078       st->reverb_level = (*(float*)ptr);
1079       break;
1080    case SPEEX_PREPROCESS_GET_DEREVERB_LEVEL:
1081       (*(float*)ptr) = st->reverb_level;
1082       break;
1083    
1084    case SPEEX_PREPROCESS_SET_DEREVERB_DECAY:
1085       st->reverb_decay = (*(float*)ptr);
1086       break;
1087    case SPEEX_PREPROCESS_GET_DEREVERB_DECAY:
1088       (*(float*)ptr) = st->reverb_decay;
1089       break;
1090
1091    case SPEEX_PREPROCESS_SET_PROB_START:
1092       *(spx_int32_t*)ptr = MIN32(Q15_ONE,MAX32(0, *(spx_int32_t*)ptr));
1093       st->speech_prob_start = DIV32_16(MULT16_16(32767,*(spx_int32_t*)ptr), 100);
1094       break;
1095    case SPEEX_PREPROCESS_GET_PROB_START:
1096       (*(spx_int32_t*)ptr) = MULT16_16_Q15(st->speech_prob_start, 100);
1097       break;
1098
1099    case SPEEX_PREPROCESS_SET_PROB_CONTINUE:
1100       *(spx_int32_t*)ptr = MIN32(Q15_ONE,MAX32(0, *(spx_int32_t*)ptr));
1101       st->speech_prob_continue = DIV32_16(MULT16_16(32767,*(spx_int32_t*)ptr), 100);
1102       break;
1103    case SPEEX_PREPROCESS_GET_PROB_CONTINUE:
1104       (*(spx_int32_t*)ptr) = MULT16_16_Q15(st->speech_prob_continue, 100);
1105       break;
1106
1107    case SPEEX_PREPROCESS_SET_NOISE_SUPPRESS:
1108       st->noise_suppress = -ABS(*(spx_int32_t*)ptr);
1109       break;
1110    case SPEEX_PREPROCESS_GET_NOISE_SUPPRESS:
1111       (*(spx_int32_t*)ptr) = st->noise_suppress;
1112       break;
1113    case SPEEX_PREPROCESS_SET_ECHO_SUPPRESS:
1114       st->echo_suppress = -ABS(*(spx_int32_t*)ptr);
1115       break;
1116    case SPEEX_PREPROCESS_GET_ECHO_SUPPRESS:
1117       (*(spx_int32_t*)ptr) = st->echo_suppress;
1118       break;
1119    case SPEEX_PREPROCESS_SET_ECHO_SUPPRESS_ACTIVE:
1120       st->echo_suppress_active = -ABS(*(spx_int32_t*)ptr);
1121       break;
1122    case SPEEX_PREPROCESS_GET_ECHO_SUPPRESS_ACTIVE:
1123       (*(spx_int32_t*)ptr) = st->echo_suppress_active;
1124       break;
1125    case SPEEX_PREPROCESS_SET_ECHO_STATE:
1126       st->echo_state = (SpeexEchoState*)ptr;
1127       break;
1128    case SPEEX_PREPROCESS_GET_ECHO_STATE:
1129       ptr = (void*)st->echo_state;
1130       break;
1131
1132    default:
1133       speex_warning_int("Unknown speex_preprocess_ctl request: ", request);
1134       return -1;
1135    }
1136    return 0;
1137 }