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 }