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 }