2 * SpanDSP - a series of DSP components for telephony
4 * echo.c - A line echo canceller. This code is being developed
5 * against and partially complies with G168.
7 * Written by Steve Underwood <steveu@coppice.org>
8 * and David Rowe <david_at_rowetel_dot_com>
10 * Copyright (C) 2001, 2003 Steve Underwood, 2007 David Rowe
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
16 * All rights reserved.
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.
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.
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.
31 * $Id: echo.c,v 1.20 2006/12/01 18:00:48 steveu Exp $
36 /* Implementation Notes
40 This code started life as Steve's NLMS algorithm with a tap
41 rotation algorithm to handle divergence during double talk. I
42 added a Geigel Double Talk Detector (DTD) [2] and performed some
43 G168 tests. However I had trouble meeting the G168 requirements,
44 especially for double talk - there were always cases where my DTD
45 failed, for example where near end speech was under the 6dB
46 threshold required for declaring double talk.
48 So I tried a two path algorithm [1], which has so far given better
49 results. The original tap rotation/Geigel algorithm is available
50 in SVN http://svn.rowetel.com/software/oslec/tags/before_16bit.
51 It's probably possible to make it work if some one wants to put some
54 At present no special treatment is provided for tones, which
55 generally cause NLMS algorithms to diverge. Initial runs of a
56 subset of the G168 tests for tones (e.g ./echo_test 6) show the
57 current algorithm is passing OK, which is kind of surprising. The
58 full set of tests needs to be performed to confirm this result.
60 One other interesting change is that I have managed to get the NLMS
61 code to work with 16 bit coefficients, rather than the original 32
62 bit coefficents. This reduces the MIPs and storage required.
63 I evaulated the 16 bit port using g168_tests.sh and listening tests
64 on 4 real-world samples.
66 I also attempted the implementation of a block based NLMS update
67 [2] but although this passes g168_tests.sh it didn't converge well
68 on the real-world samples. I have no idea why, perhaps a scaling
69 problem. The block based code is also available in SVN
70 http://svn.rowetel.com/software/oslec/tags/before_16bit. If this
71 code can be debugged, it will lead to further reduction in MIPS, as
72 the block update code maps nicely onto DSP instruction sets (it's a
73 dot product) compared to the current sample-by-sample update.
75 Steve also has some nice notes on echo cancellers in echo.h
80 [1] Ochiai, Areseki, and Ogihara, "Echo Canceller with Two Echo
81 Path Models", IEEE Transactions on communications, COM-25,
84 http://www.rowetel.com/images/echo/dual_path_paper.pdf
86 [2] The classic, very useful paper that tells you how to
87 actually build a real world echo canceller:
88 Messerschmitt, Hedberg, Cole, Haoui, Winship, "Digital Voice
89 Echo Canceller with a TMS320020,
90 http://www.rowetel.com/images/echo/spra129.pdf
92 [3] I have written a series of blog posts on this work, here is
93 Part 1: http://www.rowetel.com/blog/?p=18
95 [4] The source code http://svn.rowetel.com/software/oslec/
97 [5] A nice reference on LMS filters:
98 http://en.wikipedia.org/wiki/Least_mean_squares_filter
102 Thanks to Steve Underwood, Jean-Marc Valin, and Ramakrishnan
103 Muthukrishnan for their suggestions and email discussions. Thanks
104 also to those people who collected echo samples for me such as
105 Mark, Pawel, and Pavel.
108 #include <linux/kernel.h> /* We're doing kernel work */
109 #include <linux/module.h>
110 #include <linux/kernel.h>
111 #include <linux/slab.h>
112 #define malloc(a) kmalloc((a), GFP_KERNEL)
113 #define free(a) kfree(a)
115 #include "bit_operations.h"
118 #define MIN_TX_POWER_FOR_ADAPTION 64
119 #define MIN_RX_POWER_FOR_ADAPTION 64
120 #define DTD_HANGOVER 600 /* 600 samples, or 75ms */
121 #define DC_LOG2BETA 3 /* log2() of DC filter Beta */
123 /*-----------------------------------------------------------------------*\
125 \*-----------------------------------------------------------------------*/
127 /* adapting coeffs using the traditional stochastic descent (N)LMS algorithm */
130 #ifdef __BLACKFIN_ASM__
131 static void __inline__ lms_adapt_bg(echo_can_state_t *ec, int clean, int shift)
142 factor = clean << shift;
144 factor = clean >> -shift;
146 /* Update the FIR taps */
148 offset2 = ec->curr_pos;
149 offset1 = ec->taps - offset2;
150 phist = &ec->fir_state_bg.history[offset2];
152 /* st: and en: help us locate the assembler in echo.s */
156 for (i = 0, j = offset2; i < n; i++, j++)
158 exp = *phist++ * factor;
159 ec->fir_taps16[1][i] += (int16_t) ((exp+(1<<14)) >> 15);
163 /* Note the asm for the inner loop above generated by Blackfin gcc
164 4.1.1 is pretty good (note even parallel instructions used):
175 A block based update algorithm would be much faster but the
176 above can't be improved on much. Every instruction saved in
177 the loop above is 2 MIPs/ch! The for loop above is where the
178 Blackfin spends most of it's time - about 17 MIPs/ch measured
179 with speedtest.c with 256 taps (32ms). Write-back and
180 Write-through cache gave about the same performance.
185 IDEAS for further optimisation of lms_adapt_bg():
187 1/ The rounding is quite costly. Could we keep as 32 bit coeffs
188 then make filter pluck the MS 16-bits of the coeffs when filtering?
189 However this would lower potential optimisation of filter, as I
190 think the dual-MAC architecture requires packed 16 bit coeffs.
192 2/ Block based update would be more efficient, as per comments above,
193 could use dual MAC architecture.
195 3/ Look for same sample Blackfin LMS code, see if we can get dual-MAC
198 4/ Execute the whole e/c in a block of say 20ms rather than sample
199 by sample. Processing a few samples every ms is inefficient.
203 static __inline__ void lms_adapt_bg(echo_can_state_t *ec, int clean, int shift)
213 factor = clean << shift;
215 factor = clean >> -shift;
217 /* Update the FIR taps */
219 offset2 = ec->curr_pos;
220 offset1 = ec->taps - offset2;
222 for (i = ec->taps - 1; i >= offset1; i--)
224 exp = (ec->fir_state_bg.history[i - offset1]*factor);
225 ec->fir_taps16[1][i] += (int16_t) ((exp+(1<<14)) >> 15);
229 exp = (ec->fir_state_bg.history[i + offset2]*factor);
230 ec->fir_taps16[1][i] += (int16_t) ((exp+(1<<14)) >> 15);
235 /*- End of function --------------------------------------------------------*/
237 echo_can_state_t *echo_can_create(int len, int adaption_mode)
239 echo_can_state_t *ec;
243 ec = kmalloc(sizeof(*ec), GFP_KERNEL);
246 memset(ec, 0, sizeof(*ec));
249 ec->log2taps = top_bit(len);
250 ec->curr_pos = ec->taps - 1;
252 for (i = 0; i < 2; i++)
254 if ((ec->fir_taps16[i] = (int16_t *) malloc((ec->taps)*sizeof(int16_t))) == NULL)
256 for (j = 0; j < i; j++)
257 kfree(ec->fir_taps16[j]);
261 memset(ec->fir_taps16[i], 0, (ec->taps)*sizeof(int16_t));
264 fir16_create(&ec->fir_state,
267 fir16_create(&ec->fir_state_bg,
272 ec->xvtx[i] = ec->yvtx[i] = ec->xvrx[i] = ec->yvrx[i] = 0;
275 ec->cng_level = 1000;
276 echo_can_adaption_mode(ec, adaption_mode);
278 ec->snapshot = (int16_t*)malloc(ec->taps*sizeof(int16_t));
279 memset(ec->snapshot, 0, sizeof(int16_t)*ec->taps);
283 ec->Ltxacc = ec->Lrxacc = ec->Lcleanacc = ec->Lclean_bgacc = 0;
284 ec->Ltx = ec->Lrx = ec->Lclean = ec->Lclean_bg = 0;
285 ec->tx_1 = ec->tx_2 = ec->rx_1 = ec->rx_2 = 0;
286 ec->Lbgn = ec->Lbgn_acc = 0;
287 ec->Lbgn_upper = 200;
288 ec->Lbgn_upper_acc = ec->Lbgn_upper << 13;
292 EXPORT_SYMBOL_GPL(echo_can_create);
293 /*- End of function --------------------------------------------------------*/
295 void echo_can_free(echo_can_state_t *ec)
299 fir16_free(&ec->fir_state);
300 fir16_free(&ec->fir_state_bg);
301 for (i = 0; i < 2; i++)
302 kfree(ec->fir_taps16[i]);
306 EXPORT_SYMBOL_GPL(echo_can_free);
307 /*- End of function --------------------------------------------------------*/
309 void echo_can_adaption_mode(echo_can_state_t *ec, int adaption_mode)
311 ec->adaption_mode = adaption_mode;
313 EXPORT_SYMBOL_GPL(echo_can_adaption_mode);
314 /*- End of function --------------------------------------------------------*/
316 void echo_can_flush(echo_can_state_t *ec)
320 ec->Ltxacc = ec->Lrxacc = ec->Lcleanacc = ec->Lclean_bgacc = 0;
321 ec->Ltx = ec->Lrx = ec->Lclean = ec->Lclean_bg = 0;
322 ec->tx_1 = ec->tx_2 = ec->rx_1 = ec->rx_2 = 0;
324 ec->Lbgn = ec->Lbgn_acc = 0;
325 ec->Lbgn_upper = 200;
326 ec->Lbgn_upper_acc = ec->Lbgn_upper << 13;
328 ec->nonupdate_dwell = 0;
330 fir16_flush(&ec->fir_state);
331 fir16_flush(&ec->fir_state_bg);
332 ec->fir_state.curr_pos = ec->taps - 1;
333 ec->fir_state_bg.curr_pos = ec->taps - 1;
334 for (i = 0; i < 2; i++)
335 memset(ec->fir_taps16[i], 0, ec->taps*sizeof(int16_t));
337 ec->curr_pos = ec->taps - 1;
340 EXPORT_SYMBOL_GPL(echo_can_flush);
341 /*- End of function --------------------------------------------------------*/
343 void echo_can_snapshot(echo_can_state_t *ec) {
344 memcpy(ec->snapshot, ec->fir_taps16[0], ec->taps*sizeof(int16_t));
346 EXPORT_SYMBOL_GPL(echo_can_snapshot);
347 /*- End of function --------------------------------------------------------*/
349 /* Dual Path Echo Canceller ------------------------------------------------*/
351 int16_t echo_can_update(echo_can_state_t *ec, int16_t tx, int16_t rx)
357 /* Input scaling was found be required to prevent problems when tx
358 starts clipping. Another possible way to handle this would be the
359 filter coefficent scaling. */
361 ec->tx = tx; ec->rx = rx;
366 Filter DC, 3dB point is 160Hz (I think), note 32 bit precision required
367 otherwise values do not track down to 0. Zero at DC, Pole at (1-Beta)
368 only real axis. Some chip sets (like Si labs) don't need
369 this, but something like a $10 X100P card does. Any DC really slows
372 Note: removes some low frequency from the signal, this reduces
373 the speech quality when listening to samples through headphones
374 but may not be obvious through a telephone handset.
376 Note that the 3dB frequency in radians is approx Beta, e.g. for
377 Beta = 2^(-3) = 0.125, 3dB freq is 0.125 rads = 159Hz.
380 if (ec->adaption_mode & ECHO_CAN_USE_RX_HPF) {
383 /* Make sure the gain of the HPF is 1.0. This can still saturate a little under
384 impulse conditions, and it might roll to 32768 and need clipping on sustained peak
385 level signals. However, the scale of such clipping is small, and the error due to
386 any saturation should not markedly affect the downstream processing. */
389 ec->rx_1 += -(ec->rx_1>>DC_LOG2BETA) + tmp - ec->rx_2;
391 /* hard limit filter to prevent clipping. Note that at this stage
392 rx should be limited to +/- 16383 due to right shift above */
393 tmp1 = ec->rx_1 >> 15;
394 if (tmp1 > 16383) tmp1 = 16383;
395 if (tmp1 < -16383) tmp1 = -16383;
400 /* Block average of power in the filter states. Used for
401 adaption power calculation. */
406 /* efficient "out with the old and in with the new" algorithm so
407 we don't have to recalculate over the whole block of
409 new = (int)tx * (int)tx;
410 old = (int)ec->fir_state.history[ec->fir_state.curr_pos] *
411 (int)ec->fir_state.history[ec->fir_state.curr_pos];
412 ec->Pstates += ((new - old) + (1<<ec->log2taps)) >> ec->log2taps;
413 if (ec->Pstates < 0) ec->Pstates = 0;
416 /* Calculate short term average levels using simple single pole IIRs */
418 ec->Ltxacc += abs(tx) - ec->Ltx;
419 ec->Ltx = (ec->Ltxacc + (1<<4)) >> 5;
420 ec->Lrxacc += abs(rx) - ec->Lrx;
421 ec->Lrx = (ec->Lrxacc + (1<<4)) >> 5;
423 /* Foreground filter ---------------------------------------------------*/
425 ec->fir_state.coeffs = ec->fir_taps16[0];
426 echo_value = fir16(&ec->fir_state, tx);
427 ec->clean = rx - echo_value;
428 ec->Lcleanacc += abs(ec->clean) - ec->Lclean;
429 ec->Lclean = (ec->Lcleanacc + (1<<4)) >> 5;
431 /* Background filter ---------------------------------------------------*/
433 echo_value = fir16(&ec->fir_state_bg, tx);
434 clean_bg = rx - echo_value;
435 ec->Lclean_bgacc += abs(clean_bg) - ec->Lclean_bg;
436 ec->Lclean_bg = (ec->Lclean_bgacc + (1<<4)) >> 5;
438 /* Background Filter adaption -----------------------------------------*/
440 /* Almost always adap bg filter, just simple DT and energy
441 detection to minimise adaption in cases of strong double talk.
442 However this is not critical for the dual path algorithm.
446 if ((ec->nonupdate_dwell == 0)) {
451 f = Beta * clean_bg_rx/P ------ (1)
453 where P is the total power in the filter states.
455 The Boffins have shown that if we obey (1) we converge
456 quickly and avoid instability.
458 The correct factor f must be in Q30, as this is the fixed
459 point format required by the lms_adapt_bg() function,
460 therefore the scaled version of (1) is:
462 (2^30) * f = (2^30) * Beta * clean_bg_rx/P
463 factor = (2^30) * Beta * clean_bg_rx/P ----- (2)
465 We have chosen Beta = 0.25 by experiment, so:
467 factor = (2^30) * (2^-2) * clean_bg_rx/P
470 factor = clean_bg_rx 2 ----- (3)
472 To avoid a divide we approximate log2(P) as top_bit(P),
473 which returns the position of the highest non-zero bit in
474 P. This approximation introduces an error as large as a
475 factor of 2, but the algorithm seems to handle it OK.
477 Come to think of it a divide may not be a big deal on a
478 modern DSP, so its probably worth checking out the cycles
479 for a divide versus a top_bit() implementation.
482 P = MIN_TX_POWER_FOR_ADAPTION + ec->Pstates;
483 logP = top_bit(P) + ec->log2taps;
484 shift = 30 - 2 - logP;
487 lms_adapt_bg(ec, clean_bg, shift);
490 /* very simple DTD to make sure we dont try and adapt with strong
494 if ((ec->Lrx > MIN_RX_POWER_FOR_ADAPTION) && (ec->Lrx > ec->Ltx))
495 ec->nonupdate_dwell = DTD_HANGOVER;
496 if (ec->nonupdate_dwell)
497 ec->nonupdate_dwell--;
499 /* Transfer logic ------------------------------------------------------*/
501 /* These conditions are from the dual path paper [1], I messed with
502 them a bit to improve performance. */
504 if ((ec->adaption_mode & ECHO_CAN_USE_ADAPTION) &&
505 (ec->nonupdate_dwell == 0) &&
506 (8*ec->Lclean_bg < 7*ec->Lclean) /* (ec->Lclean_bg < 0.875*ec->Lclean) */ &&
507 (8*ec->Lclean_bg < ec->Ltx) /* (ec->Lclean_bg < 0.125*ec->Ltx) */ )
509 if (ec->cond_met == 6) {
510 /* BG filter has had better results for 6 consecutive samples */
512 memcpy(ec->fir_taps16[0], ec->fir_taps16[1], ec->taps*sizeof(int16_t));
520 /* Non-Linear Processing ---------------------------------------------------*/
522 ec->clean_nlp = ec->clean;
523 if (ec->adaption_mode & ECHO_CAN_USE_NLP)
525 /* Non-linear processor - a fancy way to say "zap small signals, to avoid
526 residual echo due to (uLaw/ALaw) non-linearity in the channel.". */
528 if ((16*ec->Lclean < ec->Ltx))
530 /* Our e/c has improved echo by at least 24 dB (each factor of 2 is 6dB,
531 so 2*2*2*2=16 is the same as 6+6+6+6=24dB) */
532 if (ec->adaption_mode & ECHO_CAN_USE_CNG)
534 ec->cng_level = ec->Lbgn;
536 /* Very elementary comfort noise generation. Just random
537 numbers rolled off very vaguely Hoth-like. DR: This
538 noise doesn't sound quite right to me - I suspect there
539 are some overlfow issues in the filtering as it's too
540 "crackly". TODO: debug this, maybe just play noise at
541 high level or look at spectrum.
544 ec->cng_rndnum = 1664525U*ec->cng_rndnum + 1013904223U;
545 ec->cng_filter = ((ec->cng_rndnum & 0xFFFF) - 32768 + 5*ec->cng_filter) >> 3;
546 ec->clean_nlp = (ec->cng_filter*ec->cng_level*8) >> 14;
549 else if (ec->adaption_mode & ECHO_CAN_USE_CLIP)
551 /* This sounds much better than CNG */
552 if (ec->clean_nlp > ec->Lbgn)
553 ec->clean_nlp = ec->Lbgn;
554 if (ec->clean_nlp < -ec->Lbgn)
555 ec->clean_nlp = -ec->Lbgn;
559 /* just mute the residual, doesn't sound very good, used mainly
565 /* Background noise estimator. I tried a few algorithms
566 here without much luck. This very simple one seems to
567 work best, we just average the level using a slow (1 sec
568 time const) filter if the current level is less than a
569 (experimentally derived) constant. This means we dont
570 include high level signals like near end speech. When
571 combined with CNG or especially CLIP seems to work OK.
573 if (ec->Lclean < 40) {
574 ec->Lbgn_acc += abs(ec->clean) - ec->Lbgn;
575 ec->Lbgn = (ec->Lbgn_acc + (1<<11)) >> 12;
580 /* Roll around the taps buffer */
581 if (ec->curr_pos <= 0)
582 ec->curr_pos = ec->taps;
585 if (ec->adaption_mode & ECHO_CAN_DISABLE)
588 /* Output scaled back up again to match input scaling */
590 return (int16_t) ec->clean_nlp << 1;
592 EXPORT_SYMBOL_GPL(echo_can_update);
593 /*- End of function --------------------------------------------------------*/
595 /* This function is seperated from the echo canceller is it is usually called
596 as part of the tx process. See rx HP (DC blocking) filter above, it's
599 Some soft phones send speech signals with a lot of low frequency
600 energy, e.g. down to 20Hz. This can make the hybrid non-linear
601 which causes the echo canceller to fall over. This filter can help
602 by removing any low frequency before it gets to the tx port of the
605 It can also help by removing and DC in the tx signal. DC is bad
608 This is one of the classic DC removal filters, adjusted to provide sufficient
609 bass rolloff to meet the above requirement to protect hybrids from things that
610 upset them. The difference between successive samples produces a lousy HPF, and
611 then a suitably placed pole flattens things out. The final result is a nicely
612 rolled off bass end. The filtering is implemented with extended fractional
613 precision, which noise shapes things, giving very clean DC removal.
616 int16_t echo_can_hpf_tx(echo_can_state_t *ec, int16_t tx) {
619 if (ec->adaption_mode & ECHO_CAN_USE_TX_HPF) {
622 /* Make sure the gain of the HPF is 1.0. The first can still saturate a little under
623 impulse conditions, and it might roll to 32768 and need clipping on sustained peak
624 level signals. However, the scale of such clipping is small, and the error due to
625 any saturation should not markedly affect the downstream processing. */
628 ec->tx_1 += -(ec->tx_1>>DC_LOG2BETA) + tmp - ec->tx_2;
629 tmp1 = ec->tx_1 >> 15;
630 if (tmp1 > 32767) tmp1 = 32767;
631 if (tmp1 < -32767) tmp1 = -32767;
638 EXPORT_SYMBOL_GPL(echo_can_hpf_tx);
640 MODULE_LICENSE("GPL");
641 MODULE_AUTHOR("David Rowe");
642 MODULE_DESCRIPTION("Open Source Line Echo Canceller");
643 MODULE_VERSION("0.3.0");