drivers/staging/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  * $Id: echo.c,v 1.20 2006/12/01 18:00:48 steveu Exp $
  32  */
  33
  34 /*! \file */
  35
  36 /* Implementation Notes
  37    David Rowe
  38    April 2007
  39
  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.
  47
  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
  52    serious work into it.
  53
  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.
  59
  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.
  65
  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.
  74
  75    Steve also has some nice notes on echo cancellers in echo.h
  76
  77
  78    References:
  79
  80    [1] Ochiai, Areseki, and Ogihara, "Echo Canceller with Two Echo
  81        Path Models", IEEE Transactions on communications, COM-25,
  82        No. 6, June
  83        1977.
  84        http://www.rowetel.com/images/echo/dual_path_paper.pdf
  85
  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
  91
  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
  94
  95    [4] The source code http://svn.rowetel.com/software/oslec/
  96
  97    [5] A nice reference on LMS filters:
  98          http://en.wikipedia.org/wiki/Least_mean_squares_filter
  99
 100    Credits:
 101
 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.
 106 */
 107
 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
 113 #include "bit_operations.h"
 114 #include "echo.h"
 115
 116 #define MIN_TX_POWER_FOR_ADAPTION   64
 117 #define MIN_RX_POWER_FOR_ADAPTION   64
 118 #define DTD_HANGOVER               600     /* 600 samples, or 75ms     */
 119 #define DC_LOG2BETA                  3     /* log2() of DC filter Beta */
 120
 121 /*-----------------------------------------------------------------------*\
 122                                FUNCTIONS
 123 \*-----------------------------------------------------------------------*/
 124
 125 /* adapting coeffs using the traditional stochastic descent (N)LMS algorithm */
 126
 127
 128 #ifdef __bfin__
 129 static void __inline__ lms_adapt_bg(struct oslec_state *ec, int clean, int shift)
 130 {
 131     int i, j;
 132     int offset1;
 133     int offset2;
 134     int factor;
 135     int exp;
 136     int16_t *phist;
 137     int n;
 138
 139     if (shift > 0)
 140         factor = clean << shift;
 141     else
 142         factor = clean >> -shift;
 143
 144     /* Update the FIR taps */
 145
 146     offset2 = ec->curr_pos;
 147     offset1 = ec->taps - offset2;
 148     phist = &ec->fir_state_bg.history[offset2];
 149
 150     /* st: and en: help us locate the assembler in echo.s */
 151
 152     //asm("st:");
 153     n = ec->taps;
 154     for (i = 0, j = offset2;  i < n;  i++, j++)
 155     {
 156        exp = *phist++ * factor;
 157        ec->fir_taps16[1][i] += (int16_t) ((exp+(1<<14)) >> 15);
 158     }
 159     //asm("en:");
 160
 161     /* Note the asm for the inner loop above generated by Blackfin gcc
 162        4.1.1 is pretty good (note even parallel instructions used):
 163
 164         R0 = W [P0++] (X);
 165         R0 *= R2;
 166         R0 = R0 + R3 (NS) ||
 167         R1 = W [P1] (X) ||
 168         nop;
 169         R0 >>>= 15;
 170         R0 = R0 + R1;
 171         W [P1++] = R0;
 172
 173         A block based update algorithm would be much faster but the
 174         above can't be improved on much.  Every instruction saved in
 175         the loop above is 2 MIPs/ch!  The for loop above is where the
 176         Blackfin spends most of it's time - about 17 MIPs/ch measured
 177         with speedtest.c with 256 taps (32ms).  Write-back and
 178         Write-through cache gave about the same performance.
 179     */
 180 }
 181
 182 /*
 183    IDEAS for further optimisation of lms_adapt_bg():
 184
 185    1/ The rounding is quite costly.  Could we keep as 32 bit coeffs
 186    then make filter pluck the MS 16-bits of the coeffs when filtering?
 187    However this would lower potential optimisation of filter, as I
 188    think the dual-MAC architecture requires packed 16 bit coeffs.
 189
 190    2/ Block based update would be more efficient, as per comments above,
 191    could use dual MAC architecture.
 192
 193    3/ Look for same sample Blackfin LMS code, see if we can get dual-MAC
 194    packing.
 195
 196    4/ Execute the whole e/c in a block of say 20ms rather than sample
 197    by sample.  Processing a few samples every ms is inefficient.
 198 */
 199
 200 #else
 201 static __inline__ void lms_adapt_bg(struct oslec_state *ec, int clean, int shift)
 202 {
 203     int i;
 204
 205     int offset1;
 206     int offset2;
 207     int factor;
 208     int exp;
 209
 210     if (shift > 0)
 211         factor = clean << shift;
 212     else
 213         factor = clean >> -shift;
 214
 215     /* Update the FIR taps */
 216
 217     offset2 = ec->curr_pos;
 218     offset1 = ec->taps - offset2;
 219
 220     for (i = ec->taps - 1;  i >= offset1;  i--)
 221     {
 222        exp = (ec->fir_state_bg.history[i - offset1]*factor);
 223        ec->fir_taps16[1][i] += (int16_t) ((exp+(1<<14)) >> 15);
 224     }
 225     for (  ;  i >= 0;  i--)
 226     {
 227        exp = (ec->fir_state_bg.history[i + offset2]*factor);
 228        ec->fir_taps16[1][i] += (int16_t) ((exp+(1<<14)) >> 15);
 229     }
 230 }
 231 #endif
 232
 233 /*- End of function --------------------------------------------------------*/
 234
 235 struct oslec_state *oslec_create(int len, int adaption_mode)
 236 {
 237     struct oslec_state *ec;
 238     int i;
 239
 240     ec = kzalloc(sizeof(*ec), GFP_KERNEL);
 241     if (!ec)
 242         return NULL;
 243
 244     ec->taps = len;
 245     ec->log2taps = top_bit(len);
 246     ec->curr_pos = ec->taps - 1;
 247
 248     for (i = 0; i < 2; i++) {
 249         ec->fir_taps16[i] = kcalloc(ec->taps, sizeof(int16_t), GFP_KERNEL);
 250         if (!ec->fir_taps16[i])
 251             goto error_oom;
 252     }
 253
 254     fir16_create(&ec->fir_state,
 255                  ec->fir_taps16[0],
 256                  ec->taps);
 257     fir16_create(&ec->fir_state_bg,
 258                  ec->fir_taps16[1],
 259                  ec->taps);
 260
 261     for(i=0; i<5; i++) {
 262       ec->xvtx[i] = ec->yvtx[i] = ec->xvrx[i] = ec->yvrx[i] = 0;
 263     }
 264
 265     ec->cng_level = 1000;
 266     oslec_adaption_mode(ec, adaption_mode);
 267
 268     ec->snapshot = kcalloc(ec->taps, sizeof(int16_t), GFP_KERNEL);
 269     if (!ec->snapshot)
 270         goto error_oom;
 271
 272     ec->cond_met = 0;
 273     ec->Pstates = 0;
 274     ec->Ltxacc = ec->Lrxacc = ec->Lcleanacc = ec->Lclean_bgacc = 0;
 275     ec->Ltx = ec->Lrx = ec->Lclean = ec->Lclean_bg = 0;
 276     ec->tx_1 = ec->tx_2 = ec->rx_1 = ec->rx_2 = 0;
 277     ec->Lbgn = ec->Lbgn_acc = 0;
 278     ec->Lbgn_upper = 200;
 279     ec->Lbgn_upper_acc = ec->Lbgn_upper << 13;
 280
 281     return  ec;
 282
 283 error_oom:
 284     for (i = 0; i < 2; i++)
 285         kfree(ec->fir_taps16[i]);
 286
 287     kfree(ec);
 288     return NULL;
 289 }
 290 EXPORT_SYMBOL_GPL(oslec_create);
 291 /*- End of function --------------------------------------------------------*/
 292
 293 void oslec_free(struct oslec_state *ec)
 294 {
 295         int i;
 296
 297         fir16_free(&ec->fir_state);
 298         fir16_free(&ec->fir_state_bg);
 299         for (i = 0;  i < 2;  i++)
 300                 kfree(ec->fir_taps16[i]);
 301         kfree(ec->snapshot);
 302         kfree(ec);
 303 }
 304 EXPORT_SYMBOL_GPL(oslec_free);
 305 /*- End of function --------------------------------------------------------*/
 306
 307 void oslec_adaption_mode(struct oslec_state *ec, int adaption_mode)
 308 {
 309     ec->adaption_mode = adaption_mode;
 310 }
 311 EXPORT_SYMBOL_GPL(oslec_adaption_mode);
 312 /*- End of function --------------------------------------------------------*/
 313
 314 void oslec_flush(struct oslec_state *ec)
 315 {
 316     int i;
 317
 318     ec->Ltxacc = ec->Lrxacc = ec->Lcleanacc = ec->Lclean_bgacc = 0;
 319     ec->Ltx = ec->Lrx = ec->Lclean = ec->Lclean_bg = 0;
 320     ec->tx_1 = ec->tx_2 = ec->rx_1 = ec->rx_2 = 0;
 321
 322     ec->Lbgn = ec->Lbgn_acc = 0;
 323     ec->Lbgn_upper = 200;
 324     ec->Lbgn_upper_acc = ec->Lbgn_upper << 13;
 325
 326     ec->nonupdate_dwell = 0;
 327
 328     fir16_flush(&ec->fir_state);
 329     fir16_flush(&ec->fir_state_bg);
 330     ec->fir_state.curr_pos = ec->taps - 1;
 331     ec->fir_state_bg.curr_pos = ec->taps - 1;
 332     for (i = 0;  i < 2;  i++)
 333         memset(ec->fir_taps16[i], 0, ec->taps*sizeof(int16_t));
 334
 335     ec->curr_pos = ec->taps - 1;
 336     ec->Pstates = 0;
 337 }
 338 EXPORT_SYMBOL_GPL(oslec_flush);
 339 /*- End of function --------------------------------------------------------*/
 340
 341 void oslec_snapshot(struct oslec_state *ec) {
 342     memcpy(ec->snapshot, ec->fir_taps16[0], ec->taps*sizeof(int16_t));
 343 }
 344 EXPORT_SYMBOL_GPL(oslec_snapshot);
 345 /*- End of function --------------------------------------------------------*/
 346
 347 /* Dual Path Echo Canceller ------------------------------------------------*/
 348
 349 int16_t oslec_update(struct oslec_state *ec, int16_t tx, int16_t rx)
 350 {
 351     int32_t echo_value;
 352     int clean_bg;
 353     int tmp, tmp1;
 354
 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. */
 358
 359     ec->tx = tx; ec->rx = rx;
 360     tx >>=1;
 361     rx >>=1;
 362
 363     /*
 364        Filter DC, 3dB point is 160Hz (I think), note 32 bit precision required
 365        otherwise values do not track down to 0. Zero at DC, Pole at (1-Beta)
 366        only real axis.  Some chip sets (like Si labs) don't need
 367        this, but something like a $10 X100P card does.  Any DC really slows
 368        down convergence.
 369
 370        Note: removes some low frequency from the signal, this reduces
 371        the speech quality when listening to samples through headphones
 372        but may not be obvious through a telephone handset.
 373
 374        Note that the 3dB frequency in radians is approx Beta, e.g. for
 375        Beta = 2^(-3) = 0.125, 3dB freq is 0.125 rads = 159Hz.
 376     */
 377
 378     if (ec->adaption_mode & ECHO_CAN_USE_RX_HPF) {
 379       tmp = rx << 15;
 380 #if 1
 381         /* Make sure the gain of the HPF is 1.0. This can still saturate a little under
 382            impulse conditions, and it might roll to 32768 and need clipping on sustained peak
 383            level signals. However, the scale of such clipping is small, and the error due to
 384            any saturation should not markedly affect the downstream processing. */
 385         tmp -= (tmp >> 4);
 386 #endif
 387       ec->rx_1 += -(ec->rx_1>>DC_LOG2BETA) + tmp - ec->rx_2;
 388
 389       /* hard limit filter to prevent clipping.  Note that at this stage
 390          rx should be limited to +/- 16383 due to right shift above */
 391       tmp1 = ec->rx_1 >> 15;
 392       if (tmp1 > 16383) tmp1 = 16383;
 393       if (tmp1 < -16383) tmp1 = -16383;
 394       rx = tmp1;
 395       ec->rx_2 = tmp;
 396     }
 397
 398     /* Block average of power in the filter states.  Used for
 399        adaption power calculation. */
 400
 401     {
 402         int new, old;
 403
 404         /* efficient "out with the old and in with the new" algorithm so
 405            we don't have to recalculate over the whole block of
 406            samples. */
 407         new = (int)tx * (int)tx;
 408         old = (int)ec->fir_state.history[ec->fir_state.curr_pos] *
 409               (int)ec->fir_state.history[ec->fir_state.curr_pos];
 410         ec->Pstates += ((new - old) + (1<<ec->log2taps)) >> ec->log2taps;
 411         if (ec->Pstates < 0) ec->Pstates = 0;
 412     }
 413
 414     /* Calculate short term average levels using simple single pole IIRs */
 415
 416     ec->Ltxacc += abs(tx) - ec->Ltx;
 417     ec->Ltx = (ec->Ltxacc + (1<<4)) >> 5;
 418     ec->Lrxacc += abs(rx) - ec->Lrx;
 419     ec->Lrx = (ec->Lrxacc + (1<<4)) >> 5;
 420
 421     /* Foreground filter ---------------------------------------------------*/
 422
 423     ec->fir_state.coeffs = ec->fir_taps16[0];
 424     echo_value = fir16(&ec->fir_state, tx);
 425     ec->clean = rx - echo_value;
 426     ec->Lcleanacc += abs(ec->clean) - ec->Lclean;
 427     ec->Lclean = (ec->Lcleanacc + (1<<4)) >> 5;
 428
 429     /* Background filter ---------------------------------------------------*/
 430
 431     echo_value = fir16(&ec->fir_state_bg, tx);
 432     clean_bg = rx - echo_value;
 433     ec->Lclean_bgacc += abs(clean_bg) - ec->Lclean_bg;
 434     ec->Lclean_bg = (ec->Lclean_bgacc + (1<<4)) >> 5;
 435
 436     /* Background Filter adaption -----------------------------------------*/
 437
 438     /* Almost always adap bg filter, just simple DT and energy
 439        detection to minimise adaption in cases of strong double talk.
 440        However this is not critical for the dual path algorithm.
 441     */
 442     ec->factor = 0;
 443     ec->shift = 0;
 444     if ((ec->nonupdate_dwell == 0)) {
 445         int   P, logP, shift;
 446
 447         /* Determine:
 448
 449            f = Beta * clean_bg_rx/P ------ (1)
 450
 451            where P is the total power in the filter states.
 452
 453            The Boffins have shown that if we obey (1) we converge
 454            quickly and avoid instability.
 455
 456            The correct factor f must be in Q30, as this is the fixed
 457            point format required by the lms_adapt_bg() function,
 458            therefore the scaled version of (1) is:
 459
 460            (2^30) * f  = (2^30) * Beta * clean_bg_rx/P
 461                factor  = (2^30) * Beta * clean_bg_rx/P         ----- (2)
 462
 463            We have chosen Beta = 0.25 by experiment, so:
 464
 465                factor  = (2^30) * (2^-2) * clean_bg_rx/P
 466
 467                                        (30 - 2 - log2(P))
 468                factor  = clean_bg_rx 2                         ----- (3)
 469
 470            To avoid a divide we approximate log2(P) as top_bit(P),
 471            which returns the position of the highest non-zero bit in
 472            P.  This approximation introduces an error as large as a
 473            factor of 2, but the algorithm seems to handle it OK.
 474
 475            Come to think of it a divide may not be a big deal on a
 476            modern DSP, so its probably worth checking out the cycles
 477            for a divide versus a top_bit() implementation.
 478         */
 479
 480         P = MIN_TX_POWER_FOR_ADAPTION + ec->Pstates;
 481         logP = top_bit(P) + ec->log2taps;
 482         shift = 30 - 2 - logP;
 483         ec->shift = shift;
 484
 485         lms_adapt_bg(ec, clean_bg, shift);
 486     }
 487
 488     /* very simple DTD to make sure we dont try and adapt with strong
 489        near end speech */
 490
 491     ec->adapt = 0;
 492     if ((ec->Lrx > MIN_RX_POWER_FOR_ADAPTION) && (ec->Lrx > ec->Ltx))
 493         ec->nonupdate_dwell = DTD_HANGOVER;
 494     if (ec->nonupdate_dwell)
 495         ec->nonupdate_dwell--;
 496
 497     /* Transfer logic ------------------------------------------------------*/
 498
 499     /* These conditions are from the dual path paper [1], I messed with
 500        them a bit to improve performance. */
 501
 502     if ((ec->adaption_mode & ECHO_CAN_USE_ADAPTION) &&
 503         (ec->nonupdate_dwell == 0) &&
 504         (8*ec->Lclean_bg < 7*ec->Lclean) /* (ec->Lclean_bg < 0.875*ec->Lclean) */ &&
 505         (8*ec->Lclean_bg < ec->Ltx)      /* (ec->Lclean_bg < 0.125*ec->Ltx)    */ )
 506     {
 507         if (ec->cond_met == 6) {
 508             /* BG filter has had better results for 6 consecutive samples */
 509             ec->adapt = 1;
 510             memcpy(ec->fir_taps16[0], ec->fir_taps16[1], ec->taps*sizeof(int16_t));
 511         }
 512         else
 513             ec->cond_met++;
 514     }
 515     else
 516         ec->cond_met = 0;
 517
 518     /* Non-Linear Processing ---------------------------------------------------*/
 519
 520     ec->clean_nlp = ec->clean;
 521     if (ec->adaption_mode & ECHO_CAN_USE_NLP)
 522     {
 523         /* Non-linear processor - a fancy way to say "zap small signals, to avoid
 524            residual echo due to (uLaw/ALaw) non-linearity in the channel.". */
 525
 526       if ((16*ec->Lclean < ec->Ltx))
 527       {
 528         /* Our e/c has improved echo by at least 24 dB (each factor of 2 is 6dB,
 529            so 2*2*2*2=16 is the same as 6+6+6+6=24dB) */
 530         if (ec->adaption_mode & ECHO_CAN_USE_CNG)
 531         {
 532             ec->cng_level = ec->Lbgn;
 533
 534             /* Very elementary comfort noise generation.  Just random
 535                numbers rolled off very vaguely Hoth-like.  DR: This
 536                noise doesn't sound quite right to me - I suspect there
 537                are some overlfow issues in the filtering as it's too
 538                "crackly".  TODO: debug this, maybe just play noise at
 539                high level or look at spectrum.
 540             */
 541
 542             ec->cng_rndnum = 1664525U*ec->cng_rndnum + 1013904223U;
 543             ec->cng_filter = ((ec->cng_rndnum & 0xFFFF) - 32768 + 5*ec->cng_filter) >> 3;
 544             ec->clean_nlp = (ec->cng_filter*ec->cng_level*8) >> 14;
 545
 546         }
 547         else if (ec->adaption_mode & ECHO_CAN_USE_CLIP)
 548         {
 549             /* This sounds much better than CNG */
 550             if (ec->clean_nlp > ec->Lbgn)
 551               ec->clean_nlp = ec->Lbgn;
 552             if (ec->clean_nlp < -ec->Lbgn)
 553               ec->clean_nlp = -ec->Lbgn;
 554         }
 555         else
 556         {
 557           /* just mute the residual, doesn't sound very good, used mainly
 558              in G168 tests */
 559           ec->clean_nlp = 0;
 560         }
 561       }
 562       else {
 563           /* Background noise estimator.  I tried a few algorithms
 564              here without much luck.  This very simple one seems to
 565              work best, we just average the level using a slow (1 sec
 566              time const) filter if the current level is less than a
 567              (experimentally derived) constant.  This means we dont
 568              include high level signals like near end speech.  When
 569              combined with CNG or especially CLIP seems to work OK.
 570           */
 571           if (ec->Lclean < 40) {
 572               ec->Lbgn_acc += abs(ec->clean) - ec->Lbgn;
 573               ec->Lbgn = (ec->Lbgn_acc + (1<<11)) >> 12;
 574           }
 575        }
 576     }
 577
 578     /* Roll around the taps buffer */
 579     if (ec->curr_pos <= 0)
 580         ec->curr_pos = ec->taps;
 581     ec->curr_pos--;
 582
 583     if (ec->adaption_mode & ECHO_CAN_DISABLE)
 584       ec->clean_nlp = rx;
 585
 586     /* Output scaled back up again to match input scaling */
 587
 588     return (int16_t) ec->clean_nlp << 1;
 589 }
 590 EXPORT_SYMBOL_GPL(oslec_update);
 591 /*- End of function --------------------------------------------------------*/
 592
 593 /* This function is seperated from the echo canceller is it is usually called
 594    as part of the tx process.  See rx HP (DC blocking) filter above, it's
 595    the same design.
 596
 597    Some soft phones send speech signals with a lot of low frequency
 598    energy, e.g. down to 20Hz.  This can make the hybrid non-linear
 599    which causes the echo canceller to fall over.  This filter can help
 600    by removing any low frequency before it gets to the tx port of the
 601    hybrid.
 602
 603    It can also help by removing and DC in the tx signal.  DC is bad
 604    for LMS algorithms.
 605
 606    This is one of the classic DC removal filters, adjusted to provide sufficient
 607    bass rolloff to meet the above requirement to protect hybrids from things that
 608    upset them. The difference between successive samples produces a lousy HPF, and
 609    then a suitably placed pole flattens things out. The final result is a nicely
 610    rolled off bass end. The filtering is implemented with extended fractional
 611    precision, which noise shapes things, giving very clean DC removal.
 612 */
 613
 614 int16_t oslec_hpf_tx(struct oslec_state *ec, int16_t tx) {
 615     int tmp, tmp1;
 616
 617     if (ec->adaption_mode & ECHO_CAN_USE_TX_HPF) {
 618         tmp = tx << 15;
 619 #if 1
 620         /* Make sure the gain of the HPF is 1.0. The first can still saturate a little under
 621            impulse conditions, and it might roll to 32768 and need clipping on sustained peak
 622            level signals. However, the scale of such clipping is small, and the error due to
 623            any saturation should not markedly affect the downstream processing. */
 624         tmp -= (tmp >> 4);
 625 #endif
 626         ec->tx_1 += -(ec->tx_1>>DC_LOG2BETA) + tmp - ec->tx_2;
 627         tmp1 = ec->tx_1 >> 15;
 628         if (tmp1 > 32767) tmp1 = 32767;
 629         if (tmp1 < -32767) tmp1 = -32767;
 630         tx = tmp1;
 631         ec->tx_2 = tmp;
 632     }
 633
 634     return tx;
 635 }
 636 EXPORT_SYMBOL_GPL(oslec_hpf_tx);
 637
 638 MODULE_LICENSE("GPL");
 639 MODULE_AUTHOR("David Rowe");
 640 MODULE_DESCRIPTION("Open Source Line Echo Canceller");
 641 MODULE_VERSION("0.3.0");