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.
34 /* Implementation Notes
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.
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
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.
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.
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.
73 Steve also has some nice notes on echo cancellers in echo.h
77 [1] Ochiai, Areseki, and Ogihara, "Echo Canceller with Two Echo
78 Path Models", IEEE Transactions on communications, COM-25,
81 http://www.rowetel.com/images/echo/dual_path_paper.pdf
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
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
92 [4] The source code http://svn.rowetel.com/software/oslec/
94 [5] A nice reference on LMS filters:
95 http://en.wikipedia.org/wiki/Least_mean_squares_filter
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.
105 #include <linux/kernel.h>
106 #include <linux/module.h>
107 #include <linux/slab.h>
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 */
116 /* adapting coeffs using the traditional stochastic descent (N)LMS algorithm */
119 static inline void lms_adapt_bg(struct oslec_state
*ec
, int clean
, int shift
)
130 factor
= clean
<< shift
;
132 factor
= clean
>> -shift
;
134 /* Update the FIR taps */
136 offset2
= ec
->curr_pos
;
137 offset1
= ec
->taps
- offset2
;
138 phist
= &ec
->fir_state_bg
.history
[offset2
];
140 /* st: and en: help us locate the assembler in echo.s */
144 for (i
= 0; i
< n
; i
++) {
145 exp
= *phist
++ * factor
;
146 ec
->fir_taps16
[1][i
] += (int16_t) ((exp
+ (1 << 14)) >> 15);
150 /* Note the asm for the inner loop above generated by Blackfin gcc
151 4.1.1 is pretty good (note even parallel instructions used):
162 A block based update algorithm would be much faster but the
163 above can't be improved on much. Every instruction saved in
164 the loop above is 2 MIPs/ch! The for loop above is where the
165 Blackfin spends most of it's time - about 17 MIPs/ch measured
166 with speedtest.c with 256 taps (32ms). Write-back and
167 Write-through cache gave about the same performance.
172 IDEAS for further optimisation of lms_adapt_bg():
174 1/ The rounding is quite costly. Could we keep as 32 bit coeffs
175 then make filter pluck the MS 16-bits of the coeffs when filtering?
176 However this would lower potential optimisation of filter, as I
177 think the dual-MAC architecture requires packed 16 bit coeffs.
179 2/ Block based update would be more efficient, as per comments above,
180 could use dual MAC architecture.
182 3/ Look for same sample Blackfin LMS code, see if we can get dual-MAC
185 4/ Execute the whole e/c in a block of say 20ms rather than sample
186 by sample. Processing a few samples every ms is inefficient.
190 static inline void lms_adapt_bg(struct oslec_state
*ec
, int clean
, int shift
)
200 factor
= clean
<< shift
;
202 factor
= clean
>> -shift
;
204 /* Update the FIR taps */
206 offset2
= ec
->curr_pos
;
207 offset1
= ec
->taps
- offset2
;
209 for (i
= ec
->taps
- 1; i
>= offset1
; i
--) {
210 exp
= (ec
->fir_state_bg
.history
[i
- offset1
] * factor
);
211 ec
->fir_taps16
[1][i
] += (int16_t) ((exp
+ (1 << 14)) >> 15);
213 for (; i
>= 0; i
--) {
214 exp
= (ec
->fir_state_bg
.history
[i
+ offset2
] * factor
);
215 ec
->fir_taps16
[1][i
] += (int16_t) ((exp
+ (1 << 14)) >> 15);
220 static inline int top_bit(unsigned int bits
)
225 return (int)fls((int32_t) bits
) - 1;
228 struct oslec_state
*oslec_create(int len
, int adaption_mode
)
230 struct oslec_state
*ec
;
232 const int16_t *history
;
234 ec
= kzalloc(sizeof(*ec
), GFP_KERNEL
);
239 ec
->log2taps
= top_bit(len
);
240 ec
->curr_pos
= ec
->taps
- 1;
243 kcalloc(ec
->taps
, sizeof(int16_t), GFP_KERNEL
);
244 if (!ec
->fir_taps16
[0])
248 kcalloc(ec
->taps
, sizeof(int16_t), GFP_KERNEL
);
249 if (!ec
->fir_taps16
[1])
252 history
= fir16_create(&ec
->fir_state
, ec
->fir_taps16
[0], ec
->taps
);
255 history
= fir16_create(&ec
->fir_state_bg
, ec
->fir_taps16
[1], ec
->taps
);
259 for (i
= 0; i
< 5; i
++)
260 ec
->xvtx
[i
] = ec
->yvtx
[i
] = ec
->xvrx
[i
] = ec
->yvrx
[i
] = 0;
262 ec
->cng_level
= 1000;
263 oslec_adaption_mode(ec
, adaption_mode
);
265 ec
->snapshot
= kcalloc(ec
->taps
, sizeof(int16_t), GFP_KERNEL
);
271 ec
->Ltxacc
= ec
->Lrxacc
= ec
->Lcleanacc
= ec
->Lclean_bgacc
= 0;
272 ec
->Ltx
= ec
->Lrx
= ec
->Lclean
= ec
->Lclean_bg
= 0;
273 ec
->tx_1
= ec
->tx_2
= ec
->rx_1
= ec
->rx_2
= 0;
274 ec
->Lbgn
= ec
->Lbgn_acc
= 0;
275 ec
->Lbgn_upper
= 200;
276 ec
->Lbgn_upper_acc
= ec
->Lbgn_upper
<< 13;
281 fir16_free(&ec
->fir_state_bg
);
283 fir16_free(&ec
->fir_state
);
285 kfree(ec
->fir_taps16
[1]);
287 kfree(ec
->fir_taps16
[0]);
292 EXPORT_SYMBOL_GPL(oslec_create
);
294 void oslec_free(struct oslec_state
*ec
)
298 fir16_free(&ec
->fir_state
);
299 fir16_free(&ec
->fir_state_bg
);
300 for (i
= 0; i
< 2; i
++)
301 kfree(ec
->fir_taps16
[i
]);
305 EXPORT_SYMBOL_GPL(oslec_free
);
307 void oslec_adaption_mode(struct oslec_state
*ec
, int adaption_mode
)
309 ec
->adaption_mode
= adaption_mode
;
311 EXPORT_SYMBOL_GPL(oslec_adaption_mode
);
313 void oslec_flush(struct oslec_state
*ec
)
317 ec
->Ltxacc
= ec
->Lrxacc
= ec
->Lcleanacc
= ec
->Lclean_bgacc
= 0;
318 ec
->Ltx
= ec
->Lrx
= ec
->Lclean
= ec
->Lclean_bg
= 0;
319 ec
->tx_1
= ec
->tx_2
= ec
->rx_1
= ec
->rx_2
= 0;
321 ec
->Lbgn
= ec
->Lbgn_acc
= 0;
322 ec
->Lbgn_upper
= 200;
323 ec
->Lbgn_upper_acc
= ec
->Lbgn_upper
<< 13;
325 ec
->nonupdate_dwell
= 0;
327 fir16_flush(&ec
->fir_state
);
328 fir16_flush(&ec
->fir_state_bg
);
329 ec
->fir_state
.curr_pos
= ec
->taps
- 1;
330 ec
->fir_state_bg
.curr_pos
= ec
->taps
- 1;
331 for (i
= 0; i
< 2; i
++)
332 memset(ec
->fir_taps16
[i
], 0, ec
->taps
* sizeof(int16_t));
334 ec
->curr_pos
= ec
->taps
- 1;
337 EXPORT_SYMBOL_GPL(oslec_flush
);
339 void oslec_snapshot(struct oslec_state
*ec
)
341 memcpy(ec
->snapshot
, ec
->fir_taps16
[0], ec
->taps
* sizeof(int16_t));
343 EXPORT_SYMBOL_GPL(oslec_snapshot
);
345 /* Dual Path Echo Canceller */
347 int16_t oslec_update(struct oslec_state
*ec
, int16_t tx
, int16_t rx
)
355 * Input scaling was found be required to prevent problems when tx
356 * starts clipping. Another possible way to handle this would be the
357 * filter coefficent scaling.
366 * Filter DC, 3dB point is 160Hz (I think), note 32 bit precision
367 * required otherwise values do not track down to 0. Zero at DC, Pole
368 * at (1-Beta) on real axis. Some chip sets (like Si labs) don't
369 * need this, but something like a $10 X100P card does. Any DC really
370 * slows down convergence.
372 * Note: removes some low frequency from the signal, this reduces the
373 * speech quality when listening to samples through headphones but may
374 * not be obvious through a telephone handset.
376 * Note that the 3dB frequency in radians is approx Beta, e.g. for Beta
377 * = 2^(-3) = 0.125, 3dB freq is 0.125 rads = 159Hz.
380 if (ec
->adaption_mode
& ECHO_CAN_USE_RX_HPF
) {
384 * Make sure the gain of the HPF is 1.0. This can still
385 * saturate a little under impulse conditions, and it might
386 * roll to 32768 and need clipping on sustained peak level
387 * signals. However, the scale of such clipping is small, and
388 * the error due to any saturation should not markedly affect
389 * the downstream processing.
393 ec
->rx_1
+= -(ec
->rx_1
>> DC_LOG2BETA
) + tmp
- ec
->rx_2
;
396 * hard limit filter to prevent clipping. Note that at this
397 * stage rx should be limited to +/- 16383 due to right shift
400 tmp1
= ec
->rx_1
>> 15;
409 /* Block average of power in the filter states. Used for
410 adaption power calculation. */
415 /* efficient "out with the old and in with the new" algorithm so
416 we don't have to recalculate over the whole block of
418 new = (int)tx
* (int)tx
;
419 old
= (int)ec
->fir_state
.history
[ec
->fir_state
.curr_pos
] *
420 (int)ec
->fir_state
.history
[ec
->fir_state
.curr_pos
];
422 ((new - old
) + (1 << (ec
->log2taps
- 1))) >> ec
->log2taps
;
427 /* Calculate short term average levels using simple single pole IIRs */
429 ec
->Ltxacc
+= abs(tx
) - ec
->Ltx
;
430 ec
->Ltx
= (ec
->Ltxacc
+ (1 << 4)) >> 5;
431 ec
->Lrxacc
+= abs(rx
) - ec
->Lrx
;
432 ec
->Lrx
= (ec
->Lrxacc
+ (1 << 4)) >> 5;
434 /* Foreground filter */
436 ec
->fir_state
.coeffs
= ec
->fir_taps16
[0];
437 echo_value
= fir16(&ec
->fir_state
, tx
);
438 ec
->clean
= rx
- echo_value
;
439 ec
->Lcleanacc
+= abs(ec
->clean
) - ec
->Lclean
;
440 ec
->Lclean
= (ec
->Lcleanacc
+ (1 << 4)) >> 5;
442 /* Background filter */
444 echo_value
= fir16(&ec
->fir_state_bg
, tx
);
445 clean_bg
= rx
- echo_value
;
446 ec
->Lclean_bgacc
+= abs(clean_bg
) - ec
->Lclean_bg
;
447 ec
->Lclean_bg
= (ec
->Lclean_bgacc
+ (1 << 4)) >> 5;
449 /* Background Filter adaption */
451 /* Almost always adap bg filter, just simple DT and energy
452 detection to minimise adaption in cases of strong double talk.
453 However this is not critical for the dual path algorithm.
457 if ((ec
->nonupdate_dwell
== 0)) {
462 f = Beta * clean_bg_rx/P ------ (1)
464 where P is the total power in the filter states.
466 The Boffins have shown that if we obey (1) we converge
467 quickly and avoid instability.
469 The correct factor f must be in Q30, as this is the fixed
470 point format required by the lms_adapt_bg() function,
471 therefore the scaled version of (1) is:
473 (2^30) * f = (2^30) * Beta * clean_bg_rx/P
474 factor = (2^30) * Beta * clean_bg_rx/P ----- (2)
476 We have chosen Beta = 0.25 by experiment, so:
478 factor = (2^30) * (2^-2) * clean_bg_rx/P
481 factor = clean_bg_rx 2 ----- (3)
483 To avoid a divide we approximate log2(P) as top_bit(P),
484 which returns the position of the highest non-zero bit in
485 P. This approximation introduces an error as large as a
486 factor of 2, but the algorithm seems to handle it OK.
488 Come to think of it a divide may not be a big deal on a
489 modern DSP, so its probably worth checking out the cycles
490 for a divide versus a top_bit() implementation.
493 P
= MIN_TX_POWER_FOR_ADAPTION
+ ec
->Pstates
;
494 logP
= top_bit(P
) + ec
->log2taps
;
495 shift
= 30 - 2 - logP
;
498 lms_adapt_bg(ec
, clean_bg
, shift
);
501 /* very simple DTD to make sure we dont try and adapt with strong
505 if ((ec
->Lrx
> MIN_RX_POWER_FOR_ADAPTION
) && (ec
->Lrx
> ec
->Ltx
))
506 ec
->nonupdate_dwell
= DTD_HANGOVER
;
507 if (ec
->nonupdate_dwell
)
508 ec
->nonupdate_dwell
--;
512 /* These conditions are from the dual path paper [1], I messed with
513 them a bit to improve performance. */
515 if ((ec
->adaption_mode
& ECHO_CAN_USE_ADAPTION
) &&
516 (ec
->nonupdate_dwell
== 0) &&
517 /* (ec->Lclean_bg < 0.875*ec->Lclean) */
518 (8 * ec
->Lclean_bg
< 7 * ec
->Lclean
) &&
519 /* (ec->Lclean_bg < 0.125*ec->Ltx) */
520 (8 * ec
->Lclean_bg
< ec
->Ltx
)) {
521 if (ec
->cond_met
== 6) {
523 * BG filter has had better results for 6 consecutive
527 memcpy(ec
->fir_taps16
[0], ec
->fir_taps16
[1],
528 ec
->taps
* sizeof(int16_t));
534 /* Non-Linear Processing */
536 ec
->clean_nlp
= ec
->clean
;
537 if (ec
->adaption_mode
& ECHO_CAN_USE_NLP
) {
539 * Non-linear processor - a fancy way to say "zap small
540 * signals, to avoid residual echo due to (uLaw/ALaw)
541 * non-linearity in the channel.".
544 if ((16 * ec
->Lclean
< ec
->Ltx
)) {
546 * Our e/c has improved echo by at least 24 dB (each
547 * factor of 2 is 6dB, so 2*2*2*2=16 is the same as
550 if (ec
->adaption_mode
& ECHO_CAN_USE_CNG
) {
551 ec
->cng_level
= ec
->Lbgn
;
554 * Very elementary comfort noise generation.
555 * Just random numbers rolled off very vaguely
556 * Hoth-like. DR: This noise doesn't sound
557 * quite right to me - I suspect there are some
558 * overflow issues in the filtering as it's too
560 * TODO: debug this, maybe just play noise at
561 * high level or look at spectrum.
565 1664525U * ec
->cng_rndnum
+ 1013904223U;
567 ((ec
->cng_rndnum
& 0xFFFF) - 32768 +
568 5 * ec
->cng_filter
) >> 3;
570 (ec
->cng_filter
* ec
->cng_level
* 8) >> 14;
572 } else if (ec
->adaption_mode
& ECHO_CAN_USE_CLIP
) {
573 /* This sounds much better than CNG */
574 if (ec
->clean_nlp
> ec
->Lbgn
)
575 ec
->clean_nlp
= ec
->Lbgn
;
576 if (ec
->clean_nlp
< -ec
->Lbgn
)
577 ec
->clean_nlp
= -ec
->Lbgn
;
580 * just mute the residual, doesn't sound very
581 * good, used mainly in G168 tests
587 * Background noise estimator. I tried a few
588 * algorithms here without much luck. This very simple
589 * one seems to work best, we just average the level
590 * using a slow (1 sec time const) filter if the
591 * current level is less than a (experimentally
592 * derived) constant. This means we dont include high
593 * level signals like near end speech. When combined
594 * with CNG or especially CLIP seems to work OK.
596 if (ec
->Lclean
< 40) {
597 ec
->Lbgn_acc
+= abs(ec
->clean
) - ec
->Lbgn
;
598 ec
->Lbgn
= (ec
->Lbgn_acc
+ (1 << 11)) >> 12;
603 /* Roll around the taps buffer */
604 if (ec
->curr_pos
<= 0)
605 ec
->curr_pos
= ec
->taps
;
608 if (ec
->adaption_mode
& ECHO_CAN_DISABLE
)
611 /* Output scaled back up again to match input scaling */
613 return (int16_t) ec
->clean_nlp
<< 1;
615 EXPORT_SYMBOL_GPL(oslec_update
);
617 /* This function is separated from the echo canceller is it is usually called
618 as part of the tx process. See rx HP (DC blocking) filter above, it's
621 Some soft phones send speech signals with a lot of low frequency
622 energy, e.g. down to 20Hz. This can make the hybrid non-linear
623 which causes the echo canceller to fall over. This filter can help
624 by removing any low frequency before it gets to the tx port of the
627 It can also help by removing and DC in the tx signal. DC is bad
630 This is one of the classic DC removal filters, adjusted to provide
631 sufficient bass rolloff to meet the above requirement to protect hybrids
632 from things that upset them. The difference between successive samples
633 produces a lousy HPF, and then a suitably placed pole flattens things out.
634 The final result is a nicely rolled off bass end. The filtering is
635 implemented with extended fractional precision, which noise shapes things,
636 giving very clean DC removal.
639 int16_t oslec_hpf_tx(struct oslec_state
*ec
, int16_t tx
)
644 if (ec
->adaption_mode
& ECHO_CAN_USE_TX_HPF
) {
648 * Make sure the gain of the HPF is 1.0. The first can still
649 * saturate a little under impulse conditions, and it might
650 * roll to 32768 and need clipping on sustained peak level
651 * signals. However, the scale of such clipping is small, and
652 * the error due to any saturation should not markedly affect
653 * the downstream processing.
657 ec
->tx_1
+= -(ec
->tx_1
>> DC_LOG2BETA
) + tmp
- ec
->tx_2
;
658 tmp1
= ec
->tx_1
>> 15;
669 EXPORT_SYMBOL_GPL(oslec_hpf_tx
);
671 MODULE_LICENSE("GPL");
672 MODULE_AUTHOR("David Rowe");
673 MODULE_DESCRIPTION("Open Source Line Echo Canceller");
674 MODULE_VERSION("0.3.0");