1 /**
2 Various window types.
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.window;
8 
9 import std.math;
10 
11 import dplug.core.vec;
12 
13 // Useful reference:
14 // https://en.wikipedia.org/wiki/Window_function
15 // These window generation functions are explicit with the end-points, 
16 // else it would be very easy to get confused and ignorant of these issues.
17 
18 enum WindowType
19 {
20     /// Constant window.
21     rect,
22     
23     /// Triangular window. You probably want something better.
24     bartlett,
25 
26     /// Good default choice if you lack ideas.
27     hann,
28 
29     hamming,
30     
31     nuttall,
32     
33     blackmannNutall, 
34     
35     blackmannHarris,
36 
37     flatTopSR785,    // Flat top window
38 
39     /// Kaiser-Bessel window, this one need a parameter alpha (typical values: 1.0 to 4.0)
40     /// Note: Kaiser windows are talked about in terms of parameters alpha or beta.
41     kaiserBessel,
42 
43     // Older names
44 
45     /// Please use the new name instead.
46     RECT = rect,
47 
48     /// Please use the new name instead.
49     BARTLETT = bartlett,
50 
51     /// Please use the new name instead.
52     HANN = hann,
53 
54     /// Please use the new name instead.
55     HAMMING = hamming,
56 
57     /// Please use the new name instead.
58     NUTTALL = nuttall,
59 
60     /// Please use the new name instead.
61     BLACKMANN_NUTTALL = blackmannNutall,
62 
63     /// Please use the new name instead. 
64     BLACKMANN_HARRIS = blackmannHarris,
65 
66     /// Please use the new name instead.
67     FLATTOP_SR785 = flatTopSR785,
68 
69     /// Please use the new name instead.
70     KAISER_BESSEL = kaiserBessel,
71 }
72 
73 
74 /// How "aligned" the window is in its support.
75 /// Very important, see Issue #236.
76 /// When choosing a window, you probably want to take a hard look about `WindowSymmetry` 
77 /// because it has implications for latency and correctness.
78 enum WindowAlignment
79 {
80     /// The window is asymmetric and if it goes to zero (eg: HANN), it will have one zero 
81     /// at the FIRST coefficient output[N-1].
82     /// Its center is exactly at output[(N/2)-1] which is an integer delay for even window
83     /// lengths.
84     /// This loose one sample of latency, however it has the easiest to compute latency
85     /// in most usage.
86     right,
87 
88     /// The window is asymmetric and if it goes to zero (eg: HANN), it will have one zero 
89     /// at the LAST coefficient output[N-1].
90     /// Its center is exactly at output[(N/2)-1] which is an integer delay for even window
91     /// lengths.
92     /// This has the best latency characteristics, however this might make the latency
93     /// computation itself trickier.
94     left,    
95 
96     /// The window is symmetric and if it goes to zero (eg: HANN), it will have two zero 
97     /// coefficients. 
98     /// Its center is exactly between samples at position output[(N - 1)/2] which is NOT 
99     /// an integer delay for even window lengths.
100     /// Such a window might also break derivativeness.
101     /// However FOR HISTORICAL REASONS THIS WAS THE DEFAULT SETTINGS.
102     /// IT IS ADVISED NOT TO USE IT, FOR THE CONSEQUENCES ARE SUBTLE AND OFTEN HARMFUL.
103     symmetric,
104 }
105 
106 
107 struct WindowDesc
108 {
109 public:
110 nothrow:
111 @nogc:
112 
113     /// Construct a window description, with a foolproof constructor.
114     this(WindowType type, WindowAlignment alignment, float param = float.nan)
115     {
116         this.type = type;
117         this.alignment = alignment;
118         this.param = param;
119     }
120 
121 private:
122     WindowType type;
123 
124     // TODO: this is a bad default! You'll probably want WindowAlignment.right instead.
125     // Make sure it isn't used implicitely anymore.
126     // Then change that default
127     WindowAlignment alignment = WindowAlignment.symmetric; 
128 
129     float param;
130 }
131 
132 /// Generates a window described by `windowDesc`, with periodicity of 
133 /// `outputWindow.length`.
134 void generateWindow(T)(WindowDesc desc, T[] outputWindow) pure nothrow @nogc
135 {
136     int N = cast(int)(outputWindow.length);
137     for (int i = 0; i < N; ++i)
138     {
139         outputWindow[i] = cast(T)(evalWindow(desc, i, N));
140     }
141 }
142 
143 /// Multiplies the given slice in-place by a window described by `windowDesc`,
144 /// whose periodicity is `inoutImpulse.length`.
145 void multiplyByWindow(T)(T[] inoutImpulse, WindowDesc windowDesc) pure nothrow @nogc
146 {
147     int N = cast(int)(inoutImpulse.length);
148     for (int i = 0; i < N; ++i)
149     {
150         inoutImpulse[i] *= evalWindow(windowDesc, i, N);
151     }
152 }
153 
154 // Note: N is the number of output samples, not necessarily the periodicity
155 double evalWindow(WindowDesc desc, int n, int N) pure nothrow @nogc
156 {
157     final switch(desc.alignment) with (WindowAlignment)
158     {
159         case right:
160             return evalWindowInternal(desc, n, N);
161 
162         case left:
163             // left are just a rotation of the right-aligned window
164             if (n == 0)
165                 return evalWindowInternal(desc, N - 1, N);
166             else
167                 return evalWindowInternal(desc, n + 1, N);
168 
169         case symmetric:
170             // Symmetric is just a shorter window
171             return evalWindowInternal(desc, n, N - 1);
172     }
173 }
174 
175 
176 struct Window(T) if (is(T == float) || is(T == double))
177 {
178     void initialize(WindowDesc desc, int lengthInSamples) nothrow @nogc
179     {
180         _lengthInSamples = lengthInSamples;
181         _window.reallocBuffer(lengthInSamples);
182         generateWindow!T(desc, _window);
183     }
184 
185     ~this() nothrow @nogc
186     {
187         _window.reallocBuffer(0);
188     }
189 
190     double sumOfWindowSamples() pure const nothrow @nogc
191     {
192         double result = 0;
193         foreach(windowValue; _window)
194             result += windowValue;
195         return result;
196     }
197 
198     @disable this(this);
199 
200     T[] _window = null;
201     int _lengthInSamples;
202     alias _window this;
203 }
204 
205 
206 private:
207 
208 // Generates WindowAlignment.right window types.
209 // N is the periodicity.
210 // Does NOT look at desc.alignment, which is taken care of in `evalWindow`.
211 double evalWindowInternal(WindowDesc desc, int n, int N) pure nothrow @nogc
212 {
213     static double computeKaiserFunction(double alpha, int n, int N) pure nothrow @nogc
214     {
215         static double I0(double x) pure nothrow @nogc
216         {
217             double sum = 1;
218             double mx = x * 0.5;
219             double denom = 1;
220             double numer = 1;
221             for (int m = 1; m <= 32; m++) 
222             {
223                 numer *= mx;
224                 denom *= m;
225                 double term = numer / denom;
226                 sum += term * term;
227             }
228             return sum;
229         }
230 
231         double piAlpha = PI * alpha;
232         double C = 2.0 * n / cast(double)N - 1.0; 
233         double result = I0(piAlpha * sqrt(1.0 - C * C)) / I0(piAlpha);
234         return result;
235     }
236 
237     final switch(desc.type)
238     {
239         case WindowType.RECT:
240             return 1.0;
241 
242         case WindowType.BARTLETT:
243         {
244             double nm1 = cast(double)N / 2;
245             return 1 - abs(n - nm1) / nm1;
246         }
247 
248         case WindowType.HANN:
249         {
250             double phi = (2 * PI * n ) / N;
251             return 0.5 - 0.5 * cos(phi);
252         }
253 
254         case WindowType.HAMMING:
255         {
256             double phi = (2 * PI * n ) / N;
257             return 0.54 - 0.46 * cos(phi);
258         }
259 
260         case WindowType.NUTTALL:
261         {
262             double phi = (2 * PI * n ) / N;
263             return 0.355768 - 0.487396 * cos(phi) + 0.144232 * cos(2 * phi) - 0.012604 * cos(3 * phi);            
264         }
265 
266         case WindowType.BLACKMANN_HARRIS:
267         {
268             double phi = (2 * PI * n ) / N;
269             return 0.35875 - 0.48829 * cos(phi) + 0.14128 * cos(2 * phi) - 0.01168 * cos(3 * phi);
270         }
271 
272         case WindowType.BLACKMANN_NUTTALL:
273         {
274             double phi = (2 * PI * n ) / N;
275             return 0.3635819 - 0.4891775 * cos(phi) + 0.1365995 * cos(2 * phi) - 0.0106411 * cos(3 * phi);
276         }
277 
278         case WindowType.FLATTOP_SR785:
279         {
280             double phi = (2 * PI * n ) / cast(double)N;
281             return 1 - 1.93 * cos(phi) 
282                      + 1.29 * cos(2 * phi) 
283                      - 0.388 * cos(3 * phi) 
284                      + 0.028 * cos(4 * phi);
285         }
286 
287         case WindowType.KAISER_BESSEL:
288         {
289             return computeKaiserFunction(desc.param, n, N);
290         }
291     }
292 }
293 
294 
295 // Tests the general shape of windows.
296 unittest
297 {
298     void checkIsSymmetric(T)(T[] inp)
299     {
300         for (int i = 0; i < inp.length/2; ++i)
301         {
302             // Window generation should be precise
303             assert(isClose(inp[i], inp[$-1-i], 1e-10));
304         }
305     }
306 
307     void checkMaxInCenter(T)(T[] inp)
308     {
309         bool odd = (inp.length & 1) != 0;
310         if (odd)
311         {
312             for (int i = 0; i < inp.length/2; ++i)
313             {
314                 assert(inp[i] <= inp[inp.length/2]);
315                 assert(inp[$-1-i] <= inp[inp.length/2]);
316             }
317         }
318         else
319         {
320             for (int i = 0; i < inp.length/2-1; ++i)
321             {
322                 // take either of the centers
323                 assert(inp[i] <= inp[inp.length/2-1]);
324                 assert(inp[$-1-i] <= inp[inp.length/2-1]);
325             }
326         }
327     }
328 
329     void testAllWindows(T)(int size)
330     {
331         Window!T window;
332 
333         foreach(type; WindowType.min..WindowType.max+1)
334         {
335             foreach(alignment; WindowAlignment.min..WindowAlignment.max+1)
336             {
337                 // Note: only Kaiser-Bessel have a parameter for now, so take 2
338                 WindowDesc desc = WindowDesc(cast(WindowType)type, 
339                                              cast(WindowAlignment)alignment, 
340                                              2.0);
341 
342                 window.initialize(desc, size);
343 
344                 final switch(alignment) with (WindowAlignment)
345                 {
346                     case WindowAlignment.symmetric:
347                         checkIsSymmetric(window[0..$]);
348                         checkMaxInCenter(window[0..$]);
349                         break;
350 
351                     case WindowAlignment.right:
352                         checkIsSymmetric(window[1..$]);
353                         checkMaxInCenter(window[1..$]);
354                         break;
355 
356                     case WindowAlignment.left:
357                         checkIsSymmetric(window[0..$-1]);
358                         checkMaxInCenter(window[0..$-1]);
359                         break;
360                 }
361             }
362         }
363     }
364 
365     testAllWindows!float(64);
366     testAllWindows!double(32);
367 
368     // It should be possible to generate odd-sized window, 
369     // though in practice with WindowAlignement you don't need them.
370     testAllWindows!float(5);
371 }