1 /* 2 * Copyright (C) 2016 The Android Open Source Project 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 * 16 * This code was translated from the JSyn Java code. 17 * JSyn is Copyright 2009 Phil Burk, Mobileer Inc 18 * JSyn is licensed under the Apache License, Version 2.0 19 */ 20 21 #ifndef SYNTHMARK_BIQUAD_FILTER_H 22 #define SYNTHMARK_BIQUAD_FILTER_H 23 24 #include <cstdint> 25 #include <math.h> 26 #include "SynthTools.h" 27 #include "UnitGenerator.h" 28 29 namespace marksynth { 30 31 #define BIQUAD_MIN_FREQ (0.00001f) // REVIEW 32 #define BIQUAD_MIN_Q (0.00001f) // REVIEW 33 34 #define RECALCULATE_PER_SAMPLE 0 35 36 /** 37 * Time varying lowpass resonant filter. 38 */ 39 class BiquadFilter : public UnitGenerator 40 { 41 public: BiquadFilter()42 BiquadFilter() 43 : mQ(1.0) 44 { 45 xn1 = xn2 = yn1 = yn2 = (synth_float_t) 0; 46 a0 = a1 = a2 = b1 = b2 = (synth_float_t) 0; 47 } 48 49 virtual ~BiquadFilter() = default; 50 51 /** 52 * Resonance, typically between 1.0 and 10.0. 53 * Input will clipped at a BIQUAD_MIN_Q. 54 */ setQ(synth_float_t q)55 void setQ(synth_float_t q) { 56 if( q < BIQUAD_MIN_Q ) { 57 q = BIQUAD_MIN_Q; 58 } 59 mQ = q; 60 } 61 getQ()62 synth_float_t getQ() { 63 return mQ; 64 } 65 generate(synth_float_t * input,synth_float_t * frequencies,int32_t numSamples)66 void generate(synth_float_t *input, 67 synth_float_t *frequencies, 68 int32_t numSamples) { 69 synth_float_t xn, yn; 70 71 #if RECALCULATE_PER_SAMPLE == 0 72 calculateCoefficients(frequencies[0], mQ); 73 #endif 74 for (int i = 0; i < numSamples; i++) { 75 #if RECALCULATE_PER_SAMPLE == 1 76 calculateCoefficients(frequencies[i], mQ); 77 #endif 78 // Generate outputs by filtering inputs. 79 xn = input[i]; 80 synth_float_t finite = (a0 * xn) + (a1 * xn1) + (a2 * xn2); 81 // Use double precision for recursive portion. 82 yn = finite - (b1 * yn1) - (b2 * yn2); 83 output[i] = (synth_float_t) yn; 84 85 // Delay input and output values. 86 xn2 = xn1; 87 xn1 = xn; 88 yn2 = yn1; 89 yn1 = yn; 90 } 91 92 // Apply a small bipolar impulse to filter to prevent arithmetic underflow. 93 yn1 += (synth_float_t) 1.0E-26; 94 yn2 -= (synth_float_t) 1.0E-26; 95 } 96 97 98 private: 99 synth_float_t mQ; 100 101 synth_float_t xn1; // delay lines 102 synth_float_t xn2; 103 double yn1; 104 double yn2; 105 106 synth_float_t a0; // coefficients 107 synth_float_t a1; 108 synth_float_t a2; 109 110 synth_float_t b1; 111 synth_float_t b2; 112 113 synth_float_t cos_omega; 114 synth_float_t sin_omega; 115 synth_float_t alpha; 116 117 // Calculate coefficients common to many parametric biquad filters. calcCommon(synth_float_t ratio,synth_float_t Q)118 void calcCommon( synth_float_t ratio, synth_float_t Q ) 119 { 120 synth_float_t omega; 121 122 /* Don't let frequency get too close to Nyquist or filter will blow up. */ 123 if( ratio >= 0.499f ) ratio = 0.499f; 124 omega = 2.0f * (synth_float_t)M_PI * ratio; 125 126 #if 1 127 // This is not significantly faster on Mac or Linux. 128 cos_omega = SynthTools::fastCosine(omega); 129 sin_omega = SynthTools::fastSine(omega ); 130 #else 131 { 132 float fsin_omega; 133 float fcos_omega; 134 sincosf(omega, &fsin_omega, &fcos_omega); 135 cos_omega = (synth_float_t) fcos_omega; 136 sin_omega = (synth_float_t) fsin_omega; 137 } 138 #endif 139 alpha = sin_omega / (2.0f * Q); 140 } 141 142 // Lowpass coefficients calculateCoefficients(synth_float_t frequency,synth_float_t Q)143 void calculateCoefficients( synth_float_t frequency, synth_float_t Q ) 144 { 145 synth_float_t scalar, omc; 146 147 if( frequency < BIQUAD_MIN_FREQ ) frequency = BIQUAD_MIN_FREQ; 148 149 calcCommon( frequency * mSamplePeriod, Q ); 150 151 scalar = 1.0f / (1.0f + alpha); 152 omc = (1.0f - cos_omega); 153 154 a0 = omc * 0.5f * scalar; 155 a1 = omc * scalar; 156 a2 = a0; 157 b1 = -2.0f * cos_omega * scalar; 158 b2 = (1.0f - alpha) * scalar; 159 } 160 }; 161 162 }; 163 #endif // SYNTHMARK_BIQUAD_FILTER_H 164