drivers/misc/echo/echo.c

   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) {
 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");