4
Copyright (C) 2012 John ffitch
5
Based in part on code from Vercoe and Whittle
7
This file is part of Csound.
9
The Csound Library is free software; you can redistribute it
10
and/or modify it under the terms of the GNU Lesser General Public
11
License as published by the Free Software Foundation; either
12
version 2.1 of the License, or (at your option) any later version.
14
Csound is distributed in the hope that it will be useful,
15
but WITHOUT ANY WARRANTY; without even the implied warranty of
16
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17
GNU Lesser General Public License for more details.
19
You should have received a copy of the GNU Lesser General Public
20
License along with Csound; if not, write to the Free Software
21
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
25
#include "csoundCore.h" /* UGENS2.C */
30
#define MYFLOOR(x) (x >= FL(0.0) ? (int32)x : (int32)((double)x - 0.99999999))
32
/*****************************************************************************/
34
/* Table read code - see TABLE data structure in ugens2.h. */
36
/*************************************/
38
static int itblchk(CSOUND *csound, TABLE *p)
40
if (UNLIKELY((p->ftp = csound->FTFind(csound, p->xfn)) == NULL))
41
return csound->InitError(csound, Str("table does not exist"));
43
/* Although TABLE has an integer variable for the table number
44
* (p->pfn) we do not need to write it. We know that the k
45
* and a rate functions which will follow will not be
46
* expecting a changed table number.
48
* p->pfn exists only for checking table number changes for
49
* functions which are expecting a k rate table number. */
51
/* Set denormalisation factor to 1 or table length, depending
52
* on the state of ixmode. */
54
p->xbmul = p->ftp->flen;
58
/* Multiply the ixoff value by the xbmul denormalisation
59
* factor and then check it is between 0 and the table length.
61
* Bug fix: was p->offset * *p->ixoff */
63
if (UNLIKELY((p->offset = p->xbmul * *p->ixoff) < FL(0.0) ||
64
p->offset > p->ftp->flen)) {
65
return csound->InitError(csound, Str("Offset %f < 0 or > tablelength"),
69
p->wrap = (int)*p->iwrap;
73
int pktable(CSOUND *,TABLE*);
74
int pktabli(CSOUND *,TABLE*);
75
int pktabl3(CSOUND *,TABLE*);
77
int pitable(CSOUND *csound, TABLE *p)
79
if (LIKELY(itblchk(csound,p)==OK)) return pktable(csound,p);
83
int pitabli(CSOUND *csound, TABLE *p)
85
if (LIKELY(itblchk(csound,p)==OK)) return pktabli(csound,p);
89
int pitabl3(CSOUND *csound, TABLE *p)
91
if (LIKELY(itblchk(csound,p)==OK)) return pktabl3(csound,p);
95
/*---------------------------------------------------------------------------*/
97
/* Functions which read the table. */
99
int pktable(CSOUND *csound, TABLE *p)
106
if (UNLIKELY(ftp==NULL)) goto err1; /* RWD fix */
109
/* Multiply ndx by denormalisation factor, and add in the offset
110
* - already denormalised - by tblchk().
111
* xbmul = 1 or table length depending on state of ixmode. */
113
ndx = ( ndx * p->xbmul) + p->offset;
115
/* ndx now includes the offset and is ready to address the table.
117
* The original code was:
118
* indx = (long) (ndx + p->offset);
120
* This is a problem, causes problems with negative numbers.
123
indx = (int32) MYFLOOR((double)ndx);
125
/* Now for "limit mode" - the "non-wrap" function, depending on
128
* The following section of code limits the final index to 0 and
129
* the last location in the table.
131
* It is only used when wrap is OFF. The wrapping is achieved by
132
* code after this - when this code is not run. */
134
/* Previously this code limited the upper range of the indx to
135
* the table length - for instance 8. Then the result was ANDed
136
* with a mask (for instance 7).
138
* This meant that when the input index was 8 or above, it got
139
* reduced to 0. What we want is for it to stick at the index
140
* which reads the last value from the table - in this example
143
* So instead of limiting to the table length, we limit to
144
* (table length - 1). */
145
if (UNLIKELY(indx > length - 1))
148
/* Now limit negative values to zero. */
149
else if (UNLIKELY(indx < 0L))
152
/* The following code used to use an AND with an integer like 0000 0111
153
* to wrap the current index within the range of the table. In
154
* the original version, this code always ran, but with the new
155
* (length - 1) code above, it would have no effect, so it is now
156
* an else operation - running only when iwrap = 1. This may save
157
* half a usec or so. */
158
/* Now safe against non power-of-2 tables */
159
else if (indx>=length) indx = indx % length;
160
else if (indx<0) indx = length - (-indx)%length;
162
/* Now find the address of the start of the table, add it to the
163
* index, read the value from there and write it to the
165
*p->rslt = *(ftp->ftable + indx);
168
return csound->PerfError(csound, Str("ptable(krate): not initialised"));
173
int ptablefn(CSOUND *csound, TABLE *p)
176
MYFLT *rslt, *pxndx, *tab;
178
int n, nsmps=csound->ksmps;
179
MYFLT ndx, xbmul, offset;
183
if (UNLIKELY(ftp==NULL)) goto err1; /* RWD fix */
187
xbmul = (MYFLT)p->xbmul;
189
//mask = ftp->lenmask;
191
for (n=0; n<nsmps; n++) {
192
/* Read in the next raw index and increment the pointer ready
193
* for the next cycle.
195
* Then multiply the ndx by the denormalising factor and add in
198
ndx = (pxndx[n] * xbmul) + offset;
199
indx = (int32) MYFLOOR((double)ndx);
201
/* Limit = non-wrap. Limits to 0 and (length - 1), or does the
202
* wrap code. See notes above in ktable(). */
204
if (UNLIKELY(indx >= length))
206
else if (UNLIKELY(indx < (int32)0))
209
/* do the wrap code only if we are not doing the non-wrap code. */
210
else if (indx >= length) indx %= length;
211
else if (indx < 0) indx = length - (-indx)%length;
216
return csound->PerfError(csound, Str("table: not initialised"));
221
int pktabli(CSOUND *csound, TABLE *p)
225
MYFLT v1, v2, fract, ndx;
228
if (UNLIKELY(ftp==NULL)) goto err1;
231
/* Multiply ndx by denormalisation factor.
232
* xbmul is 1 or table length depending on state of ixmode.
233
* Add in the offset, which has already been denormalised by
236
ndx = (ndx * p->xbmul) + p->offset;
237
indx = (int32) MYFLOOR((double)ndx);
239
/* We need to generate a fraction - How much above indx is ndx?
240
* It will be between 0 and just below 1.0. */
244
if (UNLIKELY(ndx >= length)) {
248
else if (UNLIKELY(ndx < 0)) {
253
/* We are in wrap mode, so do the wrap function. */
254
else if (indx>=length) indx %= length;
255
else if (indx<0) indx = length - (-indx)%length;
257
/* Now read the value at indx and the one beyond */
258
v1 = *(ftp->ftable + indx);
260
if (indx>=length) indx -= length;
261
v2 = *(ftp->ftable + indx);
262
*p->rslt = v1 + (v2 - v1) * fract;
265
return csound->PerfError(csound, Str("ptablei(krate): not initialised"));
269
int pktabl3(CSOUND *csound, TABLE *p)
273
MYFLT v1, v2, fract, ndx;
276
if (UNLIKELY(ftp==NULL)) goto err1;
279
/* Multiply ndx by denormalisation factor.
280
* xbmul is 1 or table length depending on state of ixmode.
281
* Add in the offset, which has already been denormalised by
284
ndx = (ndx * p->xbmul) + p->offset;
285
indx = (int32) MYFLOOR((double)ndx);
287
/* We need to generate a fraction - How much above indx is ndx?
288
* It will be between 0 and just below 1.0. */
292
if (UNLIKELY(ndx >= length)) {
296
else if (UNLIKELY(ndx < 0)) {
301
/* We are in wrap mode, so do the wrap function. */
302
else if (indx>=length) indx %= length;
303
else if (indx<0) indx = length - (-indx)%length;
305
/* interpolate with cubic if we can, else linear */
306
if (UNLIKELY(indx<1 || indx==length-2 || length <4)) {
307
v1 = *(ftp->ftable + indx);
308
v2 = *(ftp->ftable + indx + 1);
309
*p->rslt = v1 + (v2 - v1) * fract;
312
MYFLT *tab = ftp->ftable;
313
MYFLT ym1 = tab[indx-1], y0 = tab[indx];
314
MYFLT y1 = tab[indx+1], y2 = tab[indx+2];
315
MYFLT frsq = fract*fract;
316
MYFLT frcu = frsq*ym1;
317
MYFLT t1 = y2 + y0+y0+y0;
318
*p->rslt = y0 + FL(0.5)*frcu
319
+ fract*(y1 - frcu/FL(6.0) - t1/FL(6.0) - ym1/FL(3.0))
320
+ frsq*fract*(t1/FL(6.0) - FL(0.5)*y1) + frsq*(FL(0.5)* y1 - y0);
324
return csound->PerfError(csound, Str("ptable3(krate): not initialised"));
327
int ptabli(CSOUND *csound, TABLE *p)
331
int n, nsmps=csound->ksmps;
332
MYFLT *rslt, *pxndx, *tab;
333
MYFLT fract, v1, v2, ndx, xbmul, offset;
336
if (UNLIKELY(ftp==NULL)) goto err1;
340
xbmul = (MYFLT)p->xbmul;
342
//mask = ftp->lenmask;
344
/* As for ktabli() code to handle non wrap mode, and wrap mode. */
346
for (n=0; n<nsmps; n++) {
347
/* Read in the next raw index and increment the pointer ready
348
* for the next cycle.
349
* Then multiply the ndx by the denormalising factor and add in
351
ndx = (pxndx[n] * xbmul) + offset;
353
if (UNLIKELY(ndx <= FL(0.0))) {
357
if (UNLIKELY(indx >= length)) {
358
rslt[n] = tab[length-1];
361
/* We need to generate a fraction - How much above indx is ndx?
362
* It will be between 0 and just below 1.0. */
364
/* As for ktabli(), read two values and interpolate between
368
if (indx>=length) indx = length-1;
370
rslt[n] = v1 + (v2 - v1)*fract;
373
else { /* wrap mode */
374
for (n=0; n<nsmps; n++) {
376
/* Read in the next raw index and increment the pointer ready
377
* for the next cycle.
378
* Then multiply the ndx by the denormalising factor and add in
380
ndx = (pxndx[n] * xbmul) + offset;
381
indx = (int32) MYFLOOR(ndx);
382
/* We need to generate a fraction - How much above indx is ndx?
383
* It will be between 0 and just below 1.0. */
385
if (indx >= length) indx %= length;
386
else if (indx<0) indx = length-(-indx)%length;
387
/* As for ktabli(), read two values and interpolate between
391
if (j >= length) j -= length;
393
rslt[n] = v1 + (v2 - v1)*fract;
398
return csound->PerfError(csound, Str("ptablei: not initialised"));
401
int ptabl3(CSOUND *csound, TABLE *p) /* Like ptabli but cubic interpolation */
405
int n, nsmps=csound->ksmps;
406
MYFLT *rslt, *pxndx, *tab;
407
MYFLT fract, v1, v2, ndx, xbmul, offset;
411
if (UNLIKELY(ftp==NULL)) goto err1;
415
xbmul = (MYFLT)p->xbmul;
418
for (n=0; n<nsmps; n++) {
419
/* Read in the next raw index and increment the pointer ready
420
* for the next cycle.
421
* Then multiply the ndx by the denormalising factor and add in
424
ndx = (pxndx[n] * xbmul) + offset;
425
indx = (int32) MYFLOOR((double)ndx);
427
/* We need to generate a fraction - How much above indx is ndx?
428
* It will be between 0 and just below 1.0. */
430
/* As for pktabli() code to handle non wrap mode, and wrap mode. */
432
if (UNLIKELY(ndx >= length)) {
436
else if (UNLIKELY(ndx < 0)) {
441
else if (UNLIKELY(indx>=length)) indx %= length;
442
else if (UNLIKELY(indx<0)) indx = length-(-indx)%length;
443
/* interpolate with cubic if we can */
444
if (UNLIKELY(indx <1 || indx == length-2 || length<4)) {
445
/* Too short or at ends */
448
rslt[n] = v1 + (v2 - v1)*fract;
451
MYFLT ym1 = tab[indx-1], y0 = tab[indx];
452
int j = (indx+1<length ? indx+1 : indx+1-length);
453
int k = (indx+2<length ? indx+2 : indx+2-length);
454
MYFLT y1 = tab[j], y2 = tab[k];
455
MYFLT frsq = fract*fract;
456
MYFLT frcu = frsq*ym1;
457
MYFLT t1 = y2 + y0+y0+y0;
458
rslt[n] = y0 + FL(0.5)*frcu +
459
fract*(y1 - frcu/FL(6.0) - t1/FL(6.0) - ym1/FL(3.0)) +
460
frsq*fract*(t1/FL(6.0) - FL(0.5)*y1) + frsq*(FL(0.5)* y1 - y0);
465
return csound->PerfError(csound, Str("ptable3: not initialised"));
468
extern int itblchkw(CSOUND *, TABLEW*);
469
int pktablew(CSOUND *, TABLEW*);
470
int pitablew(CSOUND *csound, TABLEW *p)
472
if (LIKELY(itblchkw(csound, p) == OK))
473
return pktablew(csound, p);
477
/*---------------------------------------------------------------------------*/
479
/* pktablew is called with p pointing to the TABLEW data structure -
480
* which contains the input arguments. */
482
int pktablew(CSOUND *csound, TABLEW *p)
484
/* Pointer to data structure for accessing the table we will be
489
MYFLT ndx; /* for calculating index of read. */
490
MYFLT *ptab; /* Where we will write */
492
/*-----------------------------------*/
493
/* Assume that TABLEW has been set up correctly. */
498
/* Multiply ndx by denormalisation factor. and add in the
499
* offset - already denormalised - by tblchkw().
500
* xbmul = 1 or table length depending on state of ixmode. */
502
ndx = (ndx * p->xbmul) + p->offset;
504
/* ndx now includes the offset and is ready to address the table.
505
* However we have three modes to consider:
506
* igwm = 0 Limit mode.
511
/* Limit mode - when igmode = 0.
513
* Limit the final index to 0 and the last location in the table.
515
indx = (int32) MYFLOOR(ndx); /* Limit to (table length - 1) */
516
if (UNLIKELY(indx > length - 1))
517
indx = length - 1; /* Limit the high values. */
518
else if (UNLIKELY(indx < 0L)) indx = 0L; /* limit negative values to zero. */
520
/* Wrap and guard point mode.
521
* In guard point mode only, add 0.5 to the index. */
523
if (p->iwgm == 2) ndx += FL(0.5);
524
indx = (int32) MYFLOOR(ndx);
526
/* Both wrap and guard point mode.
527
* The following code uses an AND with an integer like 0000 0111 to wrap
528
* the current index within the range of the table. */
529
if (UNLIKELY(indx>=length)) indx %= length;
530
else if (UNLIKELY(indx<0)) indx = length-(-indx)%length;
532
/* Calculate the address of where we
533
* want to write to, from indx and the
534
* starting address of the table.
536
ptab = ftp->ftable + indx;
537
*ptab = *p->xsig; /* Write the input value to the table. */
538
/* If this is guard point mode and we
539
* have just written to location 0,
540
* then also write to the guard point.
542
if ((p->iwgm == 2) && indx == 0) { /* Fix -- JPff 2000/1/5 */
549
/*---------------------------------------------------------------------------*/
551
/* tablew() is similar to ktablew() above, except that it processes
552
* two arrays of input values and indexes. These arrays are ksmps long. */
553
int ptablew(CSOUND *csound, TABLEW *p)
555
FUNC *ftp; /* Pointer to function table data structure. */
556
MYFLT *psig; /* Array of input values to be written to table. */
557
MYFLT *pxndx; /* Array of input index values */
558
MYFLT *ptab; /* Pointer to start of table we will write. */
559
MYFLT *pwrite;/* Pointer to location in table where we will write */
560
int32 indx; /* Used to read table. */
561
int32 length; /* Length of table */
562
int liwgm; /* Local copy of iwgm for speed */
563
int n, nsmps = csound->ksmps;
564
MYFLT ndx, xbmul, offset;
565
/*-----------------------------------*/
566
/* Assume that TABLEW has been set up correctly. */
574
xbmul = (MYFLT)p->xbmul;
576
/* Main loop - for the number of a samples in a k cycle. */
577
for (n=0; n<nsmps; n++) {
578
/* Read in the next raw index and increment the pointer ready for the
579
next cycle. Then multiply the ndx by the denormalising factor and
580
add in the offset. */
581
ndx = (pxndx[n] * xbmul) + offset;
582
if (liwgm == 0) { /* Limit mode - when igmode = 0. */
583
indx = (int32) MYFLOOR(ndx);
584
if (UNLIKELY(indx > length - 1)) indx = length - 1;
585
else if (UNLIKELY(indx < 0L)) indx = 0L;
588
if (liwgm == 2) ndx += FL(0.5);
589
indx = (int32) MYFLOOR(ndx);
590
/* Both wrap and guard point mode. */
591
if (UNLIKELY(indx>=length)) indx %= length;
592
else if (UNLIKELY(indx<0)) indx = length-(-indx)%length;
594
pwrite = ptab + indx;
596
/* If this is guard point mode and we
597
* have just written to location 0,
598
* then also write to the guard point.
600
if ((liwgm == 2) && indx == 0) { /* Fix -- JPff 2000/1/5 */
601
/* Note that since pwrite is a pointer
602
* to a float, adding length to it
603
* adds (4 * length) to its value since
604
* the length of a float is 4 bytes.
607
/* Decrement psig to make it point
608
* to the same input value.
609
* Write to guard point.
added
removed
OOps/ugens2a.c
