1
2 /*
3 pitch.c:
4
5 Copyright (C) 1999 John ffitch, Istvan Varga, Peter Neubäcker,
6 rasmus ekman, Phil Burk
7
8 This file is part of Csound.
9
10 The Csound Library is free software; you can redistribute it
11 and/or modify it under the terms of the GNU Lesser General Public
12 License as published by the Free Software Foundation; either
13 version 2.1 of the License, or (at your option) any later version.
14
15 Csound is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU Lesser General Public License for more details.
19
20 You should have received a copy of the GNU Lesser General Public
21 License along with Csound; if not, write to the Free Software
22 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
23 02110-1301 USA
24 */
25
26 // #include "csdl.h"
27 #include "csoundCore.h" /* PITCH.C */
28 #include <math.h>
29 #include <limits.h>
30 #include "cwindow.h"
31 #include "spectra.h"
32 #include "pitch.h"
33 #include "uggab.h"
34 #include <inttypes.h>
35
36 #define STARTING 1
37 #define PLAYING 2
38 #define LOGTWO (0.69314718055994530942)
39
MYFLOOR(MYFLT x)40 static inline int32 MYFLOOR(MYFLT x) {
41 if (x >= 0.0) {
42 return (int32) x;
43 } else {
44 return (int32) (x - FL(0.99999999));
45 }
46 }
47
48 static const MYFLT bicoefs[] = {
49 -FL(0.2674054), FL(0.7491305), FL(0.7160484), FL(0.0496285), FL(0.7160484),
50 FL(0.0505247), FL(0.3514850), FL(0.5257536), FL(0.3505025), FL(0.5257536),
51 FL(0.3661840), FL(0.0837990), FL(0.3867783), FL(0.6764264), FL(0.3867783)
52 };
53
54 #define rand_31(x) (x->Rand31(&(x->randSeed1)) - 1)
55
pitchset(CSOUND * csound,PITCH * p)56 int32_t pitchset(CSOUND *csound, PITCH *p) /* pitch - uses spectra technology */
57 {
58 double b; /* For RMS */
59 int32_t n, nocts, nfreqs, ncoefs;
60 MYFLT Q, *fltp;
61 OCTDAT *octp;
62 DOWNDAT *dwnp = &p->downsig;
63 SPECDAT *specp = &p->wsig;
64 int32 npts, nptls, nn, lobin;
65 int32_t *dstp, ptlmax;
66 MYFLT fnfreqs, rolloff, *oct0p, *flop, *fhip, *fundp, *fendp, *fp;
67 MYFLT weight, weightsum, dbthresh, ampthresh;
68
69 /* RMS of input signal */
70 b = 2.0 - cos(10.0*(double)csound->tpidsr);
71 p->c2 = b - sqrt(b * b - 1.0);
72 p->c1 = 1.0 - p->c2;
73 if (!*p->istor) p->prvq = 0.0;
74 /* End of rms */
75 /* Initialise spectrum */
76 /* for mac roundoff */
77 p->timcount = (int32_t)(CS_EKR * *p->iprd + FL(0.001));
78 nocts = (int32_t)*p->iocts; if (UNLIKELY(nocts<=0)) nocts = 6;
79 nfreqs = (int32_t)*p->ifrqs; if (UNLIKELY(nfreqs<=0)) nfreqs = 12;
80 ncoefs = nocts * nfreqs;
81 Q = *p->iq; if (UNLIKELY(Q<=FL(0.0))) Q = FL(15.0);
82
83 if (UNLIKELY(p->timcount <= 0))
84 return csound->InitError(csound, Str("illegal iprd"));
85 if (UNLIKELY(nocts > MAXOCTS))
86 return csound->InitError(csound, Str("illegal iocts"));
87 if (UNLIKELY(nfreqs > MAXFRQS))
88 return csound->InitError(csound, Str("illegal ifrqs"));
89
90 if (nocts != dwnp->nocts ||
91 nfreqs != p->nfreqs || /* if anything has changed */
92 Q != p->curq ) { /* make new tables */
93 double basfrq, curfrq, frqmlt, Qfactor;
94 double theta, a, windamp, onedws, pidws;
95 MYFLT *sinp, *cosp;
96 int32_t k, sumk, windsiz, halfsiz, *wsizp, *woffp;
97 int32 auxsiz, bufsiz;
98 int32 majr, minr, totsamps;
99 double hicps,locps,oct; /* must alloc anew */
100
101 p->nfreqs = nfreqs;
102 p->curq = Q;
103 p->ncoefs = ncoefs;
104 dwnp->srate = CS_ESR;
105 hicps = dwnp->srate * 0.375; /* top freq is 3/4 pi/2 ... */
106 oct = log(hicps / ONEPT) / LOGTWO; /* octcps() (see aops.c) */
107 dwnp->looct = (MYFLT)(oct - nocts); /* true oct val of lowest frq */
108 locps = hicps / (1L << nocts);
109 basfrq = hicps * 0.5; /* oct below retuned top */
110 frqmlt = pow(2.0,1.0/(double)nfreqs); /* nfreq interval mult */
111 Qfactor = Q * dwnp->srate;
112 curfrq = basfrq;
113 for (sumk=0,wsizp=p->winlen,woffp=p->offset,n=nfreqs; n--; ) {
114 *wsizp++ = k = (int32_t)(Qfactor/curfrq) | 01; /* calc odd wind sizes */
115 *woffp++ = (*(p->winlen) - k) / 2; /* & symmetric offsets */
116 sumk += k; /* and find total */
117 curfrq *= frqmlt;
118 }
119 windsiz = *(p->winlen);
120 auxsiz = (windsiz + 2*sumk) * sizeof(MYFLT); /* calc lcl space rqd */
121
122 csound->AuxAlloc(csound, (size_t)auxsiz, &p->auxch1); /* & alloc auxspace */
123
124 fltp = (MYFLT *) p->auxch1.auxp;
125 p->linbufp = fltp; fltp += windsiz; /* linbuf must take nsamps */
126 p->sinp = sinp = fltp; fltp += sumk;
127 p->cosp = cosp = fltp; /* cos gets rem sumk */
128 wsizp = p->winlen;
129 curfrq = basfrq * TWOPI / dwnp->srate;
130 for (n = nfreqs; n--; ) { /* now fill tables */
131 windsiz = *wsizp++; /* (odd win size) */
132 halfsiz = windsiz >> 1;
133 onedws = 1.0 / (windsiz-1);
134 pidws = PI / (windsiz-1);
135 for (k = -halfsiz; k<=halfsiz; k++) { /* with sines */
136 a = cos(k * pidws);
137 windamp = 0.08 + 0.92 * a * a; /* times hamming */
138 windamp *= onedws; /* scaled */
139 theta = k * curfrq;
140 *sinp++ = (MYFLT)(windamp * sin(theta));
141 *cosp++ = (MYFLT)(windamp * cos(theta));
142 }
143 curfrq *= frqmlt; /* step by log freq */
144 }
145 dwnp->hifrq = (MYFLT)hicps;
146 dwnp->lofrq = (MYFLT)locps;
147 dwnp->nsamps = windsiz = *(p->winlen);
148 dwnp->nocts = nocts;
149 minr = windsiz >> 1; /* sep odd windsiz into maj, min */
150 majr = windsiz - minr; /* & calc totsamps reqd */
151 totsamps = (majr*nocts) + (minr<<nocts) - minr;
152 DOWNset(csound, dwnp, totsamps); /* auxalloc in DOWNDAT struct */
153 fltp = (MYFLT *) dwnp->auxch.auxp; /* & distrib to octdata */
154 for (n=nocts,octp=dwnp->octdata+(nocts-1); n--; octp--) {
155 bufsiz = majr + minr;
156 octp->begp = fltp; fltp += bufsiz; /* (lo oct first) */
157 octp->endp = fltp; minr *= 2;
158 }
159 SPECset(csound, specp, (int32)ncoefs);/* prep the spec dspace */
160 specp->downsrcp = dwnp; /* & record its source */
161 }
162 for (octp=dwnp->octdata; nocts--; octp++) { /* reset all oct params, & */
163 octp->curp = octp->begp;
164 memset(octp->feedback, '\0', 6*sizeof(MYFLT));
165 octp->scount = 0;
166 }
167 specp->nfreqs = p->nfreqs; /* save the spec descriptors */
168 specp->dbout = 0;
169 specp->ktimstamp = 0; /* init specdata to not new */
170 specp->ktimprd = p->timcount;
171 p->scountdown = p->timcount; /* prime the spect countdown */
172 /* Start specptrk */
173 if ((npts = specp->npts) != p->winpts) { /* if size has changed */
174 SPECset(csound, &p->wfund, (int32)npts); /* realloc for wfund */
175 p->wfund.downsrcp = specp->downsrcp;
176 p->fundp = (MYFLT *) p->wfund.auxch.auxp;
177 p->winpts = npts;
178 }
179 if (UNLIKELY(*p->inptls<=FL(0.0))) nptls = 4;
180 else nptls = (int32)*p->inptls;
181 if (UNLIKELY(nptls > MAXPTL)) {
182 return csound->InitError(csound, Str("illegal no of partials"));
183 }
184 if (UNLIKELY(*p->irolloff<=FL(0.0))) p->rolloff = FL(0.6);
185 else p->rolloff = *p->irolloff;
186 p->nptls = nptls; /* number, whether all or odd */
187 ptlmax = nptls;
188 dstp = p->pdist;
189 fnfreqs = (MYFLT)specp->nfreqs;
190 for (nn = 1; nn <= ptlmax; nn++)
191 *dstp++ = (int32_t) ((LOG((MYFLT) nn) / (MYFLT)LOGTWO) * fnfreqs + FL(0.5));
192 if (UNLIKELY((rolloff = p->rolloff) == FL(0.0) ||
193 rolloff == FL(1.0) || nptls == 1)) {
194 p->rolloff = FL(0.0);
195 weightsum = (MYFLT)nptls;
196 }
197 else {
198 MYFLT *fltp = p->pmult;
199 MYFLT octdrop = (FL(1.0) - rolloff) / fnfreqs;
200 weightsum = FL(0.0);
201 for (dstp = p->pdist, nn = nptls; nn--; ) {
202 weight = FL(1.0) - octdrop * *dstp++; /* rolloff * octdistance */
203 weightsum += weight;
204 *fltp++ = weight;
205 }
206 if (UNLIKELY(*--fltp < FL(0.0))) {
207 return csound->InitError(csound, Str("per octave rolloff too steep"));
208 }
209 p->rolloff = 1;
210 }
211 lobin = (int32)(specp->downsrcp->looct * fnfreqs);
212 oct0p = p->fundp - lobin; /* virtual loc of oct 0 */
213
214 flop = oct0p + (int32_t)(*p->ilo * fnfreqs);
215 fhip = oct0p + (int32_t)(*p->ihi * fnfreqs);
216 fundp = p->fundp;
217 fendp = fundp + specp->npts;
218 if (flop < fundp) flop = fundp;
219 if (UNLIKELY(fhip > fendp)) fhip = fendp;
220 if (UNLIKELY(flop >= fhip)) { /* chk hi-lo range valid */
221 return csound->InitError(csound, Str("illegal lo-hi values"));
222 }
223 for (fp = fundp; fp < flop; )
224 *fp++ = FL(0.0); /* clear unused lo and hi range */
225 for (fp = fhip; fp < fendp; )
226 *fp++ = FL(0.0);
227
228 dbthresh = *p->idbthresh; /* thresholds: */
229 ampthresh = (MYFLT)exp((double)dbthresh * LOG10D20);
230 p->threshon = ampthresh; /* mag */
231 p->threshoff = ampthresh * FL(0.5);
232 p->threshon *= weightsum;
233 p->threshoff *= weightsum;
234 p->oct0p = oct0p; /* virtual loc of oct 0 */
235 p->confact = *p->iconf;
236 p->flop = flop;
237 p->fhip = fhip;
238 p->playing = 0;
239 p->kvalsav = (*p->istrt>=FL(0.0) ? *p->istrt : (*p->ilo+*p->ihi)*FL(0.5));
240 p->kval = p->kinc = FL(0.0);
241 p->kavl = p->kanc = FL(0.0);
242 p->jmpcount = 0;
243 return OK;
244 }
245
pitch(CSOUND * csound,PITCH * p)246 int32_t pitch(CSOUND *csound, PITCH *p)
247 {
248 MYFLT *asig;
249 double q;
250 double c1 = p->c1, c2 = p->c2;
251
252 MYFLT a, b, *dftp, SIG, yt1, yt2;
253 uint32_t offset = p->h.insdshead->ksmps_offset;
254 uint32_t early = p->h.insdshead->ksmps_no_end;
255 uint32_t n, nsmps = CS_KSMPS;
256 int32_t nocts, winlen;
257 DOWNDAT *downp = &p->downsig;
258 OCTDAT *octp;
259 SPECDAT *specp;
260 MYFLT c;
261 MYFLT kvar;
262 /* RMS */
263 q = p->prvq;
264 asig = p->asig;
265 if (UNLIKELY(offset)) memset(asig, '\0', offset*sizeof(MYFLT));
266 if (UNLIKELY(early)) {
267 nsmps -= early;
268 memset(&asig[nsmps], '\0', early*sizeof(MYFLT));
269 }
270 for (n=offset; n<nsmps; n++) {
271 MYFLT as = asig[n]*DFLT_DBFS/csound->e0dbfs; /* Normalise.... */
272 q = c1 * as * as + c2 * q;
273 SIG = as; /* for each source sample: */
274 octp = downp->octdata; /* align onto top octave */
275 nocts = downp->nocts;
276 do { /* then for each oct: */
277 const MYFLT *coefp;
278 MYFLT *ytp, *curp;
279 int32_t nfilt;
280 curp = octp->curp;
281 *curp++ = SIG; /* write samp to cur buf */
282 if (UNLIKELY(curp >= octp->endp))
283 curp = octp->begp; /* & modulo the pointer */
284 octp->curp = curp;
285 if (UNLIKELY(!(--nocts))) break; /* if lastoct, break */
286 coefp = bicoefs; ytp = octp->feedback;
287 for (nfilt = 3; nfilt--; ) { /* apply triple biquad: */
288 yt2 = *ytp++; yt1 = *ytp--; /* get prev feedback */
289 SIG -= (*coefp++ * yt1); /* apply recurs filt */
290 SIG -= (*coefp++ * yt2);
291 *ytp++ = yt1; *ytp++ = SIG; /* stor nxt feedback */
292 SIG *= *coefp++;
293 SIG += (*coefp++ * yt1); /* apply forwrd filt */
294 SIG += (*coefp++ * yt2);
295 }
296 } while (!(++octp->scount & 01) && octp++); /* send alt samps to nxtoct */
297 }
298 p->prvq = q;
299 kvar = SQRT((MYFLT)q); /* End of spectrum part */
300
301 specp = &p->wsig;
302 if (LIKELY((--p->scountdown))) goto nxt; /* if not yet time for new spec */
303 p->scountdown = p->timcount; /* else reset counter & proceed: */
304 downp = &p->downsig;
305 nocts = downp->nocts;
306 octp = downp->octdata + nocts;
307 dftp = (MYFLT *) specp->auxch.auxp;
308 winlen = *(p->winlen);
309 while (nocts--) {
310 MYFLT *bufp, *sinp, *cosp;
311 int32_t len, *lenp, *offp, nfreqs;
312 MYFLT *begp, *curp, *endp, *linbufp;
313 int32_t len2;
314 octp--; /* for each oct (low to high) */
315 begp = octp->begp;
316 curp = octp->curp;
317 endp = octp->endp;
318 if (UNLIKELY((len = endp - curp) >= winlen)) /* if no wrap */
319 linbufp = curp; /* use samples in circbuf */
320 else {
321 len2 = winlen - len;
322 linbufp = bufp = p->linbufp; /* else cp crcbuf to linbuf */
323 while (len--)
324 *bufp++ = *curp++;
325 curp = begp;
326 while (len2--)
327 *bufp++ = *curp++;
328 }
329 cosp = p->cosp; /* get start windowed sines */
330 sinp = p->sinp;
331 lenp = p->winlen;
332 offp = p->offset;
333 for (nfreqs=p->nfreqs; nfreqs--; ) { /* now for ea. frq this oct */
334 a = FL(0.0);
335 b = FL(0.0);
336 bufp = linbufp + *offp++;
337 for (len = *lenp++; len--; bufp++) { /* apply windowed sine seg */
338 a += *bufp * *cosp++;
339 b += *bufp * *sinp++;
340 }
341 c = HYPOT(a, b); /* get magnitude */
342 *dftp++ = c; /* store in out spectrum */
343 }
344 }
345 specp->ktimstamp = CS_KCNT; /* time-stamp the output */
346
347 nxt:
348 /* specptrk */
349 {
350 MYFLT *inp = (MYFLT *) specp->auxch.auxp;
351 MYFLT *endp = inp + specp->npts;
352 MYFLT *inp2, sum, *fp;
353 int32_t nn, *pdist, confirms;
354 MYFLT kval, fmax, *fmaxp, absdiff, realbin;
355 MYFLT *flop, *fhip, *ilop, *ihip, a, b, c, denom, delta;
356 int32 lobin, hibin;
357
358 if (UNLIKELY(inp==NULL)) goto err1; /* RWD fix */
359 kval = p->playing == PLAYING ? p->kval : p->kvalsav;
360 lobin = (int32)((kval - kvar) * specp->nfreqs);/* set lims of frq interest */
361 hibin = (int32)((kval + kvar) * specp->nfreqs);
362 if ((flop = p->oct0p + lobin) < p->flop) /* as fundp bin pntrs */
363 flop = p->flop;
364 if ((fhip = p->oct0p + hibin) > p->fhip) /* within hard limits */
365 fhip = p->fhip;
366 ilop = inp + (flop - p->fundp); /* similar for input bins */
367 ihip = inp + (fhip - p->fundp);
368 inp = ilop;
369 fp = flop;
370 if (p->rolloff) {
371 MYFLT *pmult;
372 do {
373 sum = *inp;
374 pdist = p->pdist + 1;
375 pmult = p->pmult + 1;
376 for (nn = p->nptls; --nn; ) {
377 if ((inp2 = inp + *pdist++) >= endp)
378 break;
379 sum += *inp2 * *pmult++;
380 }
381 *fp++ = sum;
382 } while (++inp < ihip);
383 }
384 else {
385 do {
386 sum = *inp;
387 pdist = p->pdist + 1;
388 for (nn = p->nptls; --nn; ) {
389 if ((inp2 = inp + *pdist++) >= endp)
390 break;
391 sum += *inp2;
392 }
393 *fp++ = sum;
394 } while (++inp < ihip);
395 }
396 fp = flop; /* now srch fbins for peak */
397 for (fmaxp = fp, fmax = *fp; ++fp<fhip; )
398 if (UNLIKELY(*fp > fmax)) {
399 fmax = *fp;
400 fmaxp = fp;
401 }
402 if (!p->playing) {
403 if (fmax > p->threshon) /* not playing & threshon? */
404 p->playing = STARTING; /* prepare to turn on */
405 else goto output;
406 }
407 else {
408 if (fmax < p->threshoff) { /* playing & threshoff ? */
409 if (p->playing == PLAYING)
410 p->kvalsav = p->kval; /* save val & turn off */
411 p->kval = FL(0.0);
412 p->kavl = FL(0.0);
413 p->kinc = FL(0.0);
414 p->kanc = FL(0.0);
415 p->playing = 0;
416 goto output;
417 }
418 }
419 a = fmaxp>flop ? *(fmaxp-1) : FL(0.0); /* calc a refined bin no */
420 b = fmax;
421 c = fmaxp<fhip-1 ? *(fmaxp+1) : FL(0.0);
422 if (b < FL(2.0) * (a + c))
423 denom = b + b - a - c;
424 else denom = a + b + c;
425 if (denom != FL(0.0))
426 delta = FL(0.5) * (c - a) / denom;
427 else delta = FL(0.0);
428 realbin = (fmaxp - p->oct0p) + delta; /* get modified bin number */
429 kval = realbin / specp->nfreqs; /* & cvt to true decoct */
430
431 if (p->playing == STARTING) { /* STARTING mode: */
432 absdiff = FABS(kval - p->kvalsav);// < FL(0.0)) absdiff = -absdiff;
433 confirms = (int32_t)(absdiff * p->confact); /* get interval dependency */
434 if (UNLIKELY(p->jmpcount < confirms)) {
435 p->jmpcount += 1; /* if not enough confirms, */
436 goto output; /* must wait some more */
437 } else {
438 p->playing = PLAYING; /* else switch on playing */
439 p->jmpcount = 0;
440 p->kval = kval; /* but suppress interp */
441 p->kinc = FL(0.0);
442 }
443 } else { /* PLAYING mode: */
444 absdiff = FABS(kval - p->kval);
445 confirms = (int32_t)(absdiff * p->confact); /* get interval dependency */
446 if (p->jmpcount < confirms) {
447 p->jmpcount += 1; /* if not enough confirms, */
448 p->kinc = FL(0.0); /* must wait some more */
449 } else {
450 p->jmpcount = 0; /* else OK to jump interval */
451 p->kval = kval;
452 }
453 }
454 fmax += delta * (c - a) * FL(0.25); /* get modified amp */
455 p->kavl = fmax;
456 }
457 output:
458 *p->koct = p->kval; /* output true decoct & amp */
459 *p->kamp = p->kavl * FL(4.0);
460 return OK;
461 err1:
462 return csound->PerfError(csound, &(p->h),
463 Str("pitch: not initialised"));
464 }
465
466 /* Multiply and accumulate opcodes */
467
macset(CSOUND * csound,SUM * p)468 int32_t macset(CSOUND *csound, SUM *p)
469 {
470 if (UNLIKELY((((int32_t)p->INOCOUNT)&1)==1)) {
471 return csound->PerfError(csound, &(p->h),
472 Str("Must have even number of arguments in mac\n"));
473 }
474 return OK;
475 }
476
maca(CSOUND * csound,SUM * p)477 int32_t maca(CSOUND *csound, SUM *p)
478 {
479 IGN(csound);
480 uint32_t offset = p->h.insdshead->ksmps_offset;
481 uint32_t early = p->h.insdshead->ksmps_no_end;
482 uint32_t k, nsmps = CS_KSMPS;
483 int32_t count=(int32_t) p->INOCOUNT, j;
484 MYFLT *ar = p->ar, **args = p->argums;
485
486 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
487 if (UNLIKELY(early)) {
488 nsmps -= early;
489 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
490 }
491 for (k=offset; k<nsmps; k++) {
492 MYFLT ans = FL(0.0);
493 for (j=0; j<count; j +=2)
494 ans += args[j][k] * args[j+1][k];
495 ar[k] = ans;
496 }
497 return OK;
498 }
499
mac(CSOUND * csound,SUM * p)500 int32_t mac(CSOUND *csound, SUM *p)
501 {
502 IGN(csound);
503 uint32_t offset = p->h.insdshead->ksmps_offset;
504 uint32_t early = p->h.insdshead->ksmps_no_end;
505 uint32_t k, nsmps = CS_KSMPS;
506 int32_t count=(int32_t) p->INOCOUNT, j;
507 MYFLT *ar = p->ar, **args = p->argums;
508 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
509 if (UNLIKELY(early)) {
510 nsmps -= early;
511 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
512 }
513 for (k=offset; k<nsmps; k++) {
514 MYFLT ans = FL(0.0);
515 for (j=0; j<count; j +=2)
516 ans += *args[j] * args[j+1][k];
517 ar[k] = ans;
518 }
519 return OK;
520 }
521
522 typedef struct {
523 RTCLOCK r;
524 double counters[33];
525 int32_t running[33];
526 } CPU_CLOCK;
527
initClockStruct(CSOUND * csound,void ** p)528 static void initClockStruct(CSOUND *csound, void **p)
529 {
530 *p = csound->QueryGlobalVariable(csound, "readClock::counters");
531 if (UNLIKELY(*p == NULL)) {
532 csound->CreateGlobalVariable(csound, "readClock::counters",
533 sizeof(CPU_CLOCK));
534 *p = csound->QueryGlobalVariable(csound, "readClock::counters");
535 csound->InitTimerStruct(&(((CPU_CLOCK*) (*p))->r));
536 }
537 }
538
getClockStruct(CSOUND * csound,void ** p)539 static inline CPU_CLOCK *getClockStruct(CSOUND *csound, void **p)
540 {
541 if (UNLIKELY(*p == NULL))
542 initClockStruct(csound, p);
543 return (CPU_CLOCK*) (*p);
544 }
545
clockset(CSOUND * csound,CLOCK * p)546 int32_t clockset(CSOUND *csound, CLOCK *p)
547 {
548 IGN(csound);
549 p->c = (int32_t)*p->cnt;
550 if (UNLIKELY(p->c < 0 || p->c > 31))
551 p->c = 32;
552 return OK;
553 }
554
clockon(CSOUND * csound,CLOCK * p)555 int32_t clockon(CSOUND *csound, CLOCK *p)
556 {
557 CPU_CLOCK *clk = getClockStruct(csound, &(p->clk));
558 if (LIKELY(!clk->running[p->c])) {
559 clk->running[p->c] = 1;
560 clk->counters[p->c] -= csound->GetCPUTime(&(clk->r));
561 }
562 return OK;
563 }
564
clockoff(CSOUND * csound,CLOCK * p)565 int32_t clockoff(CSOUND *csound, CLOCK *p)
566 {
567 CPU_CLOCK *clk = getClockStruct(csound, &(p->clk));
568 if (LIKELY(clk->running[p->c])) {
569 clk->running[p->c] = 0;
570 clk->counters[p->c] += csound->GetCPUTime(&(clk->r));
571 }
572 return OK;
573 }
574
clockread(CSOUND * csound,CLKRD * p)575 int32_t clockread(CSOUND *csound, CLKRD *p)
576 {
577 CPU_CLOCK *clk = getClockStruct(csound, &(p->clk));
578 int32_t cnt = (int32_t) *p->a;
579 if (UNLIKELY(cnt < 0 || cnt > 32)) cnt = 32;
580 if (UNLIKELY(clk->running[cnt]))
581 return csound->InitError(csound, Str("clockread: clock still running, "
582 "call clockoff first"));
583 /* result in ms */
584 printf("readclock%d: %g\n", cnt, clk->counters[cnt]);
585
586 *p->r = (MYFLT) (clk->counters[cnt] * 1000.0);
587 return OK;
588 }
589
scratchread(CSOUND * csound,SCRATCHPAD * p)590 int32_t scratchread(CSOUND *csound, SCRATCHPAD *p)
591 {
592 int32_t index = MYFLT2LRND(*p->index);
593 if (index<0 || index>3)
594 return csound->PerfError(csound, &(p->h),
595 Str("scratchpad index out of range"));
596 *p->val = p->h.insdshead->scratchpad[index];
597 return OK;
598 }
599
scratchwrite(CSOUND * csound,SCRATCHPAD * p)600 int32_t scratchwrite(CSOUND *csound, SCRATCHPAD *p)
601 {
602 int32_t index = MYFLT2LRND(*p->index);
603 if (index<0 || index>3)
604 return csound->PerfError(csound, &(p->h),
605 Str("scratchpad index out of range"));
606 p->h.insdshead->scratchpad[index] = *p->val;
607 return OK;
608 }
609
610
611 /* ************************************************************ */
612 /* Opcodes from Peter Neubäcker */
613 /* ************************************************************ */
614
adsyntset(CSOUND * csound,ADSYNT * p)615 int32_t adsyntset(CSOUND *csound, ADSYNT *p)
616 {
617 FUNC *ftp;
618 uint32_t count;
619 int32 *lphs;
620
621 p->inerr = 0;
622
623 if (LIKELY((ftp = csound->FTFind(csound, p->ifn)) != NULL)) {
624 p->ftp = ftp;
625 }
626 else {
627 p->inerr = 1;
628 return csound->InitError(csound, Str("adsynt: wavetable not found!"));
629 }
630
631 count = (uint32_t)*p->icnt;
632 if (UNLIKELY(count < 1))
633 count = 1;
634 p->count = count;
635
636 if (LIKELY((ftp = csound->FTnp2Find(csound, p->ifreqtbl)) != NULL)) {
637 p->freqtp = ftp;
638 }
639 else {
640 p->inerr = 1;
641 return csound->InitError(csound, Str("adsynt: freqtable not found!"));
642 }
643 if (UNLIKELY(ftp->flen < count)) {
644 p->inerr = 1;
645 return csound->InitError(csound, Str(
646 "adsynt: partial count is greater than freqtable size!"));
647 }
648
649 if (LIKELY((ftp = csound->FTnp2Find(csound, p->iamptbl)) != NULL)) {
650 p->amptp = ftp;
651 }
652 else {
653 p->inerr = 1;
654 return csound->InitError(csound, Str("adsynt: amptable not found!"));
655 }
656 if (UNLIKELY(ftp->flen < count)) {
657 p->inerr = 1;
658 return csound->InitError(csound, Str(
659 "adsynt: partial count is greater than amptable size!"));
660 }
661
662 if (p->lphs.auxp==NULL || p->lphs.size < (size_t)sizeof(int32)*count)
663 csound->AuxAlloc(csound, sizeof(int32)*count, &p->lphs);
664
665 lphs = (int32*)p->lphs.auxp;
666 if (*p->iphs > 1) {
667 do {
668 *lphs++ = ((int32) ((MYFLT) ((double) rand_31(csound) / 2147483645.0)
669 * FMAXLEN)) & PHMASK;
670 } while (--count);
671 }
672 else if (*p->iphs >= 0) {
673 do {
674 *lphs++ = ((int32) (*p->iphs * FMAXLEN)) & PHMASK;
675 } while (--count);
676 }
677
678 return OK;
679 }
680
adsynt(CSOUND * csound,ADSYNT * p)681 int32_t adsynt(CSOUND *csound, ADSYNT *p)
682 {
683 FUNC *ftp, *freqtp, *amptp;
684 MYFLT *ar, *ftbl, *freqtbl, *amptbl;
685 MYFLT amp0, amp, cps0, cps;
686 int32 phs, inc, lobits;
687 int32 *lphs;
688 uint32_t offset = p->h.insdshead->ksmps_offset;
689 uint32_t early = p->h.insdshead->ksmps_no_end;
690 uint32_t n, nsmps = CS_KSMPS;
691 int32_t c, count;
692
693 if (UNLIKELY(p->inerr)) {
694 return csound->PerfError(csound, &(p->h),
695 Str("adsynt: not initialised"));
696 }
697 ftp = p->ftp;
698 ftbl = ftp->ftable;
699 lobits = ftp->lobits;
700 freqtp = p->freqtp;
701 freqtbl = freqtp->ftable;
702 amptp = p->amptp;
703 amptbl = amptp->ftable;
704 lphs = (int32*)p->lphs.auxp;
705
706 cps0 = *p->kcps;
707 amp0 = *p->kamp;
708 count = p->count;
709
710 ar = p->sr;
711 memset(ar, 0, nsmps*sizeof(MYFLT));
712 if (UNLIKELY(early)) nsmps -= early;
713
714 for (c=0; c<count; c++) {
715 amp = amptbl[c] * amp0;
716 cps = freqtbl[c] * cps0;
717 inc = (int32) (cps * csound->sicvt);
718 phs = lphs[c];
719 for (n=offset; n<nsmps; n++) {
720 ar[n] += *(ftbl + (phs >> lobits)) * amp;
721 phs += inc;
722 phs &= PHMASK;
723 }
724 lphs[c] = phs;
725 }
726 return OK;
727 }
728
hsboscset(CSOUND * csound,HSBOSC * p)729 int32_t hsboscset(CSOUND *csound, HSBOSC *p)
730 {
731 FUNC *ftp;
732 int32_t octcnt, i;
733
734 if (LIKELY((ftp = csound->FTnp2Finde(csound, p->ifn)) != NULL)) {
735 p->ftp = ftp;
736 if (UNLIKELY(*p->ioctcnt < 2))
737 octcnt = 3;
738 else
739 octcnt = (int32_t)*p->ioctcnt;
740 if (UNLIKELY(octcnt > 10))
741 octcnt = 10;
742 p->octcnt = octcnt;
743 if (*p->iphs >= 0) {
744 for (i=0; i<octcnt; i++)
745 p->lphs[i] = ((int32)(*p->iphs * FMAXLEN)) & PHMASK;
746 }
747 }
748 else p->ftp = NULL;
749 if (LIKELY((ftp = csound->FTnp2Finde(csound, p->imixtbl)) != NULL)) {
750 p->mixtp = ftp;
751 }
752 else p->mixtp = NULL;
753 return OK;
754 }
755
hsboscil(CSOUND * csound,HSBOSC * p)756 int32_t hsboscil(CSOUND *csound, HSBOSC *p)
757 {
758 FUNC *ftp, *mixtp;
759 MYFLT fract, v1, amp0, amp, *ar, *ftab, *mtab;
760 int32 phs, inc, lobits;
761 int32 phases[10];
762 uint32_t offset = p->h.insdshead->ksmps_offset;
763 uint32_t early = p->h.insdshead->ksmps_no_end;
764 uint32_t n, nsmps = CS_KSMPS;
765 MYFLT tonal, bright, freq, ampscl;
766 int32_t octcnt = p->octcnt;
767 MYFLT octstart, octoffs, octbase;
768 int32_t octshift, i, mtablen;
769 MYFLT hesr = CS_ESR / FL(2.0);
770
771 ftp = p->ftp;
772 mixtp = p->mixtp;
773 if (UNLIKELY(ftp==NULL || mixtp==NULL)) {
774 return csound->PerfError(csound, &(p->h),
775 Str("hsboscil: not initialised"));
776 }
777
778 tonal = *p->ktona;
779 tonal -= MYFLOOR(tonal);
780 bright = *p->kbrite - tonal;
781 octstart = bright - (MYFLT)octcnt * FL(0.5);
782 octbase = MYFLOOR(MYFLOOR(octstart) + FL(1.5));
783 octoffs = octbase - octstart;
784
785 mtab = mixtp->ftable;
786 mtablen = mixtp->flen;
787 freq = *p->ibasef * POWER(FL(2.0), tonal + octbase);
788
789 ampscl = mtab[(int32_t)((1.0 / (MYFLT)octcnt) * mtablen)];
790 amp = mtab[(int32_t)((octoffs / (MYFLT)octcnt) * mtablen)];
791 if ((amp - p->prevamp) > (ampscl * FL(0.5)))
792 octshift = 1;
793 else if ((amp - p->prevamp) < (-(ampscl * FL(0.5))))
794 octshift = -1;
795 else
796 octshift = 0;
797 p->prevamp = amp;
798
799 ampscl = FL(0.0);
800 for (i=0; i<octcnt; i++) {
801 phases[i] = p->lphs[(i+octshift+100*octcnt) % octcnt];
802 ampscl += mtab[(int32_t)(((MYFLT)i / (MYFLT)octcnt) * mtablen)];
803 }
804
805 amp0 = *p->kamp / ampscl;
806 lobits = ftp->lobits;
807 ar = p->sr;
808 memset(ar, 0, nsmps*sizeof(MYFLT));
809 if (UNLIKELY(early)) nsmps -= early;
810
811 for (i=0; i<octcnt; i++) {
812 phs = phases[i];
813 amp = mtab[(int32_t)((octoffs / (MYFLT)octcnt) * mtablen)] * amp0;
814 if (UNLIKELY(freq > hesr))
815 amp = FL(0.0);
816 inc = (int32)(freq * csound->sicvt);
817 for (n=offset; n<nsmps; n++) {
818 fract = PFRAC(phs);
819 ftab = ftp->ftable + (phs >> lobits);
820 v1 = *ftab++;
821 ar[n] += (v1 + (*ftab - v1) * fract) * amp;
822 phs += inc;
823 phs &= PHMASK;
824 }
825 p->lphs[i] = phs;
826
827 octoffs += FL(1.0);
828 freq *= FL(2.0);
829 }
830 return OK;
831 }
832
pitchamdfset(CSOUND * csound,PITCHAMDF * p)833 int32_t pitchamdfset(CSOUND *csound, PITCHAMDF *p)
834 {
835 MYFLT srate, downs;
836 int32 size, minperi, maxperi, downsamp, upsamp, msize, bufsize;
837 uint32_t interval;
838 uint32_t nsmps = CS_KSMPS;
839
840 p->inerr = 0;
841
842 downs = *p->idowns;
843 if (downs < (-FL(1.9))) {
844 upsamp = (int32_t)MYFLT2LONG((-downs));
845 downsamp = 0;
846 srate = CS_ESR * (MYFLT)upsamp;
847 }
848 else {
849 downsamp = (int32_t)MYFLT2LONG(downs);
850 if (UNLIKELY(downsamp < 1))
851 downsamp = 1;
852 srate = CS_ESR / (MYFLT)downsamp;
853 upsamp = 0;
854 }
855
856 minperi = (int32)(srate / *p->imaxcps);
857 maxperi = (int32)(FL(0.5)+srate / *p->imincps);
858 if (UNLIKELY(maxperi <= minperi)) {
859 p->inerr = 1;
860 return csound->InitError(csound,
861 Str("pitchamdf: maxcps must be > mincps !"));
862 }
863
864 if (*p->iexcps < 1)
865 interval = maxperi;
866 else
867 interval = (uint32_t)(srate / *p->iexcps);
868 if (interval < nsmps) {
869 if (downsamp)
870 interval = nsmps / downsamp;
871 else
872 interval = nsmps * upsamp;
873 }
874
875 size = maxperi + interval;
876 bufsize = sizeof(MYFLT)*(size + maxperi + 2);
877
878 p->srate = srate;
879 p->downsamp = downsamp;
880 p->upsamp = upsamp;
881 p->minperi = minperi;
882 p->maxperi = maxperi;
883 p->size = size;
884 p->readp = 0;
885 p->index = 0;
886 p->lastval = FL(0.0);
887 if (*p->icps < 1)
888 p->peri = (minperi + maxperi) / 2;
889 else
890 p->peri = (int32_t)(srate / *p->icps);
891
892 if (*p->irmsmedi < 1)
893 p->rmsmedisize = 0;
894 else
895 p->rmsmedisize = ((int32_t)MYFLT2LONG(*p->irmsmedi))*2+1;
896 p->rmsmediptr = 0;
897
898 if (p->rmsmedisize) {
899 msize = p->rmsmedisize * 3 * sizeof(MYFLT);
900 if (p->rmsmedian.auxp==NULL || p->rmsmedian.size < (size_t)msize)
901 csound->AuxAlloc(csound, msize, &p->rmsmedian);
902 else {
903 memset(p->rmsmedian.auxp, 0, msize);
904 }
905 }
906
907 if (*p->imedi < 1)
908 p->medisize = 0;
909 else
910 p->medisize = (int32_t)MYFLT2LONG(*p->imedi)*2+1;
911 p->mediptr = 0;
912
913 if (p->medisize) {
914 msize = p->medisize * 3 * sizeof(MYFLT);
915 if (p->median.auxp==NULL || p->median.size < (size_t)msize)
916 csound->AuxAlloc(csound, (size_t)msize, &p->median);
917 else {
918 memset(p->median.auxp, 0, msize);
919 }
920 }
921
922 if (p->buffer.auxp==NULL || p->buffer.size < (size_t)bufsize) {
923 csound->AuxAlloc(csound, bufsize, &p->buffer);
924 }
925 else
926 memset(p->buffer.auxp, 0, bufsize);
927 return OK;
928 }
929
930 #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp
931
medianvalue(uint32 n,MYFLT * vals)932 MYFLT medianvalue(uint32 n, MYFLT *vals)
933 { /* vals must point to 1 below relevant data! */
934 uint32 i, ir, j, l, mid;
935 uint32 k = (n + 1) / 2;
936 MYFLT a, temp;
937
938 l = 1;
939 ir = n;
940 while (1) {
941 if (ir <= l+1) {
942 if (ir == l+1 && vals[ir] < vals[l]) {
943 SWAP(vals[l], vals[ir]);
944 }
945 return vals[k];
946 }
947 else {
948 mid = (l+ir) >> 1;
949 SWAP(vals[mid], vals[l+1]);
950 if (vals[l+1] > vals[ir]) {
951 SWAP(vals[l+1], vals[ir]);
952 }
953 if (vals[l] > vals[ir]) {
954 SWAP(vals[l], vals[ir]);
955 }
956 if (vals[l+1] > vals[l]) {
957 SWAP(vals[l+1], vals[l]);
958 }
959 i = l + 1;
960 j = ir;
961 a = vals[l];
962 while (1) {
963 do i++; while (vals[i] < a);
964 do j--; while (vals[j] > a);
965 if (UNLIKELY(j < i)) break;
966 SWAP(vals[i], vals[j]);
967 }
968 vals[l] = vals[j];
969 vals[j] = a;
970 if (j >= k) ir = j-1;
971 if (j <= k) l = i;
972 }
973 }
974 }
975 #undef SWAP
976
pitchamdf(CSOUND * csound,PITCHAMDF * p)977 int32_t pitchamdf(CSOUND *csound, PITCHAMDF *p)
978 {
979 MYFLT *buffer = (MYFLT*)p->buffer.auxp;
980 MYFLT *rmsmedian = (MYFLT*)p->rmsmedian.auxp;
981 int32 rmsmedisize = p->rmsmedisize;
982 int32 rmsmediptr = p->rmsmediptr;
983 MYFLT *median = (MYFLT*)p->median.auxp;
984 int32 medisize = p->medisize;
985 int32 mediptr = p->mediptr;
986 int32 size = p->size;
987 int32 index = p->index;
988 int32 minperi = p->minperi;
989 int32 maxperi = p->maxperi;
990 MYFLT *asig = p->asig;
991 MYFLT srate = p->srate;
992 int32 peri = p->peri;
993 int32 upsamp = p->upsamp;
994 MYFLT upsmp = (MYFLT)upsamp;
995 MYFLT lastval = p->lastval;
996 MYFLT newval, delta;
997 int32 readp = p->readp;
998 int32 interval = size - maxperi;
999 int32_t nsmps = CS_KSMPS;
1000 int32_t i;
1001 int32 i1, i2;
1002 MYFLT val, rms;
1003 double sum;
1004 MYFLT acc, accmin, diff;
1005
1006 if (UNLIKELY(p->inerr)) {
1007 return csound->PerfError(csound, &(p->h),
1008 Str("pitchamdf: not initialised"));
1009 }
1010
1011 if (upsamp) {
1012 while (1) {
1013 newval = asig[readp++];
1014 delta = (newval-lastval) / upsmp;
1015 lastval = newval;
1016
1017 for (i=0; i<upsamp; i++) {
1018 newval += delta;
1019 buffer[index++] = newval;
1020
1021 if (index == size) {
1022 peri = minperi;
1023 accmin = FL(0.0);
1024 for (i2 = 0; i2 < size; ++i2) {
1025 diff = buffer[i2+minperi] - buffer[i2];
1026 if (diff > 0) accmin += diff;
1027 else accmin -= diff;
1028 }
1029 for (i1 = minperi + 1; i1 <= maxperi; ++i1) {
1030 acc = FL(0.0);
1031 for (i2 = 0; i2 < size; ++i2) {
1032 diff = buffer[i1+i2] - buffer[i2];
1033 if (diff > 0) acc += diff;
1034 else acc -= diff;
1035 }
1036 if (acc < accmin) {
1037 accmin = acc;
1038 peri = i1;
1039 }
1040 }
1041
1042 for (i1 = 0; i1 < interval; i1++)
1043 buffer[i1] = buffer[i1+interval];
1044 index = maxperi;
1045
1046 if (medisize) {
1047 median[mediptr] = (MYFLT)peri;
1048 for (i1 = 0; i1 < medisize; i1++)
1049 median[medisize+i1] = median[i1];
1050
1051 median[medisize*2+mediptr] =
1052 medianvalue(medisize, &median[medisize-1]);
1053 peri = (int32)median[medisize*2 +
1054 ((mediptr+medisize/2+1) % medisize)];
1055
1056 mediptr = (mediptr + 1) % medisize;
1057 p->mediptr = mediptr;
1058 }
1059 }
1060 }
1061 if (readp >= nsmps) break;
1062 }
1063 readp = readp % nsmps;
1064 p->lastval = lastval;
1065 }
1066 else {
1067 int32 downsamp = p->downsamp;
1068 while (1) {
1069 buffer[index++] = asig[readp];
1070 readp += downsamp;
1071
1072 if (index == size) {
1073 peri = minperi;
1074 accmin = FL(0.0);
1075 for (i2 = 0; i2 < size; ++i2) {
1076 diff = buffer[i2+minperi] - buffer[i2];
1077 if (diff > FL(0.0)) accmin += diff;
1078 else accmin -= diff;
1079 }
1080 for (i1 = minperi + 1; i1 <= maxperi; ++i1) {
1081 acc = FL(0.0);
1082 for (i2 = 0; i2 < size; ++i2) {
1083 diff = buffer[i1+i2] - buffer[i2];
1084 if (diff > FL(0.0)) acc += diff;
1085 else acc -= diff;
1086 }
1087 if (acc < accmin) {
1088 accmin = acc;
1089 peri = i1;
1090 }
1091 }
1092
1093 for (i1 = 0; i1 < interval; i1++)
1094 buffer[i1] = buffer[i1+interval];
1095 index = maxperi;
1096
1097 if (medisize) {
1098 median[mediptr] = (MYFLT)peri;
1099 for (i1 = 0; i1 < medisize; i1++)
1100 median[medisize+i1] = median[i1];
1101
1102 median[medisize*2+mediptr] =
1103 medianvalue(medisize, &median[medisize-1]);
1104 peri = (int32)median[medisize*2 +
1105 ((mediptr+medisize/2+1) % medisize)];
1106
1107 mediptr = (mediptr + 1) % medisize;
1108 p->mediptr = mediptr;
1109 }
1110 }
1111
1112 if (readp >= nsmps) break;
1113 }
1114 readp = readp % nsmps;
1115 }
1116 buffer = &buffer[(index + size - peri) % size];
1117 sum = 0.0;
1118 for (i1=0; i1<peri; i1++) {
1119 val = buffer[i1];
1120 sum += (double)(val * val);
1121 }
1122 if (UNLIKELY(peri==0)) /* How xould thus happen??? */
1123 rms = FL(0.0);
1124 else
1125 rms = (MYFLT)sqrt(sum / (double)peri);
1126 if (rmsmedisize) {
1127 rmsmedian[rmsmediptr] = rms;
1128 for (i1 = 0; i1 < rmsmedisize; i1++)
1129 rmsmedian[rmsmedisize+i1] = rmsmedian[i1];
1130
1131 rmsmedian[rmsmedisize*2+rmsmediptr] =
1132 medianvalue(rmsmedisize, &rmsmedian[rmsmedisize-1]);
1133 rms = rmsmedian[rmsmedisize*2 +
1134 ((rmsmediptr+rmsmedisize/2+1) % rmsmedisize)];
1135
1136 rmsmediptr = (rmsmediptr + 1) % rmsmedisize;
1137 p->rmsmediptr = rmsmediptr;
1138 }
1139
1140 if (UNLIKELY(peri==0)) {
1141 *p->kcps = FL(0.0);
1142 }
1143 else
1144 *p->kcps = srate / (MYFLT)peri;
1145 *p->krms = rms;
1146 p->index = index;
1147 p->peri = peri;
1148 p->readp = readp;
1149 return OK;
1150 }
1151
1152 /*==================================================================*/
1153 /* phasorbnk */
1154 /*==================================================================*/
1155
phsbnkset(CSOUND * csound,PHSORBNK * p)1156 int32_t phsbnkset(CSOUND *csound, PHSORBNK *p)
1157 {
1158 double phs;
1159 int32_t n, count;
1160 double *curphs;
1161
1162 count = (int32_t)MYFLT2LONG(*p->icnt);
1163 if (UNLIKELY(count < 2))
1164 count = 2;
1165
1166 if (p->curphs.auxp==NULL || p->curphs.size < (size_t)sizeof(double)*count)
1167 csound->AuxAlloc(csound, (size_t)sizeof(double)*count, &p->curphs);
1168
1169 curphs = (double*)p->curphs.auxp;
1170 if (*p->iphs > 1) {
1171 for (n=0; n<count;n++)
1172 curphs[n] = (double) rand_31(csound) / 2147483645.0;
1173 }
1174 else if ((phs = *p->iphs) >= 0) {
1175 for (n=0; n<count;n++) curphs[n] = phs;
1176 }
1177 return OK;
1178 }
1179
kphsorbnk(CSOUND * csound,PHSORBNK * p)1180 int32_t kphsorbnk(CSOUND *csound, PHSORBNK *p)
1181 {
1182 double phs;
1183 double *curphs = (double*)p->curphs.auxp;
1184 int32_t size = p->curphs.size / sizeof(double);
1185 int32_t index = (int32_t)(*p->kindx);
1186
1187 if (UNLIKELY(curphs == NULL)) {
1188 return csound->PerfError(csound, &(p->h),
1189 Str("phasorbnk: not initialised"));
1190 }
1191
1192 if (UNLIKELY(index<0 || index>=size)) {
1193 *p->sr = FL(0.0);
1194 return NOTOK;
1195 }
1196
1197 *p->sr = (MYFLT)(phs = curphs[index]);
1198 if (UNLIKELY((phs += *p->xcps * csound->onedkr) >= 1.0))
1199 phs -= 1.0;
1200 else if (UNLIKELY(phs < 0.0)) /* patch from Matthew Scala */
1201 phs += 1.0;
1202 curphs[index] = phs;
1203 return OK;
1204 }
1205
phsorbnk(CSOUND * csound,PHSORBNK * p)1206 int32_t phsorbnk(CSOUND *csound, PHSORBNK *p)
1207 {
1208 uint32_t offset = p->h.insdshead->ksmps_offset;
1209 uint32_t early = p->h.insdshead->ksmps_no_end;
1210 uint32_t n, nsmps = CS_KSMPS;
1211 MYFLT *rs;
1212 double phase, incr;
1213 double *curphs = (double*)p->curphs.auxp;
1214 int32_t size = p->curphs.size / sizeof(double);
1215 int32_t index = (int32_t)(*p->kindx);
1216
1217 if (UNLIKELY(curphs == NULL)) {
1218 return csound->PerfError(csound, &(p->h),
1219 Str("phasorbnk: not initialised"));
1220 }
1221
1222 if (UNLIKELY(index<0 || index>=size)) {
1223 *p->sr = FL(0.0);
1224 return NOTOK;
1225 }
1226
1227 rs = p->sr;
1228 phase = curphs[index];
1229 if (UNLIKELY(offset)) memset(rs, '\0', offset*sizeof(MYFLT));
1230 if (UNLIKELY(early)) {
1231 nsmps -= early;
1232 memset(&rs[nsmps], '\0', early*sizeof(MYFLT));
1233 }
1234 if (IS_ASIG_ARG(p->xcps)) {
1235 MYFLT *cps = p->xcps;
1236 for (n=offset; n<nsmps; n++) {
1237 incr = (double)(cps[n] * csound->onedsr);
1238 rs[n] = (MYFLT)phase;
1239 phase += incr;
1240 if (UNLIKELY(phase >= 1.0))
1241 phase -= 1.0;
1242 else if (UNLIKELY(phase < 0.0))
1243 phase += 1.0;
1244 }
1245 }
1246 else {
1247 incr = (double)(*p->xcps * csound->onedsr);
1248 for (n=offset; n<nsmps; n++) {
1249 rs[n] = (MYFLT)phase;
1250 phase += incr;
1251 if (UNLIKELY(phase >= 1.0))
1252 phase -= 1.0;
1253 else if (UNLIKELY(phase < 0.0))
1254 phase += 1.0;
1255 }
1256 }
1257 curphs[index] = phase;
1258 return OK;
1259 }
1260
1261 /* Opcodes from rasmus ekman */
1262
1263 /* pinkish: Two methods for pink-type noise generation
1264 The Moore/Gardner method, coded by Phil Burke, optimised by James McCartney;
1265 Paul Kellet's -3dB/octave white->pink filter bank, "refined" version;
1266 Paul Kellet's -3dB/octave white->pink filter bank, "economy" version
1267
1268 The Moore/Gardner method output seems to have bumps in the low-mid and
1269 mid-high ranges.
1270 The Kellet method (refined) has smooth spectrum, but goes up slightly
1271 at the far high end.
1272 */
1273
1274 #define GARDNER_PINK FL(0.0)
1275 #define KELLET_PINK FL(1.0)
1276 #define KELLET_CHEAP_PINK FL(2.0)
1277
1278 int32_t GardnerPink_init(CSOUND *csound, PINKISH *p);
1279 int32_t GardnerPink_perf(CSOUND *csound, PINKISH *p);
1280
pinkset(CSOUND * csound,PINKISH * p)1281 int32_t pinkset(CSOUND *csound, PINKISH *p)
1282 {
1283 /* Check valid method */
1284 if (UNLIKELY(*p->imethod != GARDNER_PINK && *p->imethod != KELLET_PINK
1285 && *p->imethod != KELLET_CHEAP_PINK)) {
1286 return csound->InitError(csound, Str("pinkish: Invalid method code"));
1287 }
1288 /* User range scaling can be a- or k-rate for Gardner, a-rate only
1289 for filter */
1290 if (IS_ASIG_ARG(p->xin)) {
1291 p->ampinc = 1;
1292 }
1293 else {
1294 /* Cannot accept k-rate input with filter method */
1295 if (UNLIKELY(*p->imethod != FL(0.0))) {
1296 return csound->InitError(csound, Str(
1297 "pinkish: Filter method requires a-rate (noise) input"));
1298 }
1299 p->ampinc = 0;
1300 }
1301 /* Unless we're reinitializing a tied note, zero coefs */
1302 if (*p->iskip != FL(1.0)) {
1303 if (*p->imethod == GARDNER_PINK)
1304 GardnerPink_init(csound,p);
1305 else /* Filter method */
1306 p->b0 = p->b1 = p->b2 = p->b3 = p->b4 = p->b5 = p->b6 = FL(0.0);
1307 }
1308 return OK;
1309 }
1310
pinkish(CSOUND * csound,PINKISH * p)1311 int32_t pinkish(CSOUND *csound, PINKISH *p)
1312 {
1313 MYFLT *aout, *ain;
1314 double c0, c1, c2, c3, c4, c5, c6, nxtin, nxtout;
1315 uint32_t offset = p->h.insdshead->ksmps_offset;
1316 uint32_t early = p->h.insdshead->ksmps_no_end;
1317 uint32_t n, nsmps = CS_KSMPS;
1318 aout = p->aout;
1319 ain = p->xin;
1320 if (UNLIKELY(offset)) memset(aout, '\0', offset*sizeof(MYFLT));
1321 if (UNLIKELY(early)) {
1322 nsmps -= early;
1323 memset(&aout[nsmps], '\0', early*sizeof(MYFLT));
1324 }
1325 if (*p->imethod == GARDNER_PINK) { /* Gardner method (default) */
1326 GardnerPink_perf(csound,p);
1327 }
1328 else if (*p->imethod == KELLET_PINK) {
1329 /* Paul Kellet's "refined" pink filter */
1330 /* Get filter states */
1331 c0 = p->b0; c1 = p->b1; c2 = p->b2;
1332 c3 = p->b3; c4 = p->b4; c5 = p->b5; c6 = p->b6;
1333 for (n=offset;n<nsmps;n++) {
1334 nxtin = (double)ain[n];
1335 c0 = c0 * 0.99886 + nxtin * 0.0555179;
1336 c1 = c1 * 0.99332 + nxtin * 0.0750759;
1337 c2 = c2 * 0.96900 + nxtin * 0.1538520;
1338 c3 = c3 * 0.86650 + nxtin * 0.3104856;
1339 c4 = c4 * 0.55000 + nxtin * 0.5329522;
1340 c5 = c5 * -0.7616 - nxtin * 0.0168980;
1341 nxtout = c0 + c1 + c2 + c3 + c4 + c5 + c6 + nxtin * 0.5362;
1342 aout[n] = (MYFLT)(nxtout * 0.11); /* (roughly) compensate for gain */
1343 c6 = nxtin * 0.115926;
1344 }
1345 /* Store back filter coef states */
1346 p->b0 = c0; p->b1 = c1; p->b2 = c2;
1347 p->b3 = c3; p->b4 = c4; p->b5 = c5; p->b6 = c6;
1348 }
1349 else if (*p->imethod == KELLET_CHEAP_PINK) {
1350 /* Get filter states */
1351 c0 = p->b0; c1 = p->b1; c2 = p->b2;
1352
1353 for (n=offset;n<nsmps;n++) { /* Paul Kellet's "economy" pink filter */
1354 nxtin = (double)ain[n];
1355 c0 = c0 * 0.99765 + nxtin * 0.0990460;
1356 c1 = c1 * 0.96300 + nxtin * 0.2965164;
1357 c2 = c2 * 0.57000 + nxtin * 1.0526913;
1358 nxtout = c0 + c1 + c2 + nxtin * 0.1848;
1359 aout[n] = (MYFLT)(nxtout * 0.11); /* (roughly) compensate for gain */
1360 }
1361
1362 /* Store back filter coef states */
1363 p->b0 = c0; p->b1 = c1; p->b2 = c2;
1364 }
1365 return OK;
1366 }
1367
1368 /************************************************************/
1369 /*
1370 GardnerPink_init() and GardnerPink_perf()
1371
1372 Generate Pink Noise using Gardner method.
1373 Optimization suggested by James McCartney uses a tree
1374 to select which random value to replace.
1375
1376 x x x x x x x x x x x x x x x x
1377 x x x x x x x x
1378 x x x x
1379 x x
1380 x
1381
1382 Tree is generated by counting trailing zeros in an increasing index.
1383 When the index is zero, no random number is selected.
1384
1385 Author: Phil Burk, http://www.softsynth.com
1386
1387 Revision History:
1388 Csound version by rasmus ekman May 2000
1389 Several changes, some marked "(re)"
1390
1391 Copyright 1999 Phil Burk - No rights reserved.
1392 */
1393
1394 /************************************************************/
1395
1396 /* Yet another pseudo-random generator. Could probably be changed
1397 for any of the other available ones in Csound */
1398
1399 #define PINK_RANDOM_BITS (24)
1400 /* Left-shift one bit less 24 to allow negative values (re) */
1401 #define PINK_RANDOM_SHIFT (7)
1402
1403 /* Calculate pseudo-random 32 bit number based on linear congruential method. */
GenerateRandomNumber(uint32 randSeed)1404 static int32 GenerateRandomNumber(uint32 randSeed)
1405 {
1406 randSeed = ((uint32_t) randSeed * 196314165U) + 907633515UL;
1407 return (int32) ((int32_t) ((uint32_t) randSeed));
1408 }
1409
1410 /************************************************************/
1411
1412 /* Set up for user-selected number of bands of noise generators. */
GardnerPink_init(CSOUND * csound,PINKISH * p)1413 int32_t GardnerPink_init(CSOUND *csound, PINKISH *p)
1414 {
1415 int32_t i;
1416 MYFLT pmax;
1417 int32 numRows;
1418
1419 /* Set number of rows to use (default to 20) */
1420 if (*p->iparam1 >= 4 && *p->iparam1 <= GRD_MAX_RANDOM_ROWS)
1421 p->grd_NumRows = (int32)*p->iparam1;
1422 else {
1423 p->grd_NumRows = 20;
1424 /* Warn if user tried but failed to give sensible number */
1425 if (UNLIKELY(*p->iparam1 != FL(0.0)))
1426 csound->Warning(csound, Str("pinkish: Gardner method requires 4-%d bands. "
1427 "Default %"PRIi32" substituted for %d.\n"),
1428 GRD_MAX_RANDOM_ROWS, p->grd_NumRows,
1429 (int32_t) *p->iparam1);
1430 }
1431
1432 /* Seed random generator by user value or by time (default) */
1433 if (*p->iseed != FL(0.0)) {
1434 if (*p->iseed > -1.0 && *p->iseed < 1.0)
1435 p->randSeed = (uint32) (*p->iseed * (MYFLT)0x80000000);
1436 else p->randSeed = (uint32) *p->iseed;
1437 }
1438 else p->randSeed = (uint32) csound->GetRandomSeedFromTime();
1439
1440 numRows = p->grd_NumRows;
1441 p->grd_Index = 0;
1442 if (numRows == 32) p->grd_IndexMask = 0xFFFFFFFF;
1443 else p->grd_IndexMask = (1<<numRows) - 1;
1444
1445 /* Calculate reasonable maximum signed random value. */
1446 /* Tweaked to get sameish peak value over all numRows values (re) */
1447 pmax = (MYFLT)((numRows + 30) * (1<<(PINK_RANDOM_BITS-2)));
1448 p->grd_Scalar = FL(1.0) / pmax;
1449
1450 /* Warm up by filling all rows (re) (original zeroed all rows, and runningSum) */
1451 {
1452 int32 randSeed, newRandom, runningSum = 0;
1453 randSeed = p->randSeed;
1454 for (i = 0; i < numRows; i++) {
1455 randSeed = GenerateRandomNumber(randSeed);
1456 newRandom = randSeed >> PINK_RANDOM_SHIFT;
1457 runningSum += newRandom;
1458 p->grd_Rows[i] = newRandom;
1459 }
1460 p->grd_RunningSum = runningSum;
1461 p->randSeed = randSeed;
1462 }
1463 return OK;
1464 }
1465
1466 /* Generate numRows octave-spaced white bands and sum to pink noise. */
GardnerPink_perf(CSOUND * csound,PINKISH * p)1467 int32_t GardnerPink_perf(CSOUND *csound, PINKISH *p)
1468 {
1469 IGN(csound);
1470 MYFLT *aout, *amp, scalar;
1471 int32 *rows, rowIndex, indexMask, randSeed, newRandom;
1472 int32 runningSum, sum, ampinc;
1473 uint32_t offset = p->h.insdshead->ksmps_offset;
1474 uint32_t n, nsmps = CS_KSMPS - p->h.insdshead->ksmps_no_end;
1475
1476 aout = p->aout;
1477 amp = p->xin;
1478 ampinc = p->ampinc; /* Used to increment user amp if a-rate */
1479 scalar = p->grd_Scalar;
1480 rowIndex = p->grd_Index;
1481 indexMask = p->grd_IndexMask;
1482 runningSum = p->grd_RunningSum;
1483 rows = &(p->grd_Rows[0]);
1484 randSeed = p->randSeed;
1485
1486 for (n=offset; n<nsmps;n++) {
1487 /* Increment and mask index. */
1488 rowIndex = (rowIndex + 1) & indexMask;
1489
1490 /* If index is zero, do not update any random values. */
1491 if ( rowIndex != 0 ) {
1492 /* Determine how many trailing zeros in PinkIndex. */
1493 /* This algorithm will hang if n==0 so test first. */
1494 int32_t numZeros = 0;
1495 int32_t n = rowIndex;
1496 while( (n & 1) == 0 ) {
1497 n = n >> 1;
1498 numZeros++;
1499 }
1500
1501 /* Replace the indexed ROWS random value.
1502 * Subtract and add back to RunningSum instead of adding all
1503 * the random values together. Only one changes each time.
1504 */
1505 runningSum -= rows[numZeros];
1506 randSeed = GenerateRandomNumber(randSeed);
1507 newRandom = randSeed >> PINK_RANDOM_SHIFT;
1508 runningSum += newRandom;
1509 rows[numZeros] = newRandom;
1510 }
1511
1512 /* Add extra white noise value. */
1513 randSeed = GenerateRandomNumber(randSeed);
1514 newRandom = randSeed >> PINK_RANDOM_SHIFT;
1515 sum = runningSum + newRandom;
1516
1517 /* Scale to range of +/-p->xin (user-selected amp) */
1518 aout[n] = *amp * sum * scalar;
1519 amp += ampinc; /* Increment if amp is a-rate */
1520 }
1521
1522 p->grd_RunningSum = runningSum;
1523 p->grd_Index = rowIndex;
1524 p->randSeed = randSeed;
1525 return OK;
1526 }
1527
1528 /* ************************************************************ */
1529 /* A collection of clipping techniques -- JPff */
1530 /* Method 0: Bram de Jong <Bram.DeJong@rug.ac.be> */
1531 /* x > a: f(x) = a + (x-a)/(1+((x-a)/(1-a))^2) */
1532 /* x > 1: f(x) = (a+1)/2 */
1533 /* JPff scaled this to a limit and a fraction */
1534 /* Method 1: */
1535 /* |x|<limit f(x) = limit * sin(pi x/(2*limit) */
1536 /* f(x) = limit * sign(x) */
1537 /* Method 2: */
1538 /* |x|<limit f(x) = limit * tanh(x/limit)/tanh(1) */
1539 /* f(x) = limit * sign(x) */
1540 /* ************************************************************ */
1541
1542 /* Methods 0 and 2 OK, method1 broken */
1543
1544 double tanh(double);
clip_set(CSOUND * csound,CLIP * p)1545 int32_t clip_set(CSOUND *csound, CLIP *p)
1546 {
1547 IGN(csound);
1548 int32_t meth = (int32_t)MYFLT2LONG(*p->imethod);
1549 p->meth = meth;
1550 p->arg = FABS(*p->iarg);
1551 p->lim = *p->limit;
1552 switch (meth) {
1553 case 0: /* Bram de Jong method */
1554 if (p->arg > FL(1.0) || p->arg < FL(0.0)) p->arg = FL(0.999);
1555 p->arg = p->lim * p->arg;
1556 p->k1 = FL(1.0)/(p->lim - p->arg);
1557 p->k1 = p->k1 * p->k1;
1558 p->k2 = (p->lim + p->arg)*FL(0.5);
1559 break;
1560 case 1:
1561 p->k1 = PI_F/(FL(2.0) * p->lim);
1562 break;
1563 case 2:
1564 p->k1 = FL(1.0)/TANH(FL(1.0));
1565 break;
1566 default:
1567 p->meth = 0;
1568 }
1569 return OK;
1570 }
1571
clip(CSOUND * csound,CLIP * p)1572 int32_t clip(CSOUND *csound, CLIP *p)
1573 {
1574 IGN(csound);
1575 MYFLT *aout = p->aout, *ain = p->ain;
1576 uint32_t offset = p->h.insdshead->ksmps_offset;
1577 uint32_t early = p->h.insdshead->ksmps_no_end;
1578 uint32_t n, nsmps = CS_KSMPS;
1579 MYFLT a = p->arg, k1 = p->k1, k2 = p->k2;
1580 MYFLT limit = p->lim;
1581 MYFLT rlim = FL(1.0)/limit;
1582
1583 if (UNLIKELY(offset)) memset(aout, '\0', offset*sizeof(MYFLT));
1584 if (UNLIKELY(early)) {
1585 nsmps -= early;
1586 memset(&aout[nsmps], '\0', early*sizeof(MYFLT));
1587 }
1588 switch (p->meth) {
1589 case 0: /* Soft clip with division */
1590 for (n=offset;n<nsmps;n++) {
1591 MYFLT x = ain[n];
1592 if (x>=FL(0.0)) {
1593 if (UNLIKELY(x>limit)) x = k2;
1594 else if (x>a)
1595 x = a + (x-a)/(FL(1.0)+(x-a)*(x-a)*k1);
1596 }
1597 else {
1598 if (UNLIKELY(x<-limit))
1599 x = -k2;
1600 else if (-x>a)
1601 x = -a + (x+a)/(FL(1.0)+(x+a)*(x+a)*k1);
1602 }
1603 aout[n] = x;
1604 }
1605 return OK;
1606 case 1:
1607 for (n=offset;n<nsmps;n++) {
1608 MYFLT x = ain[n];
1609 if (UNLIKELY(x>=limit))
1610 x = limit;
1611 else if (UNLIKELY(x<= -limit))
1612 x = -limit;
1613 else
1614 x = limit*SIN(k1*x);
1615 aout[n] = x;
1616 }
1617 return OK;
1618 case 2:
1619 for (n=offset;n<nsmps;n++) {
1620 MYFLT x = ain[n];
1621 if (UNLIKELY(x>=limit))
1622 x = limit;
1623 else if (UNLIKELY(x<= -limit))
1624 x = -limit;
1625 else
1626 x = limit*k1*TANH(x*rlim);
1627 aout[n] = x;
1628 }
1629 return OK;
1630 }
1631 return OK;
1632 }
1633
1634 /* ********************************************************************** */
1635 /* *************** IMPULSE ********************************************** */
1636 /* ********************************************************************** */
1637
impulse_set(CSOUND * csound,IMPULSE * p)1638 int32_t impulse_set(CSOUND *csound, IMPULSE *p)
1639 {
1640 p->next = (uint32_t)MYFLT2LONG(*p->offset * CS_ESR);
1641 return OK;
1642 }
1643
impulse(CSOUND * csound,IMPULSE * p)1644 int32_t impulse(CSOUND *csound, IMPULSE *p)
1645 {
1646 uint32_t offset = p->h.insdshead->ksmps_offset;
1647 uint32_t early = p->h.insdshead->ksmps_no_end;
1648 uint32_t n, nsmps = CS_KSMPS;
1649 int32_t next = p->next;
1650 MYFLT *ar = p->ar;
1651 if (next<0) next = -next;
1652 if (UNLIKELY(next < (int32)nsmps)) { /* Impulse in this frame */
1653 MYFLT frq = *p->freq; /* Freq at k-rate */
1654 int32_t sfreq; /* Converted to samples */
1655 if (frq == FL(0.0)) sfreq = INT_MAX; /* Zero means infinite */
1656 else if (frq < FL(0.0)) sfreq = -(int32_t)frq; /* Negative cnts in sample */
1657 else sfreq = (int32_t)(frq*CS_ESR); /* Normal case */
1658 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
1659 if (UNLIKELY(early)) {
1660 nsmps -= early;
1661 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
1662 }
1663 for (n=offset;n<nsmps;n++) {
1664 if (UNLIKELY(next-- == 0)) {
1665 ar[n] = *p->amp;
1666 next = sfreq - 1; /* Note can be less than k-rate */
1667 }
1668 else ar[n] = FL(0.0);
1669 }
1670 }
1671 else { /* Nothing this time so just fill */
1672 memset(ar, 0, nsmps*sizeof(MYFLT));
1673 next -= nsmps;
1674 }
1675 p->next = next;
1676 return OK;
1677 }
1678
1679 /* ********************************************************************** */
1680 /* Version of CMUSIC trans opcode */
1681 /* creates y0 + (y1 - y0) * (1 - exp( t*alpha )) / (1 - exp(alpha)) */
1682 /* or */
1683 /* y0 + (y1 - y0) * t if alpha is zero */
1684 /* ********************************************************************** */
trnset(CSOUND * csound,TRANSEG * p)1685 int32_t trnset(CSOUND *csound, TRANSEG *p)
1686 {
1687 NSEG *segp;
1688 int32_t nsegs;
1689 MYFLT **argp, val;
1690
1691 if (UNLIKELY(p->INOCOUNT%3!=1))
1692 return csound->InitError(csound, Str("Incorrect argument count in transeg"));
1693 nsegs = p->INOCOUNT / 3; /* count segs & alloc if nec */
1694 if ((segp = (NSEG *) p->auxch.auxp) == NULL ||
1695 (uint32_t)p->auxch.size < nsegs*sizeof(NSEG)) {
1696 csound->AuxAlloc(csound, (int32)nsegs*sizeof(NSEG), &p->auxch);
1697 p->cursegp = segp = (NSEG *) p->auxch.auxp;
1698 }
1699 segp[nsegs-1].cnt = MAXPOS; /* set endcount for safety */
1700 segp[nsegs-1].acnt = MAXPOS; /* set endcount for safety */
1701 argp = p->argums;
1702 val = **argp++;
1703 if (**argp <= FL(0.0)) return OK; /* if idur1 <= 0, skip init */
1704 p->curval = val;
1705 p->curcnt = 0;
1706 p->cursegp = segp - 1; /* else setup null seg0 */
1707 p->segsrem = nsegs + 1;
1708 p->curx = FL(0.0);
1709 do { /* init each seg .. */
1710 MYFLT dur = **argp++;
1711 MYFLT alpha = **argp++;
1712 MYFLT nxtval = **argp++;
1713 MYFLT d = dur * CS_ESR;
1714 if ((segp->acnt = segp->cnt = (int32)MYFLT2LONG(d)) < 0)
1715 segp->cnt = 0;
1716 else
1717 segp->cnt = (int32)(dur * CS_EKR);
1718 segp->nxtpt = nxtval;
1719 segp->val = val;
1720 if (alpha == FL(0.0)) {
1721 segp->c1 = (nxtval-val)/d;
1722 }
1723 else {
1724 segp->c1 = (nxtval - val)/(FL(1.0) - EXP(alpha));
1725 }
1726 segp->alpha = alpha/d;
1727 val = nxtval;
1728 segp++;
1729 } while (--nsegs);
1730 p->xtra = -1;
1731 p->alpha = ((NSEG*)p->auxch.auxp)[0].alpha;
1732 p->curinc = ((NSEG*)p->auxch.auxp)[0].c1;
1733 return OK;
1734 }
1735
trnset_bkpt(CSOUND * csound,TRANSEG * p)1736 int32_t trnset_bkpt(CSOUND *csound, TRANSEG *p)
1737 {
1738 NSEG *segp;
1739 int32_t nsegs;
1740 MYFLT **argp, val;
1741 MYFLT totdur = FL(0.0);
1742
1743 if (UNLIKELY(p->INOCOUNT%3!=1))
1744 return csound->InitError(csound, Str("Incorrect argument count in transegb"));
1745 nsegs = p->INOCOUNT / 3; /* count segs & alloc if nec */
1746 if ((segp = (NSEG *) p->auxch.auxp) == NULL ||
1747 (uint32_t)p->auxch.size < nsegs*sizeof(NSEG)) {
1748 csound->AuxAlloc(csound, (int32)nsegs*sizeof(NSEG), &p->auxch);
1749 p->cursegp = segp = (NSEG *) p->auxch.auxp;
1750 }
1751 segp[nsegs-1].cnt = MAXPOS; /* set endcount for safety */
1752 argp = p->argums;
1753 val = **argp++;
1754 if (**argp <= FL(0.0)) return OK; /* if idur1 <= 0, skip init */
1755 p->curval = val;
1756 p->curcnt = 0;
1757 p->cursegp = segp - 1; /* else setup null seg0 */
1758 p->segsrem = nsegs + 1;
1759 p->curx = FL(0.0);
1760 do { /* init each seg .. */
1761 MYFLT dur = **argp++;
1762 MYFLT alpha = **argp++;
1763 MYFLT nxtval = **argp++;
1764 MYFLT d;
1765 dur -= totdur;
1766 totdur += dur;
1767 d = dur * CS_ESR;
1768 if ((segp->cnt = (int32)MYFLT2LONG(d)) < 0)
1769 segp->cnt = 0;
1770 else
1771 segp->cnt = (int32)(dur * CS_EKR);
1772 segp->nxtpt = nxtval;
1773 segp->val = val;
1774 if (alpha == FL(0.0)) {
1775 segp->c1 = (nxtval-val)/d;
1776 }
1777 else {
1778 segp->c1 = (nxtval - val)/(FL(1.0) - EXP(alpha));
1779 }
1780 segp->alpha = alpha/d;
1781 val = nxtval;
1782 segp++;
1783 } while (--nsegs);
1784 p->xtra = -1;
1785 p->alpha = ((NSEG*)p->auxch.auxp)[0].alpha;
1786 p->curinc = ((NSEG*)p->auxch.auxp)[0].c1;
1787 return OK;
1788 }
1789
ktrnseg(CSOUND * csound,TRANSEG * p)1790 int32_t ktrnseg(CSOUND *csound, TRANSEG *p)
1791 {
1792 *p->rslt = p->curval; /* put the cur value */
1793 if (UNLIKELY(p->auxch.auxp==NULL)) { /* RWD fix */
1794 csound->PerfError(csound,&(p->h),
1795 Str("Error: transeg not initialised (krate)\n"));
1796 }
1797 if (p->segsrem) { /* done if no more segs */
1798 if (--p->curcnt <= 0) { /* if done cur segment */
1799 NSEG *segp = p->cursegp;
1800 chk1:
1801 if (!(--p->segsrem)) {
1802 p->curval = segp->nxtpt; /* advance the cur val */
1803 return OK;
1804 }
1805 p->cursegp = ++segp; /* find the next */
1806 if (!(p->curcnt = segp->cnt)) { /* nonlen = discontin */
1807 p->curval = segp->nxtpt; /* poslen = new slope */
1808 goto chk1;
1809 }
1810 p->curinc = segp->c1;
1811 p->alpha = segp->alpha;
1812 p->curx = FL(0.0);
1813 }
1814 p->curx += (MYFLT)CS_KSMPS*p->alpha;
1815 if (p->alpha == FL(0.0))
1816 p->curval += p->curinc*CS_KSMPS; /* advance the cur val */
1817 else
1818 p->curval = p->cursegp->val + p->curinc *
1819 (FL(1.0) - EXP(p->curx));
1820 }
1821 return OK;
1822 }
1823
trnseg(CSOUND * csound,TRANSEG * p)1824 int32_t trnseg(CSOUND *csound, TRANSEG *p)
1825 {
1826 MYFLT val, *rs = p->rslt;
1827 uint32_t offset = p->h.insdshead->ksmps_offset;
1828 uint32_t early = p->h.insdshead->ksmps_no_end;
1829 uint32_t n, nsmps = CS_KSMPS;
1830 NSEG *segp = p->cursegp;
1831 if (UNLIKELY(p->auxch.auxp==NULL)) {
1832 return csound->PerfError(csound, &(p->h),
1833 Str("transeg: not initialised (arate)\n"));
1834 }
1835 if (UNLIKELY(offset)) memset(rs, '\0', offset*sizeof(MYFLT));
1836 if (UNLIKELY(early)) {
1837 nsmps -= early;
1838 memset(&rs[nsmps], '\0', early*sizeof(MYFLT));
1839 }
1840 val = p->curval; /* sav the cur value */
1841 for (n=offset; n<nsmps; n++) {
1842 if (p->segsrem) { /* if no more segs putk */
1843 if (--p->curcnt <= 0) { /* if done cur segment */
1844 segp = p->cursegp;
1845 chk1:
1846 if (UNLIKELY(!--p->segsrem)) { /* if none left */
1847 val = p->curval = segp->nxtpt;
1848 goto putk; /* put endval */
1849 }
1850 p->cursegp = ++segp; /* else find the next */
1851 if (!(p->curcnt = segp->acnt)) {
1852 val = p->curval = segp->nxtpt; /* nonlen = discontin */
1853 goto chk1;
1854 } /* poslen = new slope */
1855 p->curinc = segp->c1;
1856 p->alpha = segp->alpha;
1857 p->curx = FL(0.0);
1858 p->curval = val;
1859 }
1860 if (p->alpha == FL(0.0)) {
1861 rs[n] = val;
1862 val += p->curinc;
1863 }
1864 else {
1865 rs[n] = val;
1866 p->curx += p->alpha;
1867 val = segp->val + p->curinc *
1868 (FL(1.0) - EXP(p->curx));
1869 }
1870 }
1871 else{
1872 putk:
1873 rs[n] = val;
1874 }
1875 }
1876 p->curval = val;
1877 return OK;
1878 }
1879
1880 /* MIDI aware version of transeg */
trnsetr(CSOUND * csound,TRANSEG * p)1881 int32_t trnsetr(CSOUND *csound, TRANSEG *p)
1882 {
1883 int32_t relestim;
1884 NSEG *segp;
1885 int32_t nsegs;
1886 MYFLT **argp;
1887 double val;
1888
1889 if (UNLIKELY(p->INOCOUNT%3!=1))
1890 return csound->InitError(csound, Str("Incorrect argument count in transegr"));
1891 nsegs = p->INOCOUNT / 3; /* count segs & alloc if nec */
1892 if ((segp = (NSEG *) p->auxch.auxp) == NULL ||
1893 (uint32_t)p->auxch.size < nsegs*sizeof(NSEG)) {
1894 csound->AuxAlloc(csound, (int32)nsegs*sizeof(NSEG), &p->auxch);
1895 p->cursegp = segp = (NSEG *) p->auxch.auxp;
1896 }
1897 segp[nsegs-1].cnt = MAXPOS; /* set endcount for safety */
1898 segp[nsegs-1].acnt = MAXPOS; /* set endcount for safety */
1899 argp = p->argums;
1900 val = (double)**argp++;
1901 if (UNLIKELY(**argp <= FL(0.0))) return OK; /* if idur1 <= 0, skip init */
1902 p->curval = val;
1903 p->curcnt = 0;
1904 p->cursegp = segp - 1; /* else setup null seg0 */
1905 p->segsrem = nsegs + 1;
1906 p->curx = FL(0.0);
1907 do { /* init each seg .. */
1908 double dur = (double)**argp++;
1909 MYFLT alpha = **argp++;
1910 MYFLT nxtval = **argp++;
1911 MYFLT d = dur * CS_ESR;
1912 if ((segp->acnt = segp->cnt = (int32)(d + FL(0.5))) < 0)
1913 segp->cnt = 0;
1914 else
1915 segp->cnt = (int32)(dur * CS_EKR);
1916 segp->nxtpt = nxtval;
1917 segp->val = val;
1918 if (alpha == FL(0.0)) {
1919 segp->c1 = (nxtval-val)/d;
1920 //printf("alpha zero val=%f c1=%f\n", segp->val, segp->c1);
1921 }
1922 else {
1923 p->lastalpha = alpha;
1924 segp->c1 = (nxtval - val)/(FL(1.0) - EXP(alpha));
1925 }
1926 segp->alpha = alpha/d;
1927 val = nxtval;
1928 segp++;
1929 p->finalval = nxtval;
1930 } while (--nsegs);
1931 //p->xtra = -1;
1932 p->alpha = ((NSEG*)p->auxch.auxp)[0].alpha;
1933 p->curinc = ((NSEG*)p->auxch.auxp)[0].c1;
1934 relestim = (int32_t)(p->cursegp + p->segsrem - 1)->cnt;
1935 p->xtra = relestim;
1936 if (relestim > p->h.insdshead->xtratim)
1937 p->h.insdshead->xtratim = (int32_t)relestim;
1938 /* { */
1939 /* int32_t i; */
1940 /* int32_t nseg = p->INOCOUNT / 3; */
1941 /* NSEG *segp = p->cursegp; */
1942 /* for (i=0; i<nseg; i++) */
1943 /* printf("cnt=%d alpha=%f val=%f nxtpt=%f c1=%f\n", */
1944 /* segp[i].cnt, segp[i].alpha, segp[i].val, segp[i].nxtpt, segp[i].c1); */
1945 /* } */
1946 return OK;
1947 }
1948
ktrnsegr(CSOUND * csound,TRANSEG * p)1949 int32_t ktrnsegr(CSOUND *csound, TRANSEG *p)
1950 {
1951 *p->rslt = p->curval; /* put the cur value */
1952 if (UNLIKELY(p->auxch.auxp==NULL)) { /* RWD fix */
1953 csound->PerfError(csound,&(p->h),
1954 Str("Error: transeg not initialised (krate)\n"));
1955 }
1956 if (p->segsrem) { /* done if no more segs */
1957 NSEG *segp;
1958 if (p->h.insdshead->relesing && p->segsrem > 1) {
1959 //printf("releasing\n");
1960 while (p->segsrem > 1) { /* reles flag new: */
1961 segp = ++p->cursegp; /* go to last segment */
1962 p->segsrem--;
1963 } /* get univ relestim */
1964 segp->cnt = p->xtra>=0 ? p->xtra : p->h.insdshead->xtratim;
1965 if (segp->alpha == FL(0.0)) {
1966 segp->c1 = (p->finalval-p->curval)/(segp->cnt*CS_KSMPS);
1967 //printf("finalval = %f curval = %f, cnt = %d c1 = %f\n",
1968 // p->finalval, p->curval, segp->cnt, segp->c1);
1969 }
1970 else {
1971 segp->c1 = (p->finalval - p->curval)/(FL(1.0) - EXP(p->lastalpha));
1972 segp->alpha = p->lastalpha/(segp->cnt*CS_KSMPS);
1973 segp->val = p->curval;
1974 }
1975 goto newm; /* and set new curmlt */
1976 }
1977 if (--p->curcnt <= 0) { /* if done cur segment */
1978 chk1:
1979 if (p->segsrem == 2) return OK; /* seg Y rpts lastval */
1980 if (!(--p->segsrem)) return OK; /* seg Z now done all */
1981 segp = ++p->cursegp; /* find the next */
1982 newm:
1983 //printf("curcnt = %d seg/cnt = %d\n", p->curcnt, segp->cnt);
1984 if (!(p->curcnt = segp->cnt)) { /* nonlen = discontin */
1985 p->curval = segp->nxtpt; /* poslen = new slope */
1986 //printf("curval = %f\n", p->curval);
1987 goto chk1;
1988 }
1989 p->curinc = segp->c1;
1990 p->alpha = segp->alpha;
1991 p->curx = FL(0.0);
1992 }
1993 if (p->alpha == FL(0.0)) {
1994 p->curval += p->curinc *CS_KSMPS; /* advance the cur val */
1995 //printf("curval = %f\n", p->curval);
1996 }
1997 else
1998 p->curval = p->cursegp->val + (p->curinc) *
1999 (FL(1.0) - EXP(p->curx));
2000 p->curx += (MYFLT)CS_KSMPS* p->alpha;
2001 }
2002 return OK;
2003 }
2004
trnsegr(CSOUND * csound,TRANSEG * p)2005 int32_t trnsegr(CSOUND *csound, TRANSEG *p)
2006 {
2007 MYFLT val, *rs = p->rslt;
2008 uint32_t offset = p->h.insdshead->ksmps_offset;
2009 uint32_t early = p->h.insdshead->ksmps_no_end;
2010 uint32_t n, nsmps = CS_KSMPS;
2011 if (UNLIKELY(p->auxch.auxp==NULL)) {
2012 return csound->PerfError(csound, &(p->h),
2013 Str("transeg: not initialised (arate)\n"));
2014 }
2015 if (UNLIKELY(offset)) memset(rs, '\0', offset*sizeof(MYFLT));
2016 if (UNLIKELY(early)) {
2017 nsmps -= early;
2018 memset(&rs[nsmps], '\0', early*sizeof(MYFLT));
2019 }
2020 val = p->curval; /* sav the cur value */
2021 for (n=offset; n<nsmps; n++) {
2022 if (LIKELY(p->segsrem)) { /* if no more segs putk */
2023 NSEG *segp;
2024 if (p->h.insdshead->relesing && p->segsrem > 1) {
2025 while (p->segsrem > 1) { /* if release flag new */
2026 segp = ++p->cursegp; /* go to last segment */
2027 p->segsrem--;
2028 } /* get univ relestim */
2029 segp->cnt = p->xtra>=0 ? p->xtra : p->h.insdshead->xtratim;
2030 if (segp->alpha == FL(0.0)) {
2031 segp->c1 = (p->finalval-val)/segp->acnt;
2032 }
2033 else {
2034 /* this is very wrong */
2035 segp->c1 = (p->finalval - val)/(FL(1.0) - EXP(p->lastalpha));
2036 segp->alpha = p->lastalpha/segp->acnt;
2037 segp->val = val;
2038 }
2039 goto newm; /* and set new curmlt */
2040 }
2041 if (--p->curcnt <= 0) { /* if done cur segment */
2042 //segp = p->cursegp; /* overwritten later -- coverity */
2043 chk1:
2044 if (p->segsrem == 2) goto putk; /* seg Y rpts lastval */
2045 if (UNLIKELY(!--p->segsrem)) { /* if none left */
2046 //val = p->curval = segp->nxtpt;
2047 goto putk; /* put endval */
2048 }
2049 segp = ++p->cursegp; /* else find the next */
2050 newm:
2051 if (!(p->curcnt = segp->acnt)) {
2052 val = p->curval = segp->nxtpt; /* nonlen = discontin */
2053 goto chk1;
2054 } /* poslen = new slope */
2055 p->curinc = segp->c1;
2056 p->alpha = segp->alpha;
2057 p->curx = FL(0.0);
2058 p->curval = val;
2059 }
2060 if (p->alpha == FL(0.0)) {
2061 rs[n] = val;
2062 val += p->curinc;
2063 }
2064 else {
2065 segp = p->cursegp;
2066 rs[n] = val;
2067 p->curx += p->alpha;
2068 val = segp->val + p->curinc * (FL(1.0) - EXP(p->curx));
2069 }
2070 }
2071 else {
2072 putk:
2073 rs[n] = val;
2074 }
2075 }
2076 p->curval = val;
2077
2078 return OK;
2079 }
2080
2081 extern int32 randint31(int32);
2082
varicolset(CSOUND * csound,VARI * p)2083 int32_t varicolset(CSOUND *csound, VARI *p)
2084 {
2085 IGN(csound);
2086 p->last = FL(0.0);
2087 p->lastbeta = *p->beta;
2088 p->sq1mb2 = SQRT(FL(1.0)-p->lastbeta * p->lastbeta);
2089 p->ampmod = FL(0.785)/(FL(1.0)+p->lastbeta);
2090 p->ampinc = IS_ASIG_ARG(p->kamp) ? 1 : 0;
2091 return OK;
2092 }
2093
varicol(CSOUND * csound,VARI * p)2094 int32_t varicol(CSOUND *csound, VARI *p)
2095 {
2096 uint32_t offset = p->h.insdshead->ksmps_offset;
2097 uint32_t early = p->h.insdshead->ksmps_no_end;
2098 uint32_t n, nsmps = CS_KSMPS;
2099 MYFLT beta = *p->beta;
2100 MYFLT sq1mb2 = p->sq1mb2;
2101 MYFLT lastx = p->last;
2102 MYFLT ampmod = p->ampmod;
2103 MYFLT *kamp = p->kamp;
2104 int32_t ampinc = p->ampinc;
2105 MYFLT *rslt = p->rslt;
2106
2107 if (beta != p->lastbeta) {
2108 beta = p->lastbeta = *p->beta;
2109 sq1mb2 = p->sq1mb2 = SQRT(FL(1.0)-p->lastbeta * p->lastbeta);
2110 ampmod = p->ampmod = FL(0.785)/(FL(1.0)+p->lastbeta);
2111 }
2112
2113 if (UNLIKELY(offset)) memset(rslt, '\0', offset*sizeof(MYFLT));
2114 if (UNLIKELY(early)) {
2115 nsmps -= early;
2116 memset(&rslt[nsmps], '\0', early*sizeof(MYFLT));
2117 }
2118 for (n=offset; n<nsmps; n++) {
2119 MYFLT rnd = FL(2.0) * (MYFLT) rand_31(csound) / FL(2147483645) - FL(1.0);
2120 lastx = lastx * beta + sq1mb2 * rnd;
2121 rslt[n] = lastx * *kamp * ampmod;
2122 kamp += ampinc;
2123 }
2124 p->last = lastx;
2125 return OK;
2126 }
2127
2128 /* ************************************************************************ */
2129 /* ***** Josep Comajuncosas' 18dB/oct resonant 3-pole LPF with tanh dist ** */
2130 /* ***** Coded in C by John ffitch, 2000 Dec 17 *************************** */
2131 /* ************************************************************************ */
2132 #include <math.h>
2133
2134 /* This code is transcribed from a Csound macro, so no real comments */
2135
lpf18set(CSOUND * csound,LPF18 * p)2136 int32_t lpf18set(CSOUND *csound, LPF18 *p)
2137 {
2138 IGN(csound);
2139 /* Initialise delay lines */
2140 if (*p->istor==FL(0.0)) {
2141 p->ay1 = FL(0.0);
2142 p->ay2 = FL(0.0);
2143 p->aout = FL(0.0);
2144 p->lastin = FL(0.0);
2145 }
2146 return OK;
2147 }
2148
lpf18db(CSOUND * csound,LPF18 * p)2149 int32_t lpf18db(CSOUND *csound, LPF18 *p)
2150 {
2151 uint32_t offset = p->h.insdshead->ksmps_offset;
2152 uint32_t early = p->h.insdshead->ksmps_no_end;
2153 uint32_t n, nsmps = CS_KSMPS;
2154 MYFLT ay1 = p->ay1;
2155 MYFLT ay2 = p->ay2;
2156 MYFLT aout = p->aout;
2157 MYFLT *ain = p->ain;
2158 MYFLT *ar = p->ar;
2159 MYFLT lastin = p->lastin;
2160 double value = 0.0;
2161 int32_t flag = 1;
2162 MYFLT lfc=0, lrs=0, kres=0, kfcn=0, kp=0, kp1=0, kp1h=0;
2163 double lds = 0.0;
2164 MYFLT zerodb = csound->e0dbfs;
2165 int32_t asgf = IS_ASIG_ARG(p->fco), asgr = IS_ASIG_ARG(p->res),
2166 asgd = IS_ASIG_ARG(p->dist);
2167
2168 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
2169 if (UNLIKELY(early)) {
2170 nsmps -= early;
2171 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
2172 }
2173 for (n=offset;n<nsmps;n++) {
2174 MYFLT fco, res, dist;
2175 MYFLT ax1 = lastin;
2176 MYFLT ay11 = ay1;
2177 MYFLT ay31 = ay2;
2178 fco = (asgf ? p->fco[n] : *p->fco);
2179 res = (asgr ? p->res[n] : *p->res);
2180 dist = (double)(asgd ? p->dist[n] : *p->dist);
2181 if (fco != lfc || flag) {
2182 lfc = fco;
2183 kfcn = FL(2.0) * fco * csound->onedsr;
2184 kp = ((-FL(2.7528)*kfcn + FL(3.0429))*kfcn +
2185 FL(1.718))*kfcn - FL(0.9984);
2186 kp1 = kp+FL(1.0);
2187 kp1h = FL(0.5)*kp1;
2188 flag = 1;
2189 }
2190 if (res != lrs || flag) {
2191 lrs = res;
2192 kres = res * (((-FL(2.7079)*kp1 + FL(10.963))*kp1
2193 - FL(14.934))*kp1 + FL(8.4974));
2194 flag = 1;
2195 }
2196 if (dist != lds || flag) {
2197 lds = dist;
2198 value = 1.0+(dist*(1.5+2.0*(double)kres*(1.0-(double)kfcn)));
2199 }
2200 flag = 0;
2201 lastin = ain[n]/zerodb - TANH(kres*aout);
2202 ay1 = kp1h * (lastin + ax1) - kp*ay1;
2203 ay2 = kp1h * (ay1 + ay11) - kp*ay2;
2204 aout = kp1h * (ay2 + ay31) - kp*aout;
2205
2206 ar[n] = TANH(aout*value)*zerodb;
2207 }
2208 p->ay1 = ay1;
2209 p->ay2 = ay2;
2210 p->aout = aout;
2211 p->lastin = lastin;
2212 return OK;
2213 }
2214
2215 /* ************************************************** */
2216 /* **** Wishart wavesets ************************ */
2217 /* **** from Trevor and CDP ************************ */
2218 /* **** John ffitch Jan 2001 ************************ */
2219 /* ************************************************** */
2220
wavesetset(CSOUND * csound,BARRI * p)2221 int32_t wavesetset(CSOUND *csound, BARRI *p)
2222 {
2223 if (*p->len == FL(0.0))
2224 p->length = 1 + (int32_t)(p->h.insdshead->p3.value * CS_ESR * FL(0.5));
2225 else
2226 p->length = 1 + (int32_t)*p->len;
2227 if (UNLIKELY(p->length <= 1)) p->length = (int32_t)CS_ESR;
2228 csound->AuxAlloc(csound, (int32)p->length*sizeof(MYFLT), &p->auxch);
2229 p->cnt = 1;
2230 p->start = 0;
2231 p->current = 0;
2232 p->end = 0;
2233 p->direction = 1;
2234 p->lastsamp = FL(1.0);
2235 p->noinsert = 0;
2236 return OK;
2237 }
2238
waveset(CSOUND * csound,BARRI * p)2239 int32_t waveset(CSOUND *csound, BARRI *p)
2240 {
2241 IGN(csound);
2242 MYFLT *in = p->ain;
2243 MYFLT *out = p->ar;
2244 int32_t index = p->end;
2245 MYFLT *insert = (MYFLT*)(p->auxch.auxp) + index;
2246 uint32_t offset = p->h.insdshead->ksmps_offset;
2247 uint32_t early = p->h.insdshead->ksmps_no_end;
2248 uint32_t n, nsmps = CS_KSMPS;
2249 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
2250 if (UNLIKELY(early)) {
2251 nsmps -= early;
2252 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
2253 }
2254 if (p->noinsert) goto output;
2255 for (n=offset;n<nsmps;n++) { /* Deal with inputs */
2256 *insert++ = in[n];
2257 if (++index == p->start) {
2258 p->noinsert = 1;
2259 break;
2260 }
2261 if (index==p->length) { /* Input wrapping */
2262 index = 0;
2263 insert = (MYFLT*)(p->auxch.auxp);
2264 }
2265 }
2266 output:
2267 p->end = index;
2268 index = p->current;
2269 insert = (MYFLT*)(p->auxch.auxp) + index;
2270 for (n=offset;n<nsmps;n++) {
2271 MYFLT samp = *insert++;
2272 index ++;
2273 if (index==p->length) {
2274 index = 0;
2275 insert = (MYFLT*)(p->auxch.auxp);
2276 }
2277 if (samp != FL(0.0) && p->lastsamp*samp < FL(0.0)) {
2278 if (p->direction == 1)
2279 p->direction = -1; /* First cross */
2280 else { /* Second cross */
2281 p->direction = 1;
2282 if (++p->cnt > *p->rep) {
2283 p->cnt = 1;
2284 p->start = index;
2285 p->noinsert = 0;
2286 }
2287 else {
2288 index = p->start;
2289 insert = (MYFLT*)(p->auxch.auxp) + index;
2290 }
2291 }
2292 }
2293 if (samp != FL(0.0)) p->lastsamp = samp;
2294 out[n] = samp;
2295 }
2296 p->current = index;
2297 return OK;
2298 }
2299
medfiltset(CSOUND * csound,MEDFILT * p)2300 int32_t medfiltset(CSOUND *csound, MEDFILT *p)
2301 {
2302 int32_t maxwind = (int32_t)MYFLT2LONG(*p->imaxsize);
2303 int32_t auxsize = 2*sizeof(MYFLT)*maxwind;
2304 p->ind = 0;
2305 p->maxwind = maxwind;
2306
2307 if (p->b.auxp==NULL || p->b.size < (size_t)auxsize)
2308 csound->AuxAlloc(csound, (size_t)auxsize, &p->b);
2309 else
2310 if (*p->iskip!=FL(0.0)) memset(p->b.auxp, 0, auxsize);
2311 p->buff = (MYFLT*)p->b.auxp;
2312 p->med = &(p->buff[maxwind]);
2313 return OK;
2314 }
2315
medfilt(CSOUND * csound,MEDFILT * p)2316 int32_t medfilt(CSOUND *csound, MEDFILT *p)
2317 {
2318 MYFLT *aout = p->ans;
2319 MYFLT *asig = p->asig;
2320 MYFLT *buffer = p->buff;
2321 MYFLT *med = p->med;
2322 int32_t maxwind = p->maxwind;
2323 int32_t kwind = MYFLT2LONG(*p->kwind);
2324 int32_t index = p->ind;
2325 uint32_t offset = p->h.insdshead->ksmps_offset;
2326 uint32_t early = p->h.insdshead->ksmps_no_end;
2327 uint32_t n, nsmps = CS_KSMPS;
2328 if (UNLIKELY(p->b.auxp==NULL)) {
2329 return csound->PerfError(csound, &(p->h),
2330 Str("median: not initialised (arate)\n"));
2331 }
2332 if (UNLIKELY(kwind > maxwind)) {
2333 csound->Warning(csound,
2334 Str("median: window (%d)larger than maximum(%d); truncated"),
2335 kwind, maxwind);
2336 kwind = maxwind;
2337 }
2338 if (UNLIKELY(offset)) memset(aout, '\0', offset*sizeof(MYFLT));
2339 if (UNLIKELY(early)) {
2340 nsmps -= early;
2341 memset(&aout[nsmps], '\0', early*sizeof(MYFLT));
2342 }
2343 for (n=offset; n<nsmps; n++) {
2344 MYFLT x = asig[n];
2345 buffer[index++] = x;
2346
2347 if (kwind<=index) { /* all in centre */
2348 memcpy(&med[0], &buffer[index-kwind], kwind*sizeof(MYFLT));
2349 /* printf("memcpy: 0 <- %d (%d)\n", index-kwind, kwind); */
2350 }
2351 else { /* or in two parts */
2352 memcpy(&med[0], &buffer[0], index*sizeof(MYFLT));
2353 /* printf("memcpy: 0 <- 0 (%d)\n", index); */
2354 memcpy(&med[index], &buffer[maxwind+index-kwind],
2355 (kwind-index)*sizeof(MYFLT));
2356 /* printf("memcpy: %d <- %d (%d)\n",
2357 index, maxwind+index-kwind, kwind-index); */
2358 }
2359 /* { int32_t i; */
2360 /* for (i=0; i<8; i++) printf(" %f", buffer[i]); */
2361 /* printf("\n\t:"); */
2362 /* for (i=0; i<5; i++) printf(" %f", med[i]); */
2363 /* printf("\n"); */
2364 /* } */
2365 aout[n] = medianvalue(kwind, med-1); /* -1 as should point below data */
2366 /* printf("%d/$%d: %f -> %f\n", n, index-1, x, aout[n]); */
2367 if (index>=maxwind) index = 0;
2368 }
2369 p->ind = index;
2370 return OK;
2371 }
2372
kmedfilt(CSOUND * csound,MEDFILT * p)2373 int32_t kmedfilt(CSOUND *csound, MEDFILT *p)
2374 {
2375 MYFLT *buffer = p->buff;
2376 MYFLT *med = p->med;
2377 MYFLT x = *p->asig;
2378 int32_t maxwind = p->maxwind;
2379 int32_t kwind = MYFLT2LONG(*p->kwind);
2380 int32_t index = p->ind;
2381 if (UNLIKELY(p->b.auxp==NULL)) {
2382 return csound->PerfError(csound, &(p->h),
2383 Str("median: not initialised (krate)\n"));
2384 }
2385 if (UNLIKELY(kwind > maxwind)) {
2386 csound->Warning(csound,
2387 Str("median: window (%d)larger than maximum(%d); truncated"),
2388 kwind, maxwind);
2389 kwind = maxwind;
2390 }
2391 buffer[index++] = x;
2392 if (kwind<=index) { /* all in centre */
2393 memcpy(&med[0], &buffer[index-kwind], kwind*sizeof(MYFLT));
2394 }
2395 else { /* or in two parts */
2396 memcpy(&med[0], &buffer[0], index*sizeof(MYFLT));
2397 memcpy(&med[index], &buffer[maxwind+index-kwind],
2398 (kwind-index)*sizeof(MYFLT));
2399 }
2400 *p->ans = medianvalue(kwind, med-1); /* -1 as should point below data */
2401 if (index>=maxwind) index = 0;
2402 p->ind = index;
2403 return OK;
2404 }
2405