dplug.dsp.window source code

1 /**
2 * Copyright: Copyright Auburn Sounds 2015 and later.
3 * License:   $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0)
4 * Authors:   Guillaume Piolat
5 */
6 module dplug.dsp.window;
7 
8 import std.math;
9 
10 import dplug.core.alignedbuffer;
11 
12 // Useful reference:
13 // https://en.wikipedia.org/wiki/Window_function
14 
15 
16 enum WindowType
17 {
18     RECT,
19     BARTLETT,
20     HANN,
21     HAMMING,
22     NUTTALL,
23     BLACKMANN_NUTTALL, 
24     BLACKMANN_HARRIS,
25     FLATTOP_SR785,    // Flat top window
26     KAISER_BESSEL,   // Kaiser-Bessel window, this one need a parameter alpha (typical values: 1.0 to 4.0)
27 }
28 
29 struct WindowDesc
30 {
31     WindowType type;
32     float param;
33 }
34 
35 void generateWindow(T)(WindowDesc desc, T[] output) pure nothrow @nogc
36 {
37     int N = cast(int)(output.length);
38     for (int i = 0; i < N; ++i)
39     {
40         output[i] = cast(T)(evalWindow(desc, i, N));
41     }
42 }
43 
44 void generateNormalizedWindow(T)(WindowDesc desc, T[] output) pure nothrow @nogc
45 {
46     int N = cast(int)(output.length);
47     T sum = 0;
48     for (int i = 0; i < N; ++i)
49     {
50         output[i] = cast(T)(evalWindow(desc, i, N));
51         sum += output[i];
52     }
53     T invSum = 1 / sum;
54     for (int i = 0; i < N; ++i)
55         output[i] *= invSum;    
56 }
57 
58 deprecated("This function disappeared.") double secondaryLobeAttenuationInDb(WindowType type) pure nothrow @nogc
59 {
60    return double.nan;
61 }
62 
63 double evalWindow(WindowDesc desc, int n, int N) pure nothrow @nogc
64 {
65     static double computeKaiserFunction(double alpha, int n, int N) pure nothrow @nogc
66     {
67         static double I0(double x) pure nothrow @nogc
68         {
69             double sum = 1;
70             double mx = x * 0.5;
71             double denom = 1;
72             double numer = 1;
73             for (int m = 1; m <= 32; m++) 
74             {
75                 numer *= mx;
76                 denom *= m;
77                 double term = numer / denom;
78                 sum += term * term;
79             }
80             return sum;
81         }
82 
83         double piAlpha = PI * alpha;
84         double C = 2.0 * n / (N - 1.0) - 1.0; 
85         double result = I0(piAlpha * sqrt(1.0 - C * C)) / I0(piAlpha);
86         return result;
87     }
88 
89     final switch(desc.type)
90     {
91         case WindowType.RECT:
92             return 1.0;
93 
94         case WindowType.BARTLETT:
95         {
96             double nm1 = (N - 1.0)/2;
97             return 1 - abs(n - nm1) / nm1;
98         }
99 
100         case WindowType.HANN:
101         {
102             double phi = (2 * PI * n ) / (N - 1);
103             return 0.5 - 0.5 * cos(phi);
104         }
105 
106         case WindowType.HAMMING:
107         {
108             double phi = (2 * PI * n ) / (N - 1);
109             return 0.54 - 0.46 * cos(phi);
110         }
111 
112         case WindowType.NUTTALL:
113         {
114             double phi = (2 * PI * n ) / (N - 1);
115             return 0.355768 - 0.487396 * cos(phi) + 0.144232 * cos(2 * phi) - 0.012604 * cos(3 * phi);            
116         }
117 
118         case WindowType.BLACKMANN_HARRIS:
119         {
120             double phi = (2 * PI * n ) / (N - 1);
121             return 0.35875 - 0.48829 * cos(phi) + 0.14128 * cos(2 * phi) - 0.01168 * cos(3 * phi);
122         }
123 
124         case WindowType.BLACKMANN_NUTTALL:
125         {
126             double phi = (2 * PI * n ) / (N - 1);
127             return 0.3635819 - 0.4891775 * cos(phi) + 0.1365995 * cos(2 * phi) - 0.0106411 * cos(3 * phi);
128         }
129 
130         case WindowType.FLATTOP_SR785:
131         {
132             double phi = (2 * PI * n ) / (N - 1);
133             return 1 - 1.93 * cos(phi) + 1.29 * cos(2 * phi) - 0.388 * cos(3 * phi) + 0.028 * cos(4 * phi);
134         }
135 
136         case WindowType.KAISER_BESSEL:
137             return computeKaiserFunction(desc.param, n, N);
138     }
139 }
140 
141 struct Window(T) if (is(T == float) || is(T == double))
142 {
143     void initialize(WindowDesc desc, int lengthInSamples) nothrow @nogc
144     {
145         _lengthInSamples = lengthInSamples;
146         _window.reallocBuffer(lengthInSamples);
147         generateWindow!T(desc, _window);
148     }
149 
150     ~this() nothrow @nogc
151     {
152         _window.reallocBuffer(0);
153     }
154 
155     double sumOfWindowSamples() pure const nothrow @nogc
156     {
157         double result = 0;
158         foreach(windowValue; _window)
159             result += windowValue;
160         return result;
161     }
162 
163     @disable this(this);
164 
165     T[] _window = null;
166     int _lengthInSamples;
167     alias _window this;
168 }