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 }