1 /***************************************************************************
2 Kwave::.cpp - Interpolation types
3 -------------------
4 begin : Sat Feb 03 2001
5 copyright : (C) 2001 by Thomas Eschenbacher
6 email : Thomas Eschenbacher <thomas.eschenbacher@gmx.de>
7 ***************************************************************************/
8
9 /***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
18 #include "config.h"
19
20 #include <KLocalizedString>
21
22 #include "libkwave/Curve.h"
23 #include "libkwave/Interpolation.h"
24 #include "libkwave/String.h"
25 #include "libkwave/Utils.h"
26
27 //***************************************************************************
28 //***************************************************************************
29
30 #define ADD(i,n,s,d) append(i,n,_(s), _(d))
31
fill()32 void Kwave::Interpolation::InterpolationMap::fill()
33 {
34 ADD(INTPOL_LINEAR, 0, "linear", I18N_NOOP("Linear"));
35 ADD(INTPOL_SPLINE, 1, "spline", I18N_NOOP("Spline"));
36 ADD(INTPOL_NPOLYNOMIAL, 2, "n-polynom", I18N_NOOP("Polynom, nth Degree"));
37 ADD(INTPOL_POLYNOMIAL3, 3, "3-polynom", I18N_NOOP("Polynom, 3rd Degree"));
38 ADD(INTPOL_POLYNOMIAL5, 4, "5-polynom", I18N_NOOP("Polynom, 5th Degree"));
39 ADD(INTPOL_POLYNOMIAL7, 5, "7-polynom", I18N_NOOP("Polynom, 7th Degree"));
40 ADD(INTPOL_SAH, 6, "sample_hold",I18N_NOOP("Sample and Hold"));
41 }
42
43 //***************************************************************************
44 //***************************************************************************
45
46 // static initializer
47 Kwave::Interpolation::InterpolationMap
48 Kwave::Interpolation::m_interpolation_map;
49
50 //***************************************************************************
Interpolation(interpolation_t type)51 Kwave::Interpolation::Interpolation(interpolation_t type)
52 :m_curve(), m_x(), m_y(), m_der(), m_type(type)
53 {
54 }
55
56 //***************************************************************************
~Interpolation()57 Kwave::Interpolation::~Interpolation()
58 {
59 }
60
61 //***************************************************************************
descriptions(bool localized)62 QStringList Kwave::Interpolation::descriptions(bool localized)
63 {
64 QStringList list;
65 unsigned int count = m_interpolation_map.count();
66 unsigned int i;
67 for (i = 0; i < count; i++) {
68 interpolation_t index = m_interpolation_map.findFromData(i);
69 list.append(m_interpolation_map.description(index, localized));
70 }
71 return list;
72 }
73
74 //***************************************************************************
name(interpolation_t type)75 QString Kwave::Interpolation::name(interpolation_t type)
76 {
77 return m_interpolation_map.name(type);
78 }
79
80 //***************************************************************************
count()81 unsigned int Kwave::Interpolation::count()
82 {
83 return (m_curve ? m_curve->count() : 0);
84 }
85
86 //***************************************************************************
singleInterpolation(double input)87 double Kwave::Interpolation::singleInterpolation(double input)
88 {
89 if (!count()) return 0.0; // no data ?
90
91 unsigned int degree = 0;
92 unsigned int count = this->count();
93
94 if (input < 0.0) input = 0.0;
95 if (input > 1.0) input = 1.0;
96
97 switch (m_type) {
98 case INTPOL_LINEAR:
99 {
100 unsigned int i = 1;
101 while ((i < count) && (m_x[i] < input))
102 i++;
103
104 double dif1 = m_x[i] - m_x[i-1]; //!=0 per definition
105 double dif2 = input - m_x[i-1];
106
107 return (m_y[i-1] + ((m_y[i] - m_y[i-1])*dif2 / dif1));
108 }
109 case INTPOL_SPLINE:
110 {
111 double a, b, diff;
112 unsigned int j = 1;
113
114 while ((j < count) && (m_x[j] < input))
115 j++;
116
117 diff = m_x[j] - m_x[j-1];
118
119 a = (m_x[j] - input) / diff; //div should not be 0
120 b = (input - m_x[j-1]) / diff;
121
122 return (a*m_y[j-1] + b*m_y[j] + ((a*a*a - a)*m_der[j - 1] +
123 (b*b*b - b)*m_der[j])*(diff*diff) / 6);
124 }
125 case INTPOL_NPOLYNOMIAL:
126 {
127 double ny = m_y[0];
128 for (unsigned int j = 1; j < count; j++)
129 ny = ny * (input - m_x[j]) + m_y[j];
130 return ny;
131 }
132 case INTPOL_SAH: //Sample and hold
133 {
134 unsigned int i = 1;
135 while ((i < count) && (m_x[i] < input))
136 i++;
137 return m_y[i-1];
138 }
139 case INTPOL_POLYNOMIAL3:
140 degree = 3;
141 break;
142 case INTPOL_POLYNOMIAL5:
143 degree = 5;
144 break;
145 case INTPOL_POLYNOMIAL7:
146 degree = 7;
147 break;
148 }
149
150 if (degree && (degree <= 7)) {
151 Q_ASSERT(m_curve);
152 if (!m_curve) return 0;
153
154 // use polynom
155 double ny;
156 QVector<double> ax(7);
157 QVector<double> ay(7);
158
159 unsigned int i = 1;
160 while ((i < count) && (m_x[i] < input))
161 i++;
162
163 createPolynom(*m_curve, ax, ay, i - 1 - degree/2, degree);
164
165 ny = ay[0];
166 for (unsigned int j = 1; j < degree; j++)
167 ny = ny * (input - ax[j]) + ay[j];
168 return ny;
169 }
170
171 return 0;
172 }
173
174 //***************************************************************************
prepareInterpolation(const Kwave::Curve & points)175 bool Kwave::Interpolation::prepareInterpolation(const Kwave::Curve &points)
176 {
177 m_curve = &points;
178 if (!count()) return false; // no data ?
179
180 m_x = QVector<double>((count() + 1), double(0.0));
181 m_y = QVector<double>((count() + 1), double(0.0));
182 m_der = QVector<double>();
183
184 unsigned int c = 0;
185 Kwave::Curve::ConstIterator it(points);
186 while (it.hasNext()) {
187 Kwave::Curve::Point p = it.next();
188 m_x[c] = p.x();
189 m_y[c] = p.y();
190 c++;
191 }
192 m_x[c] = m_y[c] = 0.0;
193
194 switch (m_type) {
195 case INTPOL_NPOLYNOMIAL:
196 createFullPolynom(points, m_x, m_y);
197 break;
198 case INTPOL_SPLINE:
199 m_der = QVector<double>((count() + 1), double(0.0));
200 get2Derivate(m_x, m_y, m_der, count());
201 break;
202 default:
203 ;
204 }
205 return true;
206 }
207
208 //***************************************************************************
limitedInterpolation(const Kwave::Curve & points,unsigned int len)209 QVector<double> Kwave::Interpolation::limitedInterpolation(
210 const Kwave::Curve &points, unsigned int len)
211 {
212 QVector<double> y = interpolation(points, len);
213 for (unsigned int i = 0; i < len; i++) {
214 if (y[i] > 1) y[i] = 1;
215 if (y[i] < 0) y[i] = 0;
216 }
217 return y;
218 }
219
220 //***************************************************************************
interpolation(const Kwave::Curve & points,unsigned int len)221 QVector<double> Kwave::Interpolation::interpolation(
222 const Kwave::Curve &points, unsigned int len)
223 {
224 Q_ASSERT(len);
225 if (!len) return QVector<double>();
226
227 unsigned int degree = 0;
228 QVector<double> y_out(len, 0.0);
229
230 switch (m_type) {
231 case INTPOL_LINEAR:
232 {
233 Kwave::Curve::ConstIterator it(points);
234
235 if (it.hasNext()) {
236 Kwave::Curve::Point p = it.next();
237 double x0 = p.x();
238 double y0 = p.y();
239
240 while (it.hasNext()) {
241 p = it.next();
242 const double x1 = p.x();
243 const double y1 = p.y();
244
245 double dy = (y1 - y0);
246 int dx = Kwave::toInt((x1 - x0) * len);
247 int min = Kwave::toInt(x0 * len);
248
249 Q_ASSERT(x0 >= 0.0);
250 Q_ASSERT(x1 <= 1.0);
251 for (int i = Kwave::toInt(x0 * len);
252 i < Kwave::toInt(x1 * len); i++) {
253 double h = dx ? ((double(i - min)) / dx) : 0.0;
254 y_out[i] = y0 + (h * dy);
255 }
256 x0 = x1;
257 y0 = y1;
258 }
259 }
260 break;
261 }
262 case INTPOL_SPLINE:
263 {
264 int t = 1;
265 unsigned int count = points.count();
266
267 double ny = 0;
268 QVector<double> der(count + 1);
269 QVector<double> x(count + 1);
270 QVector<double> y(count + 1);
271
272 Kwave::Curve::ConstIterator it(points);
273 while (it.hasNext()) {
274 Kwave::Curve::Point p = it.next();
275 x[t] = p.x();
276 y[t] = p.y();
277 t++;
278 }
279
280 get2Derivate(x, y, der, count);
281
282 int start = Kwave::toInt(x[1] * len);
283
284 for (unsigned int j = 2; j <= count; j++) {
285 const int ent = Kwave::toInt(x[j] * len);
286 for (int i = start; i < ent; i++) {
287 const double xin = static_cast<double>(i) / len;
288 const double h = x[j] - x[j - 1];
289
290 if (!qFuzzyIsNull(h)) {
291 const double a = (x[j] - xin) / h;
292 const double b = (xin - x[j - 1]) / h;
293
294 ny = (a * y[j - 1] + b * y[j] +
295 ((a * a * a - a) * der[j - 1] +
296 (b * b * b - b) * der[j]) * (h * h) / 6.0);
297 }
298
299 y_out[i] = ny;
300 start = ent;
301 }
302 }
303 break;
304 }
305 case INTPOL_POLYNOMIAL3:
306 if (!degree) degree = 3;
307 /* FALLTHROUGH */
308 case INTPOL_POLYNOMIAL5:
309 if (!degree) degree = 5;
310 /* FALLTHROUGH */
311 case INTPOL_POLYNOMIAL7:
312 {
313 if (!degree) degree = 7;
314 const unsigned int count = points.count();
315 QVector<double> x(7);
316 QVector<double> y(7);
317
318 if (count) {
319 for (unsigned int px = 0; px < count - 1; px++) {
320 createPolynom (points, x, y, px - degree / 2, degree);
321 const double start = points[px].x();
322
323 double ent;
324 if (px >= count - degree / 2 + 1)
325 ent = 1;
326 else
327 ent = points[px + 1].x();
328
329 for (int i = Kwave::toInt(start * len);
330 i < Kwave::toInt(ent * len); i++)
331 {
332 double ny = y[0];
333 for (unsigned int j = 1; j < degree; j++)
334 ny = ny * ((static_cast<double>(i)) / len - x[j])
335 + y[j];
336
337 y_out[i] = ny;
338 }
339 }
340 }
341 break;
342 }
343 case INTPOL_NPOLYNOMIAL:
344 {
345 const int count = points.count();
346 if (count != 0) {
347 QVector<double> x(count+1);
348 QVector<double> y(count+1);
349
350 createFullPolynom(points, x, y);
351
352 for (unsigned int i = 1; i < len; i++) {
353 const double px = static_cast<double>(i) / len;
354 double ny = y[0];
355 for (int j = 1; j < count; j++)
356 ny = ny * (px - x[j]) + y[j];
357
358 y_out[i] = ny;
359 }
360 }
361 break;
362 }
363 case INTPOL_SAH:
364 {
365 Kwave::Curve::ConstIterator it(points);
366 if (it.hasNext()) {
367 Kwave::Curve::Point p = it.next();
368 double x0 = p.x();
369 double y0 = p.y();
370
371 while (it.hasNext()) {
372 p = it.next();
373 const double x1 = p.x();
374 const double y1 = p.y();
375
376 for (int i = Kwave::toInt(x0 * len);
377 i < Kwave::toInt(x1 * len); i++)
378 y_out[i] = y0;
379
380 x0 = x1;
381 y0 = y1;
382 }
383 }
384 break;
385 }
386 }
387
388 return y_out;
389 }
390
391 //***************************************************************************
createFullPolynom(const Kwave::Curve & points,QVector<double> & x,QVector<double> & y)392 void Kwave::Interpolation::createFullPolynom(const Kwave::Curve &points,
393 QVector<double> &x, QVector<double> &y)
394 {
395 Q_ASSERT(!points.isEmpty());
396 Q_ASSERT(m_curve);
397 if (points.isEmpty()) return;
398 if (!m_curve) return;
399
400 Q_ASSERT(points.count() == m_curve->count());
401 if (points.count() != m_curve->count()) return;
402
403 unsigned int count = 0;
404 Kwave::Curve::ConstIterator it(points);
405 while (it.hasNext()) {
406 Kwave::Curve::Point p = it.next();
407 x[count] = p.x();
408 y[count] = p.y();
409 count++;
410 }
411
412 for (unsigned int k = 0; k < count; k++)
413 for (unsigned int j = k; j; ) {
414 j--;
415 y[j] = (y[j] - y[j+1]) / (x[j] - x[k]);
416 }
417 }
418
419 //***************************************************************************
get2Derivate(const QVector<double> & x,const QVector<double> & y,QVector<double> & ab,unsigned int n)420 void Kwave::Interpolation::get2Derivate(const QVector<double> &x,
421 const QVector<double> &y, QVector<double> &ab, unsigned int n)
422 {
423 Q_ASSERT(n);
424 if (!n) return;
425
426 unsigned int i, k;
427 QVector<double> u(n);
428
429 ab[0] = ab[1] = 0;
430 u[0] = u[1] = 0;
431
432 for (i = 2; i < n; ++i) {
433 const double sig = (x[i] - x[i-1]) / (x[i+1] - x[i-1]);
434 const double p = sig * ab[i-1] + 2;
435 ab[i] = (sig-1) / p;
436 u[i] = (y[i+1] - y[i]) / (x[i+1] - x[i])
437 - (y[i] - y[i-1]) / (x[i] - x[i-1]);
438 u[i] = (6 * u[i] / (x[i+1] - x[i-1]) - sig * u[i-1]) / p;
439 }
440
441 const double qn = 0;
442 const double un = 0;
443 ab[n] = (un - qn * u[n - 1]) / (qn * ab[n - 1] + 1);
444
445 for (k = n - 1; k > 0; --k)
446 ab[k] = ab[k] * ab[k + 1] + u[k];
447
448 }
449
450 //***************************************************************************
createPolynom(const Kwave::Curve & points,QVector<double> & x,QVector<double> & y,int pos,unsigned int degree)451 void Kwave::Interpolation::createPolynom(const Kwave::Curve &points,
452 QVector<double> &x, QVector<double> &y,
453 int pos, unsigned int degree)
454 {
455 unsigned int count = 0;
456 Kwave::Curve::ConstIterator it(points);
457 if (!it.hasNext()) return;
458 Kwave::Curve::Point p = it.next();
459
460 if (pos < 0) {
461 switch (pos) {
462 case -3:
463 x[count] = -1.5;
464 y[count++] = p.y();
465 pos++;
466 /* FALLTHROUGH */
467 case -2:
468 x[count] = -1.0;
469 y[count++] = p.y();
470 pos++;
471 /* FALLTHROUGH */
472 case -1:
473 x[count] = -0.5;
474 y[count++] = p.y();
475 pos++;
476 break;
477 }
478 }
479
480 for (int i = 0; i < pos; i++)
481 p = it.next();
482
483 while ((count < degree) && it.hasNext()) {
484 p = it.next();
485 x[count] = p.x();
486 y[count++] = p.y();
487 }
488
489 int i = 1;
490 it.toBack();
491 p = it.previous();
492 while (count < degree) {
493 x[count] = 1.0 + 0.5 * (i++);
494 y[count++] = p.y();
495 }
496
497 // create coefficients in y[i] and x[i];
498 for (unsigned int k = 0; k < degree; k++)
499 for (int j = k - 1; j >= 0; j--)
500 y[j] = (y[j] - y[j + 1]) / (x[j] - x[k]);
501
502 }
503
504 //***************************************************************************
505 //***************************************************************************
506