1 /*
2 ugens4.c:
3
4 Copyright (C) 1991 Barry Vercoe, John ffitch
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 #include "csoundCore.h" /* UGENS4.C */
25 #include "ugens4.h"
26 #include <math.h>
27 #include <inttypes.h>
28
29 /* The branch prediction slows it down!! */
30
bzzset(CSOUND * csound,BUZZ * p)31 int32_t bzzset(CSOUND *csound, BUZZ *p)
32 {
33 FUNC *ftp;
34
35 if (LIKELY((ftp = csound->FTFind(csound, p->ifn)) != NULL)) {
36 p->ftp = ftp;
37 if (*p->iphs >= 0)
38 p->lphs = (int32_t)(*p->iphs * FL(0.5) * FMAXLEN);
39 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0;
40 p->cpscod = IS_ASIG_ARG(p->xcps) ? 1 : 0;
41 p->reported = 0; /* No errors yet */
42 return OK;
43 }
44 return NOTOK;
45 }
46
buzz(CSOUND * csound,BUZZ * p)47 int32_t buzz(CSOUND *csound, BUZZ *p)
48 {
49 FUNC *ftp;
50 MYFLT *ar, *ampp, *cpsp, *ftbl;
51 int32_t phs, inc, lobits, dwnphs, tnp1, lenmask;
52 MYFLT sicvt2, over2n, scal, num, denom;
53 uint32_t offset = p->h.insdshead->ksmps_offset;
54 uint32_t early = p->h.insdshead->ksmps_no_end;
55 uint32_t n, nsmps = CS_KSMPS;
56 int32_t nn;
57
58 ftp = p->ftp;
59 if (UNLIKELY(ftp==NULL)) goto err1; /* RWD fix */
60 ftbl = ftp->ftable;
61 sicvt2 = csound->sicvt * FL(0.5); /* for theta/2 */
62 lobits = ftp->lobits;
63 lenmask = ftp->lenmask;
64 ampp = p->xamp;
65 cpsp = p->xcps;
66 if ((nn = (int32_t)*p->knh) < 0) nn = -nn;
67 if (UNLIKELY(nn == 0)) { /* fix nn = knh */
68 nn = 1;
69 }
70 tnp1 = (nn<<1) + 1; /* calc 2n + 1 */
71 over2n = FL(0.5) / (MYFLT)nn;
72 scal = *ampp * over2n;
73 inc = (int32_t)(*cpsp * sicvt2);
74 ar = p->ar;
75 phs = p->lphs;
76 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
77 if (UNLIKELY(early)) {
78 //printf("early=%d\n", early);
79 nsmps -= early;
80 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
81 }
82 for (n=offset; n<nsmps; n++) {
83 if (p->ampcod)
84 scal = ampp[n] * over2n;
85 if (p->cpscod)
86 inc = (int32_t)(cpsp[n] * sicvt2);
87 dwnphs = phs >> lobits;
88 denom = ftbl[dwnphs];
89 if (denom > FL(0.0002) || denom < -FL(0.0002)) {
90 num = ftbl[dwnphs * tnp1 & lenmask];
91 ar[n] = (num / denom - FL(1.0)) * scal;
92 }
93 else ar[n] = *ampp;
94 phs += inc;
95 phs &= PHMASK;
96 }
97 p->lphs = phs;
98 return OK;
99 err1:
100 return csound->PerfError(csound, &(p->h), Str("buzz: not initialised"));
101 }
102
gbzset(CSOUND * csound,GBUZZ * p)103 int32_t gbzset(CSOUND *csound, GBUZZ *p)
104 {
105 FUNC *ftp;
106
107 if (LIKELY((ftp = csound->FTFind(csound, p->ifn)) != NULL)) {
108 p->ftp = ftp;
109 if (*p->iphs >= 0) {
110 p->lphs = (int32_t)(*p->iphs * FMAXLEN);
111 p->prvr = FL(0.0);
112 }
113 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0;
114 p->cpscod = IS_ASIG_ARG(p->xcps) ? 1 : 0;
115 p->reported = 0;
116 p->last = FL(1.0);
117 return OK;
118 }
119 return NOTOK;
120 }
121
intpow1(MYFLT x,int32_t n)122 static inline MYFLT intpow1(MYFLT x, int32_t n) /* Binary positive power function */
123 {
124 MYFLT ans = FL(1.0);
125 while (n!=0) {
126 if (n&1) ans = ans * x;
127 n >>= 1;
128 x = x*x;
129 }
130 return ans;
131 }
132
intpow(MYFLT x,int32_t n)133 MYFLT intpow(MYFLT x, int32_t n) /* Binary power function */
134 {
135 if (n<0) {
136 n = -n;
137 x = FL(1.0)/x;
138 }
139 return intpow1(x, n);
140 }
141
gbuzz(CSOUND * csound,GBUZZ * p)142 int32_t gbuzz(CSOUND *csound, GBUZZ *p)
143 {
144 FUNC *ftp;
145 MYFLT *ar, *ampp, *cpsp, *ftbl;
146 int32_t phs, inc, lobits, lenmask, k, km1, kpn, kpnm1;
147 uint32_t offset = p->h.insdshead->ksmps_offset;
148 uint32_t early = p->h.insdshead->ksmps_no_end;
149 uint32_t n, nsmps = CS_KSMPS;
150 MYFLT r, absr, num, denom, scal, last = p->last;
151 int32_t nn, lphs = p->lphs;
152
153 ftp = p->ftp;
154 if (UNLIKELY(ftp==NULL)) goto err1;
155 ftbl = ftp->ftable;
156 lobits = ftp->lobits;
157 lenmask = ftp->lenmask;
158 ampp = p->xamp;
159 cpsp = p->xcps;
160 k = (int32_t)*p->kk; /* fix k and n */
161 if ((nn = (int32_t)*p->kn)<0) nn = -nn;
162 if (UNLIKELY(nn == 0)) { /* n must be > 0 */
163 nn = 1;
164 }
165 km1 = k - 1;
166 kpn = k + nn;
167 kpnm1 = kpn - 1;
168 if ((r = *p->kr) != p->prvr || nn != p->prvn) {
169 p->twor = r + r;
170 p->rsqp1 = r * r + FL(1.0);
171 p->rtn = intpow1(r, nn);
172 p->rtnp1 = p->rtn * r;
173 if ((absr = FABS(r)) > FL(0.999) && absr < FL(1.001))
174 p->rsumr = FL(1.0) / nn;
175 else p->rsumr = (FL(1.0) - absr) / (FL(1.0) - FABS(p->rtn));
176 p->prvr = r;
177 p->prvn = (int16)nn;
178 }
179 scal = *ampp * p->rsumr;
180 inc = (int32_t)(*cpsp * csound->sicvt);
181 ar = p->ar;
182 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
183 if (UNLIKELY(early)) {
184 nsmps -= early;
185 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
186 }
187 for (n=offset; n<nsmps; n++) {
188 if (p->ampcod)
189 scal = p->rsumr * ampp[n];
190 if (p->cpscod)
191 inc = (int32_t)(cpsp[n] * csound->sicvt);
192 phs = lphs >>lobits;
193 denom = p->rsqp1 - p->twor * ftbl[phs];
194 num = ftbl[phs * k & lenmask]
195 - r * ftbl[phs * km1 & lenmask]
196 - p->rtn * ftbl[phs * kpn & lenmask]
197 + p->rtnp1 * ftbl[phs * kpnm1 & lenmask];
198 if (LIKELY(denom > FL(0.0002) || denom < -FL(0.0002))) {
199 ar[n] = last = num / denom * scal;
200 }
201 else if (last<0)
202 ar[n] = last = - (p->ampcod ? ampp[n] : *ampp);
203 else
204 ar[n] = last = (p->ampcod ? ampp[n] : *ampp);
205 lphs += inc;
206 lphs &= PHMASK;
207 }
208 p->last = last;
209 p->lphs = lphs;
210 return OK;
211 err1:
212 return csound->PerfError(csound, &(p->h), Str("gbuzz: not initialised"));
213 }
214
215 #define PLUKMIN 64
216
217 static int16 rand15(CSOUND *);
218 static int16 rand16(CSOUND *);
219
plukset(CSOUND * csound,PLUCK * p)220 int plukset(CSOUND *csound, PLUCK *p)
221 {
222 int32_t n;
223 int32_t npts, iphs;
224 char *auxp;
225 FUNC *ftp;
226 MYFLT *ap, *fp;
227 MYFLT phs, phsinc;
228
229 if (UNLIKELY((npts = (int32_t)(csound->esr / *p->icps)) < PLUKMIN)) {
230 /* npts is wavelen in sampls */
231 npts = PLUKMIN; /* (but at least min size) */
232 }
233 if ((auxp = p->auxch.auxp) == NULL ||
234 npts > p->maxpts) { /* get newspace */
235 csound->AuxAlloc(csound, (npts+1)*sizeof(MYFLT),&p->auxch);
236 auxp = p->auxch.auxp;
237 p->maxpts = npts; /* if reqd */
238 }
239 ap = (MYFLT *)auxp; /* as MYFLT array */
240 if (*p->ifn == 0.0)
241 for (n=npts; n--; ) /* f0: fill w. rands */
242 *ap++ = (MYFLT) rand16(csound) * DV32768;
243 else if ((ftp = csound->FTnp2Finde(csound, p->ifn)) != NULL) {
244 fp = ftp->ftable; /* else from ftable */
245 phs = FL(0.0);
246 phsinc = (MYFLT)(ftp->flen/npts);
247 for (n=npts; n--; phs += phsinc) {
248 iphs = (int32_t)phs;
249 *ap++ = fp[iphs];
250 }
251 }
252 *ap = *(MYFLT *)auxp; /* last = copy of 1st */
253 p->npts = npts;
254 /* tuned pitch convt */
255 p->sicps = (npts * FL(256.0) + FL(128.0)) * csound->onedsr;
256 p->phs256 = 0;
257 p->method = (int16)*p->imeth;
258 p->param1 = *p->ipar1;
259 p->param2 = *p->ipar2;
260 switch(p->method) {
261 case 1: /* ignore any given parameters */
262 break;
263 case 2: /* stretch factor: param1 >= 1 */
264 if (UNLIKELY(p->param1 < FL(1.0)))
265 return csound->InitError(csound,
266 Str("illegal stretch factor(param1) value"));
267 else p->thresh1 = (int16)(FL(32768.0) / p->param1);
268 break;
269 case 3: /* roughness factor: 0 <= param1 <= 1 */
270 if (UNLIKELY(p->param1 < FL(0.0) || p->param1 > FL(1.0)))
271 return csound->InitError(csound,
272 Str("illegal roughness factor(param1) value"));
273 else
274 p->thresh1 = (int16)(FL(32768.0) * p->param1);
275 break;
276 case 4: /* rough and stretch factor: 0 <= param1 <= 1, param2 >= 1 */
277 if (UNLIKELY(p->param1 < FL(0.0) || p->param1 > FL(1.0)))
278 return csound->InitError(csound,
279 Str("illegal roughness factor(param1) value"));
280 else p->thresh1 = (int16)(FL(32768.0) * p->param1);
281 if (UNLIKELY(p->param2 < FL(1.0)))
282 return csound->InitError(csound,
283 Str("illegal stretch factor(param2) value"));
284 else p->thresh2 = (int16)(FL(32768.0) / p->param2);
285 break;
286 case 5: /* weighting coeff's: param1 + param2 <= 1 */
287 if (UNLIKELY(p->param1 + p->param2 > 1))
288 return csound->InitError(csound, Str("coefficients too large "
289 "(param1 + param2)"));
290 break;
291 case 6: /* ignore any given parameters */
292 break;
293
294 default:
295 return csound->InitError(csound, Str("unknown method code"));
296 }
297 return OK;
298 }
299
pluck(CSOUND * csound,PLUCK * p)300 int32_t pluck(CSOUND *csound, PLUCK *p)
301 {
302 MYFLT *ar, *fp;
303 int32_t phs256, phsinc, ltwopi, offset;
304 uint32_t koffset = p->h.insdshead->ksmps_offset;
305 uint32_t early = p->h.insdshead->ksmps_no_end;
306 uint32_t n, nsmps = CS_KSMPS;
307 MYFLT frac, diff;
308
309 if (UNLIKELY(p->auxch.auxp==NULL)) goto err1; /* RWD FIX */
310 ar = p->ar;
311 phsinc = (int32_t)(*p->kcps * p->sicps);
312 phs256 = p->phs256;
313 ltwopi = p->npts << 8;
314 if (UNLIKELY(phsinc > ltwopi)) goto err2;
315 if (UNLIKELY(koffset)) memset(ar, '\0', koffset*sizeof(MYFLT));
316 if (UNLIKELY(early)) {
317 nsmps -= early;
318 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
319 }
320 for (n=koffset; n<nsmps; n++) {
321 offset = phs256 >> 8;
322 fp = (MYFLT *)p->auxch.auxp + offset; /* lookup position */
323 diff = fp[1] - fp[0];
324 frac = (MYFLT)(phs256 & 255) / FL(256.0); /* w. interpolation */
325 ar[n] = (fp[0] + diff*frac) * *p->kamp; /* gives output val */
326 if ((phs256 += phsinc) >= ltwopi) {
327 int32_t nn;
328 MYFLT newval, preval;
329 phs256 -= ltwopi; /* at phase wrap, */
330 fp=(MYFLT *)p->auxch.auxp;
331 preval = fp[0]; /* copy last pnt */
332 fp[0] = fp[p->npts]; /* to first, */
333 fp++; /* apply smoothing */
334 nn = p->npts; /* up to npts+1 */
335 switch(p->method) {
336 case 1:
337 do { /* simple averaging */
338 newval = (*fp + preval) * FL(0.5);
339 preval = *fp;
340 *fp++ = newval;
341 } while (--nn);
342 break;
343 case 2:
344 do { /* stretched avrging */
345 if (rand15(csound) < p->thresh1) {
346 newval = (*fp + preval) * FL(0.5);
347 preval = *fp;
348 *fp++ = newval;
349 }
350 else preval = *fp++;
351 } while (--nn);
352 break;
353 case 3:
354 do { /* simple drum */
355 if (rand15(csound) < p->thresh1)
356 newval = -(*fp + preval) * FL(0.5);
357 else
358 newval = (*fp + preval) * FL(0.5);
359 preval = *fp;
360 *fp++ = newval;
361 } while (--nn);
362 break;
363 case 4:
364 do { /* stretched drum */
365 if (rand15(csound) < p->thresh2) {
366 if (rand15(csound) < p->thresh1)
367 newval = -(*fp + preval) * FL(0.5);
368 else
369 newval = (*fp + preval) * FL(0.5);
370 preval = *fp;
371 *fp++ = newval;
372 }
373 else preval = *fp++;
374 } while (--nn);
375 break;
376 case 5:
377 do { /* weighted avraging */
378 newval = p->param1 * *fp + p->param2 * preval;
379 preval = *fp;
380 *fp++ = newval;
381 } while (--nn);
382 break;
383 case 6:
384 do { /* 1st order recursive filter*/
385 preval = (*fp + preval) * FL(0.5);
386 *fp++ = preval;
387 } while (--nn);
388 break;
389 }
390 }
391 }
392 p->phs256 = phs256;
393 return OK;
394 err1:
395 return csound->PerfError(csound, &(p->h), Str("pluck: not initialised"));
396 err2:
397 return csound->PerfError(csound, &(p->h),
398 Str("pluck: kcps more than sample rate"));
399 }
400
401 #define RNDMUL 15625
402 #define MASK16 0xFFFF
403 #define MASK15 0x7FFF
404
405 /* quick generate a random int16 between -32768 and 32767 */
406
rand16(CSOUND * csound)407 static int16 rand16(CSOUND *csound)
408 {
409 csound->ugens4_rand_16 = (csound->ugens4_rand_16 * RNDMUL + 1) & MASK16;
410 return (int16) ((int16_t) csound->ugens4_rand_16);
411 }
412
413 /* quick generate a random int16 between 0 and 32767 */
414
rand15(CSOUND * csound)415 static int16 rand15(CSOUND *csound)
416 {
417 csound->ugens4_rand_15 = (csound->ugens4_rand_15 * RNDMUL + 1) & MASK15;
418 return (int16) csound->ugens4_rand_15;
419 }
420
421 /*=========================================================================
422 *
423 * randint31()
424 *
425 * 31 bit Park Miller PRNG using Linus Schrage's method for doing it all
426 * with 32 bit variables.
427 *
428 * Code adapted from Ray Garder's public domain code of July 1997 at:
429 * http://www.snippets.org/RG_RAND.C Thanks!
430 *
431 * Based on "Random Number Generators: Good Ones Are Hard to Find",
432 * S.K. Park and K.W. Miller, Communications of the ACM 31:10 (Oct 1988),
433 * and "Two Fast Implementations of the 'Minimal Standard' Random
434 * Number Generator", David G. Carta, Comm. ACM 33, 1 (Jan 1990), p. 87-88
435 *
436 * Linear congruential generator: f(z) = (16807 * z) mod (2 ** 31 - 1)
437 *
438 * Uses L. Schrage's method to avoid overflow problems.
439 */
440
441 #define BIPOLAR 0x7FFFFFFF /* Constant to make bipolar */
442 #define RIA (16807) /* multiplier */
443 #define RIM 0x7FFFFFFFL /* 2**31 - 1 */
444
445 #define dv2_31 (FL(4.656612873077392578125e-10))
446
randint31(int32 seed31)447 int32 randint31(int32 seed31)
448 {
449 uint32 rilo, rihi;
450
451 rilo = RIA * (int32_t)(seed31 & 0xFFFF);
452 rihi = RIA * (int32_t)((uint32)seed31 >> 16);
453 rilo += (rihi & 0x7FFF) << 16;
454 if (rilo > RIM) {
455 rilo &= RIM;
456 ++rilo;
457 }
458 rilo += rihi >> 15;
459 if (rilo > RIM) {
460 rilo &= RIM;
461 ++rilo;
462 }
463 return (int32_t)rilo;
464 }
465
rndset(CSOUND * csound,RAND * p)466 int32_t rndset(CSOUND *csound, RAND *p)
467 {
468 p->new = (*p->sel!=FL(0.0));
469 if (*p->iseed >= FL(0.0)) {
470 if (*p->iseed > FL(1.0)) { /* As manual suggest seed in range [0,1] */
471 uint32 seed; /* I reinterpret >1 as a time seed */
472 seed = csound->GetRandomSeedFromTime();
473 csound->Warning(csound, Str("Seeding from current time %"PRIu32"\n"), seed);
474 if (!p->new) {
475 p->rand = (int32_t) (seed & 0xFFFFUL);
476 }
477 else {
478 p->rand = (int32_t) (seed % 0x7FFFFFFEUL) + 1L;
479 }
480 }
481 else {
482 if (p->new) {
483 MYFLT seed = *p->iseed;
484 if (seed==FL(0.0)) seed = FL(0.5);
485 p->rand = (int32_t) (seed * FL(2147483648.0));
486 p->rand = randint31(p->rand);
487 p->rand = randint31(p->rand);
488 }
489 else
490 p->rand = ((int16)(*p->iseed * FL(32768.0)))&0xffff;
491 }
492 }
493 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0; /* (not used by krand) */
494
495 return OK;
496 }
497
krand(CSOUND * csound,RAND * p)498 int32_t krand(CSOUND *csound, RAND *p)
499 {
500 IGN(csound);
501 if (p->new) {
502 int32_t r = randint31(p->rand); /* result is a 31-bit value */
503 p->rand = r;
504 *p->ar = *p->base +
505 dv2_31 * (MYFLT)((int32_t)((uint32_t)r<<1)-BIPOLAR) * *p->xamp;
506 }
507 else {
508 int16 rand = (int16)p->rand;
509 rand *= RNDMUL;
510 rand += 1;
511 /* IV - Jul 11 2002 */
512 *p->ar = *p->base + (MYFLT)rand * *p->xamp * DV32768;
513 p->rand = rand;
514 }
515 return OK;
516 }
517
arand(CSOUND * csound,RAND * p)518 int32_t arand(CSOUND *csound, RAND *p)
519 {
520 IGN(csound);
521 MYFLT *ar;
522 int16 rndmul = RNDMUL;
523 uint32_t offset = p->h.insdshead->ksmps_offset;
524 uint32_t early = p->h.insdshead->ksmps_no_end;
525 uint32_t n, nsmps = CS_KSMPS;
526 MYFLT ampscl;
527 MYFLT base = *p->base;
528
529 ar = p->ar;
530 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
531 if (UNLIKELY(early)) {
532 nsmps -= early;
533 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
534 }
535 if (!p->new) {
536 int16 rand = p->rand;
537 if (!(p->ampcod)) {
538 ampscl = *p->xamp * DV32768; /* IV - Jul 11 2002 */
539 for (n=offset;n<nsmps;n++) {
540 rand *= rndmul;
541 rand += 1;
542 ar[n] = base + (MYFLT)rand * ampscl; /* IV - Jul 11 2002 */
543 }
544 }
545 else {
546 MYFLT *xamp = p->xamp;
547 for (n=offset;n<nsmps;n++) {
548 rand *= rndmul;
549 rand += 1;
550 /* IV - Jul 11 2002 */
551 ar[n] = base + (MYFLT)rand * xamp[n] * DV32768;
552 }
553 }
554 p->rand = rand; /* save current rand */
555 }
556 else {
557 int32_t rand = p->rand;
558 if (!(p->ampcod)) {
559 ampscl = *p->xamp * dv2_31;
560 for (n=offset;n<nsmps;n++) {
561 rand = randint31(rand);
562 ar[n] = base + (MYFLT)((int32_t)((uint32_t)rand<<1)-BIPOLAR) * ampscl;
563 }
564 }
565 else {
566 MYFLT *xamp = p->xamp;
567 for (n=offset;n<nsmps;n++) {
568 rand = randint31(rand);
569 ar[n] = base +
570 dv2_31 * (MYFLT)((int32_t)((uint32_t)rand<<1)-BIPOLAR) * xamp[n];
571 }
572 }
573 p->rand = rand; /* save current rand */
574 }
575
576 return OK;
577 }
578
rhset(CSOUND * csound,RANDH * p)579 int32_t rhset(CSOUND *csound, RANDH *p)
580 {
581 p->new = (*p->sel!=FL(0.0));
582 if (*p->iseed >= FL(0.0)) { /* new seed: */
583 if (*p->iseed > FL(1.0)) { /* As manual suggest sseed in range [0,1] */
584 uint32 seed; /* I reinterpret >1 as a time seed */
585 seed = csound->GetRandomSeedFromTime();
586 csound->Warning(csound, Str("Seeding from current time %"PRIu32"\n"), seed);
587 if (!p->new) {
588 p->rand = (int32_t) (seed & 0xFFFFUL);
589 p->num1 = (MYFLT) ((int16) p->rand) * DV32768;
590 }
591 else {
592 p->rand = (int32_t) (seed % 0x7FFFFFFEUL) + 1L;
593 p->num1 = (MYFLT) ((int32_t) ((uint32_t)p->rand<<1) - BIPOLAR) * dv2_31;
594 }
595 }
596 else if (!p->new) {
597 p->rand = 0xffff&(int16)(*p->iseed * 32768L); /* init rand integ */
598 p->num1 = *p->iseed; /* store fnum */
599 }
600 else {
601 MYFLT ss = *p->iseed;
602 if (ss>FL(1.0)) p->rand = (int32_t) ss;
603 else p->rand = (int32_t) (*p->iseed * FL(2147483648.0));
604 p->rand = randint31(p->rand);
605 p->rand = randint31(p->rand);
606 p->num1 = (MYFLT) ((int32_t) ((uint32_t)p->rand<<1) - BIPOLAR) * dv2_31;
607 }
608 p->phs = 0; /* & phs */
609 }
610 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0; /* (not used by krandh) */
611 p->cpscod = IS_ASIG_ARG(p->xcps) ? 1 : 0;
612 return OK;
613 }
614
krandh(CSOUND * csound,RANDH * p)615 int32_t krandh(CSOUND *csound, RANDH *p)
616 {
617 IGN(csound);
618 *p->ar = *p->base + p->num1 * *p->xamp; /* rslt = num * amp */
619 p->phs += (int64_t)(*p->xcps * CS_KICVT); /* phs += inc */
620
621 if (p->phs >= MAXLEN) { /* when phs overflows, */
622 p->phs &= PHMASK; /* mod the phs */
623 if (!p->new) {
624 int16 rand = (int16)p->rand;
625 rand *= RNDMUL; /* & recalc number */
626 rand += 1;
627 p->num1 = (MYFLT)rand * DV32768; /* IV - Jul 11 2002 */
628 p->rand = rand;
629 }
630 else {
631 int32_t r = randint31(p->rand); /* & recalc number */
632 p->rand = r;
633 p->num1 = (MYFLT)((int32_t)((uint32_t)r<<1)-BIPOLAR) * dv2_31;
634 }
635 }
636 return OK;
637 }
638
randh(CSOUND * csound,RANDH * p)639 int32_t randh(CSOUND *csound, RANDH *p)
640 {
641 int64_t phs = p->phs, inc;
642 uint32_t offset = p->h.insdshead->ksmps_offset;
643 uint32_t early = p->h.insdshead->ksmps_no_end;
644 uint32_t n, nsmps = CS_KSMPS;
645 MYFLT *ar, *ampp, *cpsp;
646 MYFLT base = *p->base;
647
648 cpsp = p->xcps;
649 ampp = p->xamp;
650 ar = p->ar;
651 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
652 if (UNLIKELY(early)) {
653 nsmps -= early;
654 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
655 }
656 inc = (int64_t)(*cpsp++ * csound->sicvt);
657 for (n=offset;n<nsmps;n++) {
658 /* IV - Jul 11 2002 */
659 ar[n] = base + p->num1 * *ampp; /* rslt = num * amp */
660 if (p->ampcod)
661 ampp++;
662 phs += inc; /* phs += inc */
663 if (p->cpscod)
664 inc = (int64_t)(*cpsp++ * csound->sicvt);
665 if (phs >= MAXLEN) { /* when phs o'flows, */
666 phs &= PHMASK;
667 if (!p->new) {
668 int16 rand = p->rand;
669 rand *= RNDMUL; /* calc new number */
670 rand += 1;
671 p->num1 = (MYFLT)rand * DV32768; /* IV - Jul 11 2002 */
672 p->rand = rand;
673 }
674 else {
675 int32_t r = randint31(p->rand); /* calc new number */
676 p->rand = r;
677 p->num1 = (MYFLT)((int32_t)((uint32_t)r<<1)-BIPOLAR) * dv2_31;
678 }
679 }
680 }
681 p->phs = phs;
682 return OK;
683 }
684
riset(CSOUND * csound,RANDI * p)685 int32_t riset(CSOUND *csound, RANDI *p)
686 {
687 p->new = (*p->sel!=FL(0.0));
688 if (*p->iseed >= FL(0.0)) { /* new seed: */
689 if (*p->iseed > FL(1.0)) { /* As manual suggest seed in range [0,1] */
690 uint32 seed; /* I reinterpret >1 as a time seed */
691 seed = csound->GetRandomSeedFromTime();
692 csound->Warning(csound, Str("Seeding from current time %"PRIu32"\n"), seed);
693 if (!p->new) {
694 int16 rand = (int16)seed;
695 /* int16 ss = rand; */
696 /* IV - Jul 11 2002 */
697 p->num1 = (MYFLT)(rand) * DV32768; /* store num1,2 */
698 rand *= RNDMUL; /* recalc random */
699 rand += 1;
700 /* IV - Jul 11 2002 */
701 p->num2 = (MYFLT)(p->rand=rand) * DV32768;
702 /* printf("seed, rand, num1, num2 = %d(%x), %d*%x), %f, %f\n", */
703 /* ss,ss,p->rand, p->rand, p->num1, p->num2); */
704 }
705 else {
706 p->rand = randint31((int32_t) (seed % 0x7FFFFFFEUL) + 1L);
707 p->rand = randint31(p->rand);
708 p->num1 = (MYFLT)(p->rand<<1) * dv2_31; /* store num1,2 */
709 p->rand = randint31(p->rand);
710 p->num2 = (MYFLT)(p->rand<<1) * dv2_31;
711 }
712 }
713 else if (!p->new) {
714 int16 rand = (int16)(*p->iseed * FL(32768.0)); /* init rand integ */
715 rand *= RNDMUL; /* to 2nd value */
716 rand += 1;
717 p->num1 = *p->iseed; /* store num1,2 */
718 p->num2 = (MYFLT)rand * DV32768; /* IV - Jul 11 2002 */
719 p->rand = rand;
720 }
721 else {
722 MYFLT ss = *p->iseed;
723 if (ss>FL(1.0)) p->rand = (int32_t) ss;
724 else if (ss==FL(0.0)) p->rand = (int32_t) (FL(0.5) * FL(2147483648.0));
725 else p->rand = (int32_t) (ss * FL(2147483648.0));
726 p->rand = randint31(p->rand);
727 p->rand = randint31(p->rand);
728 p->num1 = (MYFLT)(p->rand<1) * dv2_31; /* store num1,2 */
729 p->rand = randint31(p->rand);
730 p->num2 = (MYFLT)(p->rand<<1) * dv2_31;
731 }
732 p->phs = 0; /* & clear phs */
733 p->dfdmax = (p->num2 - p->num1) / FMAXLEN; /* & diff */
734 }
735 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0; /* (not used by krandi) */
736 p->cpscod = IS_ASIG_ARG(p->xcps) ? 1 : 0;
737 return OK;
738 }
739
krandi(CSOUND * csound,RANDI * p)740 int32_t krandi(CSOUND *csound, RANDI *p)
741 { /* rslt = (num1 + diff*phs) * amp */
742 IGN(csound);
743 *p->ar = *p->base + (p->num1 + (MYFLT)p->phs * p->dfdmax) * *p->xamp;
744 p->phs += (int64_t)(*p->xcps * CS_KICVT); /* phs += inc */
745 if (p->phs >= MAXLEN) { /* when phs overflows, */
746 p->phs &= PHMASK; /* mod the phs */
747 if (!p->new) {
748 int16 rand = p->rand;
749 rand *= RNDMUL; /* recalc random */
750 rand += 1;
751 p->num1 = p->num2; /* & new num vals */
752 p->num2 = (MYFLT)rand * DV32768; /* IV - Jul 11 2002 */
753 p->rand = rand;
754 }
755 else {
756 int32_t r = randint31(p->rand); /* recalc random */
757 p->rand = r;
758 p->num1 = p->num2; /* & new num vals */
759 p->num2 = (MYFLT)((int32_t)((uint32_t)r<<1)-BIPOLAR) * dv2_31;
760 }
761 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
762 }
763 return OK;
764 }
765
randi(CSOUND * csound,RANDI * p)766 int32_t randi(CSOUND *csound, RANDI *p)
767 {
768 int64_t phs = p->phs, inc;
769 uint32_t offset = p->h.insdshead->ksmps_offset;
770 uint32_t early = p->h.insdshead->ksmps_no_end;
771 uint32_t n, nsmps = CS_KSMPS;
772 MYFLT *ar, *ampp, *cpsp;
773 MYFLT base = *p->base;
774
775 cpsp = p->xcps;
776 ampp = p->xamp;
777 ar = p->ar;
778 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
779 if (UNLIKELY(early)) {
780 nsmps -= early;
781 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
782 }
783 inc = (int64_t)(*cpsp++ * csound->sicvt);
784 for (n=offset;n<nsmps;n++) {
785 /* IV - Jul 11 2002 */
786 ar[n] = base + (p->num1 + (MYFLT)phs * p->dfdmax) * *ampp;
787 if (p->ampcod)
788 ampp++;
789 phs += inc; /* phs += inc */
790 if (p->cpscod)
791 inc = (int64_t)(*cpsp++ * csound->sicvt); /* (nxt inc) */
792 if (phs >= MAXLEN) { /* when phs o'flows, */
793 phs &= PHMASK;
794 if (!p->new) {
795 int16 rand = p->rand;
796 rand *= RNDMUL; /* calc new numbers*/
797 rand += 1;
798 p->num1 = p->num2;
799 p->num2 = (MYFLT)rand * DV32768; /* IV - Jul 11 2002 */
800 p->rand = rand;
801 }
802 else {
803 int32_t r = randint31(p->rand); /* calc new numbers*/
804 //printf("r = %x\n", r);
805 p->rand = r;
806 p->num1 = p->num2;
807 p->num2 = (MYFLT)((int32_t)((uint32_t)r<<1)-BIPOLAR) * dv2_31;
808 }
809 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
810 }
811 }
812 p->phs = phs;
813 return OK;
814 }
815
816 /* Cubic interpolation between random values -- JPff 2019 */
817
rcset(CSOUND * csound,RANDC * p)818 int32_t rcset(CSOUND *csound, RANDC *p)
819 {
820 p->new = (*p->sel!=FL(0.0));
821 if (*p->iseed >= FL(0.0)) { /* new seed: */
822 if (*p->iseed > FL(1.0)) { /* As manual suggest sseed in range [0,1] */
823 uint32 seed; /* I reinterpret >1 as a time seed */
824 seed = csound->GetRandomSeedFromTime();
825 csound->Warning(csound, Str("Seeding from current time %"PRIu32"\n"), seed);
826 if (!p->new) {
827 int16 rand = (int16)seed;
828 /* int16 ss = rand; */
829 /* IV - Jul 11 2002 */
830 p->num1 = (MYFLT)(rand) * DV32768; /* store num1,2 */
831 rand *= RNDMUL; /* recalc random */
832 rand += 1;
833 p->num2 = (MYFLT)(rand) * DV32768;
834 rand *= RNDMUL; /* recalc random */
835 rand += 1;
836 p->num3 = (MYFLT)(rand) * DV32768;
837 rand *= RNDMUL; /* recalc random */
838 rand += 1;
839 p->num4 = (MYFLT)(p->rand=rand) * DV32768;
840 }
841 else {
842 p->rand = randint31((int32_t) (seed % 0x7FFFFFFEUL) + 1L);
843 p->rand = randint31(p->rand);
844 p->num1 = (MYFLT)(p->rand<<1) * dv2_31; /* store num1,2 */
845 p->rand = randint31(p->rand);
846 p->num2 = (MYFLT)(p->rand<<1) * dv2_31;
847 p->rand = randint31(p->rand);
848 p->num3 = (MYFLT)(p->rand<<1) * dv2_31;
849 p->rand = randint31(p->rand);
850 p->num4 = (MYFLT)(p->rand<<1) * dv2_31;
851 }
852 }
853 else if (!p->new) {
854 int16 rand = (int16)(*p->iseed * FL(32768.0)); /* init rand integ */
855 rand *= RNDMUL; /* to 2nd value */
856 rand += 1;
857 p->num1 = *p->iseed; /* store num1,2 */
858 p->num2 = (MYFLT)rand * DV32768; /* IV - Jul 11 2002 */
859 rand *= RNDMUL; /* to 2nd value */
860 rand += 1;
861 p->num3 = (MYFLT)rand * DV32768;
862 rand *= RNDMUL; /* to 2nd value */
863 rand += 1;
864 p->num4 = (MYFLT)rand * DV32768;
865 p->rand = rand;
866 }
867 else {
868 MYFLT ss = *p->iseed;
869 if (ss>FL(1.0)) p->rand = (int32_t) ss;
870 else p->rand = (int32_t) (*p->iseed * FL(2147483648.0));
871 p->rand = randint31(p->rand);
872 p->rand = randint31(p->rand);
873 p->num1 = (MYFLT)(p->rand<1) * dv2_31; /* store num1,2 */
874 p->rand = randint31(p->rand);
875 p->num2 = (MYFLT)(p->rand<<1) * dv2_31;
876 p->rand = randint31(p->rand);
877 p->num3 = (MYFLT)(p->rand<<1) * dv2_31;
878 p->rand = randint31(p->rand);
879 p->num4 = (MYFLT)(p->rand<<1) * dv2_31;
880 }
881 }
882 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0; /* (not used by krandi) */
883 p->cpscod = IS_ASIG_ARG(p->xcps) ? 1 : 0;
884 p->phs = 0;
885 return OK;
886 }
887
888
krandc(CSOUND * csound,RANDC * p)889 int32_t krandc(CSOUND *csound, RANDC *p)
890 { /* rslt = (num1 + diff*phs) * amp */
891 IGN(csound);
892 MYFLT a0 = p->num4 - p->num3 - p->num1 + p->num2;
893 MYFLT a1 = p->num1 - p->num2 - a0;
894 MYFLT a2 = p->num3 - p->num1;
895 MYFLT a3 = p->num2;
896 MYFLT mu = (MYFLT)p->phs / (MYFLT)MAXLEN;
897 *p->ar = *p->base + (((a0 * mu +a1) * mu+a2) * mu + a3) * *p->xamp;
898 p->phs += (int64_t)(*p->xcps * CS_KICVT); /* phs += inc */
899 if (p->phs >= MAXLEN) { /* when phs overflows, */
900 p->phs &= PHMASK; /* mod the phs */
901 if (!p->new) {
902 int16 rand = p->rand;
903 rand *= RNDMUL; /* recalc random */
904 rand += 1;
905 p->num1 = p->num2; /* & new num vals */
906 p->num2 = p->num3;
907 p->num3 = p->num4;
908 p->num4 = (MYFLT)rand * DV32768;
909 p->rand = rand;
910 }
911 else {
912 int32_t r = randint31(p->rand); /* recalc random */
913 p->rand = r;
914 p->num1 = p->num2; /* & new num vals */
915 p->num2 = p->num3;
916 p->num3 = p->num4;
917 p->num4 = (MYFLT)((int32_t)((uint32_t)r<<1)-BIPOLAR) * dv2_31;
918 }
919 }
920 return OK;
921 }
922
randc(CSOUND * csound,RANDC * p)923 int32_t randc(CSOUND *csound, RANDC *p)
924 {
925 int64_t phs = p->phs, inc;
926 uint32_t offset = p->h.insdshead->ksmps_offset;
927 uint32_t early = p->h.insdshead->ksmps_no_end;
928 uint32_t n, nsmps = CS_KSMPS;
929 MYFLT *ar, *ampp, *cpsp;
930 MYFLT mu;
931 MYFLT base = *p->base;
932 MYFLT a0 = p->num4 - p->num3 - p->num1 + p->num2;
933 MYFLT a1 = p->num1 - p->num2 - a0;
934 MYFLT a2 = p->num3 - p->num1;
935 MYFLT a3 = p->num2;
936 cpsp = p->xcps;
937 ampp = p->xamp;
938 inc = (int64_t)(*cpsp++ * csound->sicvt);
939 ar = p->ar;
940 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
941 if (UNLIKELY(early)) {
942 nsmps -= early;
943 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
944 }
945
946 for (n=offset;n<nsmps;n++) {
947
948 mu = (MYFLT)phs/(MYFLT)MAXLEN;
949
950 ar[n] = base + (((a0 * mu +a1) * mu+a2) * mu + a3) * *ampp;
951
952 if (p->ampcod)
953 ampp++;
954 phs += inc;
955 //printf("mu = %g phs, inc, MAXLEN = %ld, %ld, %d\n", mu, phs, inc, MAXLEN);
956 if (p->cpscod)
957 inc = (int64_t)(*cpsp++ * csound->sicvt); /* (nxt inc) */
958 if (phs >= MAXLEN) { /* when phs o'flows, */
959 phs &= PHMASK;
960 if (!p->new) {
961 int16 rand = p->rand;
962 rand *= RNDMUL; /* calc new numbers*/
963 rand += 1;
964 p->num1 = p->num2; /* & new num vals */
965 p->num2 = p->num3;
966 p->num3 = p->num4;
967 p->num4 = (MYFLT)rand * DV32768;
968 p->rand = rand;
969 }
970 else {
971 int32_t r = randint31(p->rand); /* calc new numbers*/
972 p->rand = r;
973 p->num1 = p->num2;
974 p->num2 = p->num3;
975 p->num3 = p->num4;
976 p->num4 = (MYFLT)((int32_t)((uint32_t)r<<1)-BIPOLAR) * dv2_31;
977 }
978 a0 = p->num4 - p->num3 - p->num1 + p->num2;
979 a1 = p->num1 - p->num2 - a0;
980 a2 = p->num3 - p->num1;
981 a3 = p->num2;
982
983 }
984 }
985 p->phs = phs;
986 return OK;
987 }
988
989