OpenShot Audio Library | OpenShotAudio  0.3.0
juce_LagrangeInterpolator.cpp
1 /*
2  ==============================================================================
3 
4  This file is part of the JUCE library.
5  Copyright (c) 2017 - ROLI Ltd.
6 
7  JUCE is an open source library subject to commercial or open-source
8  licensing.
9 
10  The code included in this file is provided under the terms of the ISC license
11  http://www.isc.org/downloads/software-support-policy/isc-license. Permission
12  To use, copy, modify, and/or distribute this software for any purpose with or
13  without fee is hereby granted provided that the above copyright notice and
14  this permission notice appear in all copies.
15 
16  JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
17  EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
18  DISCLAIMED.
19 
20  ==============================================================================
21 */
22 
23 namespace juce
24 {
25 
26 namespace
27 {
28  static forcedinline void pushInterpolationSample (float* lastInputSamples, float newValue) noexcept
29  {
30  lastInputSamples[4] = lastInputSamples[3];
31  lastInputSamples[3] = lastInputSamples[2];
32  lastInputSamples[2] = lastInputSamples[1];
33  lastInputSamples[1] = lastInputSamples[0];
34  lastInputSamples[0] = newValue;
35  }
36 
37  static forcedinline void pushInterpolationSamples (float* lastInputSamples, const float* input, int numOut) noexcept
38  {
39  if (numOut >= 5)
40  {
41  for (int i = 0; i < 5; ++i)
42  lastInputSamples[i] = input[--numOut];
43  }
44  else
45  {
46  for (int i = 0; i < numOut; ++i)
47  pushInterpolationSample (lastInputSamples, input[i]);
48  }
49  }
50 
51  static forcedinline void pushInterpolationSamples (float* lastInputSamples, const float* input,
52  int numOut, int available, int wrapAround) noexcept
53  {
54  if (numOut >= 5)
55  {
56  if (available >= 5)
57  {
58  for (int i = 0; i < 5; ++i)
59  lastInputSamples[i] = input[--numOut];
60  }
61  else
62  {
63  for (int i = 0; i < available; ++i)
64  lastInputSamples[i] = input[--numOut];
65 
66  if (wrapAround > 0)
67  {
68  numOut -= wrapAround;
69 
70  for (int i = available; i < 5; ++i)
71  lastInputSamples[i] = input[--numOut];
72  }
73  else
74  {
75  for (int i = available; i < 5; ++i)
76  lastInputSamples[i] = 0.0f;
77  }
78  }
79  }
80  else
81  {
82  if (numOut > available)
83  {
84  for (int i = 0; i < available; ++i)
85  pushInterpolationSample (lastInputSamples, input[i]);
86 
87  if (wrapAround > 0)
88  {
89  for (int i = 0; i < numOut - available; ++i)
90  pushInterpolationSample (lastInputSamples, input[i + available - wrapAround]);
91  }
92  else
93  {
94  for (int i = 0; i < numOut - available; ++i)
95  pushInterpolationSample (lastInputSamples, 0);
96  }
97  }
98  else
99  {
100  for (int i = 0; i < numOut; ++i)
101  pushInterpolationSample (lastInputSamples, input[i]);
102  }
103  }
104  }
105 
106  template <typename InterpolatorType>
107  static int interpolate (float* lastInputSamples, double& subSamplePos, double actualRatio,
108  const float* in, float* out, int numOut) noexcept
109  {
110  auto pos = subSamplePos;
111 
112  if (actualRatio == 1.0 && pos == 1.0)
113  {
114  memcpy (out, in, (size_t) numOut * sizeof (float));
115  pushInterpolationSamples (lastInputSamples, in, numOut);
116  return numOut;
117  }
118 
119  int numUsed = 0;
120 
121  while (numOut > 0)
122  {
123  while (pos >= 1.0)
124  {
125  pushInterpolationSample (lastInputSamples, in[numUsed++]);
126  pos -= 1.0;
127  }
128 
129  *out++ = InterpolatorType::valueAtOffset (lastInputSamples, (float) pos);
130  pos += actualRatio;
131  --numOut;
132  }
133 
134  subSamplePos = pos;
135  return numUsed;
136  }
137 
138  template <typename InterpolatorType>
139  static int interpolate (float* lastInputSamples, double& subSamplePos, double actualRatio,
140  const float* in, float* out, int numOut, int available, int wrap) noexcept
141  {
142  if (actualRatio == 1.0)
143  {
144  if (available >= numOut)
145  {
146  memcpy (out, in, (size_t) numOut * sizeof (float));
147  pushInterpolationSamples (lastInputSamples, in, numOut, available, wrap);
148  }
149  else
150  {
151  memcpy (out, in, (size_t) available * sizeof (float));
152  pushInterpolationSamples (lastInputSamples, in, numOut, available, wrap);
153 
154  if (wrap > 0)
155  {
156  memcpy (out + available, in + available - wrap, (size_t) (numOut - available) * sizeof (float));
157  pushInterpolationSamples (lastInputSamples, in, numOut, available, wrap);
158  }
159  else
160  {
161  for (int i = 0; i < numOut - available; ++i)
162  pushInterpolationSample (lastInputSamples, 0);
163  }
164  }
165 
166  return numOut;
167  }
168 
169  auto originalIn = in;
170  auto pos = subSamplePos;
171  bool exceeded = false;
172 
173  if (actualRatio < 1.0)
174  {
175  for (int i = numOut; --i >= 0;)
176  {
177  if (pos >= 1.0)
178  {
179  if (exceeded)
180  {
181  pushInterpolationSample (lastInputSamples, 0);
182  }
183  else
184  {
185  pushInterpolationSample (lastInputSamples, *in++);
186 
187  if (--available <= 0)
188  {
189  if (wrap > 0)
190  {
191  in -= wrap;
192  available += wrap;
193  }
194  else
195  {
196  exceeded = true;
197  }
198  }
199  }
200 
201  pos -= 1.0;
202  }
203 
204  *out++ = InterpolatorType::valueAtOffset (lastInputSamples, (float) pos);
205  pos += actualRatio;
206  }
207  }
208  else
209  {
210  for (int i = numOut; --i >= 0;)
211  {
212  while (pos < actualRatio)
213  {
214  if (exceeded)
215  {
216  pushInterpolationSample (lastInputSamples, 0);
217  }
218  else
219  {
220  pushInterpolationSample (lastInputSamples, *in++);
221 
222  if (--available <= 0)
223  {
224  if (wrap > 0)
225  {
226  in -= wrap;
227  available += wrap;
228  }
229  else
230  {
231  exceeded = true;
232  }
233  }
234  }
235 
236  pos += 1.0;
237  }
238 
239  pos -= actualRatio;
240  *out++ = InterpolatorType::valueAtOffset (lastInputSamples, jmax (0.0f, 1.0f - (float) pos));
241  }
242  }
243 
244  subSamplePos = pos;
245 
246  if (wrap == 0)
247  return (int) (in - originalIn);
248 
249  return ((int) (in - originalIn) + wrap) % wrap;
250  }
251 
252  template <typename InterpolatorType>
253  static int interpolateAdding (float* lastInputSamples, double& subSamplePos, double actualRatio,
254  const float* in, float* out, int numOut,
255  int available, int wrap, float gain) noexcept
256  {
257  if (actualRatio == 1.0)
258  {
259  if (available >= numOut)
260  {
261  FloatVectorOperations::addWithMultiply (out, in, gain, numOut);
262  pushInterpolationSamples (lastInputSamples, in, numOut, available, wrap);
263  }
264  else
265  {
266  FloatVectorOperations::addWithMultiply (out, in, gain, available);
267  pushInterpolationSamples (lastInputSamples, in, available, available, wrap);
268 
269  if (wrap > 0)
270  {
271  FloatVectorOperations::addWithMultiply (out, in - wrap, gain, numOut - available);
272  pushInterpolationSamples (lastInputSamples, in - wrap, numOut - available, available, wrap);
273  }
274  else
275  {
276  for (int i = 0; i < numOut-available; ++i)
277  pushInterpolationSample (lastInputSamples, 0.0);
278  }
279  }
280 
281  return numOut;
282  }
283 
284  auto originalIn = in;
285  auto pos = subSamplePos;
286  bool exceeded = false;
287 
288  if (actualRatio < 1.0)
289  {
290  for (int i = numOut; --i >= 0;)
291  {
292  if (pos >= 1.0)
293  {
294  if (exceeded)
295  {
296  pushInterpolationSample (lastInputSamples, 0.0);
297  }
298  else
299  {
300  pushInterpolationSample (lastInputSamples, *in++);
301 
302  if (--available <= 0)
303  {
304  if (wrap > 0)
305  {
306  in -= wrap;
307  available += wrap;
308  }
309  else
310  {
311  exceeded = true;
312  }
313  }
314  }
315 
316  pos -= 1.0;
317  }
318 
319  *out++ += gain * InterpolatorType::valueAtOffset (lastInputSamples, (float) pos);
320  pos += actualRatio;
321  }
322  }
323  else
324  {
325  for (int i = numOut; --i >= 0;)
326  {
327  while (pos < actualRatio)
328  {
329  if (exceeded)
330  {
331  pushInterpolationSample (lastInputSamples, 0.0);
332  }
333  else
334  {
335  pushInterpolationSample (lastInputSamples, *in++);
336 
337  if (--available <= 0)
338  {
339  if (wrap > 0)
340  {
341  in -= wrap;
342  available += wrap;
343  }
344  else
345  {
346  exceeded = true;
347  }
348  }
349  }
350 
351  pos += 1.0;
352  }
353 
354  pos -= actualRatio;
355  *out++ += gain * InterpolatorType::valueAtOffset (lastInputSamples, jmax (0.0f, 1.0f - (float) pos));
356  }
357  }
358 
359  subSamplePos = pos;
360 
361  if (wrap == 0)
362  return (int) (in - originalIn);
363 
364  return ((int) (in - originalIn) + wrap) % wrap;
365  }
366 
367  template <typename InterpolatorType>
368  static int interpolateAdding (float* lastInputSamples, double& subSamplePos, double actualRatio,
369  const float* in, float* out, int numOut, float gain) noexcept
370  {
371  auto pos = subSamplePos;
372 
373  if (actualRatio == 1.0 && pos == 1.0)
374  {
375  FloatVectorOperations::addWithMultiply (out, in, gain, numOut);
376  pushInterpolationSamples (lastInputSamples, in, numOut);
377  return numOut;
378  }
379 
380  int numUsed = 0;
381 
382  while (numOut > 0)
383  {
384  while (pos >= 1.0)
385  {
386  pushInterpolationSample (lastInputSamples, in[numUsed++]);
387  pos -= 1.0;
388  }
389 
390  *out++ += gain * InterpolatorType::valueAtOffset (lastInputSamples, (float) pos);
391  pos += actualRatio;
392  --numOut;
393  }
394 
395  subSamplePos = pos;
396  return numUsed;
397  }
398 }
399 
400 //==============================================================================
401 template <int k>
402 struct LagrangeResampleHelper
403 {
404  static forcedinline void calc (float& a, float b) noexcept { a *= b * (1.0f / k); }
405 };
406 
407 template<>
408 struct LagrangeResampleHelper<0>
409 {
410  static forcedinline void calc (float&, float) noexcept {}
411 };
412 
413 struct LagrangeAlgorithm
414 {
415  static forcedinline float valueAtOffset (const float* inputs, float offset) noexcept
416  {
417  return calcCoefficient<0> (inputs[4], offset)
418  + calcCoefficient<1> (inputs[3], offset)
419  + calcCoefficient<2> (inputs[2], offset)
420  + calcCoefficient<3> (inputs[1], offset)
421  + calcCoefficient<4> (inputs[0], offset);
422  }
423 
424  template <int k>
425  static forcedinline float calcCoefficient (float input, float offset) noexcept
426  {
427  LagrangeResampleHelper<0 - k>::calc (input, -2.0f - offset);
428  LagrangeResampleHelper<1 - k>::calc (input, -1.0f - offset);
429  LagrangeResampleHelper<2 - k>::calc (input, 0.0f - offset);
430  LagrangeResampleHelper<3 - k>::calc (input, 1.0f - offset);
431  LagrangeResampleHelper<4 - k>::calc (input, 2.0f - offset);
432  return input;
433  }
434 };
435 
436 LagrangeInterpolator::LagrangeInterpolator() noexcept { reset(); }
437 LagrangeInterpolator::~LagrangeInterpolator() noexcept {}
438 
440 {
441  subSamplePos = 1.0;
442 
443  for (auto& s : lastInputSamples)
444  s = 0;
445 }
446 
447 int LagrangeInterpolator::process (double actualRatio, const float* in, float* out, int numOut, int available, int wrap) noexcept
448 {
449  return interpolate<LagrangeAlgorithm> (lastInputSamples, subSamplePos, actualRatio, in, out, numOut, available, wrap);
450 }
451 
452 int LagrangeInterpolator::process (double actualRatio, const float* in, float* out, int numOut) noexcept
453 {
454  return interpolate<LagrangeAlgorithm> (lastInputSamples, subSamplePos, actualRatio, in, out, numOut);
455 }
456 
457 int LagrangeInterpolator::processAdding (double actualRatio, const float* in, float* out, int numOut, int available, int wrap, float gain) noexcept
458 {
459  return interpolateAdding<LagrangeAlgorithm> (lastInputSamples, subSamplePos, actualRatio, in, out, numOut, available, wrap, gain);
460 }
461 
462 int LagrangeInterpolator::processAdding (double actualRatio, const float* in, float* out, int numOut, float gain) noexcept
463 {
464  return interpolateAdding<LagrangeAlgorithm> (lastInputSamples, subSamplePos, actualRatio, in, out, numOut, gain);
465 }
466 
467 } // namespace juce
static void JUCE_CALLTYPE addWithMultiply(float *dest, const float *src, float multiplier, int numValues) noexcept
int processAdding(double speedRatio, const float *inputSamples, float *outputSamples, int numOutputSamplesToProduce, float gain) noexcept
int process(double speedRatio, const float *inputSamples, float *outputSamples, int numOutputSamplesToProduce) noexcept