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 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 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 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 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 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 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 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 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 == 0)) { 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 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