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);
44 assert (phase.size () > channel);
46 _speed = fabs (_speed);
50 /* z[0] not set. Two possibilities
52 * 1) we have just been constructed or ::reset()
54 * 2) we were only given 1 sample after construction or
55 * ::reset, and stored it in z[1]
60 /* first call after construction or after ::reset */
62 switch (input_samples) {
64 /* store one sample for use next time. We don't
65 * have enough points to interpolate or even
66 * compute the first z[0] value, but keep z[1]
69 z[1] = input[0]; validate (1);
73 /* store two samples for use next time, and
74 * compute a value for z[0] that will maintain
75 * the slope of the first actual segment. We
76 * still don't have enough samples to interpolate.
78 z[0] = input[0] - (input[1] - input[0]); validate (0);
79 z[1] = input[0]; validate (1);
80 z[2] = input[1]; validate (2);
84 /* We have enough samples to interpolate this time,
85 * but don't have a valid z[0] value because this is the
86 * first call after construction or ::reset.
88 * First point is based on a requirement to maintain
89 * the slope of the first actual segment
91 z[0] = input[0] - (input[1] - input[0]); validate (0);
96 /* at least one call since construction or
97 * after::reset, since we have z[1] set
99 * we can now compute z[0] as required
102 z[0] = z[1] - (input[0] - z[1]); validate (0);
104 /* we'll check the number of samples we've been given
105 in the next switch() statement below, and either
106 just save some more samples or actual interpolate
110 assert (is_valid (0));
113 switch (input_samples) {
115 /* one more sample of input. find the right vX to store
116 it in, and decide if we're ready to interpolate
119 z[1] = input[0]; validate (1);
120 /* still not ready to interpolate */
123 } else if (invalid (2)) {
124 /* still not ready to interpolate */
125 z[2] = input[0]; validate (2);
128 } else if (invalid (3)) {
129 z[3] = input[0]; validate (3);
130 /* ready to interpolate */
134 /* two more samples of input. find the right vX to store
135 them in, and decide if we're ready to interpolate
138 z[1] = input[0]; validate (1);
139 z[2] = input[1]; validate (2);
140 /* still not ready to interpolate */
143 } else if (invalid (2)) {
144 z[2] = input[0]; validate (2);
145 z[3] = input[1]; validate (3);
146 /* ready to interpolate */
147 } else if (invalid (3)) {
148 z[3] = input[0]; validate (3);
149 /* ready to interpolate */
154 /* caller has given us at least enough samples to interpolate a
157 z[1] = input[0]; validate (1);
158 z[2] = input[1]; validate (2);
159 z[3] = input[2]; validate (3);
162 /* ready to interpolate using z[0], z[1], z[2] and z[3] */
164 assert (is_valid (0));
165 assert (is_valid (1));
166 assert (is_valid (2));
167 assert (is_valid (3));
169 /* we can use up to (input_samples - 2) of the input, so compute the
170 * maximum number of output samples that represents.
172 * Remember that the expected common case here is to be given
173 * input_samples that is substantially larger than output_samples,
174 * thus allowing us to always compute output_samples in one call.
177 const samplecnt_t output_from_input = floor ((input_samples - 2) / _speed);
179 /* limit output to either the caller's requested number or the number
180 * determined by the input size.
183 const samplecnt_t limit = std::min (output_samples, output_from_input);
185 samplecnt_t outsample = 0;
186 double distance = phase[channel];
187 samplecnt_t used = floor (distance);
190 while (outsample < limit) {
192 i = floor (distance);
194 /* this call may stop the loop from being vectorized */
195 float fractional_phase_part = fmod (distance, 1.0);
197 /* Cubically interpolate into the output buffer */
198 output[outsample++] = z[1] + 0.5f * fractional_phase_part *
199 (z[2] - z[0] + fractional_phase_part * (4.0f * z[2] + 2.0f * z[0] - 5.0f * z[1] - z[3] +
200 fractional_phase_part * (3.0f * (z[1] - z[2]) - z[0] + z[3])));
210 output_samples = outsample;
211 phase[channel] = fmod (distance, 1.0);
216 CubicInterpolation::reset ()
218 Interpolation::reset ();
223 CubicInterpolation::distance (samplecnt_t nsamples)
225 assert (phase.size () > 0);
226 return floor (floor (phase[0]) + (_speed * nsamples));