1 /*
2 oscils.c:
3
4 Copyright (C) 2002 Istvan Varga
5
6 This file is part of Csound.
7
8 The Csound Library is free software; you can redistribute it
9 and/or modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2.1 of the License, or (at your option) any later version.
12
13 Csound 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 Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with Csound; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
21 02110-1301 USA
22 */
23
24 /* ------ oscils, lphasor, and tablexkt by Istvan Varga (Jan 5 2002) ------ */
25
26 #include "csoundCore.h"
27 #include <math.h>
28 #define CSOUND_OSCILS_C 1
29 #include "oscils.h"
30
31 /* ------------- set up fast sine generator ------------- */
32 /* Input args: */
33 /* a: amplitude */
34 /* f: frequency (-PI - PI) */
35 /* p: initial phase (0 - PI/2) */
36 /* Output args: */
37 /* *x: first output sample */
38 /* *c, *v: coefficients for calculating next sample as */
39 /* shown below: */
40 /* v = v + c * x */
41 /* x = x + v */
42 /* These values are calculated by: */
43 /* x = y[0] */
44 /* c = 2.0 * cos(f) - 2.0 */
45 /* v = y[1] - (c + 1.0) * y[0] */
46 /* where y[0], and y[1] are the first, and */
47 /* second sample of the sine wave to be */
48 /* generated, respectively. */
49 /* -------- written by Istvan Varga, Jan 28 2002 -------- */
50
init_sine_gen(double a,double f,double p,double * x,double * c,double * v)51 static void init_sine_gen(double a, double f, double p,
52 double *x, double *c, double *v)
53 {
54 double y0, y1; /* these should be doubles */
55
56 y0 = sin(p);
57 y1 = sin(p + f);
58 *x = y0;
59 *c = 2.0 * cos(f) - 2.0;
60 *v = y1 - *c * y0 - y0;
61 /* amp. scale */
62 *x *= a; *v *= a;
63 }
64
65 /* -------- oscils set-up -------- */
66
oscils_set(CSOUND * csound,OSCILS * p)67 int32_t oscils_set(CSOUND *csound, OSCILS *p)
68 {
69 int32_t iflg;
70
71 iflg = (int32_t) (*(p->iflg) + FL(0.5)) & 0x07; /* check flags */
72 if (UNLIKELY(iflg & 1)) return OK; /* skip init, nothing to do */
73 p->use_double = (iflg & 2 ? 1 : 0); /* use doubles internally */
74 init_sine_gen((double)*(p->iamp), (double)(*(p->icps) * csound->tpidsr),
75 (double)(*(p->iphs) * TWOPI_F),
76 &(p->xd), &(p->cd), &(p->vd));
77 if (!(p->use_double)) {
78 p->x = (MYFLT) p->xd; /* use floats */
79 p->c = (MYFLT) p->cd;
80 p->v = (MYFLT) p->vd;
81 }
82 return OK;
83 }
84
85 /* -------- oscils performance -------- */
86
oscils(CSOUND * csound,OSCILS * p)87 int32_t oscils(CSOUND *csound, OSCILS *p)
88 {
89 IGN(csound);
90 uint32_t offset = p->h.insdshead->ksmps_offset;
91 uint32_t early = p->h.insdshead->ksmps_no_end;
92 uint32_t n, nsmps = CS_KSMPS;
93 MYFLT *ar, x, c, v;
94 double xd, cd, vd;
95
96 /* copy object data to local variables */
97 ar = p->ar;
98
99 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
100 if (UNLIKELY(early)) {
101 nsmps -= early;
102 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
103 }
104 if (p->use_double) { /* use doubles */
105 xd = p->xd; cd = p->cd; vd = p->vd;
106 for (n=offset; n<nsmps; n++) {
107 ar[n] = (MYFLT) xd;
108 vd += cd * xd;
109 xd += vd;
110 }
111 p->xd = xd; p->vd = vd;
112 }
113 else { /* use floats */
114 x = p->x; c = p->c; v = p->v;
115 for (n=offset; n<nsmps; n++) {
116 ar[n] = x;
117 v += c * x;
118 x += v;
119 }
120 p->x = x; p->v = v;
121 }
122 return OK;
123 }
124
125 /* -------- lphasor set-up -------- */
126
lphasor_set(CSOUND * csound,LPHASOR * p)127 int32_t lphasor_set(CSOUND *csound, LPHASOR *p)
128 {
129 IGN(csound);
130 if (UNLIKELY(*(p->istor) != FL(0.0))) return OK; /* nothing to do */
131
132 p->phs = (double)*(p->istrt); /* start phase */
133 p->lps = (double)*(p->ilps); /* loop start */
134 p->lpe = (double)*(p->ilpe); /* loop end */
135 p->loop_mode = (int32_t) (*(p->imode) + FL(0.5)) & 0x03; /* loop mode */
136 if (p->lpe <= p->lps) p->loop_mode = 0; /* disable loop */
137 p->dir = 1; /* direction */
138 return OK;
139 }
140
141 /* -------- lphasor performance -------- */
142
lphasor(CSOUND * csound,LPHASOR * p)143 int32_t lphasor(CSOUND *csound, LPHASOR *p)
144 {
145 IGN(csound);
146 uint32_t offset = p->h.insdshead->ksmps_offset;
147 uint32_t early = p->h.insdshead->ksmps_no_end;
148 uint32_t n, nsmps = CS_KSMPS;
149 int32_t loop_mode, dir;
150 MYFLT *ar, *xtrns;
151 double trns, phs, lps, lpe, lpt;
152 int32_t assxtr = IS_ASIG_ARG(p->xtrns);
153
154 /* copy object data to local variables */
155 ar = p->ar; xtrns = p->xtrns;
156 phs = p->phs; lps = p->lps; lpe = p->lpe;
157 lpt = lpe - lps;
158 loop_mode = p->loop_mode;
159 trns = (double)*xtrns;
160
161 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
162 if (UNLIKELY(early)) {
163 nsmps -= early;
164 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
165 }
166 for (n=offset; n<nsmps; n++) {
167 if (assxtr) trns = (double)xtrns[n];
168 ar[n] = (MYFLT) phs;
169 phs += (p->dir ? trns : -trns);
170 if (loop_mode) {
171 dir = (trns < 0.0 ? !(p->dir) : p->dir);
172 if (dir && (phs >= lpe)) {
173 phs += lpt * (double)((int32_t)((lps - phs) / lpt));
174 if (loop_mode & 2) {
175 phs = lps + lpe - phs; /* reverse direction */
176 p->dir = !(p->dir);
177 }
178 }
179 else if (!dir && (phs <= lps)) {
180 phs += lpt * (double)((int32_t)((lpe - phs) / lpt));
181 if (loop_mode & 1) {
182 phs = lps + lpe - phs; /* reverse direction */
183 p->dir = !(p->dir);
184 }
185 }
186 }
187 }
188 /* store phase */
189 p->phs = phs;
190 return OK;
191 }
192
193 /* -------- tablexkt set-up -------- */
194
tablexkt_set(CSOUND * csound,TABLEXKT * p)195 int32_t tablexkt_set(CSOUND *csound, TABLEXKT *p)
196 {
197 IGN(csound);
198 p->wsize = (int32_t)(*(p->iwsize) + 0.5); /* window size */
199 if (UNLIKELY(p->wsize < 3)) {
200 p->wsize = 2;
201 }
202 else {
203 p->wsize = ((p->wsize + 2) >> 2) << 2; /* round to nearest */
204 if (p->wsize > 1024) p->wsize = 1024; /* integer multiply of 4 */
205 }
206 /* constant for window calculation */
207 p->win_fact = (FL(1.0) - POWER(p->wsize * FL(0.85172), -FL(0.89624)))
208 / ((MYFLT)((p->wsize * p->wsize) >> 2));
209
210 p->ndx_scl = (*(p->ixmode) == FL(0.0) ? 0 : 1); /* index mode */
211 p->wrap_ndx = (*(p->iwrap) == FL(0.0) ? 0 : 1); /* wrap index */
212 /* use raw index values without scale / offset */
213 if ((*(p->ixoff) != FL(0.0)) || p->ndx_scl) p->raw_ndx = 0;
214 else p->raw_ndx = 1;
215 return OK;
216 }
217
218 /* -------- tablexkt opcode -------- */
219
tablexkt(CSOUND * csound,TABLEXKT * p)220 int32_t tablexkt(CSOUND *csound, TABLEXKT *p)
221 {
222 uint32_t offset = p->h.insdshead->ksmps_offset;
223 uint32_t early = p->h.insdshead->ksmps_no_end;
224 uint32_t n, nsmps = CS_KSMPS;
225 int32_t i, wsize, wsized2, wrap_ndx, warp;
226 double ndx, d, x, c, v, flen_d, onedpi_d, pidwarp_d;
227 int32_t ndx_i=0, flen;
228 MYFLT *ar, *xndx, ndx_f, a0, a1, a2, a3, v0, v1, v2, v3, *ftable;
229 MYFLT onedwarp, win_fact;
230 FUNC *ftp;
231 int32_t asgx = IS_ASIG_ARG(p->xndx);
232
233 /* window size */
234 wsize = p->wsize;
235 if (UNLIKELY((wsize < 2) || (wsize > 1024))) {
236 return csound->PerfError(csound, &(p->h),
237 Str("tablexkt: not initialised"));
238 }
239 wsized2 = wsize >> 1;
240
241 /* check ftable */
242 if (UNLIKELY((ftp = csound->FTnp2Finde(csound, p->kfn)) == NULL))
243 return NOTOK; /* invalid table */
244 if (UNLIKELY((ftable = ftp->ftable) == NULL)) return NOTOK;
245 flen = ftp->flen; /* table length */
246 flen_d = (double)flen;
247
248 /* copy object data to local variables */
249 ar = p->ar;
250 xndx = p->xndx;
251 wrap_ndx = p->wrap_ndx;
252 if ((wsize > 4) && UNLIKELY((*(p->kwarp) > FL(1.001)))) {
253 warp = 1; /* enable warp */
254 onedwarp = FL(1.0) / *(p->kwarp);
255 pidwarp_d = PI / (double)*(p->kwarp);
256 /* correct window for kwarp */
257 x = v = (double)wsized2; x *= x; x = 1.0 / x;
258 v *= (double)onedwarp; v -= (double)((int32_t)v) + 0.5; v *= 4.0 * v;
259 win_fact = (MYFLT)(((double)p->win_fact - x) * v + x);
260 }
261 else {
262 warp = 0; onedwarp = FL(0.0); pidwarp_d = 0.0;
263 win_fact = p->win_fact;
264 }
265 onedpi_d = 1.0 / PI;
266
267 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
268 if (UNLIKELY(early)) {
269 nsmps -= early;
270 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
271 }
272 for (n=offset; n<nsmps; n++) {
273 ndx = (double)*xndx;
274 if (asgx) xndx++;
275 /* calculate table index */
276 if (!(p->raw_ndx)) {
277 ndx += (double)*(p->ixoff);
278 if (p->ndx_scl) ndx *= flen_d;
279 }
280 /* integer and fractional part of table index */
281 ndx_i = (int32_t)ndx;
282 ndx_f = (MYFLT) (ndx - (double)ndx_i);
283 if (ndx_f < FL(0.0)) {
284 ndx_f++; ndx_i--;
285 }
286 /* wrap or limit to allowed range */
287 if (wrap_ndx) {
288 while (ndx_i >= flen) ndx_i -= flen;
289 while (ndx_i < 0) ndx_i += flen;
290 }
291 else { /* limit */
292 if (UNLIKELY(ndx_i < 0)) {
293 ndx_i = 0; ndx_f = FL(0.0);
294 }
295 else if (UNLIKELY(ndx_i >= flen)) {
296 ndx_i = flen; ndx_f = FL(0.0);
297 }
298 }
299 switch (wsize) {
300 case 2: /* ---- linear interpolation ---- */
301 ar[n] = ftable[ndx_i];
302 if (++ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
303 ar[n] += (ftable[ndx_i] - ar[n]) * ndx_f;
304 break;
305 case 4: /* ---- cubic interpolation ---- */
306 /* sample 0 */
307 v1 = ftable[ndx_i];
308 /* sample -1 */
309 if (ndx_i) {
310 v0 = ftable[ndx_i - 1L];
311 } else {
312 v0 = ftable[(wrap_ndx ? flen - 1L : 0L)];
313 }
314 /* sample +1 */
315 if (++ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
316 v2 = ftable[ndx_i];
317 /* sample +2 */
318 if (++ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
319 v3 = ftable[ndx_i];
320 a3 = ndx_f * ndx_f; a3--; a3 *= FL(0.1666666667);
321 a2 = ndx_f; a2++; a0 = (a2 *= FL(0.5)); a0--;
322 a1 = FL(3.0) * a3; a2 -= a1; a0 -= a3; a1 -= ndx_f;
323 ar[n] = (a0 * v0 + a1 * v1 + a2 * v2 + a3 * v3) * ndx_f + v1;
324 break;
325 default: /* ---- sinc interpolation ---- */
326 ar[n] = FL(0.0); /* clear output */
327 ndx = (double)ndx_f;
328 ndx_i += (int32_t)(1 - wsized2);
329 d = (double)(1 - wsized2) - ndx;
330 if (warp) { /* ---- warp enabled ---- */
331 init_sine_gen(onedpi_d, pidwarp_d, pidwarp_d * d, &x, &c, &v);
332 /* samples -(window size / 2 - 1) to -1 */
333 i = wsized2 - 1;
334 do {
335 a1 = (MYFLT) d; a1 = FL(1.0) - a1 * a1 * win_fact;
336 a1 = a1 * a1 / (MYFLT) d;
337 ar[n] += ftable[(ndx_i < 0L ? (wrap_ndx ? ndx_i + flen : 0L)
338 : ndx_i)] * (MYFLT) x * a1;
339 ndx_i++;
340 d++; v += c * x; x += v;
341 } while (--i);
342 /* sample 0 */
343 /* avoid division by zero */
344 if (UNLIKELY(ndx < 0.00003)) ar[n] += onedwarp * ftable[ndx_i];
345 else {
346 a1 = (MYFLT) d; a1 = FL(1.0) - a1 * a1 * win_fact;
347 a1 = a1 * a1 / (MYFLT) d;
348 ar[n] += (MYFLT) x * a1 * ftable[ndx_i];
349 }
350 d++; v += c * x; x += v;
351 if (++ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
352 /* sample 1 */
353 /* avoid division by zero */
354 if (ndx > 0.99997) ar[n] += onedwarp * ftable[ndx_i];
355 else {
356 a1 = (MYFLT) d; a1 = FL(1.0) - a1 * a1 * win_fact;
357 a1 = a1 * a1 / (MYFLT) d;
358 ar[n] += (MYFLT) x * a1 * ftable[ndx_i];
359 }
360 d++; v += c * x; x += v;
361 if (++ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
362 /* samples 2 to (window size / 2) */
363 i = wsized2 - 1;
364 do {
365 a1 = (MYFLT) d; a1 = FL(1.0) - a1 * a1 * win_fact;
366 a1 = a1 * a1 / (MYFLT) d;
367 ar[n] += (MYFLT) x * a1 * ftable[ndx_i];
368 d++; v += c * x; x += v;
369 if (++ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
370 } while (--i);
371 }
372 else { /* ---- warp disabled ---- */
373 /* avoid division by zero */
374 if (UNLIKELY(ndx < 0.00001)) {
375 ndx_i += (int32_t) (wsized2 - 1); /* no need to check here */
376 ar[n] = ftable[ndx_i];
377 }
378 else if (ndx > 0.99999) {
379 ndx_i += (int32_t) wsized2; /* does need range checking */
380 if (ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
381 ar[n] = ftable[ndx_i];
382 }
383 else {
384 /* samples -(window size / 2 - 1) to 0 */
385 i = wsized2 >> 1;
386 do {
387 a1 = (MYFLT) d; a1 = FL(1.0) - a1 * a1 * win_fact;
388 a1 = a1 * a1 / (MYFLT) d;
389 ar[n] += ftable[(ndx_i < 0L ? (wrap_ndx ? ndx_i + flen : 0L)
390 : ndx_i)] * a1;
391 d+=1.0; ndx_i++;
392 a1 = (MYFLT) d; a1 = FL(1.0) - a1 * a1 * win_fact;
393 a1 = a1 * a1 / (MYFLT) d;
394 ar[n] -= ftable[(ndx_i < 0L ? (wrap_ndx ? ndx_i + flen : 0L)
395 : ndx_i)] * a1;
396 d+=1.0; ndx_i++;
397
398 } while (--i);
399 /* samples 1 to (window size / 2) */
400 i = wsized2 >> 1;
401 do {
402 a1 = (MYFLT) d; a1 = FL(1.0) - a1 * a1 * win_fact;
403 a1 = a1 * a1 / (MYFLT) d;
404 ar[n] += a1 * ftable[ndx_i];
405 d+=1.0;
406 if (++ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
407 a1 = (MYFLT) d; a1 = FL(1.0) - a1 * a1 * win_fact;
408 a1 = a1 * a1 / (MYFLT) d;
409 ar[n] -= a1 * ftable[ndx_i];
410 d+=1.0;
411 if (++ndx_i >= flen) ndx_i = (wrap_ndx ? ndx_i - flen : flen);
412 } while (--i);
413 ar[n] *= SIN(PI_F * ndx) / PI_F;
414 }
415 }
416 break;
417 }
418 }
419 return OK;
420 }
421