1 /* Copyright (C) 2006-2008 CSIRO, Jean-Marc Valin, Xiph.Org Foundation
4 Shaped comb-allpass filter for channel decorrelation
6 Redistribution and use in source and binary forms, with or without
7 modification, are permitted provided that the following conditions are
10 1. Redistributions of source code must retain the above copyright notice,
11 this list of conditions and the following disclaimer.
13 2. Redistributions in binary form must reproduce the above copyright
14 notice, this list of conditions and the following disclaimer in the
15 documentation and/or other materials provided with the distribution.
17 3. The name of the author may not be used to endorse or promote products
18 derived from this software without specific prior written permission.
20 THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
21 IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
22 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
23 DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
24 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
25 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
26 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
27 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
28 STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
29 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30 POSSIBILITY OF SUCH DAMAGE.
34 The algorithm implemented here is described in:
36 * J.-M. Valin, Perceptually-Motivated Nonlinear Channel Decorrelation For
37 Stereo Acoustic Echo Cancellation, Accepted for Joint Workshop on
38 Handsfree Speech Communication and Microphone Arrays (HSCMA), 2008.
39 http://people.xiph.org/~jm/papers/valin_hscma2008.pdf
47 #include "speex/speex_echo.h"
48 #include "vorbis_psy.h"
50 #include "os_support.h"
56 #define M_PI 3.14159261
59 #define ALLPASS_ORDER 20
61 struct SpeexDecorrState_ {
67 struct drft_lookup lookup;
75 /* Per-channel stuff */
77 float (*ring)[ALLPASS_ORDER];
85 EXPORT SpeexDecorrState *speex_decorrelate_new(int rate, int channels, int frame_size)
88 SpeexDecorrState *st = speex_alloc(sizeof(SpeexDecorrState));
90 st->channels = channels;
91 st->frame_size = frame_size;
93 st->psy = vorbis_psy_init(rate, 2*frame_size);
94 spx_drft_init(&st->lookup, 2*frame_size);
95 st->wola_mem = speex_alloc(frame_size*sizeof(float));
96 st->curve = speex_alloc(frame_size*sizeof(float));
98 st->y = speex_alloc(frame_size*sizeof(float));
100 st->buff = speex_alloc(channels*2*frame_size*sizeof(float));
101 st->ringID = speex_alloc(channels*sizeof(int));
102 st->order = speex_alloc(channels*sizeof(int));
103 st->alpha = speex_alloc(channels*sizeof(float));
104 st->ring = speex_alloc(channels*ALLPASS_ORDER*sizeof(float));
106 /*FIXME: The +20 is there only as a kludge for ALL_PASS_OLA*/
107 st->vorbis_win = speex_alloc((2*frame_size+20)*sizeof(float));
108 for (i=0;i<2*frame_size;i++)
109 st->vorbis_win[i] = sin(.5*M_PI* sin(M_PI*i/(2*frame_size))*sin(M_PI*i/(2*frame_size)) );
112 for (ch=0;ch<channels;ch++)
114 for (i=0;i<ALLPASS_ORDER;i++)
123 static float uni_rand(int *seed)
125 const unsigned int jflone = 0x3f800000;
126 const unsigned int jflmsk = 0x007fffff;
127 union {int i; float f;} ran;
128 *seed = 1664525 * *seed + 1013904223;
129 ran.i = jflone | (jflmsk & *seed);
134 static unsigned int irand(int *seed)
136 *seed = 1664525 * *seed + 1013904223;
137 return ((unsigned int)*seed)>>16;
141 EXPORT void speex_decorrelate(SpeexDecorrState *st, const spx_int16_t *in, spx_int16_t *out, int strength)
151 amount = .01*strength;
152 for (ch=0;ch<st->channels;ch++)
155 int N=2*st->frame_size;
166 buff = st->buff+ch*2*st->frame_size;
168 ringID = st->ringID[ch];
169 order = st->order[ch];
170 alpha = st->alpha[ch];
172 for (i=0;i<st->frame_size;i++)
173 buff[i] = buff[i+st->frame_size];
174 for (i=0;i<st->frame_size;i++)
175 buff[i+st->frame_size] = in[i*st->channels+ch];
177 x = buff+st->frame_size;
178 beta = 1.-.3*amount*amount;
180 beta = 1-sqrt(.4*amount);
182 beta = 1-0.63246*amount;
187 for (i=0;i<st->frame_size;i++)
189 st->y[i] = alpha*(x[i-ALLPASS_ORDER+order]-beta*x[i-ALLPASS_ORDER+order-1])*st->vorbis_win[st->frame_size+i+order]
190 + x[i-ALLPASS_ORDER]*st->vorbis_win[st->frame_size+i]
191 - alpha*(ring[ringID]
192 - beta*ring[ringID+1>=order?0:ringID+1]);
193 ring[ringID++]=st->y[i];
194 st->y[i] *= st->vorbis_win[st->frame_size+i];
198 order = order+(irand(&st->seed)%3)-1;
203 /*order = 5+(irand(&st->seed)%6);*/
204 max_alpha = pow(.96+.04*(amount-1),order);
205 if (max_alpha > .98/(1.+beta2))
206 max_alpha = .98/(1.+beta2);
208 alpha = alpha + .4*uni_rand(&st->seed);
209 if (alpha > max_alpha)
211 if (alpha < -max_alpha)
213 for (i=0;i<ALLPASS_ORDER;i++)
216 for (i=0;i<st->frame_size;i++)
218 float tmp = alpha*(x[i-ALLPASS_ORDER+order]-beta*x[i-ALLPASS_ORDER+order-1])*st->vorbis_win[i+order]
219 + x[i-ALLPASS_ORDER]*st->vorbis_win[i]
220 - alpha*(ring[ringID]
221 - beta*ring[ringID+1>=order?0:ringID+1]);
223 tmp *= st->vorbis_win[i];
232 for (i=0;i<2*st->frame_size;i++)
234 //float coef = .5*0.78130;
235 float coef = M_PI*0.075063 * 0.93763 * amount * .8 * 0.707;
236 compute_curve(st->psy, buff, st->curve);
237 for (i=1;i<st->frame_size;i++)
242 x1 = uni_rand(&st->seed);
243 x2 = uni_rand(&st->seed);
244 } while (x1*x1+x2*x2 > 1.);
245 gain = coef*sqrt(.1+st->curve[i]);
246 frame[2*i-1] = gain*x1;
247 frame[2*i] = gain*x2;
249 frame[0] = coef*uni_rand(&st->seed)*sqrt(.1+st->curve[0]);
250 frame[2*st->frame_size-1] = coef*uni_rand(&st->seed)*sqrt(.1+st->curve[st->frame_size-1]);
251 spx_drft_backward(&st->lookup,frame);
252 for (i=0;i<2*st->frame_size;i++)
253 frame[i] *= st->vorbis_win[i];
256 for (i=0;i<st->frame_size;i++)
259 float tmp = st->y[i] + frame[i] + st->wola_mem[i];
260 st->wola_mem[i] = frame[i+st->frame_size];
262 float tmp = st->y[i];
268 out[i*st->channels+ch] = tmp;
271 st->ringID[ch] = ringID;
272 st->order[ch] = order;
273 st->alpha[ch] = alpha;
278 EXPORT void speex_decorrelate_destroy(SpeexDecorrState *st)
281 vorbis_psy_destroy(st->psy);
282 speex_free(st->wola_mem);
283 speex_free(st->curve);
285 speex_free(st->buff);
286 speex_free(st->ring);
287 speex_free(st->ringID);
288 speex_free(st->alpha);
289 speex_free(st->vorbis_win);
290 speex_free(st->order);