1 /** 2 Basic IIR 1-pole and 2-pole filters through biquads. 3 4 Copyright: Guillaume Piolat 2015-2017. 5 License: $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0) 6 */ 7 module dplug.dsp.iir; 8 9 import std.math; 10 import dplug.core.math; 11 import inteli.emmintrin; 12 13 // DMD with a 32-bit target uses the FPU 14 version(X86) 15 { 16 version(DigitalMars) 17 { 18 version = killDenormals; 19 } 20 } 21 22 // TODO: function to make biquads from z-plane poles and zeroes 23 24 public 25 { 26 27 /// Type which hold the biquad coefficients. 28 /// Important: Coefficients are considered always normalized by a0. 29 /// Note: coeff[0] is b0, 30 /// coeff[1] is b1, 31 /// coeff[2] is b2, 32 /// coeff[3] is a1, 33 /// coeff[4] is a2 in the litterature. 34 alias BiquadCoeff = double[5]; 35 36 /// Maintain state for a biquad state. 37 /// A biquad is a realization that can model two poles and two zeros. 38 struct BiquadDelay 39 { 40 enum order = 2; 41 42 double _x0; 43 double _x1; 44 double _y0; 45 double _y1; 46 47 void initialize() nothrow @nogc 48 { 49 _x0 = 0; 50 _x1 = 0; 51 _y0 = 0; 52 _y1 = 0; 53 } 54 55 static if (order == 2) 56 { 57 /// Process a single sample through the biquad filter. 58 /// This is rather inefficient, in general you'll prefer to use the `nextBuffer` 59 /// functions. 60 float nextSample(float input, const(BiquadCoeff) coeff) nothrow @nogc 61 { 62 double x1 = _x0, 63 x2 = _x1, 64 y1 = _y0, 65 y2 = _y1; 66 67 double a0 = coeff[0], 68 a1 = coeff[1], 69 a2 = coeff[2], 70 a3 = coeff[3], 71 a4 = coeff[4]; 72 73 double current = a0 * input + a1 * x1 + a2 * x2 - a3 * y1 - a4 * y2; 74 75 // kill denormals,and double values that would be converted 76 // to float denormals 77 version(killDenormals) 78 { 79 current += 1e-18f; 80 current -= 1e-18f; 81 } 82 83 _x0 = input; 84 _x1 = x1; 85 _y0 = current; 86 _y1 = y1; 87 return current; 88 } 89 ///ditto 90 double nextSample(double input, const(BiquadCoeff) coeff) nothrow @nogc 91 { 92 double x1 = _x0, 93 x2 = _x1, 94 y1 = _y0, 95 y2 = _y1; 96 97 double a0 = coeff[0], 98 a1 = coeff[1], 99 a2 = coeff[2], 100 a3 = coeff[3], 101 a4 = coeff[4]; 102 103 double current = a0 * input + a1 * x1 + a2 * x2 - a3 * y1 - a4 * y2; 104 105 // kill denormals, and double values that would be converted 106 // to float denormals 107 version(killDenormals) 108 { 109 current += 1e-18f; 110 current -= 1e-18f; 111 } 112 113 _x0 = input; 114 _x1 = x1; 115 _y0 = current; 116 _y1 = y1; 117 return current; 118 } 119 120 void nextBuffer(const(float)* input, float* output, int frames, const(BiquadCoeff) coeff) nothrow @nogc 121 { 122 static if (D_InlineAsm_Any) 123 { 124 double x0 = _x0, 125 x1 = _x1, 126 y0 = _y0, 127 y1 = _y1; 128 129 double a0 = coeff[0], 130 a1 = coeff[1], 131 a2 = coeff[2], 132 a3 = coeff[3], 133 a4 = coeff[4]; 134 135 version(LDC) 136 { 137 import inteli.emmintrin; 138 139 __m128d XMM0 = _mm_set_pd(x1, x0); 140 __m128d XMM1 = _mm_set_pd(y1, y0); 141 __m128d XMM2 = _mm_set_pd(a2, a1); 142 __m128d XMM3 = _mm_set_pd(a4, a3); 143 __m128d XMM4 = _mm_set_sd(a0); 144 145 __m128d XMM6 = _mm_undefined_pd(); 146 __m128d XMM7 = _mm_undefined_pd(); 147 __m128d XMM5 = _mm_undefined_pd(); 148 for (int n = 0; n < frames; ++n) 149 { 150 __m128 INPUT = _mm_load_ss(input + n); 151 XMM5 = _mm_setzero_pd(); 152 XMM5 = _mm_cvtss_sd(XMM5,INPUT); 153 154 XMM6 = XMM0; 155 XMM7 = XMM1; 156 XMM5 = _mm_mul_pd(XMM5, XMM4); // input[i]*a0 157 XMM6 = _mm_mul_pd(XMM6, XMM2); // x1*a2 x0*a1 158 XMM7 = _mm_mul_pd(XMM7, XMM3); // y1*a4 y0*a3 159 XMM5 = _mm_add_pd(XMM5, XMM6); 160 XMM5 = _mm_sub_pd(XMM5, XMM7); // x1*a2 - y1*a4 | input[i]*a0 + x0*a1 - y0*a3 161 XMM6 = XMM5; 162 163 XMM0 = cast(double2) _mm_slli_si128!8(cast(__m128i) XMM0); 164 XMM6 = cast(double2) _mm_srli_si128!8(cast(__m128i) XMM6); 165 166 XMM0 = _mm_cvtss_sd(XMM0, INPUT); 167 XMM5 = _mm_add_pd(XMM5, XMM6); 168 XMM7 = cast(double2) _mm_cvtsd_ss(_mm_undefined_ps(), XMM5); 169 XMM5 = _mm_unpacklo_pd(XMM5, XMM1); 170 XMM1 = XMM5; 171 _mm_store_ss(output + n, cast(__m128) XMM7); 172 } 173 _mm_storel_pd(&x0, XMM0); 174 _mm_storeh_pd(&x1, XMM0); 175 _mm_storel_pd(&y0, XMM1); 176 _mm_storeh_pd(&y1, XMM1); 177 } 178 else version(D_InlineAsm_X86) 179 { 180 asm nothrow @nogc 181 { 182 mov EAX, input; 183 mov EDX, output; 184 mov ECX, frames; 185 186 movlpd XMM0, qword ptr x0; // XMM0 = x1 x0 187 movhpd XMM0, qword ptr x1; 188 movlpd XMM1, qword ptr y0; // XMM1 = y1 y0 189 movhpd XMM1, qword ptr y1; 190 191 movlpd XMM2, qword ptr a1; // XMM2 = a2 a1 192 movhpd XMM2, qword ptr a2; 193 movlpd XMM3, qword ptr a3; // XMM3 = a4 a3 194 movhpd XMM3, qword ptr a4; 195 196 movq XMM4, qword ptr a0; // XMM4 = 0 a0 197 198 align 16; 199 loop: 200 pxor XMM5, XMM5; 201 cvtss2sd XMM5, dword ptr [EAX]; 202 203 movapd XMM6, XMM0; 204 movapd XMM7, XMM1; 205 206 mulpd XMM5, XMM4; // input[i]*a0 207 mulpd XMM6, XMM2; // x1*a2 x0*a1 208 mulpd XMM7, XMM3; // y1*a4 y0*a3 209 210 addpd XMM5, XMM6; 211 subpd XMM5, XMM7; // x1*a2 - y1*a4 | input[i]*a0 + x0*a1 - y0*a3 212 213 movapd XMM6, XMM5; 214 pslldq XMM0, 8; 215 psrldq XMM6, 8; 216 217 cvtss2sd XMM0, dword ptr [EAX]; // XMM0 <- x0 input[i] 218 addpd XMM5, XMM6; // garbage | input[i]*a0 + x0*a1 - y0*a3 + x1*a2 - y1*a4 219 220 cvtsd2ss XMM7, XMM5; 221 punpcklqdq XMM5, XMM1; // XMM5 <- y0 current 222 add EAX, 4; 223 movd dword ptr [EDX], XMM7; 224 add EDX, 4; 225 movapd XMM1, XMM5; 226 227 sub ECX, 1; 228 jnz short loop; 229 230 movlpd qword ptr x0, XMM0; 231 movhpd qword ptr x1, XMM0; 232 movlpd qword ptr y0, XMM1; 233 movhpd qword ptr y1, XMM1; 234 } 235 } 236 else version(D_InlineAsm_X86_64) 237 { 238 ubyte[16] storage0; 239 ubyte[16] storage1; 240 asm nothrow @nogc 241 { 242 movups storage0, XMM6; 243 movups storage1, XMM7; 244 245 mov RAX, input; 246 mov RDX, output; 247 mov ECX, frames; 248 249 movlpd XMM0, qword ptr x0; // XMM0 = x1 x0 250 movhpd XMM0, qword ptr x1; 251 movlpd XMM1, qword ptr y0; // XMM1 = y1 y0 252 movhpd XMM1, qword ptr y1; 253 254 movlpd XMM2, qword ptr a1; // XMM2 = a2 a1 255 movhpd XMM2, qword ptr a2; 256 movlpd XMM3, qword ptr a3; // XMM3 = a4 a3 257 movhpd XMM3, qword ptr a4; 258 259 movq XMM4, qword ptr a0; // XMM4 = 0 a0 260 261 align 16; 262 loop: 263 pxor XMM5, XMM5; 264 cvtss2sd XMM5, dword ptr [RAX]; 265 266 movapd XMM6, XMM0; 267 movapd XMM7, XMM1; 268 269 mulpd XMM5, XMM4; // input[i]*a0 270 mulpd XMM6, XMM2; // x1*a2 x0*a1 271 mulpd XMM7, XMM3; // y1*a4 y0*a3 272 273 addpd XMM5, XMM6; 274 subpd XMM5, XMM7; // x1*a2 - y1*a4 | input[i]*a0 + x0*a1 - y0*a3 275 276 movapd XMM6, XMM5; 277 pslldq XMM0, 8; 278 psrldq XMM6, 8; 279 280 addpd XMM5, XMM6; // garbage | input[i]*a0 + x0*a1 - y0*a3 + x1*a2 - y1*a4 281 cvtss2sd XMM0, dword ptr [RAX]; // XMM0 <- x0 input[i] 282 cvtsd2ss XMM7, XMM5; 283 punpcklqdq XMM5, XMM1; // XMM5 <- y0 current 284 add RAX, 4; 285 movd dword ptr [RDX], XMM7; 286 add RDX, 4; 287 movapd XMM1, XMM5; 288 289 sub ECX, 1; 290 jnz short loop; 291 292 movlpd qword ptr x0, XMM0; 293 movhpd qword ptr x1, XMM0; 294 movlpd qword ptr y0, XMM1; 295 movhpd qword ptr y1, XMM1; 296 297 movups XMM6, storage0; // XMMx with x >= 6 registers need to be preserved 298 movups XMM7, storage1; 299 } 300 } 301 else 302 static assert(false, "Not implemented for this platform."); 303 304 // Kill small signals that can cause denormals (no precision loss was measurable) 305 version(killDenormals) 306 { 307 x0 += 1e-10; 308 x0 -= 1e-10; 309 x1 += 1e-10; 310 x1 -= 1e-10; 311 y0 += 1e-10; 312 y0 -= 1e-10; 313 y1 += 1e-10; 314 y1 -= 1e-10; 315 } 316 317 _x0 = x0; 318 _x1 = x1; 319 _y0 = y0; 320 _y1 = y1; 321 } 322 else 323 { 324 double x0 = _x0, 325 x1 = _x1, 326 y0 = _y0, 327 y1 = _y1; 328 329 double a0 = coeff[0], 330 a1 = coeff[1], 331 a2 = coeff[2], 332 a3 = coeff[3], 333 a4 = coeff[4]; 334 335 for(int i = 0; i < frames; ++i) 336 { 337 double current = a0 * input[i] + a1 * x0 + a2 * x1 - a3 * y0 - a4 * y1; 338 339 // kill denormals,and double values that would be converted 340 // to float denormals 341 version(killDenormals) 342 { 343 current += 1e-18f; 344 current -= 1e-18f; 345 } 346 347 x1 = x0; 348 x0 = input[i]; 349 y1 = y0; 350 y0 = current; 351 output[i] = current; 352 } 353 354 _x0 = x0; 355 _x1 = x1; 356 _y0 = y0; 357 _y1 = y1; 358 } 359 } 360 ///ditto 361 void nextBuffer(const(double)* input, double* output, int frames, const(BiquadCoeff) coeff) nothrow @nogc 362 { 363 double x0 = _x0, 364 x1 = _x1, 365 y0 = _y0, 366 y1 = _y1; 367 368 double a0 = coeff[0], 369 a1 = coeff[1], 370 a2 = coeff[2], 371 a3 = coeff[3], 372 a4 = coeff[4]; 373 374 __m128d XMM0 = _mm_set_pd(x1, x0); 375 __m128d XMM1 = _mm_set_pd(y1, y0); 376 __m128d XMM2 = _mm_set_pd(a2, a1); 377 __m128d XMM3 = _mm_set_pd(a4, a3); 378 __m128d XMM4 = _mm_set_sd(a0); 379 380 __m128d XMM6 = _mm_undefined_pd(); 381 __m128d XMM7 = _mm_undefined_pd(); 382 __m128d XMM5 = _mm_undefined_pd(); 383 for (int n = 0; n < frames; ++n) 384 { 385 __m128d INPUT = _mm_load_sd(input + n); 386 XMM5 = INPUT; 387 388 XMM6 = XMM0; 389 XMM7 = XMM1; 390 XMM5 = _mm_mul_pd(XMM5, XMM4); // input[i]*a0 391 XMM6 = _mm_mul_pd(XMM6, XMM2); // x1*a2 x0*a1 392 XMM7 = _mm_mul_pd(XMM7, XMM3); // y1*a4 y0*a3 393 XMM5 = _mm_add_pd(XMM5, XMM6); 394 XMM5 = _mm_sub_pd(XMM5, XMM7); // x1*a2 - y1*a4 | input[i]*a0 + x0*a1 - y0*a3 395 XMM6 = XMM5; 396 397 XMM0 = cast(double2) _mm_slli_si128!8(cast(__m128i) XMM0); 398 XMM6 = cast(double2) _mm_srli_si128!8(cast(__m128i) XMM6); 399 XMM0.ptr[0] = INPUT.array[0]; 400 XMM5 = _mm_add_pd(XMM5, XMM6); 401 XMM7 = XMM5; 402 XMM5 = _mm_unpacklo_pd(XMM5, XMM1); 403 XMM1 = XMM5; 404 _mm_store_sd(output + n, XMM7); 405 } 406 _mm_storel_pd(&x0, XMM0); 407 _mm_storeh_pd(&x1, XMM0); 408 _mm_storel_pd(&y0, XMM1); 409 _mm_storeh_pd(&y1, XMM1); 410 } 411 412 /// Special version of biquad processing, for a constant DC input. 413 void nextBuffer(float input, float* output, int frames, const(BiquadCoeff) coeff) nothrow @nogc 414 { 415 // Note: this naive version performs better than an intel-intrinsics one 416 double x0 = _x0, 417 x1 = _x1, 418 y0 = _y0, 419 y1 = _y1; 420 421 double a0 = coeff[0], 422 a1 = coeff[1], 423 a2 = coeff[2], 424 a3 = coeff[3], 425 a4 = coeff[4]; 426 427 for(int i = 0; i < frames; ++i) 428 { 429 double current = a0 * input + a1 * x0 + a2 * x1 - a3 * y0 - a4 * y1; 430 431 // kill denormals,and double values that would be converted 432 // to float denormals 433 version(killDenormals) 434 { 435 current += 1e-18f; 436 current -= 1e-18f; 437 } 438 439 x1 = x0; 440 x0 = input; 441 y1 = y0; 442 y0 = current; 443 output[i] = current; 444 } 445 446 _x0 = x0; 447 _x1 = x1; 448 _y0 = y0; 449 _y1 = y1; 450 } 451 ///ditto 452 void nextBuffer(double input, double* output, int frames, const(BiquadCoeff) coeff) nothrow @nogc 453 { 454 double x0 = _x0, 455 x1 = _x1, 456 y0 = _y0, 457 y1 = _y1; 458 459 double a0 = coeff[0], 460 a1 = coeff[1], 461 a2 = coeff[2], 462 a3 = coeff[3], 463 a4 = coeff[4]; 464 465 for(int i = 0; i < frames; ++i) 466 { 467 double current = a0 * input + a1 * x0 + a2 * x1 - a3 * y0 - a4 * y1; 468 469 // kill denormals,and double values that would be converted 470 // to float denormals 471 version(killDenormals) 472 { 473 current += 1e-18f; 474 current -= 1e-18f; 475 } 476 477 x1 = x0; 478 x0 = input; 479 y1 = y0; 480 y0 = current; 481 output[i] = current; 482 } 483 484 _x0 = x0; 485 _x1 = x1; 486 _y0 = y0; 487 _y1 = y1; 488 } 489 } 490 } 491 492 /// 1-pole low-pass filter. 493 /// Note: the cutoff frequency can be >= nyquist, in which case it asymptotically approaches a bypass. 494 /// the cutoff frequency can be below 0 Hz, in which case it is equal to zero. 495 /// This filter is normalized on DC. 496 /// You always have -3 dB at cutoff in the valid range. 497 /// 498 /// See_also: http://www.earlevel.com/main/2012/12/15/a-one-pole-filter/ 499 BiquadCoeff biquadOnePoleLowPass(double frequency, double sampleRate) nothrow @nogc 500 { 501 double fc = frequency / sampleRate; 502 if (fc < 0.0f) 503 fc = 0.0f; 504 double t2 = fast_exp(-2.0 * PI * fc); 505 BiquadCoeff result; 506 result[0] = 1 - t2; 507 result[1] = 0; 508 result[2] = 0; 509 result[3] = -t2; 510 result[4] = 0; 511 return result; 512 } 513 514 /// 1-pole high-pass filter. 515 /// Note: Like the corresponding one-pole lowpass, this is normalized for DC. 516 /// The cutoff frequency can be <= 0 Hz, in which case it is a bypass. 517 /// Going in very high frequency does NOT give zero. 518 /// You always have -3 dB at cutoff in the valid range. 519 /// 520 /// See_also: https://www.dspguide.com/ch19/2.html 521 BiquadCoeff biquadOnePoleHighPass(double frequency, double sampleRate) nothrow @nogc 522 { 523 double fc = frequency / sampleRate; 524 if (fc < 0.0f) 525 fc = 0.0f; 526 double t2 = fast_exp(-2.0 * PI * fc); 527 BiquadCoeff result; 528 result[0] = (1 + t2) * 0.5; 529 result[1] = -(1 + t2) * 0.5; 530 result[2] = 0; 531 result[3] = -t2; 532 result[4] = 0; 533 return result; 534 } 535 536 /// 1-pole low-pass filter, frequency mapping is not precise. 537 /// Not accurate across sample rates, but coefficient computation is cheap. Not advised. 538 BiquadCoeff biquadOnePoleLowPassImprecise(double frequency, double samplerate) nothrow @nogc 539 { 540 double t0 = frequency / samplerate; 541 if (t0 > 0.5f) 542 t0 = 0.5f; 543 544 double t1 = (1 - 2 * t0); 545 double t2 = t1 * t1; 546 547 BiquadCoeff result; 548 result[0] = cast(float)(1 - t2); 549 result[1] = 0; 550 result[2] = 0; 551 result[3] = cast(float)(-t2); 552 result[4] = 0; 553 return result; 554 } 555 556 /// 1-pole high-pass filter, frequency mapping is not precise. 557 /// Not accurate across sample rates, but coefficient computation is cheap. Not advised. 558 BiquadCoeff biquadOnePoleHighPassImprecise(double frequency, double samplerate) nothrow @nogc 559 { 560 double t0 = frequency / samplerate; 561 if (t0 > 0.5f) 562 t0 = 0.5f; 563 564 double t1 = (2 * t0); 565 double t2 = t1 * t1; 566 567 BiquadCoeff result; 568 result[0] = cast(float)(1 - t2); 569 result[1] = 0; 570 result[2] = 0; 571 result[3] = cast(float)(t2); 572 result[4] = 0; 573 return result; 574 } 575 576 // Cookbook formulae for audio EQ biquad filter coefficients 577 // by Robert Bristow-Johnson 578 579 /// Low-pass filter 12 dB/oct as described by Robert Bristow-Johnson. 580 BiquadCoeff biquadRBJLowPass(double frequency, double samplerate, double Q = SQRT1_2) nothrow @nogc 581 { 582 return generateBiquad(BiquadType.LOW_PASS_FILTER, frequency, samplerate, 0, Q); 583 } 584 585 /// High-pass filter 12 dB/oct as described by Robert Bristow-Johnson. 586 BiquadCoeff biquadRBJHighPass(double frequency, double samplerate, double Q = SQRT1_2) nothrow @nogc 587 { 588 return generateBiquad(BiquadType.HIGH_PASS_FILTER, frequency, samplerate, 0, Q); 589 } 590 591 /// Band-pass filter as described by Robert Bristow-Johnson. 592 BiquadCoeff biquadRBJBandPass(double frequency, double samplerate, double Q = SQRT1_2) nothrow @nogc 593 { 594 return generateBiquad(BiquadType.BAND_PASS_FILTER, frequency, samplerate, 0, Q); 595 } 596 597 /// Notch filter as described by Robert Bristow-Johnson. 598 BiquadCoeff biquadRBJNotch(double frequency, double samplerate, double Q = SQRT1_2) nothrow @nogc 599 { 600 return generateBiquad(BiquadType.NOTCH_FILTER, frequency, samplerate, 0, Q); 601 } 602 603 /// Peak filter as described by Robert Bristow-Johnson. 604 BiquadCoeff biquadRBJPeak(double frequency, double samplerate, double gain, double Q = SQRT1_2) nothrow @nogc 605 { 606 return generateBiquad(BiquadType.PEAK_FILTER, frequency, samplerate, gain, Q); 607 } 608 609 /// Low-shelf filter as described by Robert Bristow-Johnson. 610 BiquadCoeff biquadRBJLowShelf(double frequency, double samplerate, double gain, double Q = SQRT1_2) nothrow @nogc 611 { 612 return generateBiquad(BiquadType.LOW_SHELF, frequency, samplerate, gain, Q); 613 } 614 615 /// High-shelf filter as described by Robert Bristow-Johnson. 616 BiquadCoeff biquadRBJHighShelf(double frequency, double samplerate, double gain, double Q = SQRT1_2) nothrow @nogc 617 { 618 return generateBiquad(BiquadType.HIGH_SHELF, frequency, samplerate, gain, Q); 619 } 620 621 /// 2nd order All-pass filter as described by Robert Bristow-Johnson. 622 /// This is helpful to introduce the exact same phase response as the RBJ low-pass, but doesn't affect magnitude. 623 BiquadCoeff biquadRBJAllPass(double frequency, double samplerate, double Q = SQRT1_2) nothrow @nogc 624 { 625 return generateBiquad(BiquadType.ALL_PASS_FILTER, frequency, samplerate, 0, Q); 626 } 627 628 /// Identity biquad, pass signal unchanged. 629 BiquadCoeff biquadBypass() nothrow @nogc 630 { 631 BiquadCoeff result; 632 result[0] = 1; 633 result[1] = 0; 634 result[2] = 0; 635 result[3] = 0; 636 result[4] = 0; 637 return result; 638 } 639 640 /// Zero biquad, gives zero output. 641 BiquadCoeff biquadZero() nothrow @nogc 642 { 643 BiquadCoeff result; 644 result[0] = 0; 645 result[1] = 0; 646 result[2] = 0; 647 result[3] = 0; 648 result[4] = 0; 649 return result; 650 } 651 } 652 653 private 654 { 655 enum BiquadType 656 { 657 LOW_PASS_FILTER, 658 HIGH_PASS_FILTER, 659 BAND_PASS_FILTER, 660 NOTCH_FILTER, 661 PEAK_FILTER, 662 LOW_SHELF, 663 HIGH_SHELF, 664 ALL_PASS_FILTER 665 } 666 667 // generates RBJ biquad coefficients 668 BiquadCoeff generateBiquad(BiquadType type, double frequency, double samplerate, double gaindB, double Q) nothrow @nogc 669 { 670 // regardless of the output precision, always compute coefficients in double precision 671 672 double A = pow(10.0, gaindB / 40.0); 673 double w0 = (2.0 * PI) * frequency / samplerate; 674 double sin_w0 = sin(w0); 675 double cos_w0 = cos(w0); 676 677 double alpha = sin_w0 / (2 * Q); 678 679 // = sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0) ) (case: BW) 680 // = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 ) (case: S) 681 682 double b0, b1, b2, a0, a1, a2; 683 684 final switch(type) 685 { 686 case BiquadType.LOW_PASS_FILTER: 687 b0 = (1 - cos_w0) / 2; 688 b1 = 1 - cos_w0; 689 b2 = (1 - cos_w0) / 2; 690 a0 = 1 + alpha; 691 a1 = -2 * cos_w0; 692 a2 = 1 - alpha; 693 break; 694 695 case BiquadType.HIGH_PASS_FILTER: 696 b0 = (1 + cos_w0) / 2; 697 b1 = -(1 + cos_w0); 698 b2 = (1 + cos_w0) / 2; 699 a0 = 1 + alpha; 700 a1 = -2 * cos_w0; 701 a2 = 1 - alpha; 702 break; 703 704 case BiquadType.BAND_PASS_FILTER: 705 b0 = alpha; 706 b1 = 0; 707 b2 = -alpha; 708 a0 = 1 + alpha; 709 a1 = -2 * cos_w0; 710 a2 = 1 - alpha; 711 break; 712 713 case BiquadType.NOTCH_FILTER: 714 b0 = 1; 715 b1 = -2*cos_w0; 716 b2 = 1; 717 a0 = 1 + alpha; 718 a1 = -2*cos(w0); 719 a2 = 1 - alpha; 720 break; 721 722 case BiquadType.PEAK_FILTER: 723 b0 = 1 + alpha * A; 724 b1 = -2 * cos_w0; 725 b2 = 1 - alpha * A; 726 a0 = 1 + alpha / A; 727 a1 = -2 * cos_w0; 728 a2 = 1 - alpha / A; 729 break; 730 731 case BiquadType.LOW_SHELF: 732 { 733 double ap1 = A + 1; 734 double am1 = A - 1; 735 double M = 2 * sqrt(A) * alpha; 736 b0 = A * (ap1 - am1 * cos_w0 + M); 737 b1 = 2 * A * (am1 - ap1 * cos_w0); 738 b2 = A * (ap1 - am1 * cos_w0 - M); 739 a0 = ap1 + am1 * cos_w0 + M; 740 a1 = -2 * (am1 + ap1 * cos_w0); 741 a2 = ap1 + am1 * cos_w0 - M; 742 } 743 break; 744 745 case BiquadType.HIGH_SHELF: 746 { 747 double ap1 = A + 1; 748 double am1 = A - 1; 749 double M = 2 * sqrt(A) * alpha; 750 b0 = A * (ap1 + am1 * cos_w0 + M); 751 b1 = -2 * A * (am1 + ap1 * cos_w0); 752 b2 = A * (ap1 + am1 * cos_w0 - M); 753 a0 = ap1 - am1 * cos_w0 + M; 754 a1 = 2 * (am1 - ap1 * cos_w0); 755 a2 = ap1 - am1 * cos_w0 - M; 756 } 757 break; 758 759 case BiquadType.ALL_PASS_FILTER: 760 { 761 b0 = 1 - alpha; 762 b1 = -2 * cos_w0; 763 b2 = 1 + alpha; 764 a0 = 1 + alpha; 765 a1 = -2 * cos_w0; 766 a2 = 1 - alpha; 767 } 768 break; 769 } 770 771 BiquadCoeff result; 772 result[0] = cast(float)(b0 / a0); 773 result[1] = cast(float)(b1 / a0); 774 result[2] = cast(float)(b2 / a0); 775 result[3] = cast(float)(a1 / a0); 776 result[4] = cast(float)(a2 / a0); 777 return result; 778 } 779 } 780 781 782 version(D_InlineAsm_X86) 783 private enum D_InlineAsm_Any = true; 784 else version(D_InlineAsm_X86_64) 785 private enum D_InlineAsm_Any = true; 786 else 787 private enum D_InlineAsm_Any = false;