1 /** 2 Naive FIR implementation. 3 4 Copyright: Guillaume Piolat 2015. 5 License: $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0) 6 */ 7 module dplug.dsp.fir; 8 9 import core.stdc.complex; 10 11 import std.math: PI, sin; 12 import std.complex; 13 14 import dplug.core.math; 15 import dplug.core.vec; 16 import dplug.dsp.delayline; 17 import dplug.dsp.window; 18 19 // FUTURE: should probably be removed, not good enough. 20 // Need a real convolver library. 21 22 // Basic sinc impulse functions 23 24 // Generates a sinc lowpass impulse, centered on floor(output.length / 2). 25 deprecated("This will be removed in Dplug v16") 26 void generateLowpassImpulse(T)(T[] output, double cutoff, double samplerate) nothrow @nogc 27 { 28 checkFilterParams(output.length, cutoff, samplerate); 29 double cutoffNormalized = cutoff / samplerate; 30 double wc = cutoffNormalized * 2.0 * PI; 31 32 int len = cast(int)(output.length); 33 for (int i = 0; i < len; ++i) 34 { 35 int x = i - (len / 2); 36 if (x == 0) 37 output[i] = cutoffNormalized * 2.0; 38 else 39 output[i] = sin(wc * x) / cast(double)(PI * x); 40 } 41 } 42 43 // Generates a sinc highpass impulse, centered on floor(output.length / 2). 44 // When convolved with, preserve amplitude of the pass-band. 45 deprecated("This will be removed in Dplug v16") 46 void generateHighpassImpulse(T)(T[] output, double cutoff, double samplerate) nothrow @nogc 47 { 48 checkFilterParams(output.length, cutoff, samplerate); 49 double cutoffNormalized = cutoff / samplerate; 50 double wc = cutoffNormalized * 2.0 * PI; 51 52 int len = cast(int)(output.length); 53 for (int i = 0; i < len; ++i) 54 { 55 int x = i - (len / 2); 56 if (x == 0) 57 output[i] = 1.0 - 2 * cutoffNormalized; 58 else 59 output[i] = sinc(PI * x) / cast(double)(PI * x) - 2.0 * cutoffNormalized * sin(wc * x) / (wc * x); 60 } 61 } 62 63 // Generates a hilbert transformer impulse, centered on floor(output.length / 2). 64 deprecated("This will be removed in Dplug v16") 65 void generateHilbertTransformer(T)(T[] outImpulse, WindowDesc windowDesc, double samplerate) nothrow @nogc 66 { 67 static bool isOdd(int i) pure nothrow @nogc @safe 68 { 69 return (i & 1) != 0; 70 } 71 72 static bool isEven(int i) pure nothrow @nogc @safe 73 { 74 return (i & 1) == 0; 75 } 76 77 int size = cast(int)(outImpulse.length); 78 assert(isOdd(size)); 79 int center = size / 2; 80 81 for (int i = 0; i < center; ++i) 82 { 83 int xi = i - center; 84 double x = cast(double)xi; 85 if (isEven(xi)) 86 { 87 outImpulse[i] = 0; 88 outImpulse[size - 1 - i] = 0; 89 } 90 else 91 { 92 double y = x * cast(double)PI / 2; 93 double sine = sin(y); 94 T value = cast(T)(-sine*sine / y); 95 value *= evalWindow(windowDesc, i, size); 96 outImpulse[i] = value; 97 outImpulse[size - 1 - i] = -value; 98 } 99 } 100 outImpulse[center] = 0; 101 } 102 103 private static void checkFilterParams(size_t length, double cutoff, double sampleRate) nothrow @nogc 104 { 105 assert((length & 1) == 0, "FIR impulse length must be even"); 106 assert(cutoff * 2 < sampleRate, "2x the cutoff exceed sampling rate, Nyquist disapproving"); 107 } 108 109 // Composed of a delay-line, and an inpulse. It's a bad convolver. 110 deprecated("This will be removed in Dplug v16") struct FIR(T) 111 { 112 /// Initializes the FIR filter. It's up to you to fill the impulse with something worthwhile. 113 void initialize(int sizeOfImpulse) nothrow @nogc 114 { 115 assert(sizeOfImpulse > 0); 116 _delayline.initialize(sizeOfImpulse); 117 _impulse.reallocBuffer(sizeOfImpulse); 118 _impulse[] = T.nan; 119 } 120 121 ~this() nothrow @nogc 122 { 123 _impulse.reallocBuffer(0); 124 _windowBuffer.reallocBuffer(0); 125 } 126 127 @disable this(this); 128 129 /// Returns: Filtered input sample, naive convolution. 130 T nextSample(T input) nothrow @nogc 131 { 132 _delayline.feedSample(input); 133 T sum = 0; 134 int N = cast(int)impulse.length; 135 for (int i = 0; i < N; ++i) 136 sum += _impulse.ptr[i] * _delayline.sampleFull(i); 137 return sum; 138 } 139 140 void nextBuffer(T* input, T* output, int samples) nothrow @nogc 141 { 142 for (int i = 0; i < samples; ++i) 143 { 144 output[i] = nextSample(input[i]); 145 } 146 } 147 148 /// Returns: Impulse response. If you write it, you can call makeMinimumPhase() next. 149 inout(T)[] impulse() inout nothrow @nogc 150 { 151 return _impulse; 152 } 153 154 void applyWindow(WindowDesc windowDesc) nothrow @nogc 155 { 156 _windowBuffer.reallocBuffer(_impulse.length); 157 generateWindow(windowDesc, _windowBuffer); 158 foreach(i; 0.._impulse.length) 159 { 160 _impulse[i] *= _windowBuffer[i]; 161 } 162 } 163 164 private: 165 T[] _impulse; 166 167 Delayline!T _delayline; 168 T[] _windowBuffer; 169 }