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; 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 deprecated and 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 void generateLowpassImpulse(T)(T[] output, double cutoff, double samplerate) nothrow @nogc 26 { 27 checkFilterParams(output.length, cutoff, samplerate); 28 double cutoffNormalized = cutoff / samplerate; 29 double wc = cutoffNormalized * 2.0 * PI; 30 31 int len = cast(int)(output.length); 32 for (int i = 0; i < len; ++i) 33 { 34 int x = i - (len / 2); 35 if (x == 0) 36 output[i] = cutoffNormalized * 2.0; 37 else 38 output[i] = sin(wc * x) / cast(double)(PI * x); 39 } 40 } 41 42 /// Generates a sinc highpass impulse, centered on floor(output.length / 2). 43 /// When convolved with, preserve amplitude of the pass-band. 44 void generateHighpassImpulse(T)(T[] output, double cutoff, double samplerate) nothrow @nogc 45 { 46 checkFilterParams(output.length, cutoff, samplerate); 47 double cutoffNormalized = cutoff / samplerate; 48 double wc = cutoffNormalized * 2.0 * PI; 49 50 int len = cast(int)(output.length); 51 for (int i = 0; i < len; ++i) 52 { 53 int x = i - (len / 2); 54 if (x == 0) 55 output[i] = 1.0 - 2 * cutoffNormalized; 56 else 57 output[i] = sinc(PI * x) / cast(double)(PI * x) - 2.0 * cutoffNormalized * sin(wc * x) / (wc * x); 58 } 59 } 60 61 /// Generates a hilbert transformer impulse, centered on floor(output.length / 2). 62 void generateHilbertTransformer(T)(T[] outImpulse, WindowDesc windowDesc, double samplerate) nothrow @nogc 63 { 64 int size = cast(int)(outImpulse.length); 65 assert(isOdd(size)); 66 int center = size / 2; 67 68 for (int i = 0; i < center; ++i) 69 { 70 int xi = i - center; 71 double x = cast(double)xi; 72 if (isEven(xi)) 73 { 74 outImpulse[i] = 0; 75 outImpulse[size - 1 - i] = 0; 76 } 77 else 78 { 79 double y = x * cast(double)PI / 2; 80 double sine = sin(y); 81 T value = cast(T)(-sine*sine / y); 82 value *= evalWindow(windowDesc, i, size); 83 outImpulse[i] = value; 84 outImpulse[size - 1 - i] = -value; 85 } 86 } 87 outImpulse[center] = 0; 88 } 89 90 private static void checkFilterParams(size_t length, double cutoff, double sampleRate) nothrow @nogc 91 { 92 assert((length & 1) == 0, "FIR impulse length must be even"); 93 assert(cutoff * 2 < sampleRate, "2x the cutoff exceed sampling rate, Nyquist disapproving"); 94 } 95 96 unittest 97 { 98 double[256] lp_impulse; 99 double[256] hp_impulse; 100 generateLowpassImpulse(lp_impulse[], 40.0, 44100.0); 101 generateHighpassImpulse(hp_impulse[], 40.0, 44100.0); 102 } 103 104 unittest 105 { 106 double[256] lp_impulse; 107 generateHilbertTransformer(lp_impulse[0..$-1], 108 WindowDesc(WindowType.blackmannHarris, 109 WindowAlignment.right), 44100.0); 110 } 111 112 113 // Composed of a delay-line, and an inpulse. 114 struct FIR(T) 115 { 116 /// Initializes the FIR filter. It's up to you to fill the impulse with something worthwhile. 117 void initialize(int sizeOfImpulse) nothrow @nogc 118 { 119 assert(sizeOfImpulse > 0); 120 _delayline.initialize(sizeOfImpulse); 121 _impulse.reallocBuffer(sizeOfImpulse); 122 _impulse[] = T.nan; 123 } 124 125 ~this() nothrow @nogc 126 { 127 _impulse.reallocBuffer(0); 128 _tempBuffer.reallocBuffer(0); 129 _windowBuffer.reallocBuffer(0); 130 } 131 132 @disable this(this); 133 134 /// Returns: Filtered input sample, naive convolution. 135 T nextSample(T input) nothrow @nogc 136 { 137 _delayline.feedSample(input); 138 T sum = 0; 139 int N = cast(int)impulse.length; 140 for (int i = 0; i < N; ++i) 141 sum += _impulse.ptr[i] * _delayline.sampleFull(i); 142 return sum; 143 } 144 145 void nextBuffer(T* input, T* output, int samples) nothrow @nogc 146 { 147 for (int i = 0; i < samples; ++i) 148 { 149 output[i] = nextSample(input[i]); 150 } 151 } 152 153 /// Returns: Impulse response. If you write it, you can call makeMinimumPhase() next. 154 inout(T)[] impulse() inout nothrow @nogc 155 { 156 return _impulse; 157 } 158 159 void applyWindow(WindowDesc windowDesc) nothrow @nogc 160 { 161 _windowBuffer.reallocBuffer(_impulse.length); 162 generateWindow(windowDesc, _windowBuffer); 163 foreach(i; 0.._impulse.length) 164 { 165 _impulse[i] *= _windowBuffer[i]; 166 } 167 } 168 169 private: 170 T[] _impulse; 171 172 Delayline!T _delayline; 173 Complex!T[] _tempBuffer; 174 T[] _windowBuffer; 175 } 176 177 unittest 178 { 179 FIR!double fir; 180 fir.initialize(32); 181 generateLowpassImpulse(fir.impulse(), 40.0, 44100.0); 182 fir.applyWindow(WindowDesc(WindowType.hann, WindowAlignment.right)); 183 }