1 /** 2 Various DSP smoothers. 3 4 Copyright: Guillaume Piolats 2015-2022. 5 License: $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0) 6 */ 7 module dplug.dsp.smooth; 8 9 import std.algorithm.comparison; 10 import std.math; 11 12 import dplug.core.math; 13 import dplug.core.ringbuf; 14 import dplug.core.nogc; 15 import dplug.core.vec; 16 17 /// Smooth values exponentially with a 1-pole lowpass. 18 /// This is usually sufficient for most parameter smoothing. 19 deprecated("ExpSmoother will be removed as of Dplug v14") 20 struct ExpSmoother(T) if (is(T == float) || is(T == double)) 21 { 22 public: 23 /// time: the time constant of the smoother. 24 /// threshold: absolute difference below which we consider current value and target equal 25 void initialize(float samplerate, float timeAttackRelease, T initialValue) nothrow @nogc 26 { 27 assert(isFinite(initialValue)); 28 29 _current = cast(T)(initialValue); 30 _sampleRate = samplerate; 31 32 setAttackReleaseTime(timeAttackRelease); 33 34 assert(isFinite(_expFactor)); 35 } 36 37 /// Changes attack and release time (given in seconds). 38 void setAttackReleaseTime(float timeAttackRelease) nothrow @nogc 39 { 40 _expFactor = cast(T)(expDecayFactor(timeAttackRelease, _sampleRate)); 41 } 42 43 /// Advance smoothing and return the next smoothed sample with respect 44 /// to tau time and samplerate. 45 T nextSample(T target) nothrow @nogc 46 { 47 T diff = target - _current; 48 if (diff != 0) 49 { 50 if (fast_fabs(diff) < 1e-10f) // to avoid subnormal, and excess churn 51 { 52 _current = target; 53 } 54 else 55 { 56 double temp = _current + diff * _expFactor; // Is double-precision really needed here? 57 T newCurrent = cast(T)(temp); 58 _current = newCurrent; 59 } 60 } 61 return _current; 62 } 63 64 bool hasConverged(T target) nothrow @nogc 65 { 66 return target == _current; 67 } 68 69 void nextBuffer(const(T)* input, T* output, int frames) nothrow @nogc 70 { 71 for (int i = 0; i < frames; ++i) 72 { 73 output[i] = nextSample(input[i]); 74 } 75 } 76 77 void nextBuffer(T input, T* output, int frames) nothrow @nogc 78 { 79 for (int i = 0; i < frames; ++i) 80 { 81 output[i] = nextSample(input); 82 } 83 } 84 85 private: 86 T _current; 87 T _expFactor; 88 float _sampleRate; 89 } 90 91 92 /// Same as ExpSmoother but have different attack and release decay factors. 93 deprecated("AttackReleaseSmoother will be removed as of Dplug v14") 94 struct AttackReleaseSmoother(T) if (is(T == float) || is(T == double)) 95 { 96 public: 97 /// time: the time constant of the smoother. 98 /// threshold: absolute difference below which we consider current value and target equal 99 void initialize(float sampleRate, float timeAttackSecs, float timeReleaseSecs, T initialValue) nothrow @nogc 100 { 101 assert(isFinite(initialValue)); 102 _sampleRate = sampleRate; 103 _current = cast(T)(initialValue); 104 setAttackTime(timeAttackSecs); 105 setReleaseTime(timeReleaseSecs); 106 } 107 108 /// Changes attack time (given in seconds). 109 void setAttackTime(float timeAttackSecs) nothrow @nogc 110 { 111 _expFactorAttack = cast(T)(expDecayFactor(timeAttackSecs, _sampleRate)); 112 } 113 114 /// Changes release time (given in seconds). 115 void setReleaseTime(float timeReleaseSecs) nothrow @nogc 116 { 117 _expFactorRelease = cast(T)(expDecayFactor(timeReleaseSecs, _sampleRate)); 118 } 119 120 /// Advance smoothing and return the next smoothed sample with respect 121 /// to tau time and samplerate. 122 T nextSample(T target) nothrow @nogc 123 { 124 T diff = target - _current; 125 if (diff != 0) 126 { 127 if (fast_fabs(diff) < 1e-10f) // to avoid subnormal, and excess churn 128 { 129 _current = target; 130 } 131 else 132 { 133 double expFactor = (diff > 0) ? _expFactorAttack : _expFactorRelease; 134 double temp = _current + diff * expFactor; // Is double-precision really needed here? 135 T newCurrent = cast(T)(temp); 136 _current = newCurrent; 137 } 138 } 139 return _current; 140 } 141 142 void nextBuffer(const(T)* input, T* output, int frames) 143 { 144 for (int i = 0; i < frames; ++i) 145 { 146 output[i] = nextSample(input[i]); 147 } 148 } 149 150 void nextBuffer(T input, T* output, int frames) 151 { 152 for (int i = 0; i < frames; ++i) 153 { 154 output[i] = nextSample(input); 155 } 156 } 157 158 private: 159 T _current; 160 T _expFactorAttack; 161 T _expFactorRelease; 162 float _sampleRate; 163 } 164 165 166 /// Non-linear smoother using absolute difference. 167 /// Designed to have a nice phase response. 168 /// Warning: samplerate-dependent. 169 deprecated("AbsSmoother will be removed as of Dplug v14") 170 struct AbsSmoother(T) if (is(T == float) || is(T == double)) 171 { 172 public: 173 174 /// Initialize the AbsSmoother. 175 /// maxAbsDiff: maximum difference between filtered consecutive samples 176 void initialize(T initialValue, T maxAbsDiff) nothrow @nogc 177 { 178 assert(isFinite(initialValue)); 179 _maxAbsDiff = maxAbsDiff; 180 _current = initialValue; 181 } 182 183 T nextSample(T input) nothrow @nogc 184 { 185 T absDiff = abs(input - _current); 186 if (absDiff <= _maxAbsDiff) 187 _current = input; 188 else 189 _current = _current + absDiff * (input > _current ? 1 : -1); 190 return _current; 191 } 192 193 void nextBuffer(const(T)* input, T* output, int frames) nothrow @nogc 194 { 195 for(int i = 0; i < frames; ++i) 196 output[i] = nextSample(input[i]); 197 } 198 199 private: 200 T _current; 201 T _maxAbsDiff; 202 } 203 204 /// Smooth values over time with a linear slope. 205 /// This can be useful for some smoothing needs. 206 /// Intermediate between fast phase and actual smoothing. 207 deprecated("LinearSmoother will be removed as of Dplug v14") 208 struct LinearSmoother(T) if (is(T == float) || is(T == double)) 209 { 210 public: 211 212 /// Initialize the LinearSmoother. 213 void initialize(T initialValue, float periodSecs, float sampleRate) nothrow @nogc 214 { 215 _period = periodSecs; 216 _periodInv = 1 / periodSecs; 217 _sampleRateInv = 1 / sampleRate; 218 219 // clear state 220 _current = initialValue; 221 _phase = 0; 222 _firstNextAfterInit = true; 223 } 224 225 /// Set the target value and return the next sample. 226 T nextSample(T input) nothrow @nogc 227 { 228 _phase += _sampleRateInv; 229 if (_firstNextAfterInit || _phase > _period) 230 { 231 _phase -= _period; 232 _increment = (input - _current) * (_sampleRateInv * _periodInv); 233 _firstNextAfterInit = false; 234 } 235 _current += _increment; 236 return _current; 237 } 238 239 void nextBuffer(const(T)* input, T* output, int frames) nothrow @nogc 240 { 241 for(int i = 0; i < frames; ++i) 242 output[i] = nextSample(input[i]); 243 } 244 245 private: 246 T _current; 247 T _increment; 248 float _period; 249 float _periodInv; 250 float _sampleRateInv; 251 float _phase; 252 bool _firstNextAfterInit; 253 } 254 255 /// Can be very useful when filtering values with outliers. 256 /// For what it's meant to do, excellent phase response. 257 deprecated("MedianFilter will be removed as of Dplug v14") 258 struct MedianFilter(T) if (is(T == float) || is(T == double)) 259 { 260 public: 261 262 void initialize(T initialValue, int samples) nothrow @nogc 263 { 264 assert(samples >= 2, "N must be >= 2"); 265 assert(samples % 2 == 1, "N must be odd"); 266 267 _delay.reallocBuffer(samples - 1); 268 _delay[] = initialValue; 269 270 _arr.reallocBuffer(samples); 271 _N = samples; 272 273 _tempbuf = makeVec!T(); 274 } 275 276 ~this() 277 { 278 _delay.reallocBuffer(0); 279 _arr.reallocBuffer(0); 280 } 281 282 T nextSample(T input) nothrow @nogc 283 { 284 // dramatically inefficient 285 286 _arr[0] = input; 287 for (int i = 0; i < _N - 1; ++i) 288 _arr[i + 1] = _delay[i]; 289 290 // sort in place 291 // this sort doesn't even need to be stable, but well. 292 timSort!T(_arr[], 293 _tempbuf, 294 (a, b) nothrow @nogc 295 { 296 if (a > b) return 1; 297 else if (a < b) return -1; 298 else return 0; 299 } 300 ); 301 302 T median = _arr[_N/2]; 303 304 for (int i = _N - 3; i >= 0; --i) 305 _delay[i + 1] = _delay[i]; 306 _delay[0] = input; 307 return median; 308 } 309 310 void nextBuffer(const(T)* input, T* output, int frames) nothrow @nogc 311 { 312 for(int i = 0; i < frames; ++i) 313 output[i] = nextSample(input[i]); 314 } 315 316 private: 317 T[] _delay; 318 T[] _arr; 319 int _N; 320 Vec!T _tempbuf; 321 } 322 323 /// Simple FIR to smooth things cheaply. 324 /// Introduces (samples - 1) / 2 latency. 325 /// This one doesn't convert to integers internally so it may 326 /// loose precision over time. Meants for finite signals. 327 deprecated("UnstableMeanFilter will be removed as of Dplug v14") 328 struct UnstableMeanFilter(T) if (is(T == float) || is(T == double)) 329 { 330 public: 331 /// Initialize mean filter with given number of samples. 332 void initialize(T initialValue, int samples) nothrow @nogc 333 { 334 _delay = RingBufferNoGC!T(samples); 335 336 _invNFactor = cast(T)1 / samples; 337 338 while(!_delay.isFull()) 339 _delay.pushBack(initialValue); 340 341 _sum = _delay.length * initialValue; 342 } 343 344 /// Initialize with with cutoff frequency and samplerate. 345 void initialize(T initialValue, double cutoffHz, double samplerate) nothrow @nogc 346 { 347 int nSamples = cast(int)(0.5 + samplerate / (2 * cutoffHz)); 348 349 if (nSamples < 1) 350 nSamples = 1; 351 352 initialize(initialValue, nSamples); 353 } 354 355 // process next sample 356 T nextSample(T x) nothrow @nogc 357 { 358 _sum = _sum + x; 359 _sum = _sum - _delay.popFront(); 360 _delay.pushBack(x); 361 return _sum * _invNFactor; 362 } 363 364 void nextBuffer(const(T)* input, T* output, int frames) nothrow @nogc 365 { 366 for(int i = 0; i < frames; ++i) 367 output[i] = nextSample(input[i]); 368 } 369 370 private: 371 RingBufferNoGC!T _delay; 372 double _sum; // should be approximately the sum of samples in delay 373 T _invNFactor; 374 T _factor; 375 } 376 377 /// Simple FIR to smooth things cheaply. 378 /// Introduces (samples - 1) / 2 latency. 379 /// Converts everything to long for stability purpose. 380 /// So this may run forever as long as the input is below some threshold. 381 deprecated("MeanFilter will be removed as of Dplug v14") 382 struct MeanFilter(T) if (is(T == float) || is(T == double)) 383 { 384 public: 385 /// Initialize mean filter with given number of samples. 386 void initialize(T initialValue, int samples, T maxExpectedValue) nothrow @nogc 387 { 388 _delay = RingBufferNoGC!long(samples); 389 390 _factor = cast(T)(2147483648.0 / maxExpectedValue); 391 _invNFactor = cast(T)1 / (_factor * samples); 392 393 // clear state 394 // round to integer 395 long ivInt = toIntDomain(initialValue); 396 397 while(!_delay.isFull()) 398 _delay.pushBack(ivInt); 399 400 _sum = cast(int)(_delay.length) * ivInt; 401 } 402 403 /// Initialize with with cutoff frequency and samplerate. 404 void initialize(T initialValue, double cutoffHz, double samplerate, T maxExpectedValue) nothrow @nogc 405 { 406 int nSamples = cast(int)(0.5 + samplerate / (2 * cutoffHz)); 407 408 if (nSamples < 1) 409 nSamples = 1; 410 411 initialize(initialValue, nSamples, maxExpectedValue); 412 } 413 414 // process next sample 415 T nextSample(T x) nothrow @nogc 416 { 417 // round to integer 418 long input = cast(long)(cast(T)0.5 + x * _factor); 419 _sum = _sum + input; 420 _sum = _sum - _delay.popFront(); 421 _delay.pushBack(input); 422 return cast(T)_sum * _invNFactor; 423 } 424 425 void nextBuffer(const(T)* input, T* output, int frames) nothrow @nogc 426 { 427 for(int i = 0; i < frames; ++i) 428 output[i] = nextSample(input[i]); 429 } 430 431 private: 432 433 long toIntDomain(T x) pure const nothrow @nogc 434 { 435 return cast(long)(cast(T)0.5 + x * _factor); 436 } 437 438 RingBufferNoGC!long _delay; 439 long _sum; // should always be the sum of samples in delay 440 T _invNFactor; 441 T _factor; 442 } 443