1/* aec.h 2 * 3 * Copyright (C) DFS Deutsche Flugsicherung (2004, 2005). 4 * All Rights Reserved. 5 * Author: Andre Adrian 6 * 7 * Acoustic Echo Cancellation Leaky NLMS-pw algorithm 8 * 9 * Version 0.3 filter created with www.dsptutor.freeuk.com 10 * Version 0.3.1 Allow change of stability parameter delta 11 * Version 0.4 Leaky Normalized LMS - pre whitening algorithm 12 */ 13 14#ifndef _AEC_H /* include only once */ 15 16#ifdef HAVE_CONFIG_H 17#include <config.h> 18#endif 19 20#include <pulse/gccmacro.h> 21#include <pulse/xmalloc.h> 22 23#include <pulsecore/macro.h> 24 25#define WIDEB 2 26 27// use double if your CPU does software-emulation of float 28#define REAL float 29 30/* dB Values */ 31#define M0dB 1.0f 32#define M3dB 0.71f 33#define M6dB 0.50f 34#define M9dB 0.35f 35#define M12dB 0.25f 36#define M18dB 0.125f 37#define M24dB 0.063f 38 39/* dB values for 16bit PCM */ 40/* MxdB_PCM = 32767 * 10 ^(x / 20) */ 41#define M10dB_PCM 10362.0f 42#define M20dB_PCM 3277.0f 43#define M25dB_PCM 1843.0f 44#define M30dB_PCM 1026.0f 45#define M35dB_PCM 583.0f 46#define M40dB_PCM 328.0f 47#define M45dB_PCM 184.0f 48#define M50dB_PCM 104.0f 49#define M55dB_PCM 58.0f 50#define M60dB_PCM 33.0f 51#define M65dB_PCM 18.0f 52#define M70dB_PCM 10.0f 53#define M75dB_PCM 6.0f 54#define M80dB_PCM 3.0f 55#define M85dB_PCM 2.0f 56#define M90dB_PCM 1.0f 57 58#define MAXPCM 32767.0f 59 60/* Design constants (Change to fine tune the algorithms */ 61 62/* The following values are for hardware AEC and studio quality 63 * microphone */ 64 65/* NLMS filter length in taps (samples). A longer filter length gives 66 * better Echo Cancellation, but maybe slower convergence speed and 67 * needs more CPU power (Order of NLMS is linear) */ 68#define NLMS_LEN (100*WIDEB*8) 69 70/* Vector w visualization length in taps (samples). 71 * Must match argv value for wdisplay.tcl */ 72#define DUMP_LEN (40*WIDEB*8) 73 74/* minimum energy in xf. Range: M70dB_PCM to M50dB_PCM. Should be equal 75 * to microphone ambient Noise level */ 76#define NoiseFloor M55dB_PCM 77 78/* Leaky hangover in taps. 79 */ 80#define Thold (60 * WIDEB * 8) 81 82// Adrian soft decision DTD 83// left point. X is ratio, Y is stepsize 84#define STEPX1 1.0 85#define STEPY1 1.0 86// right point. STEPX2=2.0 is good double talk, 3.0 is good single talk. 87#define STEPX2 2.5 88#define STEPY2 0 89#define ALPHAFAST (1.0f / 100.0f) 90#define ALPHASLOW (1.0f / 20000.0f) 91 92 93 94/* Ageing multiplier for LMS memory vector w */ 95#define Leaky 0.9999f 96 97/* Double Talk Detector Speaker/Microphone Threshold. Range <=1 98 * Large value (M0dB) is good for Single-Talk Echo cancellation, 99 * small value (M12dB) is good for Double-Talk AEC */ 100#define GeigelThreshold M6dB 101 102/* for Non Linear Processor. Range >0 to 1. Large value (M0dB) is good 103 * for Double-Talk, small value (M12dB) is good for Single-Talk */ 104#define NLPAttenuation M12dB 105 106/* Below this line there are no more design constants */ 107 108typedef struct IIR_HP IIR_HP; 109 110/* Exponential Smoothing or IIR Infinite Impulse Response Filter */ 111struct IIR_HP { 112 REAL x; 113}; 114 115static IIR_HP* IIR_HP_init(void) { 116 IIR_HP *i = pa_xnew(IIR_HP, 1); 117 i->x = 0.0f; 118 return i; 119 } 120 121static REAL IIR_HP_highpass(IIR_HP *i, REAL in) { 122 const REAL a0 = 0.01f; /* controls Transfer Frequency */ 123 /* Highpass = Signal - Lowpass. Lowpass = Exponential Smoothing */ 124 i->x += a0 * (in - i->x); 125 return in - i->x; 126 } 127 128typedef struct FIR_HP_300Hz FIR_HP_300Hz; 129 130#if WIDEB==1 131/* 17 taps FIR Finite Impulse Response filter 132 * Coefficients calculated with 133 * www.dsptutor.freeuk.com/KaiserFilterDesign/KaiserFilterDesign.html 134 */ 135class FIR_HP_300Hz { 136 REAL z[18]; 137 138public: 139 FIR_HP_300Hz() { 140 memset(this, 0, sizeof(FIR_HP_300Hz)); 141 } 142 143 REAL highpass(REAL in) { 144 const REAL a[18] = { 145 // Kaiser Window FIR Filter, Filter type: High pass 146 // Passband: 300.0 - 4000.0 Hz, Order: 16 147 // Transition band: 75.0 Hz, Stopband attenuation: 10.0 dB 148 -0.034870606, -0.039650206, -0.044063766, -0.04800318, 149 -0.051370874, -0.054082647, -0.056070227, -0.057283327, 150 0.8214126, -0.057283327, -0.056070227, -0.054082647, 151 -0.051370874, -0.04800318, -0.044063766, -0.039650206, 152 -0.034870606, 0.0 153 }; 154 memmove(z + 1, z, 17 * sizeof(REAL)); 155 z[0] = in; 156 REAL sum0 = 0.0, sum1 = 0.0; 157 int j; 158 159 for (j = 0; j < 18; j += 2) { 160 // optimize: partial loop unrolling 161 sum0 += a[j] * z[j]; 162 sum1 += a[j + 1] * z[j + 1]; 163 } 164 return sum0 + sum1; 165 } 166}; 167 168#else 169 170/* 35 taps FIR Finite Impulse Response filter 171 * Passband 150Hz to 4kHz for 8kHz sample rate, 300Hz to 8kHz for 16kHz 172 * sample rate. 173 * Coefficients calculated with 174 * www.dsptutor.freeuk.com/KaiserFilterDesign/KaiserFilterDesign.html 175 */ 176struct FIR_HP_300Hz { 177 REAL z[36]; 178}; 179 180static FIR_HP_300Hz* FIR_HP_300Hz_init(void) { 181 FIR_HP_300Hz *ret = pa_xnew(FIR_HP_300Hz, 1); 182 memset(ret, 0, sizeof(FIR_HP_300Hz)); 183 return ret; 184 } 185 186static REAL FIR_HP_300Hz_highpass(FIR_HP_300Hz *f, REAL in) { 187 REAL sum0 = 0.0, sum1 = 0.0; 188 int j; 189 const REAL a[36] = { 190 // Kaiser Window FIR Filter, Filter type: High pass 191 // Passband: 150.0 - 4000.0 Hz, Order: 34 192 // Transition band: 34.0 Hz, Stopband attenuation: 10.0 dB 193 -0.016165324, -0.017454365, -0.01871232, -0.019931411, 194 -0.021104068, -0.022222936, -0.02328091, -0.024271343, 195 -0.025187887, -0.02602462, -0.026776174, -0.027437767, 196 -0.028004972, -0.028474221, -0.028842418, -0.029107114, 197 -0.02926664, 0.8524841, -0.02926664, -0.029107114, 198 -0.028842418, -0.028474221, -0.028004972, -0.027437767, 199 -0.026776174, -0.02602462, -0.025187887, -0.024271343, 200 -0.02328091, -0.022222936, -0.021104068, -0.019931411, 201 -0.01871232, -0.017454365, -0.016165324, 0.0 202 }; 203 memmove(f->z + 1, f->z, 35 * sizeof(REAL)); 204 f->z[0] = in; 205 206 for (j = 0; j < 36; j += 2) { 207 // optimize: partial loop unrolling 208 sum0 += a[j] * f->z[j]; 209 sum1 += a[j + 1] * f->z[j + 1]; 210 } 211 return sum0 + sum1; 212 } 213#endif 214 215typedef struct IIR1 IIR1; 216 217/* Recursive single pole IIR Infinite Impulse response High-pass filter 218 * 219 * Reference: The Scientist and Engineer's Guide to Digital Processing 220 * 221 * output[N] = A0 * input[N] + A1 * input[N-1] + B1 * output[N-1] 222 * 223 * X = exp(-2.0 * pi * Fc) 224 * A0 = (1 + X) / 2 225 * A1 = -(1 + X) / 2 226 * B1 = X 227 * Fc = cutoff freq / sample rate 228 */ 229struct IIR1 { 230 REAL in0, out0; 231 REAL a0, a1, b1; 232}; 233 234#if 0 235 IIR1() { 236 memset(this, 0, sizeof(IIR1)); 237 } 238#endif 239 240static IIR1* IIR1_init(REAL Fc) { 241 IIR1 *i = pa_xnew(IIR1, 1); 242 i->b1 = expf(-2.0f * M_PI * Fc); 243 i->a0 = (1.0f + i->b1) / 2.0f; 244 i->a1 = -(i->a0); 245 i->in0 = 0.0f; 246 i->out0 = 0.0f; 247 return i; 248 } 249 250static REAL IIR1_highpass(IIR1 *i, REAL in) { 251 REAL out = i->a0 * in + i->a1 * i->in0 + i->b1 * i->out0; 252 i->in0 = in; 253 i->out0 = out; 254 return out; 255 } 256 257 258#if 0 259/* Recursive two pole IIR Infinite Impulse Response filter 260 * Coefficients calculated with 261 * http://www.dsptutor.freeuk.com/IIRFilterDesign/IIRFiltDes102.html 262 */ 263class IIR2 { 264 REAL x[2], y[2]; 265 266public: 267 IIR2() { 268 memset(this, 0, sizeof(IIR2)); 269 } 270 271 REAL highpass(REAL in) { 272 // Butterworth IIR filter, Filter type: HP 273 // Passband: 2000 - 4000.0 Hz, Order: 2 274 const REAL a[] = { 0.29289323f, -0.58578646f, 0.29289323f }; 275 const REAL b[] = { 1.3007072E-16f, 0.17157288f }; 276 REAL out = 277 a[0] * in + a[1] * x[0] + a[2] * x[1] - b[0] * y[0] - b[1] * y[1]; 278 279 x[1] = x[0]; 280 x[0] = in; 281 y[1] = y[0]; 282 y[0] = out; 283 return out; 284 } 285}; 286#endif 287 288 289// Extension in taps to reduce mem copies 290#define NLMS_EXT (10*8) 291 292// block size in taps to optimize DTD calculation 293#define DTD_LEN 16 294 295typedef struct AEC AEC; 296 297struct AEC { 298 // Time domain Filters 299 IIR_HP *acMic, *acSpk; // DC-level remove Highpass) 300 FIR_HP_300Hz *cutoff; // 150Hz cut-off Highpass 301 REAL gain; // Mic signal amplify 302 IIR1 *Fx, *Fe; // pre-whitening Highpass for x, e 303 304 // Adrian soft decision DTD (Double Talk Detector) 305 REAL dfast, xfast; 306 REAL dslow, xslow; 307 308 // NLMS-pw 309 REAL x[NLMS_LEN + NLMS_EXT]; // tap delayed loudspeaker signal 310 REAL xf[NLMS_LEN + NLMS_EXT]; // pre-whitening tap delayed signal 311 REAL w_arr[NLMS_LEN + (16 / sizeof(REAL))]; // tap weights 312 REAL *w; // this will be a 16-byte aligned pointer into w_arr 313 int j; // optimize: less memory copies 314 double dotp_xf_xf; // double to avoid loss of precision 315 float delta; // noise floor to stabilize NLMS 316 317 // AES 318 float aes_y2; // not in use! 319 320 // w vector visualization 321 REAL ws[DUMP_LEN]; // tap weights sums 322 int fdwdisplay; // TCP file descriptor 323 int dumpcnt; // wdisplay output counter 324 325 // variables are public for visualization 326 int hangover; 327 float stepsize; 328 329 // vfuncs that are picked based on processor features available 330 REAL (*dotp) (REAL[], REAL[]); 331}; 332 333/* Double-Talk Detector 334 * 335 * in d: microphone sample (PCM as REALing point value) 336 * in x: loudspeaker sample (PCM as REALing point value) 337 * return: from 0 for doubletalk to 1.0 for single talk 338 */ 339static float AEC_dtd(AEC *a, REAL d, REAL x); 340 341static void AEC_leaky(AEC *a); 342 343/* Normalized Least Mean Square Algorithm pre-whitening (NLMS-pw) 344 * The LMS algorithm was developed by Bernard Widrow 345 * book: Haykin, Adaptive Filter Theory, 4. edition, Prentice Hall, 2002 346 * 347 * in d: microphone sample (16bit PCM value) 348 * in x_: loudspeaker sample (16bit PCM value) 349 * in stepsize: NLMS adaptation variable 350 * return: echo cancelled microphone sample 351 */ 352static REAL AEC_nlms_pw(AEC *a, REAL d, REAL x_, float stepsize); 353 354AEC* AEC_init(int RATE, int have_vector); 355void AEC_done(AEC *a); 356 357/* Acoustic Echo Cancellation and Suppression of one sample 358 * in d: microphone signal with echo 359 * in x: loudspeaker signal 360 * return: echo cancelled microphone signal 361 */ 362 int AEC_doAEC(AEC *a, int d_, int x_); 363 364PA_GCC_UNUSED static float AEC_getambient(AEC *a) { 365 return a->dfast; 366 } 367static void AEC_setambient(AEC *a, float Min_xf) { 368 a->dotp_xf_xf -= a->delta; // subtract old delta 369 a->delta = (NLMS_LEN-1) * Min_xf * Min_xf; 370 a->dotp_xf_xf += a->delta; // add new delta 371 } 372PA_GCC_UNUSED static void AEC_setgain(AEC *a, float gain_) { 373 a->gain = gain_; 374 } 375#if 0 376 void AEC_openwdisplay(AEC *a); 377#endif 378PA_GCC_UNUSED static void AEC_setaes(AEC *a, float aes_y2_) { 379 a->aes_y2 = aes_y2_; 380 } 381 382#define _AEC_H 383#endif 384