trying some ideas for soft-decision DTD based on residual-to-signal ratio
[speexdsp.git] / libspeex / filters.c
1 /* Copyright (C) 2002 Jean-Marc Valin 
2    File: filters.c
3    Various analysis/synthesis filters
4
5    Redistribution and use in source and binary forms, with or without
6    modification, are permitted provided that the following conditions
7    are met:
8    
9    - Redistributions of source code must retain the above copyright
10    notice, this list of conditions and the following disclaimer.
11    
12    - Redistributions in binary form must reproduce the above copyright
13    notice, this list of conditions and the following disclaimer in the
14    documentation and/or other materials provided with the distribution.
15    
16    - Neither the name of the Xiph.org Foundation nor the names of its
17    contributors may be used to endorse or promote products derived from
18    this software without specific prior written permission.
19    
20    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
24    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
25    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
26    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 #ifdef HAVE_CONFIG_H
34 #include "config.h"
35 #endif
36
37 #include "filters.h"
38 #include "stack_alloc.h"
39 #include "misc.h"
40 #include "math_approx.h"
41 #include "ltp.h"
42 #include <math.h>
43
44 void bw_lpc(spx_word16_t gamma, const spx_coef_t *lpc_in, spx_coef_t *lpc_out, int order)
45 {
46    int i;
47    spx_word16_t tmp=gamma;
48    lpc_out[0] = lpc_in[0];
49    for (i=1;i<order+1;i++)
50    {
51       lpc_out[i] = MULT16_16_P15(tmp,lpc_in[i]);
52       tmp = MULT16_16_P15(tmp, gamma);
53    }
54 }
55
56
57 #ifdef FIXED_POINT
58
59 /* FIXME: These functions are ugly and probably introduce too much error */
60 void signal_mul(const spx_sig_t *x, spx_sig_t *y, spx_word32_t scale, int len)
61 {
62    int i;
63    for (i=0;i<len;i++)
64    {
65       y[i] = SHL32(MULT16_32_Q14(EXTRACT16(SHR32(x[i],7)),scale),7);
66    }
67 }
68
69 void signal_div(const spx_sig_t *x, spx_sig_t *y, spx_word32_t scale, int len)
70 {
71    int i;
72    spx_word16_t scale_1;
73
74    scale = PSHR32(scale, SIG_SHIFT);
75    if (scale<2)
76       scale_1 = 32767;
77    else
78       scale_1 = EXTRACT16(DIV32(32767,scale));
79    for (i=0;i<len;i++)
80    {
81       y[i] = MULT16_32_Q15(scale_1,x[i]);
82    }
83 }
84
85 #else
86
87 void signal_mul(const spx_sig_t *x, spx_sig_t *y, spx_word32_t scale, int len)
88 {
89    int i;
90    for (i=0;i<len;i++)
91       y[i] = scale*x[i];
92 }
93
94 void signal_div(const spx_sig_t *x, spx_sig_t *y, spx_word32_t scale, int len)
95 {
96    int i;
97    float scale_1 = 1/scale;
98    for (i=0;i<len;i++)
99       y[i] = scale_1*x[i];
100 }
101 #endif
102
103
104
105 #ifdef FIXED_POINT
106
107
108
109 spx_word16_t compute_rms(const spx_sig_t *x, int len)
110 {
111    int i;
112    spx_word32_t sum=0;
113    spx_sig_t max_val=1;
114    int sig_shift;
115
116    for (i=0;i<len;i++)
117    {
118       spx_sig_t tmp = x[i];
119       if (tmp<0)
120          tmp = -tmp;
121       if (tmp > max_val)
122          max_val = tmp;
123    }
124
125    sig_shift=0;
126    while (max_val>16383)
127    {
128       sig_shift++;
129       max_val >>= 1;
130    }
131
132    for (i=0;i<len;i+=4)
133    {
134       spx_word32_t sum2=0;
135       spx_word16_t tmp;
136       tmp = EXTRACT16(SHR32(x[i],sig_shift));
137       sum2 = MAC16_16(sum2,tmp,tmp);
138       tmp = EXTRACT16(SHR32(x[i+1],sig_shift));
139       sum2 = MAC16_16(sum2,tmp,tmp);
140       tmp = EXTRACT16(SHR32(x[i+2],sig_shift));
141       sum2 = MAC16_16(sum2,tmp,tmp);
142       tmp = EXTRACT16(SHR32(x[i+3],sig_shift));
143       sum2 = MAC16_16(sum2,tmp,tmp);
144       sum = ADD32(sum,SHR32(sum2,6));
145    }
146    
147    return EXTRACT16(SHR32(SHL32(EXTEND32(spx_sqrt(1+DIV32(sum,len))),(sig_shift+3)),SIG_SHIFT));
148 }
149
150 #if defined(ARM4_ASM) || defined(ARM5E_ASM)
151 #include "filters_arm4.h"
152 #else
153
154
155 int normalize16(const spx_sig_t *x, spx_word16_t *y, spx_sig_t max_scale, int len)
156 {
157    int i;
158    spx_sig_t max_val=1;
159    int sig_shift;
160    
161    for (i=0;i<len;i++)
162    {
163       spx_sig_t tmp = x[i];
164       if (tmp<0)
165          tmp = NEG32(tmp);
166       if (tmp >= max_val)
167          max_val = tmp;
168    }
169
170    sig_shift=0;
171    while (max_val>max_scale)
172    {
173       sig_shift++;
174       max_val >>= 1;
175    }
176
177    for (i=0;i<len;i++)
178       y[i] = EXTRACT16(SHR32(x[i], sig_shift));
179    
180    return sig_shift;
181 }
182
183 #ifdef PRECISION16
184 void filter_mem2(const spx_sig_t *x, const spx_coef_t *num, const spx_coef_t *den, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
185 {
186    int i,j;
187    spx_word16_t xi,yi,nyi;
188
189    for (i=0;i<N;i++)
190    {
191       xi= EXTRACT16(PSHR32(SATURATE(x[i],536870911),SIG_SHIFT));
192       yi = EXTRACT16(PSHR32(SATURATE(ADD32(x[i], SHL32(mem[0],1)),536870911),SIG_SHIFT));
193       nyi = NEG16(yi);
194       for (j=0;j<ord-1;j++)
195       {
196          mem[j] = MAC16_16(MAC16_16(mem[j+1], num[j+1],xi), den[j+1],nyi);
197       }
198       mem[ord-1] = ADD32(MULT16_16(num[ord],xi), MULT16_16(den[ord],nyi));
199       y[i] = SHL32(EXTEND32(yi),SIG_SHIFT);
200    }
201 }
202 #else
203 void filter_mem2(const spx_sig_t *x, const spx_coef_t *num, const spx_coef_t *den, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
204 {
205    int i,j;
206    spx_sig_t xi,yi,nyi;
207
208    for (i=0;i<N;i++)
209    {
210       xi=SATURATE(x[i],805306368);
211       yi = SATURATE(ADD32(xi, SHL32(mem[0],2)),805306368);
212       nyi = NEG32(yi);
213       for (j=0;j<ord-1;j++)
214       {
215          mem[j] = MAC16_32_Q15(MAC16_32_Q15(mem[j+1], num[j+1],xi), den[j+1],nyi);
216       }
217       mem[ord-1] = SUB32(MULT16_32_Q15(num[ord],xi), MULT16_32_Q15(den[ord],yi));
218       y[i] = yi;
219    }
220 }
221 #endif
222
223 #ifdef PRECISION16
224 void iir_mem2(const spx_sig_t *x, const spx_coef_t *den, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
225 {
226    int i,j;
227    spx_word16_t yi,nyi;
228
229    for (i=0;i<N;i++)
230    {
231       yi = EXTRACT16(PSHR32(SATURATE(x[i] + SHL32(mem[0],1),536870911),SIG_SHIFT));
232       nyi = NEG16(yi);
233       for (j=0;j<ord-1;j++)
234       {
235          mem[j] = MAC16_16(mem[j+1],den[j+1],nyi);
236       }
237       mem[ord-1] = MULT16_16(den[ord],nyi);
238       y[i] = SHL32(EXTEND32(yi),SIG_SHIFT);
239    }
240 }
241 #else
242 void iir_mem2(const spx_sig_t *x, const spx_coef_t *den, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
243 {
244    int i,j;
245    spx_word32_t xi,yi,nyi;
246
247    for (i=0;i<N;i++)
248    {
249       xi=SATURATE(x[i],805306368);
250       yi = SATURATE(xi + SHL32(mem[0],2),805306368);
251       nyi = NEG32(yi);
252       for (j=0;j<ord-1;j++)
253       {
254          mem[j] = MAC16_32_Q15(mem[j+1],den[j+1],nyi);
255       }
256       mem[ord-1] = MULT16_32_Q15(den[ord],nyi);
257       y[i] = yi;
258    }
259 }
260 #endif
261
262 #endif
263
264
265 #ifdef PRECISION16
266 void fir_mem2(const spx_sig_t *x, const spx_coef_t *num, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
267 {
268    int i,j;
269    spx_word16_t xi,yi;
270
271    for (i=0;i<N;i++)
272    {
273       xi= EXTRACT16(PSHR32(SATURATE(x[i],536870911),SIG_SHIFT));
274       yi = EXTRACT16(PSHR32(SATURATE(x[i] + SHL32(mem[0],1),536870911),SIG_SHIFT));
275       for (j=0;j<ord-1;j++)
276       {
277          mem[j] = MAC16_16(mem[j+1], num[j+1],xi);
278       }
279       mem[ord-1] = MULT16_16(num[ord],xi);
280       y[i] = SHL32(EXTEND32(yi),SIG_SHIFT);
281    }
282 }
283 #else
284 void fir_mem2(const spx_sig_t *x, const spx_coef_t *num, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
285 {
286    int i,j;
287    spx_word32_t xi,yi;
288
289    for (i=0;i<N;i++)
290    {
291       xi=SATURATE(x[i],805306368);
292       yi = xi + SHL32(mem[0],2);
293       for (j=0;j<ord-1;j++)
294       {
295          mem[j] = MAC16_32_Q15(mem[j+1], num[j+1],xi);
296       }
297       mem[ord-1] = MULT16_32_Q15(num[ord],xi);
298       y[i] = SATURATE(yi,805306368);
299    }
300 }
301 #endif
302
303 #else
304
305
306
307 spx_word16_t compute_rms(const spx_sig_t *x, int len)
308 {
309    int i;
310    float sum=0;
311    for (i=0;i<len;i++)
312    {
313       sum += x[i]*x[i];
314    }
315    return sqrt(.1+sum/len);
316 }
317
318 #ifdef _USE_SSE
319 #include "filters_sse.h"
320 #else
321
322
323 void filter_mem2(const spx_sig_t *x, const spx_coef_t *num, const spx_coef_t *den, spx_sig_t *y, int N, int ord,  spx_mem_t *mem)
324 {
325    int i,j;
326    float xi,yi;
327    for (i=0;i<N;i++)
328    {
329       xi=x[i];
330       y[i] = num[0]*xi + mem[0];
331       yi=y[i];
332       for (j=0;j<ord-1;j++)
333       {
334          mem[j] = mem[j+1] + num[j+1]*xi - den[j+1]*yi;
335       }
336       mem[ord-1] = num[ord]*xi - den[ord]*yi;
337    }
338 }
339
340
341 void iir_mem2(const spx_sig_t *x, const spx_coef_t *den, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
342 {
343    int i,j;
344    for (i=0;i<N;i++)
345    {
346       y[i] = x[i] + mem[0];
347       for (j=0;j<ord-1;j++)
348       {
349          mem[j] = mem[j+1] - den[j+1]*y[i];
350       }
351       mem[ord-1] = - den[ord]*y[i];
352    }
353 }
354
355
356 void fir_mem2(const spx_sig_t *x, const spx_coef_t *num, spx_sig_t *y, int N, int ord, spx_mem_t *mem)
357 {
358    int i,j;
359    float xi;
360    for (i=0;i<N;i++)
361    {
362       xi=x[i];
363       y[i] = num[0]*xi + mem[0];
364       for (j=0;j<ord-1;j++)
365       {
366          mem[j] = mem[j+1] + num[j+1]*xi;
367       }
368       mem[ord-1] = num[ord]*xi;
369    }
370 }
371 #endif
372
373 #endif
374
375
376 void syn_percep_zero(const spx_sig_t *xx, const spx_coef_t *ak, const spx_coef_t *awk1, const spx_coef_t *awk2, spx_sig_t *y, int N, int ord, char *stack)
377 {
378    int i;
379    VARDECL(spx_mem_t *mem);
380    ALLOC(mem, ord, spx_mem_t);
381    for (i=0;i<ord;i++)
382      mem[i]=0;
383    iir_mem2(xx, ak, y, N, ord, mem);
384    for (i=0;i<ord;i++)
385       mem[i]=0;
386    filter_mem2(y, awk1, awk2, y, N, ord, mem);
387 }
388
389 void residue_percep_zero(const spx_sig_t *xx, const spx_coef_t *ak, const spx_coef_t *awk1, const spx_coef_t *awk2, spx_sig_t *y, int N, int ord, char *stack)
390 {
391    int i;
392    VARDECL(spx_mem_t *mem);
393    ALLOC(mem, ord, spx_mem_t);
394    for (i=0;i<ord;i++)
395       mem[i]=0;
396    filter_mem2(xx, ak, awk1, y, N, ord, mem);
397    for (i=0;i<ord;i++)
398      mem[i]=0;
399    fir_mem2(y, awk2, y, N, ord, mem);
400 }
401
402 void compute_impulse_response(const spx_coef_t *ak, const spx_coef_t *awk1, const spx_coef_t *awk2, spx_word16_t *y, int N, int ord, char *stack)
403 {
404    int i,j;
405    spx_word16_t y1, ny1i, ny2i;
406    VARDECL(spx_mem_t *mem1);
407    VARDECL(spx_mem_t *mem2);
408    ALLOC(mem1, ord, spx_mem_t);
409    ALLOC(mem2, ord, spx_mem_t);
410    
411    for (i=0;i<ord+1;i++)
412       y[i] = awk1[i];
413    for (;i<N;i++)
414       y[i] = VERY_SMALL;
415    
416    for (i=0;i<ord;i++)
417       mem1[i] = mem2[i] = 0;
418    for (i=0;i<N;i++)
419    {
420       y1 = ADD16(y[i], EXTRACT16(PSHR32(mem1[0],LPC_SHIFT)));
421       ny1i = NEG16(y1);
422       y[i] = ADD16(SHL16(y1,1), EXTRACT16(PSHR32(mem2[0],LPC_SHIFT)));
423       ny2i = NEG16(y[i]);
424       for (j=0;j<ord-1;j++)
425       {
426          mem1[j] = MAC16_16(mem1[j+1], awk2[j+1],ny1i);
427          mem2[j] = MAC16_16(mem2[j+1], ak[j+1],ny2i);
428       }
429       mem1[ord-1] = MULT16_16(awk2[ord],ny1i);
430       mem2[ord-1] = MULT16_16(ak[ord],ny2i);
431    }
432 }
433
434 void qmf_decomp(const spx_word16_t *xx, const spx_word16_t *aa, spx_sig_t *y1, spx_sig_t *y2, int N, int M, spx_word16_t *mem, char *stack)
435 {
436    int i,j,k,M2;
437    VARDECL(spx_word16_t *a);
438    VARDECL(spx_word16_t *x);
439    spx_word16_t *x2;
440    
441    ALLOC(a, M, spx_word16_t);
442    ALLOC(x, N+M-1, spx_word16_t);
443    x2=x+M-1;
444    M2=M>>1;
445    for (i=0;i<M;i++)
446       a[M-i-1]= aa[i];
447
448    for (i=0;i<M-1;i++)
449       x[i]=mem[M-i-2];
450    for (i=0;i<N;i++)
451       x[i+M-1]=SATURATE(PSHR(xx[i],1),16383);
452    for (i=0,k=0;i<N;i+=2,k++)
453    {
454       y1[k]=0;
455       y2[k]=0;
456       for (j=0;j<M2;j++)
457       {
458          y1[k]=ADD32(y1[k],SHR(MULT16_16(a[j],ADD16(x[i+j],x2[i-j])),1));
459          y2[k]=SUB32(y2[k],SHR(MULT16_16(a[j],SUB16(x[i+j],x2[i-j])),1));
460          j++;
461          y1[k]=ADD32(y1[k],SHR(MULT16_16(a[j],ADD16(x[i+j],x2[i-j])),1));
462          y2[k]=ADD32(y2[k],SHR(MULT16_16(a[j],SUB16(x[i+j],x2[i-j])),1));
463       }
464    }
465    for (i=0;i<M-1;i++)
466      mem[i]=SATURATE(PSHR(xx[N-i-1],1),16383);
467 }
468
469
470 /* By segher */
471 void fir_mem_up(const spx_sig_t *x, const spx_word16_t *a, spx_sig_t *y, int N, int M, spx_word32_t *mem, char *stack)
472    /* assumptions:
473       all odd x[i] are zero -- well, actually they are left out of the array now
474       N and M are multiples of 4 */
475 {
476    int i, j;
477    VARDECL(spx_word16_t *xx);
478    
479    ALLOC(xx, M+N-1, spx_word16_t);
480
481    for (i = 0; i < N/2; i++)
482       xx[2*i] = SHR(x[N/2-1-i],SIG_SHIFT+1);
483    for (i = 0; i < M - 1; i += 2)
484       xx[N+i] = mem[i+1];
485
486    for (i = 0; i < N; i += 4) {
487       spx_sig_t y0, y1, y2, y3;
488       spx_word16_t x0;
489
490       y0 = y1 = y2 = y3 = 0;
491       x0 = xx[N-4-i];
492
493       for (j = 0; j < M; j += 4) {
494          spx_word16_t x1;
495          spx_word16_t a0, a1;
496
497          a0 = a[j];
498          a1 = a[j+1];
499          x1 = xx[N-2+j-i];
500
501          y0 = ADD32(y0,SHR(MULT16_16(a0, x1),1));
502          y1 = ADD32(y1,SHR(MULT16_16(a1, x1),1));
503          y2 = ADD32(y2,SHR(MULT16_16(a0, x0),1));
504          y3 = ADD32(y3,SHR(MULT16_16(a1, x0),1));
505
506          a0 = a[j+2];
507          a1 = a[j+3];
508          x0 = xx[N+j-i];
509
510          y0 = ADD32(y0,SHR(MULT16_16(a0, x0),1));
511          y1 = ADD32(y1,SHR(MULT16_16(a1, x0),1));
512          y2 = ADD32(y2,SHR(MULT16_16(a0, x1),1));
513          y3 = ADD32(y3,SHR(MULT16_16(a1, x1),1));
514       }
515       y[i] = y0;
516       y[i+1] = y1;
517       y[i+2] = y2;
518       y[i+3] = y3;
519    }
520
521    for (i = 0; i < M - 1; i += 2)
522       mem[i+1] = xx[i];
523 }
524
525
526
527 void comb_filter_mem_init (CombFilterMem *mem)
528 {
529    mem->last_pitch=0;
530    mem->last_pitch_gain[0]=mem->last_pitch_gain[1]=mem->last_pitch_gain[2]=0;
531    mem->smooth_gain=1;
532 }
533
534 #ifdef FIXED_POINT
535 #define COMB_STEP 32767
536 #else
537 #define COMB_STEP 1.0
538 #endif
539
540 void comb_filter(
541 spx_sig_t *exc,          /*decoded excitation*/
542 spx_sig_t *new_exc,      /*enhanced excitation*/
543 spx_coef_t *ak,           /*LPC filter coefs*/
544 int p,               /*LPC order*/
545 int nsf,             /*sub-frame size*/
546 int pitch,           /*pitch period*/
547 spx_word16_t *pitch_gain,   /*pitch gain (3-tap)*/
548 spx_word16_t  comb_gain,    /*gain of comb filter*/
549 CombFilterMem *mem
550 )
551 {
552    int i;
553    spx_word16_t exc_energy=0, new_exc_energy=0;
554    spx_word16_t gain;
555    spx_word16_t step;
556    spx_word16_t fact;
557
558    /*Compute excitation amplitude prior to enhancement*/
559    exc_energy = compute_rms(exc, nsf);
560    /*for (i=0;i<nsf;i++)
561      exc_energy+=((float)exc[i])*exc[i];*/
562
563    /*Some gain adjustment if pitch is too high or if unvoiced*/
564 #ifdef FIXED_POINT
565    {
566       spx_word16_t g = gain_3tap_to_1tap(pitch_gain)+gain_3tap_to_1tap(mem->last_pitch_gain);
567       if (g > 166)
568          comb_gain = MULT16_16_Q15(DIV32_16(SHL(165,15),g), comb_gain);
569       if (g < 64)
570          comb_gain = MULT16_16_Q15(SHL(g, 9), comb_gain);
571    }
572 #else
573    {
574       float g=0;
575       g = GAIN_SCALING_1*.5*(gain_3tap_to_1tap(pitch_gain)+gain_3tap_to_1tap(mem->last_pitch_gain));
576       if (g>1.3)
577          comb_gain*=1.3/g;
578       if (g<.5)
579          comb_gain*=2.*g;
580    }
581 #endif
582    step = DIV32(COMB_STEP, nsf);
583    fact=0;
584
585    /*Apply pitch comb-filter (filter out noise between pitch harmonics)*/
586    for (i=0;i<nsf;i++)
587    {
588       spx_word32_t exc1, exc2;
589
590       fact = ADD16(fact,step);
591       
592       exc1 = SHL32(MULT16_32_Q15(SHL16(pitch_gain[0],7),exc[i-pitch+1]) +
593                  MULT16_32_Q15(SHL16(pitch_gain[1],7),exc[i-pitch]) +
594                  MULT16_32_Q15(SHL16(pitch_gain[2],7),exc[i-pitch-1]) , 2);
595       exc2 = SHL32(MULT16_32_Q15(SHL16(mem->last_pitch_gain[0],7),exc[i-mem->last_pitch+1]) +
596                  MULT16_32_Q15(SHL16(mem->last_pitch_gain[1],7),exc[i-mem->last_pitch]) +
597                  MULT16_32_Q15(SHL16(mem->last_pitch_gain[2],7),exc[i-mem->last_pitch-1]),2);
598
599       new_exc[i] = exc[i] + MULT16_32_Q15(comb_gain, ADD32(MULT16_32_Q15(fact,exc1), MULT16_32_Q15(SUB16(COMB_STEP,fact), exc2)));
600    }
601
602    mem->last_pitch_gain[0] = pitch_gain[0];
603    mem->last_pitch_gain[1] = pitch_gain[1];
604    mem->last_pitch_gain[2] = pitch_gain[2];
605    mem->last_pitch = pitch;
606
607    /*Amplitude after enhancement*/
608    new_exc_energy = compute_rms(new_exc, nsf);
609
610    if (exc_energy > new_exc_energy)
611       exc_energy = new_exc_energy;
612    
613    gain = DIV32_16(SHL32(EXTEND32(exc_energy),15),ADD16(1,new_exc_energy));
614
615 #ifdef FIXED_POINT
616    if (gain < 16384)
617       gain = 16384;
618 #else
619    if (gain < .5)
620       gain=.5;
621 #endif
622
623 #ifdef FIXED_POINT
624    for (i=0;i<nsf;i++)
625    {
626       mem->smooth_gain = ADD16(MULT16_16_Q15(31457,mem->smooth_gain), MULT16_16_Q15(1311,gain));
627       new_exc[i] = MULT16_32_Q15(mem->smooth_gain, new_exc[i]);
628    }
629 #else
630    for (i=0;i<nsf;i++)
631    {
632       mem->smooth_gain = .96*mem->smooth_gain + .04*gain;
633       new_exc[i] *= mem->smooth_gain;
634    }
635 #endif
636 }