2 Copyright (C) 2012 Paul Davis
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
25 #include "ardour/interpolation.h"
26 #include "ardour/midi_buffer.h"
28 using namespace ARDOUR;
32 CubicInterpolation::CubicInterpolation ()
38 CubicInterpolation::interpolate (int channel, samplecnt_t input_samples, Sample *input, samplecnt_t & output_samples, Sample *output)
40 assert (input_samples > 0);
41 assert (output_samples > 0);
45 _speed = fabs (_speed);
49 /* z[0] not set. Two possibilities
51 * 1) we have just been constructed or ::reset()
53 * 2) we were only given 1 sample after construction or
54 * ::reset, and stored it in z[1]
59 /* first call after construction or after ::reset */
61 switch (input_samples) {
63 /* store one sample for use next time. We don't
64 * have enough points to interpolate or even
65 * compute the first z[0] value, but keep z[1]
68 z[1] = input[0]; validate (1);
72 /* store two samples for use next time, and
73 * compute a value for z[0] that will maintain
74 * the slope of the first actual segment. We
75 * still don't have enough samples to interpolate.
77 z[0] = input[0] - (input[1] - input[0]); validate (0);
78 z[1] = input[0]; validate (1);
79 z[2] = input[1]; validate (2);
83 /* We have enough samples to interpolate this time,
84 * but don't have a valid z[0] value because this is the
85 * first call after construction or ::reset.
87 * First point is based on a requirement to maintain
88 * the slope of the first actual segment
90 z[0] = input[0] - (input[1] - input[0]); validate (0);
95 /* at least one call since construction or
96 * after::reset, since we have z[1] set
98 * we can now compute z[0] as required
101 z[0] = z[1] - (input[0] - z[1]); validate (0);
103 /* we'll check the number of samples we've been given
104 in the next switch() statement below, and either
105 just save some more samples or actual interpolate
109 assert (is_valid (0));
112 switch (input_samples) {
114 /* one more sample of input. find the right vX to store
115 it in, and decide if we're ready to interpolate
118 z[1] = input[0]; validate (1);
119 /* still not ready to interpolate */
122 } else if (invalid (2)) {
123 /* still not ready to interpolate */
124 z[2] = input[0]; validate (2);
127 } else if (invalid (3)) {
128 z[3] = input[0]; validate (3);
129 /* ready to interpolate */
133 /* two more samples of input. find the right vX to store
134 them in, and decide if we're ready to interpolate
137 z[1] = input[0]; validate (1);
138 z[2] = input[1]; validate (2);
139 /* still not ready to interpolate */
142 } else if (invalid (2)) {
143 z[2] = input[0]; validate (2);
144 z[3] = input[1]; validate (3);
145 /* ready to interpolate */
146 } else if (invalid (3)) {
147 z[3] = input[0]; validate (3);
148 /* ready to interpolate */
153 /* caller has given us at least enough samples to interpolate a
156 z[1] = input[0]; validate (1);
157 z[2] = input[1]; validate (2);
158 z[3] = input[2]; validate (3);
161 /* ready to interpolate using z[0], z[1], z[2] and z[3] */
163 assert (is_valid (0));
164 assert (is_valid (1));
165 assert (is_valid (2));
166 assert (is_valid (3));
168 /* we can use up to (input_samples - 2) of the input, so compute the
169 * maximum number of output samples that represents.
171 * Remember that the expected common case here is to be given
172 * input_samples that is substantially larger than output_samples,
173 * thus allowing us to always compute output_samples in one call.
176 const samplecnt_t output_from_input = floor ((input_samples - 2) / _speed);
178 /* limit output to either the caller's requested number or the number
179 * determined by the input size.
182 const samplecnt_t limit = std::min (output_samples, output_from_input);
184 samplecnt_t outsample = 0;
185 double distance = phase[channel];
186 samplecnt_t used = floor (distance);
189 while (outsample < limit) {
191 i = floor (distance);
193 /* this call may stop the loop from being vectorized */
194 float fractional_phase_part = fmod (distance, 1.0);
196 /* Cubically interpolate into the output buffer */
197 output[outsample++] = z[1] + 0.5f * fractional_phase_part *
198 (z[2] - z[0] + fractional_phase_part * (4.0f * z[2] + 2.0f * z[0] - 5.0f * z[1] - z[3] +
199 fractional_phase_part * (3.0f * (z[1] - z[2]) - z[0] + z[3])));
209 output_samples = outsample;
210 phase[channel] = fmod (distance, 1.0);
215 CubicInterpolation::reset ()
217 Interpolation::reset ();
222 CubicInterpolation::distance (samplecnt_t nsamples)
224 return floor (floor (phase[0]) + (_speed * nsamples));