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