1 /*
2 disprep.c:
3
4 Copyright (C) 1991 Barry Vercoe
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" /* DISPREP.C */
25 #include <math.h>
26 #include "cwindow.h"
27 #include "disprep.h"
28
29
30 #ifdef MSVC /* Thanks to Richard Dobson */
31 # define hypot _hypot
32 #endif
33
printv(CSOUND * csound,PRINTV * p)34 int32_t printv(CSOUND *csound, PRINTV *p)
35 {
36 int32_t nargs = p->INOCOUNT;
37 char **txtp = p->h.optext->t.inlist->arg;
38 MYFLT **valp = p->iargs;
39
40 csound->MessageS(csound, CSOUNDMSG_ORCH,
41 "instr %d:", (int32_t) p->h.insdshead->p1.value);
42 while (nargs--) {
43 csound->MessageS(csound, CSOUNDMSG_ORCH,
44 " %s = %5.3f", *txtp++, **valp++);
45 }
46 csound->MessageS(csound, CSOUNDMSG_ORCH, "\n");
47 return OK;
48 }
49
fdspset(CSOUND * csound,FSIGDISP * p)50 int32_t fdspset(CSOUND *csound, FSIGDISP *p){
51 char strmsg[256];
52 p->size = p->fin->N/2 + 1;
53 if ((*p->points != (MYFLT) 0) && (p->size > (int32_t) *p->points)) {
54 p->size = *p->points;
55 }
56 if ((p->fdata.auxp == NULL) ||
57 (p->fdata.size < (uint32_t) (p->size*sizeof(MYFLT)))) {
58 csound->AuxAlloc(csound, p->size*sizeof(MYFLT), &p->fdata);
59 }
60 snprintf(strmsg, 256, Str("instr %d, pvs-signal %s:"),
61 (int32_t) p->h.insdshead->p1.value, p->h.optext->t.inlist->arg[0]);
62 dispset(csound, &p->dwindow, (MYFLT*) p->fdata.auxp, p->size, strmsg,
63 (int32_t) *p->flag, Str("display"));
64 p->lastframe = 0;
65 return OK;
66
67 }
68
fdsplay(CSOUND * csound,FSIGDISP * p)69 int32_t fdsplay(CSOUND *csound, FSIGDISP *p)
70 {
71 float *fin = p->fin->frame.auxp;
72 MYFLT *pdata = p->fdata.auxp;
73 int32_t i,k, end = p->size;
74
75 if (p->lastframe < p->fin->framecount) {
76 for (i=0, k=0; k < end; i+=2,k++) pdata[k] = fin[i];
77 display(csound, &p->dwindow);
78 p->lastframe = p->fin->framecount;
79 }
80 return OK;
81 }
82
dspset(CSOUND * csound,DSPLAY * p)83 int32_t dspset(CSOUND *csound, DSPLAY *p)
84 {
85 int32_t npts, nprds, bufpts, totpts;
86 char *auxp;
87 char strmsg[256];
88
89 if (p->h.optext->t.intype == 'k')
90 npts = (int32_t)(*p->iprd * CS_EKR);
91 else npts = (int32_t)(*p->iprd * csound->esr);
92 if (UNLIKELY(npts <= 0)) {
93 return csound->InitError(csound, Str("illegal iprd in display"));
94
95 }
96 if ((nprds = (int32_t)*p->inprds) <= 1) {
97 nprds = 0;
98 bufpts = npts;
99 totpts = npts;
100 }
101 else {
102 bufpts = npts * nprds;
103 totpts = bufpts * 2;
104 }
105 if ((auxp = p->auxch.auxp) == NULL || totpts != p->totpts) {
106 csound->AuxAlloc(csound, totpts * sizeof(MYFLT), &p->auxch);
107 auxp = p->auxch.auxp;
108 p->begp = (MYFLT *) auxp;
109 p->endp = p->begp + bufpts;
110 p->npts = npts;
111 p->nprds = nprds;
112 p->bufpts = bufpts;
113 p->totpts = totpts;
114 }
115 p->nxtp = (MYFLT *) auxp;
116 p->pntcnt = npts;
117 snprintf(strmsg, 256, Str("instr %d, signal %s:"),
118 (int32_t) p->h.insdshead->p1.value, p->h.optext->t.inlist->arg[0]);
119 dispset(csound, &p->dwindow, (MYFLT*) auxp, bufpts, strmsg,
120 (int32_t) *p->iwtflg, Str("display"));
121 return OK;
122 }
123
kdsplay(CSOUND * csound,DSPLAY * p)124 int32_t kdsplay(CSOUND *csound, DSPLAY *p)
125 {
126 MYFLT *fp = p->nxtp;
127
128 if (UNLIKELY(p->auxch.auxp==NULL)) goto err1; /* RWD fix */
129 if (!p->nprds) {
130 *fp++ = *p->signal;
131 if (fp >= p->endp) {
132 fp = p->begp;
133 display(csound, &p->dwindow);
134 }
135 }
136 else {
137 MYFLT *fp2 = fp + p->bufpts;
138 *fp++ = *p->signal;
139 *fp2 = *p->signal;
140 if (!(--p->pntcnt)) {
141 p->pntcnt = p->npts;
142 if (fp >= p->endp) {
143 fp = p->begp;
144 fp2 = fp + p->bufpts;
145 }
146 p->dwindow.fdata = fp; /* display from fp */
147 display(csound, &p->dwindow);
148 }
149 }
150 p->nxtp = fp;
151 return OK;
152 err1:
153 return csound->PerfError(csound, &(p->h),
154 Str("display: not initialised"));
155 }
156
dsplay(CSOUND * csound,DSPLAY * p)157 int32_t dsplay(CSOUND *csound, DSPLAY *p)
158 {
159 MYFLT *fp = p->nxtp, *sp = p->signal, *endp = p->endp;
160 uint32_t offset = p->h.insdshead->ksmps_offset;
161 uint32_t early = p->h.insdshead->ksmps_no_end;
162 uint32_t n, nsmps = CS_KSMPS;
163
164 nsmps -= early;
165 if (!p->nprds) {
166 for (n=offset; n<nsmps; n++) {
167 *fp++ = sp[n];
168 if (fp >= endp) {
169 fp = p->begp;
170 display(csound, &p->dwindow);
171 }
172 }
173 }
174 else {
175 MYFLT *fp2 = fp + p->bufpts;
176 for (n=offset; n<nsmps; n++) {
177 *fp++ = sp[n];
178 *fp2++ = sp[n];
179 if (!(--p->pntcnt)) {
180 p->pntcnt = p->npts;
181 if (fp >= endp) {
182 fp = p->begp;
183 fp2 = fp + p->bufpts;
184 }
185 p->dwindow.fdata = fp; /* display from fp */
186 display(csound, &p->dwindow);
187 }
188 }
189 }
190 p->nxtp = fp;
191 return OK;
192 }
193
194 /* Write window coefs into buffer, don't malloc */
FillHalfWin(MYFLT * wBuf,int32_t size,MYFLT max,int32_t hannq)195 static void FillHalfWin(MYFLT *wBuf, int32_t size, MYFLT max, int32_t hannq)
196 /* 1 => hanning window else hamming */
197 {
198 MYFLT a,b;
199 int32_t i;
200
201 if (hannq) {
202 a = FL(0.50);
203 b = FL(0.50);
204 }
205 else {
206 a = FL(0.54);
207 b = FL(0.46);
208 }
209
210 if (wBuf!= NULL) { /* NB: size/2 + 1 long - just indep terms */
211 MYFLT tmp;
212 size /= 2; /* to fix scaling */
213 tmp = PI_F/(MYFLT)size; /* optimisatuon?? */
214 for (i=0; i<=size;++i)
215 wBuf[i] = max * (a-b*COS(tmp*(MYFLT)i) );
216 }
217 return;
218 }
219
ApplyHalfWin(MYFLT * buf,MYFLT * win,int32_t len)220 static void ApplyHalfWin(MYFLT *buf, MYFLT *win, int32_t len)
221 { /* Window only store 1st half, is symmetric */
222 int32_t j;
223 int32_t lenOn2 = (len/2);
224
225 for (j=lenOn2+1; j--; )
226 *buf++ *= *win++;
227 for (j =len-lenOn2-1, win--; j--; )
228 *buf++ *= *--win;
229 }
230
fftset(CSOUND * csound,DSPFFT * p)231 int32_t fftset(CSOUND *csound, DSPFFT *p) /* fftset, dspfft -- calc Fast Fourier */
232 /* Transform of collected samples and */
233 /* displays coefficients (mag or db) */
234 {
235 int32_t window_size, step_size;
236 int32_t hanning;
237 char strmsg[256];
238 int32_t minbin, maxbin;
239 minbin = *p->imin;
240 maxbin = *p->imax;
241
242 if (p->smpbuf.auxp == NULL)
243 csound->AuxAlloc(csound, sizeof(MYFLT)*WINDMAX, &(p->smpbuf));
244
245 p->sampbuf = (MYFLT *) p->smpbuf.auxp;
246
247 window_size = (int32_t)*p->inpts;
248 if (UNLIKELY(window_size > WINDMAX)) {
249 return csound->InitError(csound, Str("too many points requested (%d)"), window_size);
250 }
251 if (UNLIKELY(window_size < WINDMIN)) {
252 return csound->InitError(csound, Str("too few points requested (%d), minimum is %d"),
253 window_size, WINDMIN);
254 }
255 if (UNLIKELY(window_size < 1L || (window_size & (window_size - 1L)) != 0L)) {
256 return csound->InitError(csound, Str("window size must be power of two"));
257 }
258 if (p->h.optext->t.intype == 'k')
259 step_size = (int32_t)(*p->iprd * CS_EKR);
260 else step_size = (int32_t)(*p->iprd * csound->esr);
261 if (UNLIKELY(step_size <= 0)) {
262 return csound->InitError(csound, Str("illegal iprd in ffy display"));
263 }
264 hanning = (int32_t)*p->ihann;
265 p->dbout = (int32_t)*p->idbout;
266 p->overlap = window_size - step_size;
267
268
269
270 if ( (maxbin - minbin) != p->npts ||
271 minbin != p->start ||
272 window_size != p->windsize ||
273 hanning != p->hanning) { /* if windowing has changed: */
274 int32_t auxsiz;
275 MYFLT *hWin;
276 p->windsize = window_size; /* set new parameter values */
277 p->hanning = hanning;
278 p->bufp = p->sampbuf;
279 p->endp = p->bufp + window_size;
280 p->overN = FL(1.0)/(*p->inpts);
281 p->ncoefs = window_size >>1;
282 auxsiz = (window_size/2 + 1) * sizeof(MYFLT); /* size for half window */
283 csound->AuxAlloc(csound, (int32_t)auxsiz, &p->auxch); /* alloc or realloc */
284 hWin = (MYFLT *) p->auxch.auxp;
285 FillHalfWin(hWin, window_size,
286 FL(1.0), hanning); /* fill with proper values */
287 if (csound->disprep_fftcoefs == NULL) {
288 /* room for WINDMAX*2 floats (fft size) */
289 csound->disprep_fftcoefs = (MYFLT*) csound->Malloc(csound, WINDMAX * 2
290 * sizeof(MYFLT));
291 }
292 snprintf(strmsg, 256, Str("instr %d, signal %s, fft (%s):"),
293 (int32_t) p->h.insdshead->p1.value, p->h.optext->t.inlist->arg[0],
294 p->dbout ? Str("db") : Str("mag"));
295 if (maxbin == 0) maxbin = p->ncoefs;
296 if (minbin > maxbin) minbin = 0;
297 p->npts = maxbin - minbin;
298 p->start = minbin;
299 dispset(csound, &p->dwindow,
300 csound->disprep_fftcoefs+p->start, p->npts, strmsg,
301 (int32_t) *p->iwtflg, Str("fft"));
302 }
303
304 return OK;
305 }
306
307 /* pack re,im,re,im into re,re */
PackReals(MYFLT * buffer,int32_t size)308 static void PackReals(MYFLT *buffer, int32_t size)
309 {
310 MYFLT *b2 = buffer;
311
312 ++size;
313 while (--size) {
314 *b2++ = *buffer++;
315 ++buffer;
316 }
317 }
318
319 /* Convert Real & Imaginary spectra into Amplitude & Phase */
Rect2Polar(MYFLT * buffer,int32_t size,MYFLT scal)320 static void Rect2Polar(MYFLT *buffer, int32_t size, MYFLT scal)
321 {
322 int32_t i;
323 MYFLT *real,*imag;
324 MYFLT re,im;
325 MYFLT mag;
326
327 real = buffer;
328 imag = buffer+1;
329 for (i = 0; i < size; i++) {
330 re = real[i+i]*scal;
331 im = imag[i+i]*scal;
332 real[2L*i] = mag = HYPOT(re,im);
333 if (mag == FL(0.0))
334 imag[i+i] = FL(0.0);
335 else
336 imag[i+i] = ATAN2(im,re);
337 }
338 }
339
340 /* packed buffer ie. reals, not complex */
Lin2DB(MYFLT * buffer,int32_t size)341 static void Lin2DB(MYFLT *buffer, int32_t size)
342 {
343 while (size--) {
344 if (*buffer > 0.0)
345 *buffer = /* FL(20.0)*log10 */ FL(8.68589)*LOG(*buffer);
346 buffer++;
347 }
348 }
349
d_fft(CSOUND * csound,MYFLT * sce,MYFLT * dst,int32_t size,MYFLT * hWin,int32_t dbq,MYFLT scal)350 static void d_fft( /* perform an FFT as reqd below */
351 CSOUND *csound,
352 MYFLT *sce, /* input array - pure packed real */
353 MYFLT *dst, /* output array - packed magnitude, only half-length */
354 int32_t size, /* number of points in input */
355 MYFLT *hWin, /* hanning window lookup table */
356 int32_t dbq, MYFLT scal) /* flag: 1-> convert output into db */
357 {
358 memcpy(dst, sce, sizeof(MYFLT) * size); /* copy into scratch buffer */
359 ApplyHalfWin(dst, hWin, size);
360 csound->RealFFT(csound, dst, (int32_t) size); /* perform the FFT */
361 dst[size] = dst[1];
362 dst[1] = dst[size + 1L] = FL(0.0);
363 Rect2Polar(dst, (size >> 1) + 1, scal);
364 PackReals(dst, (size >> 1) + 1);
365 if (dbq)
366 Lin2DB(dst, (size >> 1) + 1);
367 }
368
kdspfft(CSOUND * csound,DSPFFT * p)369 int32_t kdspfft(CSOUND *csound, DSPFFT *p)
370 {
371 MYFLT *bufp = p->bufp, *endp = p->endp;
372
373 if (p->dbout) p->dwindow.polarity = NEGPOL;
374 else p->dwindow.polarity = POSPOL;
375
376 if (UNLIKELY(p->auxch.auxp==NULL)) goto err1; /* RWD fix */
377 if (bufp < p->sampbuf) /* skip any spare samples */
378 bufp++;
379 else { /* then start collecting */
380 *bufp++ = *p->signal;
381 if (bufp >= endp) { /* when full, do fft: */
382 MYFLT *tp;
383 //MYFLT *tplim;
384 MYFLT *hWin = (MYFLT *) p->auxch.auxp;
385 d_fft(csound, p->sampbuf, csound->disprep_fftcoefs,
386 p->windsize, hWin, p->dbout, p->overN);
387 tp = csound->disprep_fftcoefs;
388 //tplim = tp + p->ncoefs;
389 //do {
390 // *tp *= p->overN; /* scale 1/N */
391 //} while (++tp < tplim);
392 display(csound, &p->dwindow); /* & display */
393 if (p->overlap > 0) {
394 bufp = p->sampbuf;
395 tp = endp - p->overlap;
396 do {
397 *bufp++ = *tp++;
398 } while (tp < endp);
399 }
400 else bufp = p->sampbuf + p->overlap;
401 }
402 }
403 p->bufp = bufp;
404 return OK;
405 err1:
406 return csound->PerfError(csound, &(p->h),
407 Str("dispfft: not initialised"));
408 }
409
dspfft(CSOUND * csound,DSPFFT * p)410 int32_t dspfft(CSOUND *csound, DSPFFT *p)
411 {
412 MYFLT *sigp = p->signal, *bufp = p->bufp, *endp = p->endp;
413 uint32_t offset = p->h.insdshead->ksmps_offset;
414 uint32_t early = p->h.insdshead->ksmps_no_end;
415 uint32_t n, nsmps = CS_KSMPS;
416
417 if (p->dbout) {
418 p->dwindow.polarity = NEGPOL;
419 p->dwindow.absflag = 1;
420 }
421 else p->dwindow.polarity = POSPOL;
422
423 if (UNLIKELY(p->auxch.auxp==NULL)) goto err1;
424 nsmps -= early;
425 for (n=offset; n<nsmps; n++) {
426 if (bufp < p->sampbuf) { /* skip any spare samples */
427 bufp++; sigp++;
428 }
429 else { /* then start collecting */
430 *bufp++ = *sigp++;
431 if (bufp >= endp) { /* when full, do fft: */
432 MYFLT *tp;
433 //MYFLT *tplim;
434 MYFLT *hWin = (MYFLT *) p->auxch.auxp;
435 d_fft(csound, p->sampbuf, csound->disprep_fftcoefs,
436 p->windsize, hWin, p->dbout, p->overN);
437 tp = csound->disprep_fftcoefs;
438 //tplim = tp + p->ncoefs;
439 //do {
440 // *tp *= p->overN; /* scale 1/N */
441 //} while (++tp < tplim);
442
443 display(csound, &p->dwindow); /* & display */
444 if (p->overlap > 0) {
445 bufp = p->sampbuf;
446 tp = endp - p->overlap;
447 do {
448 *bufp++ = *tp++;
449 } while (tp < endp);
450 }
451 else bufp = p->sampbuf + p->overlap;
452 }
453 }
454 }
455 p->bufp = bufp;
456 return OK;
457 err1:
458 return csound->PerfError(csound, &(p->h),
459 Str("dispfft: not initialised"));
460 }
461
462 #define NTERMS 4
463 #define NCROSS (NTERMS * (NTERMS-1))
464
tempeset(CSOUND * csound,TEMPEST * p)465 int32_t tempeset(CSOUND *csound, TEMPEST *p)
466 {
467 int32_t npts = 0, nptsm1, minlam = 0, maxlam, lamspan, auxsiz;
468 MYFLT *fltp;
469 FUNC *ftp;
470 MYFLT b, iperiod = *p->iprd;
471 char strmsg[256];
472
473 if (UNLIKELY((p->timcount = (int32_t)(CS_EKR * iperiod)) <= 0))
474 return csound->InitError(csound, Str("illegal iperiod"));
475 if (UNLIKELY((p->dtimcnt = (int32_t)(CS_EKR * *p->idisprd)) < 0))
476 return csound->InitError(csound, Str("illegal idisprd"));
477 if (UNLIKELY((p->tweek = *p->itweek) <= 0))
478 return csound->InitError(csound, Str("illegal itweek"));
479 if (iperiod != FL(0.0)) {
480 if (UNLIKELY((minlam = (int32_t)(*p->imindur/iperiod)) <= 0))
481 return csound->InitError(csound, Str("illegal imindur"));
482 if (UNLIKELY((npts = (int32_t)(*p->imemdur / iperiod)) <= 0))
483 return csound->InitError(csound, Str("illegal imemdur"));
484 }
485 if (UNLIKELY(*p->ihtim <= FL(0.0)))
486 return csound->InitError(csound, Str("illegal ihtim"));
487 if (UNLIKELY(*p->istartempo <= FL(0.0)))
488 return csound->InitError(csound, Str("illegal startempo"));
489 ftp = csound->FTFind(csound, p->ifn);
490 if (UNLIKELY(ftp != NULL && *ftp->ftable == FL(0.0)))
491 return csound->InitError(csound, Str("ifn table begins with zero"));
492 if (UNLIKELY(ftp==NULL)) return NOTOK;
493
494 if (npts==0) return NOTOK;
495
496 nptsm1 = npts - 1;
497 if (npts != p->npts || minlam != p->minlam) {
498 p->npts = npts;
499 p->minlam = minlam;
500 p->maxlam = maxlam = nptsm1/(NTERMS-1);
501 lamspan = maxlam - minlam + 1; /* alloc 8 bufs: 2 circ, 6 lin */
502 auxsiz = (npts * 5 + lamspan * 3) * sizeof(MYFLT);
503 csound->AuxAlloc(csound, (int32_t)auxsiz, &p->auxch);
504 fltp = (MYFLT *) p->auxch.auxp;
505 p->hbeg = fltp; fltp += npts;
506 p->hend = fltp;
507 p->xbeg = fltp; fltp += npts;
508 p->xend = fltp;
509 p->stmemp = fltp; fltp += npts;
510 p->linexp = fltp; fltp += npts;
511 p->ftable = fltp; fltp += npts;
512 p->xscale = fltp; fltp += lamspan;
513 p->lmults = fltp; fltp += lamspan;
514 p->lambdas = (int16 *) fltp;
515 p->stmemnow = p->stmemp + nptsm1;
516 }
517 if (p->dtimcnt && !(p->dwindow.windid)) { /* init to display stmem & exp */
518 snprintf(strmsg, 256, "instr %d tempest:", (int32_t) p->h.insdshead->p1.value);
519 dispset(csound, &p->dwindow, p->stmemp, (int32_t)npts * 2, strmsg, 0,
520 Str("tempest"));
521 p->dwindow.danflag = 1; /* for mid-scale axis */
522 }
523 {
524 MYFLT *funp = ftp->ftable;
525 int32_t phs = 0;
526 int32_t inc = (int32_t)PHMASK / npts;
527 int32_t nn, lobits = ftp->lobits;
528 for (fltp=p->hbeg, nn=npts*4; nn--; ) /* clr 2 circ & 1st 2 lin bufs */
529 *fltp++ = FL(0.0);
530 for (fltp=p->ftable+npts, nn=npts; nn--; ) { /* now sample the ftable */
531 *--fltp = *(funp + (phs >> lobits)); /* backwards into tbl buf */
532 phs += inc;
533 }
534 }
535 {
536 MYFLT *tblp, sumraw, sumsqr; /* calc the CROSS prod scalers */
537 int32_t terms;
538 int32_t lambda, maxlam;
539 MYFLT crossprods, RMS, *endtable = p->ftable + nptsm1;
540 /* MYFLT coef, log001 = -6.9078; */
541 MYFLT *xscale = p->xscale;
542
543 p->ncross = (MYFLT) NCROSS;
544 for (lambda=p->minlam,maxlam=p->maxlam; lambda <= maxlam; lambda++) {
545 tblp = endtable;
546 sumraw = *tblp;
547 sumsqr = *tblp * *tblp;
548 terms = NTERMS - 1;
549 do {
550 tblp -= lambda;
551 sumraw += *tblp;
552 sumsqr += *tblp * *tblp;
553 } while (--terms);
554 crossprods = sumraw * sumraw - sumsqr;
555 RMS = SQRT(crossprods / p->ncross);
556 /* coef = exp(log001 * lambda / npts);
557 *xscale++ = coef / RMS / (NTERMS - 1); */
558 *xscale++ = FL(0.05)/ RMS / lambda;
559 }
560 }
561 /* calc input lo-pass filter coefs */
562 b = FL(2.0) - COS((*p->ihp * 6.28318 * CS_ONEDKR));
563 p->coef1 = b - SQRT(b * b - FL(1.0));
564 p->coef0 = FL(1.0) - p->coef1;
565 p->yt1 = FL(0.0);
566 p->fwdcoef = POWER(FL(0.5), p->timcount*CS_ONEDKR/(*p->ihtim));
567 p->fwdmask = FL(0.0);
568 #ifdef DEBUG
569 csound->Message(csound,
570 Str("kin lopass coef1 %6.4f, fwd mask coef1 %6.4f\n"),
571 p->coef1, p->fwdcoef);
572 #endif
573 p->thresh = *p->ithresh; /* record incoming loudness threshold */
574 p->xfdbak = *p->ixfdbak; /* & expectation feedback fraction */
575 p->tempscal = FL(60.0) * CS_EKR / p->timcount;
576 p->avglam = p->tempscal / *p->istartempo; /* init the tempo factors */
577 p->tempo = FL(0.0);
578 p->hcur = p->hbeg; /* init the circular ptrs */
579 p->xcur = p->xbeg;
580 p->countdown = p->timcount; /* & prime the countdowns */
581 p->dcntdown = p->dtimcnt;
582 return OK;
583 }
584
585 #define NMULTS 5
586
587 static const MYFLT lenmults[NMULTS] = {
588 FL(3.0), FL(2.0), FL(1.0), FL(0.5), FL(0.333)
589 };
590
591 static const MYFLT lenfracs[NMULTS*2] = {
592 FL(0.30), FL(0.3667), FL(0.45), FL(0.55), FL(0.92), FL(1.08),
593 FL(1.88), FL(2.12), FL(2.85), FL(3.15)
594 };
595
tempest(CSOUND * csound,TEMPEST * p)596 int32_t tempest(CSOUND *csound, TEMPEST *p)
597 {
598 p->yt1 = p->coef0 * *p->kin + p->coef1 * p->yt1; /* get lo-pass of kinput */
599
600 if (UNLIKELY(p->auxch.auxp==NULL)) goto err1; /* RWD fix */
601 if (!(--p->countdown)) { /* then on countdown: */
602 MYFLT *memp;
603 MYFLT kin, expect, *xcur = p->xcur; /* xcur from prv pass */
604 MYFLT lamtot = FL(0.0), weightot = FL(0.0);
605
606 p->countdown = p->timcount; /* reset the countdown */
607 expect = *xcur; /* get expected val from prv calc */
608 *xcur++ = FL(0.0); /* & clear the loc it occupied */
609 if (xcur >= p->xend) xcur = p->xbeg; /* xcur now points to cur xarray */
610 p->xcur = xcur;
611 #ifdef DEBUG
612 csound->Message(csound, "**kin -> %f (%f,%f)\n",
613 *p->kin - p->yt1, *p->kin, p->yt1);
614 #endif
615 if ((kin = *p->kin - p->yt1) < FL(0.0))
616 kin = FL(0.0);
617 { /* ignore input below lopass */
618 MYFLT *hcur = p->hcur;
619 MYFLT *hend = p->hend;
620 MYFLT *tblp = p->ftable;
621 int32_t wrap;
622 *hcur++ = kin + expect * p->xfdbak; /* join insample & expect val */
623 if (hcur < hend) p->hcur = hcur; /* stor pntr for next insamp */
624 else p->hcur = p->hbeg;
625 wrap = hcur - p->hbeg;
626 memp = p->stmemp;
627 while (hcur < hend) /* now lineariz & envlp hbuf */
628 *memp++ = *hcur++ * *tblp++; /* into st_mem buf */
629 for (hcur=p->hbeg; wrap--; )
630 *memp++ = *hcur++ * *tblp++;
631 }
632 if (p->yt1 > p->thresh /* if lo-pass of kinput now significant */
633 && kin > p->fwdmask) { /* & kin > masking due to prev kin */
634 MYFLT sumraw, sumsqr;
635 int32_t lambda, minlam, maxlam;
636 int32_t terms, nn, npts = p->npts;
637 MYFLT mult, crossprods, RMScross, RMStot, unilam, rd;
638 MYFLT *xend = p->xend;
639 /* MYFLT *xscale = p->xscale; */
640 const MYFLT *mults, *fracs;
641 MYFLT *mulp;
642 int16 minlen, maxlen, *lenp, *endlens;
643
644 for (memp=p->stmemp,nn=npts,sumsqr=FL(0.0); nn--; memp++)
645 sumsqr += *memp * *memp;
646 RMStot = (MYFLT)sqrt(sumsqr/npts);
647 #ifdef DEBUG
648 csound->Message(csound, "RMStot = %6.1f\n", RMStot);
649 #endif
650 mults = lenmults; /* use the static lentables */
651 fracs = lenfracs;
652 mulp = p->lmults;
653 lenp = p->lambdas;
654 minlam = p->minlam;
655 maxlam = p->maxlam;
656 nn = NMULTS;
657 do {
658 mult = *mults++;
659 minlen = (int16)(p->avglam * *fracs++); /* & the current avglam */
660 maxlen = (int16)(p->avglam * *fracs++);
661 if (minlen >= minlam && maxlen <= maxlam)
662 do {
663 *lenp++ = minlen++; /* creat lst of lambda lens */
664 *mulp++ = mult; /* & their unit multipliers */
665 } while (minlen <= maxlen);
666 } while (--nn);
667 endlens = lenp; /* now for these lambda lens: */
668 for (lenp=p->lambdas,mulp=p->lmults; lenp < endlens; ) {
669 lambda = *lenp++;
670 mult = *mulp++;
671 memp = p->stmemnow;
672 sumraw = *memp;
673 sumsqr = *memp * *memp; /* autocorrelate the st_mem buf */
674 terms = NTERMS - 1;
675 do {
676 memp -= lambda;
677 sumraw += *memp;
678 sumsqr += *memp * *memp;
679 } while (--terms);
680 crossprods = sumraw * sumraw - sumsqr;
681 if (crossprods >= 0)
682 RMScross = SQRT(crossprods / p->ncross);
683 else
684 RMScross = FL(0.0);
685 if (RMScross < FL(1.4) * RMStot) /* if RMScross significant: */
686 continue;
687 #ifdef DEBUG
688 csound->Message(csound, "RMScross = %6.1f, lambda = %ld\n",
689 RMScross, lambda);
690 #endif
691 /* RMS *= *xscale++; */
692 unilam = lambda * mult; /* get unit lambda implied */
693 lamtot += unilam * RMScross; /* & add weighted to total */
694 weightot += RMScross;
695 RMScross /= FL(5.0);
696 memp = xcur - 1; /* multiply project into expect buf */
697 for (terms=1; terms < NTERMS; ++terms) {
698 if ((memp += (lambda-terms+1)) >= xend)
699 memp -= npts;
700 for (nn=terms,rd=RMScross/terms; nn--; ) {
701 *memp++ += rd;
702 if (memp >= xend)
703 memp -= npts;
704 }
705 }
706 }
707 }
708 if (weightot) { /* if accumed weights, */
709 p->avglam = (p->avglam + lamtot/weightot)/FL(2.0); /* update the avglam */
710 p->avglam /= p->tweek;
711 p->tempo = p->tempscal / p->avglam; /* & cvt to tempo */
712 #ifdef DEBUG
713 csound->Message(csound, "lamtot %6.2f, weightot %6.2f, "
714 "newavglam %6.2f, tempo %6.2f\n",
715 lamtot, weightot, p->avglam, p->tempo);
716 #endif
717 }
718 else {
719 if (kin < -(p->fwdmask)) {
720 p->tempo = FL(0.0); /* else tempo is 0 */
721 }
722 }
723 p->fwdmask = p->fwdmask * p->fwdcoef + kin;
724 }
725 if (!(--p->dcntdown)) { /* on display countdown */
726 MYFLT *linp = p->linexp;
727 MYFLT *xcur = p->xcur;
728 MYFLT *xend = p->xend;
729 int32_t wrap = xcur - p->xbeg;
730 while (xcur < xend) /* lineariz the circ xbuf */
731 *linp++ = *xcur++; /* into linexp buf */
732 for (xcur=p->xbeg; wrap--; )
733 *linp++ = *xcur++;
734 display(csound, &p->dwindow); /* display double window */
735 p->dcntdown = p->dtimcnt; /* & reset the counter */
736 }
737 /* if (p->tempo != 0.0) */
738 *p->kout = p->tempo; /* put current tempo */
739 return OK;
740 err1:
741 return csound->PerfError(csound, &(p->h),
742 Str("tempest: not initialised"));
743 }
744