2 * Copyright (C) 2009-2011 David Robillard <d@drobilla.net>
3 * Copyright (C) 2009-2012 Carl Hetherington <carl@carlh.net>
4 * Copyright (C) 2009-2017 Paul Davis <paul@linuxaudiosystems.com>
5 * Copyright (C) 2009 Hans Baier <hansfbaier@googlemail.com>
6 * Copyright (C) 2013-2017 Robin Gareus <robin@gareus.org>
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License along
19 * with this program; if not, write to the Free Software Foundation, Inc.,
20 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
28 #include "ardour/interpolation.h"
29 #include "ardour/midi_buffer.h"
31 using namespace ARDOUR;
35 CubicInterpolation::CubicInterpolation ()
41 CubicInterpolation::interpolate (int channel, samplecnt_t input_samples, Sample *input, samplecnt_t & output_samples, Sample *output)
43 assert (input_samples > 0);
44 assert (output_samples > 0);
47 assert (phase.size () > channel);
49 _speed = fabs (_speed);
53 /* z[0] not set. Two possibilities
55 * 1) we have just been constructed or ::reset()
57 * 2) we were only given 1 sample after construction or
58 * ::reset, and stored it in z[1]
63 /* first call after construction or after ::reset */
65 switch (input_samples) {
67 /* store one sample for use next time. We don't
68 * have enough points to interpolate or even
69 * compute the first z[0] value, but keep z[1]
72 z[1] = input[0]; validate (1);
76 /* store two samples for use next time, and
77 * compute a value for z[0] that will maintain
78 * the slope of the first actual segment. We
79 * still don't have enough samples to interpolate.
81 z[0] = input[0] - (input[1] - input[0]); validate (0);
82 z[1] = input[0]; validate (1);
83 z[2] = input[1]; validate (2);
87 /* We have enough samples to interpolate this time,
88 * but don't have a valid z[0] value because this is the
89 * first call after construction or ::reset.
91 * First point is based on a requirement to maintain
92 * the slope of the first actual segment
94 z[0] = input[0] - (input[1] - input[0]); validate (0);
99 /* at least one call since construction or
100 * after::reset, since we have z[1] set
102 * we can now compute z[0] as required
105 z[0] = z[1] - (input[0] - z[1]); validate (0);
107 /* we'll check the number of samples we've been given
108 in the next switch() statement below, and either
109 just save some more samples or actual interpolate
113 assert (is_valid (0));
116 switch (input_samples) {
118 /* one more sample of input. find the right vX to store
119 it in, and decide if we're ready to interpolate
122 z[1] = input[0]; validate (1);
123 /* still not ready to interpolate */
126 } else if (invalid (2)) {
127 /* still not ready to interpolate */
128 z[2] = input[0]; validate (2);
131 } else if (invalid (3)) {
132 z[3] = input[0]; validate (3);
133 /* ready to interpolate */
137 /* two more samples of input. find the right vX to store
138 them in, and decide if we're ready to interpolate
141 z[1] = input[0]; validate (1);
142 z[2] = input[1]; validate (2);
143 /* still not ready to interpolate */
146 } else if (invalid (2)) {
147 z[2] = input[0]; validate (2);
148 z[3] = input[1]; validate (3);
149 /* ready to interpolate */
150 } else if (invalid (3)) {
151 z[3] = input[0]; validate (3);
152 /* ready to interpolate */
157 /* caller has given us at least enough samples to interpolate a
160 z[1] = input[0]; validate (1);
161 z[2] = input[1]; validate (2);
162 z[3] = input[2]; validate (3);
165 /* ready to interpolate using z[0], z[1], z[2] and z[3] */
167 assert (is_valid (0));
168 assert (is_valid (1));
169 assert (is_valid (2));
170 assert (is_valid (3));
172 /* we can use up to (input_samples - 2) of the input, so compute the
173 * maximum number of output samples that represents.
175 * Remember that the expected common case here is to be given
176 * input_samples that is substantially larger than output_samples,
177 * thus allowing us to always compute output_samples in one call.
180 const samplecnt_t output_from_input = floor ((input_samples - 2) / _speed);
182 /* limit output to either the caller's requested number or the number
183 * determined by the input size.
186 const samplecnt_t limit = std::min (output_samples, output_from_input);
188 samplecnt_t outsample = 0;
189 double distance = phase[channel];
190 samplecnt_t used = floor (distance);
193 while (outsample < limit) {
195 i = floor (distance);
197 /* this call may stop the loop from being vectorized */
198 float fractional_phase_part = fmod (distance, 1.0);
200 /* Cubically interpolate into the output buffer */
201 output[outsample++] = z[1] + 0.5f * fractional_phase_part *
202 (z[2] - z[0] + fractional_phase_part * (4.0f * z[2] + 2.0f * z[0] - 5.0f * z[1] - z[3] +
203 fractional_phase_part * (3.0f * (z[1] - z[2]) - z[0] + z[3])));
213 output_samples = outsample;
214 phase[channel] = fmod (distance, 1.0);
219 CubicInterpolation::reset ()
221 Interpolation::reset ();
226 CubicInterpolation::distance (samplecnt_t nsamples)
228 assert (phase.size () > 0);
229 return floor (floor (phase[0]) + (_speed * nsamples));