1 /*
2  * SpanDSP - a series of DSP components for telephony
3  *
4  * echo.c - A line echo canceller.  This code is being developed
5  *          against and partially complies with G168.
6  *
7  * Written by Steve Underwood <steveu@coppice.org>
8  *         and David Rowe <david_at_rowetel_dot_com>
9  *
10  * Copyright (C) 2001, 2003 Steve Underwood, 2007 David Rowe
11  *
12  * Based on a bit from here, a bit from there, eye of toad, ear of
13  * bat, 15 years of failed attempts by David and a few fried brain
14  * cells.
15  *
16  * All rights reserved.
17  *
18  * This program is free software; you can redistribute it and/or modify
19  * it under the terms of the GNU General Public License version 2, as
20  * published by the Free Software Foundation.
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
30  */
31 
32 /*! \file */
33 
34 /* Implementation Notes
35    David Rowe
36    April 2007
37 
38    This code started life as Steve's NLMS algorithm with a tap
39    rotation algorithm to handle divergence during double talk.  I
40    added a Geigel Double Talk Detector (DTD) [2] and performed some
41    G168 tests.  However I had trouble meeting the G168 requirements,
42    especially for double talk - there were always cases where my DTD
43    failed, for example where near end speech was under the 6dB
44    threshold required for declaring double talk.
45 
46    So I tried a two path algorithm [1], which has so far given better
47    results.  The original tap rotation/Geigel algorithm is available
48    in SVN http://svn.rowetel.com/software/oslec/tags/before_16bit.
49    It's probably possible to make it work if some one wants to put some
50    serious work into it.
51 
52    At present no special treatment is provided for tones, which
53    generally cause NLMS algorithms to diverge.  Initial runs of a
54    subset of the G168 tests for tones (e.g ./echo_test 6) show the
55    current algorithm is passing OK, which is kind of surprising.  The
56    full set of tests needs to be performed to confirm this result.
57 
58    One other interesting change is that I have managed to get the NLMS
59    code to work with 16 bit coefficients, rather than the original 32
60    bit coefficents.  This reduces the MIPs and storage required.
61    I evaulated the 16 bit port using g168_tests.sh and listening tests
62    on 4 real-world samples.
63 
64    I also attempted the implementation of a block based NLMS update
65    [2] but although this passes g168_tests.sh it didn't converge well
66    on the real-world samples.  I have no idea why, perhaps a scaling
67    problem.  The block based code is also available in SVN
68    http://svn.rowetel.com/software/oslec/tags/before_16bit.  If this
69    code can be debugged, it will lead to further reduction in MIPS, as
70    the block update code maps nicely onto DSP instruction sets (it's a
71    dot product) compared to the current sample-by-sample update.
72 
73    Steve also has some nice notes on echo cancellers in echo.h
74 
75    References:
76 
77    [1] Ochiai, Areseki, and Ogihara, "Echo Canceller with Two Echo
78        Path Models", IEEE Transactions on communications, COM-25,
79        No. 6, June
80        1977.
81        http://www.rowetel.com/images/echo/dual_path_paper.pdf
82 
83    [2] The classic, very useful paper that tells you how to
84        actually build a real world echo canceller:
85 	 Messerschmitt, Hedberg, Cole, Haoui, Winship, "Digital Voice
86 	 Echo Canceller with a TMS320020,
87 	 http://www.rowetel.com/images/echo/spra129.pdf
88 
89    [3] I have written a series of blog posts on this work, here is
90        Part 1: http://www.rowetel.com/blog/?p=18
91 
92    [4] The source code http://svn.rowetel.com/software/oslec/
93 
94    [5] A nice reference on LMS filters:
95 	 http://en.wikipedia.org/wiki/Least_mean_squares_filter
96 
97    Credits:
98 
99    Thanks to Steve Underwood, Jean-Marc Valin, and Ramakrishnan
100    Muthukrishnan for their suggestions and email discussions.  Thanks
101    also to those people who collected echo samples for me such as
102    Mark, Pawel, and Pavel.
103 */
104 
105 #include <linux/kernel.h>
106 #include <linux/module.h>
107 #include <linux/slab.h>
108 
109 #include "echo.h"
110 
111 #define MIN_TX_POWER_FOR_ADAPTION	64
112 #define MIN_RX_POWER_FOR_ADAPTION	64
113 #define DTD_HANGOVER			600	/* 600 samples, or 75ms     */
114 #define DC_LOG2BETA			3	/* log2() of DC filter Beta */
115 
116 /* adapting coeffs using the traditional stochastic descent (N)LMS algorithm */
117 
lms_adapt_bg(struct oslec_state * ec,int clean,int shift)118 static inline void lms_adapt_bg(struct oslec_state *ec, int clean, int shift)
119 {
120 	int i;
121 
122 	int offset1;
123 	int offset2;
124 	int factor;
125 	int exp;
126 
127 	if (shift > 0)
128 		factor = clean << shift;
129 	else
130 		factor = clean >> -shift;
131 
132 	/* Update the FIR taps */
133 
134 	offset2 = ec->curr_pos;
135 	offset1 = ec->taps - offset2;
136 
137 	for (i = ec->taps - 1; i >= offset1; i--) {
138 		exp = (ec->fir_state_bg.history[i - offset1] * factor);
139 		ec->fir_taps16[1][i] += (int16_t) ((exp + (1 << 14)) >> 15);
140 	}
141 	for (; i >= 0; i--) {
142 		exp = (ec->fir_state_bg.history[i + offset2] * factor);
143 		ec->fir_taps16[1][i] += (int16_t) ((exp + (1 << 14)) >> 15);
144 	}
145 }
146 
top_bit(unsigned int bits)147 static inline int top_bit(unsigned int bits)
148 {
149 	if (bits == 0)
150 		return -1;
151 	else
152 		return (int)fls((int32_t) bits) - 1;
153 }
154 
oslec_create(int len,int adaption_mode)155 struct oslec_state *oslec_create(int len, int adaption_mode)
156 {
157 	struct oslec_state *ec;
158 	int i;
159 	const int16_t *history;
160 
161 	ec = kzalloc(sizeof(*ec), GFP_KERNEL);
162 	if (!ec)
163 		return NULL;
164 
165 	ec->taps = len;
166 	ec->log2taps = top_bit(len);
167 	ec->curr_pos = ec->taps - 1;
168 
169 	ec->fir_taps16[0] =
170 	    kcalloc(ec->taps, sizeof(int16_t), GFP_KERNEL);
171 	if (!ec->fir_taps16[0])
172 		goto error_oom_0;
173 
174 	ec->fir_taps16[1] =
175 	    kcalloc(ec->taps, sizeof(int16_t), GFP_KERNEL);
176 	if (!ec->fir_taps16[1])
177 		goto error_oom_1;
178 
179 	history = fir16_create(&ec->fir_state, ec->fir_taps16[0], ec->taps);
180 	if (!history)
181 		goto error_state;
182 	history = fir16_create(&ec->fir_state_bg, ec->fir_taps16[1], ec->taps);
183 	if (!history)
184 		goto error_state_bg;
185 
186 	for (i = 0; i < 5; i++)
187 		ec->xvtx[i] = ec->yvtx[i] = ec->xvrx[i] = ec->yvrx[i] = 0;
188 
189 	ec->cng_level = 1000;
190 	oslec_adaption_mode(ec, adaption_mode);
191 
192 	ec->snapshot = kcalloc(ec->taps, sizeof(int16_t), GFP_KERNEL);
193 	if (!ec->snapshot)
194 		goto error_snap;
195 
196 	ec->cond_met = 0;
197 	ec->pstates = 0;
198 	ec->ltxacc = ec->lrxacc = ec->lcleanacc = ec->lclean_bgacc = 0;
199 	ec->ltx = ec->lrx = ec->lclean = ec->lclean_bg = 0;
200 	ec->tx_1 = ec->tx_2 = ec->rx_1 = ec->rx_2 = 0;
201 	ec->lbgn = ec->lbgn_acc = 0;
202 	ec->lbgn_upper = 200;
203 	ec->lbgn_upper_acc = ec->lbgn_upper << 13;
204 
205 	return ec;
206 
207 error_snap:
208 	fir16_free(&ec->fir_state_bg);
209 error_state_bg:
210 	fir16_free(&ec->fir_state);
211 error_state:
212 	kfree(ec->fir_taps16[1]);
213 error_oom_1:
214 	kfree(ec->fir_taps16[0]);
215 error_oom_0:
216 	kfree(ec);
217 	return NULL;
218 }
219 EXPORT_SYMBOL_GPL(oslec_create);
220 
oslec_free(struct oslec_state * ec)221 void oslec_free(struct oslec_state *ec)
222 {
223 	int i;
224 
225 	fir16_free(&ec->fir_state);
226 	fir16_free(&ec->fir_state_bg);
227 	for (i = 0; i < 2; i++)
228 		kfree(ec->fir_taps16[i]);
229 	kfree(ec->snapshot);
230 	kfree(ec);
231 }
232 EXPORT_SYMBOL_GPL(oslec_free);
233 
oslec_adaption_mode(struct oslec_state * ec,int adaption_mode)234 void oslec_adaption_mode(struct oslec_state *ec, int adaption_mode)
235 {
236 	ec->adaption_mode = adaption_mode;
237 }
238 EXPORT_SYMBOL_GPL(oslec_adaption_mode);
239 
oslec_flush(struct oslec_state * ec)240 void oslec_flush(struct oslec_state *ec)
241 {
242 	int i;
243 
244 	ec->ltxacc = ec->lrxacc = ec->lcleanacc = ec->lclean_bgacc = 0;
245 	ec->ltx = ec->lrx = ec->lclean = ec->lclean_bg = 0;
246 	ec->tx_1 = ec->tx_2 = ec->rx_1 = ec->rx_2 = 0;
247 
248 	ec->lbgn = ec->lbgn_acc = 0;
249 	ec->lbgn_upper = 200;
250 	ec->lbgn_upper_acc = ec->lbgn_upper << 13;
251 
252 	ec->nonupdate_dwell = 0;
253 
254 	fir16_flush(&ec->fir_state);
255 	fir16_flush(&ec->fir_state_bg);
256 	ec->fir_state.curr_pos = ec->taps - 1;
257 	ec->fir_state_bg.curr_pos = ec->taps - 1;
258 	for (i = 0; i < 2; i++)
259 		memset(ec->fir_taps16[i], 0, ec->taps * sizeof(int16_t));
260 
261 	ec->curr_pos = ec->taps - 1;
262 	ec->pstates = 0;
263 }
264 EXPORT_SYMBOL_GPL(oslec_flush);
265 
oslec_snapshot(struct oslec_state * ec)266 void oslec_snapshot(struct oslec_state *ec)
267 {
268 	memcpy(ec->snapshot, ec->fir_taps16[0], ec->taps * sizeof(int16_t));
269 }
270 EXPORT_SYMBOL_GPL(oslec_snapshot);
271 
272 /* Dual Path Echo Canceller */
273 
oslec_update(struct oslec_state * ec,int16_t tx,int16_t rx)274 int16_t oslec_update(struct oslec_state *ec, int16_t tx, int16_t rx)
275 {
276 	int32_t echo_value;
277 	int clean_bg;
278 	int tmp;
279 	int tmp1;
280 
281 	/*
282 	 * Input scaling was found be required to prevent problems when tx
283 	 * starts clipping.  Another possible way to handle this would be the
284 	 * filter coefficent scaling.
285 	 */
286 
287 	ec->tx = tx;
288 	ec->rx = rx;
289 	tx >>= 1;
290 	rx >>= 1;
291 
292 	/*
293 	 * Filter DC, 3dB point is 160Hz (I think), note 32 bit precision
294 	 * required otherwise values do not track down to 0. Zero at DC, Pole
295 	 * at (1-Beta) on real axis.  Some chip sets (like Si labs) don't
296 	 * need this, but something like a $10 X100P card does.  Any DC really
297 	 * slows down convergence.
298 	 *
299 	 * Note: removes some low frequency from the signal, this reduces the
300 	 * speech quality when listening to samples through headphones but may
301 	 * not be obvious through a telephone handset.
302 	 *
303 	 * Note that the 3dB frequency in radians is approx Beta, e.g. for Beta
304 	 * = 2^(-3) = 0.125, 3dB freq is 0.125 rads = 159Hz.
305 	 */
306 
307 	if (ec->adaption_mode & ECHO_CAN_USE_RX_HPF) {
308 		tmp = rx << 15;
309 
310 		/*
311 		 * Make sure the gain of the HPF is 1.0. This can still
312 		 * saturate a little under impulse conditions, and it might
313 		 * roll to 32768 and need clipping on sustained peak level
314 		 * signals. However, the scale of such clipping is small, and
315 		 * the error due to any saturation should not markedly affect
316 		 * the downstream processing.
317 		 */
318 		tmp -= (tmp >> 4);
319 
320 		ec->rx_1 += -(ec->rx_1 >> DC_LOG2BETA) + tmp - ec->rx_2;
321 
322 		/*
323 		 * hard limit filter to prevent clipping.  Note that at this
324 		 * stage rx should be limited to +/- 16383 due to right shift
325 		 * above
326 		 */
327 		tmp1 = ec->rx_1 >> 15;
328 		if (tmp1 > 16383)
329 			tmp1 = 16383;
330 		if (tmp1 < -16383)
331 			tmp1 = -16383;
332 		rx = tmp1;
333 		ec->rx_2 = tmp;
334 	}
335 
336 	/* Block average of power in the filter states.  Used for
337 	   adaption power calculation. */
338 
339 	{
340 		int new, old;
341 
342 		/* efficient "out with the old and in with the new" algorithm so
343 		   we don't have to recalculate over the whole block of
344 		   samples. */
345 		new = (int)tx * (int)tx;
346 		old = (int)ec->fir_state.history[ec->fir_state.curr_pos] *
347 		    (int)ec->fir_state.history[ec->fir_state.curr_pos];
348 		ec->pstates +=
349 		    ((new - old) + (1 << (ec->log2taps - 1))) >> ec->log2taps;
350 		if (ec->pstates < 0)
351 			ec->pstates = 0;
352 	}
353 
354 	/* Calculate short term average levels using simple single pole IIRs */
355 
356 	ec->ltxacc += abs(tx) - ec->ltx;
357 	ec->ltx = (ec->ltxacc + (1 << 4)) >> 5;
358 	ec->lrxacc += abs(rx) - ec->lrx;
359 	ec->lrx = (ec->lrxacc + (1 << 4)) >> 5;
360 
361 	/* Foreground filter */
362 
363 	ec->fir_state.coeffs = ec->fir_taps16[0];
364 	echo_value = fir16(&ec->fir_state, tx);
365 	ec->clean = rx - echo_value;
366 	ec->lcleanacc += abs(ec->clean) - ec->lclean;
367 	ec->lclean = (ec->lcleanacc + (1 << 4)) >> 5;
368 
369 	/* Background filter */
370 
371 	echo_value = fir16(&ec->fir_state_bg, tx);
372 	clean_bg = rx - echo_value;
373 	ec->lclean_bgacc += abs(clean_bg) - ec->lclean_bg;
374 	ec->lclean_bg = (ec->lclean_bgacc + (1 << 4)) >> 5;
375 
376 	/* Background Filter adaption */
377 
378 	/* Almost always adap bg filter, just simple DT and energy
379 	   detection to minimise adaption in cases of strong double talk.
380 	   However this is not critical for the dual path algorithm.
381 	 */
382 	ec->factor = 0;
383 	ec->shift = 0;
384 	if (!ec->nonupdate_dwell) {
385 		int p, logp, shift;
386 
387 		/* Determine:
388 
389 		   f = Beta * clean_bg_rx/P ------ (1)
390 
391 		   where P is the total power in the filter states.
392 
393 		   The Boffins have shown that if we obey (1) we converge
394 		   quickly and avoid instability.
395 
396 		   The correct factor f must be in Q30, as this is the fixed
397 		   point format required by the lms_adapt_bg() function,
398 		   therefore the scaled version of (1) is:
399 
400 		   (2^30) * f  = (2^30) * Beta * clean_bg_rx/P
401 		   factor      = (2^30) * Beta * clean_bg_rx/P     ----- (2)
402 
403 		   We have chosen Beta = 0.25 by experiment, so:
404 
405 		   factor      = (2^30) * (2^-2) * clean_bg_rx/P
406 
407 		   (30 - 2 - log2(P))
408 		   factor      = clean_bg_rx 2                     ----- (3)
409 
410 		   To avoid a divide we approximate log2(P) as top_bit(P),
411 		   which returns the position of the highest non-zero bit in
412 		   P.  This approximation introduces an error as large as a
413 		   factor of 2, but the algorithm seems to handle it OK.
414 
415 		   Come to think of it a divide may not be a big deal on a
416 		   modern DSP, so its probably worth checking out the cycles
417 		   for a divide versus a top_bit() implementation.
418 		 */
419 
420 		p = MIN_TX_POWER_FOR_ADAPTION + ec->pstates;
421 		logp = top_bit(p) + ec->log2taps;
422 		shift = 30 - 2 - logp;
423 		ec->shift = shift;
424 
425 		lms_adapt_bg(ec, clean_bg, shift);
426 	}
427 
428 	/* very simple DTD to make sure we dont try and adapt with strong
429 	   near end speech */
430 
431 	ec->adapt = 0;
432 	if ((ec->lrx > MIN_RX_POWER_FOR_ADAPTION) && (ec->lrx > ec->ltx))
433 		ec->nonupdate_dwell = DTD_HANGOVER;
434 	if (ec->nonupdate_dwell)
435 		ec->nonupdate_dwell--;
436 
437 	/* Transfer logic */
438 
439 	/* These conditions are from the dual path paper [1], I messed with
440 	   them a bit to improve performance. */
441 
442 	if ((ec->adaption_mode & ECHO_CAN_USE_ADAPTION) &&
443 	    (ec->nonupdate_dwell == 0) &&
444 	    /* (ec->Lclean_bg < 0.875*ec->Lclean) */
445 	    (8 * ec->lclean_bg < 7 * ec->lclean) &&
446 	    /* (ec->Lclean_bg < 0.125*ec->Ltx) */
447 	    (8 * ec->lclean_bg < ec->ltx)) {
448 		if (ec->cond_met == 6) {
449 			/*
450 			 * BG filter has had better results for 6 consecutive
451 			 * samples
452 			 */
453 			ec->adapt = 1;
454 			memcpy(ec->fir_taps16[0], ec->fir_taps16[1],
455 			       ec->taps * sizeof(int16_t));
456 		} else
457 			ec->cond_met++;
458 	} else
459 		ec->cond_met = 0;
460 
461 	/* Non-Linear Processing */
462 
463 	ec->clean_nlp = ec->clean;
464 	if (ec->adaption_mode & ECHO_CAN_USE_NLP) {
465 		/*
466 		 * Non-linear processor - a fancy way to say "zap small
467 		 * signals, to avoid residual echo due to (uLaw/ALaw)
468 		 * non-linearity in the channel.".
469 		 */
470 
471 		if ((16 * ec->lclean < ec->ltx)) {
472 			/*
473 			 * Our e/c has improved echo by at least 24 dB (each
474 			 * factor of 2 is 6dB, so 2*2*2*2=16 is the same as
475 			 * 6+6+6+6=24dB)
476 			 */
477 			if (ec->adaption_mode & ECHO_CAN_USE_CNG) {
478 				ec->cng_level = ec->lbgn;
479 
480 				/*
481 				 * Very elementary comfort noise generation.
482 				 * Just random numbers rolled off very vaguely
483 				 * Hoth-like.  DR: This noise doesn't sound
484 				 * quite right to me - I suspect there are some
485 				 * overflow issues in the filtering as it's too
486 				 * "crackly".
487 				 * TODO: debug this, maybe just play noise at
488 				 * high level or look at spectrum.
489 				 */
490 
491 				ec->cng_rndnum =
492 				    1664525U * ec->cng_rndnum + 1013904223U;
493 				ec->cng_filter =
494 				    ((ec->cng_rndnum & 0xFFFF) - 32768 +
495 				     5 * ec->cng_filter) >> 3;
496 				ec->clean_nlp =
497 				    (ec->cng_filter * ec->cng_level * 8) >> 14;
498 
499 			} else if (ec->adaption_mode & ECHO_CAN_USE_CLIP) {
500 				/* This sounds much better than CNG */
501 				if (ec->clean_nlp > ec->lbgn)
502 					ec->clean_nlp = ec->lbgn;
503 				if (ec->clean_nlp < -ec->lbgn)
504 					ec->clean_nlp = -ec->lbgn;
505 			} else {
506 				/*
507 				 * just mute the residual, doesn't sound very
508 				 * good, used mainly in G168 tests
509 				 */
510 				ec->clean_nlp = 0;
511 			}
512 		} else {
513 			/*
514 			 * Background noise estimator.  I tried a few
515 			 * algorithms here without much luck.  This very simple
516 			 * one seems to work best, we just average the level
517 			 * using a slow (1 sec time const) filter if the
518 			 * current level is less than a (experimentally
519 			 * derived) constant.  This means we dont include high
520 			 * level signals like near end speech.  When combined
521 			 * with CNG or especially CLIP seems to work OK.
522 			 */
523 			if (ec->lclean < 40) {
524 				ec->lbgn_acc += abs(ec->clean) - ec->lbgn;
525 				ec->lbgn = (ec->lbgn_acc + (1 << 11)) >> 12;
526 			}
527 		}
528 	}
529 
530 	/* Roll around the taps buffer */
531 	if (ec->curr_pos <= 0)
532 		ec->curr_pos = ec->taps;
533 	ec->curr_pos--;
534 
535 	if (ec->adaption_mode & ECHO_CAN_DISABLE)
536 		ec->clean_nlp = rx;
537 
538 	/* Output scaled back up again to match input scaling */
539 
540 	return (int16_t) ec->clean_nlp << 1;
541 }
542 EXPORT_SYMBOL_GPL(oslec_update);
543 
544 /* This function is separated from the echo canceller is it is usually called
545    as part of the tx process.  See rx HP (DC blocking) filter above, it's
546    the same design.
547 
548    Some soft phones send speech signals with a lot of low frequency
549    energy, e.g. down to 20Hz.  This can make the hybrid non-linear
550    which causes the echo canceller to fall over.  This filter can help
551    by removing any low frequency before it gets to the tx port of the
552    hybrid.
553 
554    It can also help by removing and DC in the tx signal.  DC is bad
555    for LMS algorithms.
556 
557    This is one of the classic DC removal filters, adjusted to provide
558    sufficient bass rolloff to meet the above requirement to protect hybrids
559    from things that upset them. The difference between successive samples
560    produces a lousy HPF, and then a suitably placed pole flattens things out.
561    The final result is a nicely rolled off bass end. The filtering is
562    implemented with extended fractional precision, which noise shapes things,
563    giving very clean DC removal.
564 */
565 
oslec_hpf_tx(struct oslec_state * ec,int16_t tx)566 int16_t oslec_hpf_tx(struct oslec_state *ec, int16_t tx)
567 {
568 	int tmp;
569 	int tmp1;
570 
571 	if (ec->adaption_mode & ECHO_CAN_USE_TX_HPF) {
572 		tmp = tx << 15;
573 
574 		/*
575 		 * Make sure the gain of the HPF is 1.0. The first can still
576 		 * saturate a little under impulse conditions, and it might
577 		 * roll to 32768 and need clipping on sustained peak level
578 		 * signals. However, the scale of such clipping is small, and
579 		 * the error due to any saturation should not markedly affect
580 		 * the downstream processing.
581 		 */
582 		tmp -= (tmp >> 4);
583 
584 		ec->tx_1 += -(ec->tx_1 >> DC_LOG2BETA) + tmp - ec->tx_2;
585 		tmp1 = ec->tx_1 >> 15;
586 		if (tmp1 > 32767)
587 			tmp1 = 32767;
588 		if (tmp1 < -32767)
589 			tmp1 = -32767;
590 		tx = tmp1;
591 		ec->tx_2 = tmp;
592 	}
593 
594 	return tx;
595 }
596 EXPORT_SYMBOL_GPL(oslec_hpf_tx);
597 
598 MODULE_LICENSE("GPL");
599 MODULE_AUTHOR("David Rowe");
600 MODULE_DESCRIPTION("Open Source Line Echo Canceller");
601 MODULE_VERSION("0.3.0");
602