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