1 /*
2 spectra.c:
3
4 Copyright (C) 1995 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 "csdl.h"
25 #include "csoundCore.h"
26 #include "interlocks.h"
27 #include <math.h>
28 #include "cwindow.h"
29 #include "spectra.h"
30 #include "pitch.h"
31 #include "uggab.h"
32 #include <inttypes.h>
33
34 #define LOGTWO (0.69314718056)
35
DOWNset(CSOUND * p,DOWNDAT * downdp,int32_t npts)36 void DOWNset(CSOUND *p, DOWNDAT *downdp, int32_t npts)
37 {
38 int32_t nbytes = npts * sizeof(MYFLT);
39
40 if (downdp->auxch.auxp == NULL || downdp->auxch.size != (uint32_t)nbytes)
41 p->AuxAlloc(p, nbytes, &downdp->auxch);
42 downdp->npts = npts;
43 }
44
SPECset(CSOUND * p,SPECDAT * specdp,int32_t npts)45 void SPECset(CSOUND *p, SPECDAT *specdp, int32_t npts)
46 {
47 int32_t nbytes = npts * sizeof(MYFLT);
48
49 if (specdp->auxch.auxp == NULL || (uint32_t)nbytes != specdp->auxch.size)
50 p->AuxAlloc(p, nbytes, &specdp->auxch);
51 specdp->npts = npts;
52 }
53
54 static const char *outstring[] = {"mag", "db", "mag sqrd", "root mag"};
55
spectset(CSOUND * csound,SPECTRUM * p)56 int32_t spectset(CSOUND *csound, SPECTRUM *p)
57 /* spectrum - calcs disc Fourier transform of */
58 /* oct-downsampled data outputs coefs (mag, */
59 /* db or mag2) of log freq within each octave */
60 {
61 int32_t n, nocts, nfreqs, ncoefs, hanning;
62 MYFLT Q, *fltp;
63 OCTDAT *octp;
64 DOWNDAT *dwnp = &p->downsig;
65 SPECDAT *specp = p->wsig;
66
67 /* for mac roundoff */
68 p->timcount = (int32_t)(CS_EKR * *p->iprd + FL(0.001));
69 nocts = (int32_t)*p->iocts;
70 nfreqs = (int32_t)*p->ifrqs;
71 ncoefs = nocts * nfreqs;
72 Q = *p->iq;
73 hanning = (*p->ihann) ? 1 : 0;
74 p->dbout = (int32_t)*p->idbout;
75 if (UNLIKELY((p->disprd = (int32_t)(CS_EKR * *p->idisprd)) < 0))
76 p->disprd = 0;
77
78 if (UNLIKELY(p->timcount <= 0))
79 return csound->InitError(csound, Str("illegal iprd"));
80 if (UNLIKELY(nocts <= 0 || nocts > MAXOCTS))
81 return csound->InitError(csound, Str("illegal iocts"));
82 if (UNLIKELY(nfreqs <= 0 || nfreqs > MAXFRQS))
83 return csound->InitError(csound, Str("illegal ifrqs"));
84 if (UNLIKELY(Q <= FL(0.0)))
85 return csound->InitError(csound, Str("illegal Q value"));
86 if (UNLIKELY(p->dbout < 0 || p->dbout > 3))
87 return csound->InitError(csound, Str("unknown dbout code"));
88
89 if (nocts != dwnp->nocts ||
90 nfreqs != p->nfreqs || /* if anything has changed */
91 Q != p->curq ||
92 (p->disprd && !p->octwindow.windid) ||
93 hanning != p->hanning) { /* make new tables */
94 double basfrq, curfrq, frqmlt, Qfactor;
95 double theta, a, windamp, onedws, pidws;
96 MYFLT *sinp, *cosp;
97 int32_t k, sumk, windsiz, halfsiz, *wsizp, *woffp;
98 int32_t auxsiz, bufsiz = 0;
99 int32_t majr, minr, totsamps, totsize;
100 double hicps,locps,oct; /* must alloc anew */
101
102 p->nfreqs = nfreqs;
103 p->curq = Q;
104 p->hanning = hanning;
105 p->ncoefs = ncoefs;
106 csound->Warning(csound,
107 Str("spectrum: %s window, %s out, making tables ...\n"),
108 (hanning) ? "hanning":"hamming", outstring[p->dbout]);
109
110 if (p->h.optext->t.intype == 'k') {
111 dwnp->srate = CS_EKR; /* define the srate */
112 p->nsmps = 1;
113 }
114 else {
115 dwnp->srate = CS_ESR;
116 p->nsmps = CS_KSMPS;
117 }
118 hicps = dwnp->srate * 0.375; /* top freq is 3/4 pi/2 ... */
119 oct = log(hicps / ONEPT) / LOGTWO; /* octcps() (see aops.c) */
120 if (p->h.optext->t.intype != 'k') { /* for sr sampling: */
121 oct = ((int32_t)(oct*12.0 + 0.5)) / 12.0; /* semitone round to A440 */
122 hicps = pow(2.0, oct) * ONEPT; /* cpsoct() */
123 }
124 dwnp->looct = (MYFLT)(oct - nocts); /* true oct val of lowest frq */
125 locps = hicps / (1L << nocts);
126 csound->Warning(csound, Str("\thigh cps %7.1f\n\t low cps %7.1f\n"),
127 hicps, locps);
128
129 basfrq = hicps/2.0; /* oct below retuned top */
130 frqmlt = pow(2.0,(double)1.0/nfreqs); /* nfreq interval mult */
131 Qfactor = Q * dwnp->srate;
132 curfrq = basfrq;
133 for (sumk=0,wsizp=p->winlen,woffp=p->offset,n=nfreqs; n--; ) {
134 *wsizp++ = k = (int32_t)(Qfactor/curfrq) | 01; /* calc odd wind sizes */
135 *woffp++ = (*(p->winlen) - k) / 2; /* & symmetric offsets */
136 sumk += k; /* and find total */
137 curfrq *= frqmlt;
138 }
139 windsiz = *(p->winlen);
140 csound->Warning(csound,
141 Str("\tQ %4.1f uses a %d sample window each octdown\n"),
142 Q, windsiz);
143 auxsiz = (windsiz + 2*sumk) * sizeof(MYFLT); /* calc lcl space rqd */
144
145 csound->AuxAlloc(csound, (size_t)auxsiz, &p->auxch1); /* & alloc auxspace */
146
147 fltp = (MYFLT *) p->auxch1.auxp;
148 p->linbufp = fltp; fltp += windsiz; /* linbuf must take nsamps */
149 p->sinp = sinp = fltp; fltp += sumk;
150 p->cosp = cosp = fltp; /* cos gets rem sumk */
151 wsizp = p->winlen;
152 curfrq = basfrq * TWOPI / dwnp->srate;
153 for (n = nfreqs; n--; ) { /* now fill tables */
154 windsiz = *wsizp++; /* (odd win size) */
155 halfsiz = windsiz >> 1;
156 onedws = 1.0 / (windsiz-1);
157 pidws = PI / (windsiz-1);
158 for (k = -halfsiz; k<=halfsiz; k++) { /* with sines */
159 a = cos(k * pidws);
160 windamp = a * a; /* times hanning */
161 if (!hanning)
162 windamp = 0.08 + 0.92 * windamp; /* or hamming */
163 windamp *= onedws; /* scaled */
164 theta = k * curfrq;
165 *sinp++ = (MYFLT)(windamp * sin(theta));
166 *cosp++ = (MYFLT)(windamp * cos(theta));
167 }
168 curfrq *= frqmlt; /* step by log freq */
169 }
170 if (*p->idsines != FL(0.0)) { /* if reqd, dsply windowed sines now! */
171 csound->dispset(csound, &p->sinwindow, p->sinp, (int32_t) sumk,
172 Str("spectrum windowed sines:"), 0, "spectrum");
173 csound->display(csound, &p->sinwindow);
174 }
175
176 dwnp->hifrq = (MYFLT)hicps;
177 dwnp->lofrq = (MYFLT)locps;
178 dwnp->nsamps = windsiz = *(p->winlen);
179 dwnp->nocts = nocts;
180 minr = windsiz >> 1; /* sep odd windsiz into maj, min */
181 majr = windsiz - minr; /* & calc totsamps reqd */
182 totsamps = (majr*nocts) + (minr<<nocts) - minr;
183 DOWNset(csound, dwnp, totsamps); /* auxalloc in DOWNDAT struct */
184 fltp = (MYFLT *) dwnp->auxch.auxp; /* & distrib to octdata */
185 for (n=nocts,octp=dwnp->octdata+(nocts-1); n--; octp--) {
186 bufsiz = majr + minr;
187 octp->begp = fltp; fltp += bufsiz; /* (lo oct first) */
188 octp->endp = fltp; minr *= 2;
189 }
190 csound->Warning(csound, Str("\t%d oct analysis window "
191 "delay = %"PRIi32" samples (%d msecs)\n"),
192 nocts, bufsiz, (int32_t)(bufsiz*1000/dwnp->srate));
193 if (p->disprd) { /* if display requested, */
194 totsize = totsamps * sizeof(MYFLT); /* alloc an equiv local */
195 csound->AuxAlloc(csound,
196 (size_t)totsize, &p->auxch2);/* linear output window */
197 csound->dispset(csound, &p->octwindow, (MYFLT *)p->auxch2.auxp,
198 (int32_t)totsamps, Str("octdown buffers:"), 0, "spectrum");
199 }
200 SPECset(csound, specp, (int32_t)ncoefs); /* prep the spec dspace */
201 specp->downsrcp = dwnp; /* & record its source */
202 }
203 for (octp=dwnp->octdata; nocts--; octp++) { /* reset all oct params, & */
204 octp->curp = octp->begp;
205 for (fltp=octp->feedback,n=6; n--; )
206 *fltp++ = FL(0.0);
207 octp->scount = 0;
208 }
209 specp->nfreqs = p->nfreqs; /* save the spec descriptors */
210 specp->dbout = p->dbout;
211 specp->ktimstamp = 0; /* init specdata to not new */
212 specp->ktimprd = p->timcount;
213 p->scountdown = p->timcount; /* prime the spect countdown */
214 p->dcountdown = p->disprd; /* & the display countdown */
215 return OK;
216 }
217
linocts(DOWNDAT * dwnp,MYFLT * bufp)218 static void linocts(DOWNDAT *dwnp, MYFLT *bufp)
219 /* linearize octdown dat to 1 buf */
220 { /* presumes correct buffer alloc'd in set */
221 MYFLT *curp, *endp;
222 int32_t wrap;
223 OCTDAT *octp;
224 int32_t nocts;
225 MYFLT *begp;
226
227 nocts = dwnp->nocts;
228 octp = dwnp->octdata + nocts;
229 while (nocts--) {
230 octp--; /* for each octave (low to high) */
231 begp = octp->begp;
232 curp = octp->curp;
233 endp = octp->endp;
234 wrap = curp - begp;
235 while (curp < endp) /* copy circbuf to linbuf */
236 *bufp++ = *curp++;
237 for (curp=begp; wrap--; )
238 *bufp++ = *curp++;
239 }
240 }
241
242 static const MYFLT bicoefs[] = {
243 -FL(0.2674054), FL(0.7491305), FL(0.7160484), FL(0.0496285), FL(0.7160484),
244 FL(0.0505247), FL(0.3514850), FL(0.5257536), FL(0.3505025), FL(0.5257536),
245 FL(0.3661840), FL(0.0837990), FL(0.3867783), FL(0.6764264), FL(0.3867783)
246 };
247
spectrum(CSOUND * csound,SPECTRUM * p)248 int32_t spectrum(CSOUND *csound, SPECTRUM *p)
249 {
250 MYFLT a, b, *dftp, *sigp = p->signal, SIG, yt1, yt2;
251 int32_t nocts, nsmps = p->nsmps, winlen;
252 uint32_t offset = p->h.insdshead->ksmps_offset;
253 uint32_t early = p->h.insdshead->ksmps_no_end;
254 DOWNDAT *downp = &p->downsig;
255 OCTDAT *octp;
256 SPECDAT *specp;
257 double c;
258
259 if (UNLIKELY(early)) nsmps -= early;
260 do {
261 SIG = *sigp++; /* for each source sample: */
262 if (offset--) SIG = FL(0.0); /* for sample accuracy */
263 octp = downp->octdata; /* align onto top octave */
264 nocts = downp->nocts;
265 do { /* then for each oct: */
266 const MYFLT *coefp;
267 MYFLT *ytp, *curp;
268 int32_t nfilt;
269 curp = octp->curp;
270 *curp++ = SIG; /* write samp to cur buf */
271 if (curp >= octp->endp)
272 curp = octp->begp; /* & modulo the pointer */
273 octp->curp = curp;
274 if (!(--nocts)) break; /* if lastoct, break */
275 coefp = bicoefs; ytp = octp->feedback;
276 for (nfilt = 3; nfilt--; ) { /* apply triple biquad: */
277 yt2 = *ytp++; yt1 = *ytp--; /* get prev feedback */
278 SIG -= (*coefp++ * yt1); /* apply recurs filt */
279 SIG -= (*coefp++ * yt2);
280 *ytp++ = yt1; *ytp++ = SIG; /* stor nxt feedback */
281 SIG *= *coefp++;
282 SIG += (*coefp++ * yt1); /* apply forwrd filt */
283 SIG += (*coefp++ * yt2);
284 }
285 } while (!(++octp->scount & 01) && octp++); /* send alt samps to nxtoct */
286 } while (--nsmps);
287
288 if (p->disprd) /* if displays requested, */
289 if (!(--p->dcountdown)) { /* on countdown */
290 linocts(downp, (MYFLT *)p->auxch2.auxp); /* linearize the oct bufs */
291 csound->display(csound, &p->octwindow); /* & display */
292 p->dcountdown = p->disprd;
293 }
294
295 if ((--p->scountdown)) return OK;/* if not yet time for new spec, return */
296 p->scountdown = p->timcount; /* else reset counter & proceed: */
297 downp = &p->downsig;
298 specp = p->wsig;
299 nocts = downp->nocts;
300 octp = downp->octdata + nocts;
301 dftp = (MYFLT *) specp->auxch.auxp;
302 winlen = *(p->winlen);
303 while (nocts--) {
304 MYFLT *bufp, *sinp, *cosp;
305 int32_t len, *lenp, *offp, nfreqs;
306 MYFLT *begp, *curp, *endp, *linbufp;
307 int32_t len2;
308 octp--; /* for each oct (low to high) */
309 begp = octp->begp;
310 curp = octp->curp;
311 endp = octp->endp;
312 if ((len = endp - curp) >= winlen) /* if no wrap */
313 linbufp = curp; /* use samples in circbuf */
314 else {
315 len2 = winlen - len;
316 linbufp = bufp = p->linbufp; /* else cp crcbuf to linbuf */
317 while (len--)
318 *bufp++ = *curp++;
319 curp = begp;
320 while (len2--)
321 *bufp++ = *curp++;
322 }
323 cosp = p->cosp; /* get start windowed sines */
324 sinp = p->sinp;
325 lenp = p->winlen;
326 offp = p->offset;
327 for (nfreqs=p->nfreqs; nfreqs--; ) { /* now for ea. frq this oct */
328 a = FL(0.0);
329 b = FL(0.0);
330 bufp = linbufp + *offp++;
331 for (len = *lenp++; len--; bufp++) { /* apply windowed sine seg */
332 a += *bufp * *cosp++;
333 b += *bufp * *sinp++;
334 }
335 c = a*a + b*b; /* get magnitude squared */
336 switch (p->dbout) {
337 case 1:
338 if (c < .001) c = .001; /* and convert to db */
339 c = 10.0 * log10(c);
340 break;
341 case 3:
342 c = sqrt(c); /* or root mag */
343 /* FALLTHRU */
344 case 0:
345 c = sqrt(c); /* or mag */
346 /* FALLTHRU */
347 case 2:
348 break; /* or leave mag sqrd */
349 }
350 *dftp++ = (MYFLT)c; /* store in out spectrum */
351 }
352 }
353 specp->ktimstamp = CS_KCNT; /* time-stamp the output */
354 return OK;
355 }
356
357 /* int32_t nocdfset(CSOUND *csound, NOCTDFT *p) */
358 /* /\* noctdft - calcs disc Fourier transform of oct-downsampled data *\/ */
359 /* /\* outputs coefs (mag, db or mag2) of log freq within each octave *\/ */
360 /* { */
361 /* int32_t nfreqs, hanning, nocts, ncoefs; */
362 /* MYFLT Q, *fltp; */
363 /* DOWNDAT *downp = p->dsig; */
364 /* SPECDAT *specp = p->wsig; */
365
366 /* p->timcount = CS_EKR * *p->iprd; */
367 /* nfreqs = *p->ifrqs; */
368 /* Q = *p->iq; */
369 /* hanning = (*p->ihann) ? 1 : 0; */
370 /* if ((p->dbout = *p->idbout) && p->dbout != 1 && p->dbout != 2) { */
371 /* return csound->InitError(csound,
372 Str("noctdft: unknown dbout code of %d"), */
373 /* p->dbout); */
374 /* } */
375 /* nocts = downp->nocts; */
376 /* ncoefs = nocts * nfreqs; */
377 /* if (nfreqs != p->nfreqs || Q != p->curq /\* if anything changed *\/ */
378 /* || p->timcount <= 0 || Q <= 0. */
379 /* || hanning != p->hanning */
380 /* || ncoefs != p->ncoefs) { /\* make new tables *\/ */
381 /* double basfrq, curfrq, frqmlt, Qfactor; */
382 /* double theta, a, windamp, onedws, pidws; */
383 /* MYFLT *sinp, *cosp; */
384 /* int32_t n, k, sumk, windsiz, *wsizp, nsamps; */
385 /* int64_t auxsiz; */
386
387 /* csound->Message(csound, */
388 /* Str("noctdft: %s window, %s out, making tables ...\n"), */
389 /* (hanning) ? "hanning":"hamming", outstring[p->dbout]); */
390 /* if (p->timcount <= 0) */
391 /* return csound->InitError(csound, Str("illegal iprd")); */
392 /* if (nfreqs <= 0 || nfreqs > MAXFRQS) */
393 /* return csound->InitError(csound, Str("illegal ifrqs")); */
394 /* if (Q <= FL(0.0)) */
395 /* return csound->InitError(csound, Str("illegal Q value")); */
396 /* nsamps = downp->nsamps; */
397 /* p->nfreqs = nfreqs; */
398 /* p->curq = Q; */
399 /* p->hanning = hanning; */
400 /* p->ncoefs = ncoefs; */
401 /* basfrq = downp->hifrq/2.0 * TWOPI/downp->srate;
402 /\* oct below retuned top *\/ */
403 /* frqmlt = pow(2.0,1.0/(double)nfreqs); /\* nfreq interval mult *\/ */
404 /* Qfactor = TWOPI * Q; /\* Was incorrect value for 2pi?? *\/ */
405 /* curfrq = basfrq; */
406 /* for (sumk=0,wsizp=p->winlen,n=nfreqs; n--; ) { */
407 /* *wsizp++ = k = Qfactor/curfrq + 0.5; /\* calc window sizes *\/ */
408 /* sumk += k; /\* and find total *\/ */
409 /* curfrq *= frqmlt; */
410 /* } */
411 /* if ((windsiz = *(p->winlen)) > nsamps) {/\* chk longest windsiz *\/ */
412 /* return csound->InitError(csound, Str("Q %4.1f needs %d samples, " */
413 /* "octdown has just %d"), */
414 /* Q, windsiz, nsamps); */
415 /* } */
416 /* else csound->Message(csound, Str("noctdft: Q %4.1f uses %d of " */
417 /* "%d samps per octdown\n"), */
418 /* Q, windsiz, nsamps); */
419 /* auxsiz = (nsamps + 2*sumk) * sizeof(MYFLT);/\* calc local space reqd *\/ */
420 /* csound->AuxAlloc(csound, (size_t)auxsiz, &p->auxch);
421 /\* & alloc auxspace *\/ */
422 /* fltp = (MYFLT *) p->auxch.auxp; */
423 /* p->linbufp = fltp; fltp += nsamps;
424 /\* linbuf must handle nsamps *\/ */
425 /* p->sinp = sinp = fltp; fltp += sumk; */
426 /* p->cosp = cosp = fltp; /\* cos gets rem sumk *\/ */
427 /* wsizp = p->winlen; */
428 /* for (curfrq=basfrq,n=nfreqs; n--; ) { /\* now fill tables *\/ */
429 /* windsiz = *wsizp++; */
430 /* onedws = 1.0 / windsiz; */
431 /* pidws = PI / windsiz; */
432 /* for (k=0; k<windsiz; k++) { /\* with sines *\/ */
433 /* a = sin(k * pidws); */
434 /* windamp = a * a; /\* times hanning *\/ */
435 /* if (!hanning) */
436 /* windamp = 0.08 + 0.92 * windamp; /\* or hamming *\/ */
437 /* windamp *= onedws; /\* scaled *\/ */
438 /* theta = k * curfrq; */
439 /* *sinp++ = windamp * sin(theta); */
440 /* *cosp++ = windamp * cos(theta); */
441 /* } */
442 /* curfrq *= frqmlt; /\* step by log freq *\/ */
443 /* } */
444 /* if (*p->idsines != FL(0.0)) { */
445 /* /\* if reqd, display windowed sines immediately *\/ */
446 /* csound->dispset(csound, &p->dwindow, p->sinp, (int32_t) sumk, */
447 /* Str("octdft windowed sines:"), 0, "octdft"); */
448 /* csound->display(csound, &p->dwindow); */
449 /* } */
450 /* SPECset(csound, */
451 /* specp, (int64_t)ncoefs); /\* prep the spec dspace *\/ */
452 /* specp->downsrcp = downp; /\* & record its source *\/ */
453 /* } */
454 /* specp->nfreqs = p->nfreqs; \* save the spec descriptors *\/ */
455 /* specp->dbout = p->dbout; */
456 /* specp->ktimstamp = 0; \* init specdata to not new *\/ */
457 /* specp->ktimprd = p->timcount; */
458 /* p->countdown = p->timcount; \* & prime the countdown *\/ */
459 /* return OK; */
460 /* } */
461
462 /* int32_t noctdft(CSOUND *csound, NOCTDFT *p) */
463 /* { */
464 /* DOWNDAT *downp; */
465 /* SPECDAT *specp; */
466 /* OCTDAT *octp; */
467 /* MYFLT *dftp; */
468 /* int32_t nocts, wrap; */
469 /* MYFLT a, b; */
470 /* double c; */
471
472 /* if ((--p->countdown)) return;
473 /\* if not yet time for new spec, return *\/ */
474 /* if (p->auxch.auxp==NULL) { /\* RWD fix *\/ */
475 /* return csound->PerfError(csound, &(p->h), */
476 /* Str("noctdft: not initialised")); */
477 /* } */
478 /* p->countdown = p->timcount; /\* else reset counter & proceed: *\/ */
479 /* downp = p->dsig; */
480 /* specp = p->wsig; */
481 /* nocts = downp->nocts; */
482 /* octp = downp->octdata + nocts; */
483 /* dftp = (MYFLT *) specp->auxch.auxp; */
484 /* while (nocts--) { */
485 /* MYFLT *bufp, *sinp, *cosp; */
486 /* int32_t len, *lenp, nfreqs; */
487 /* MYFLT *begp, *curp, *endp; */
488 /* octp--; /\* for each octave (low to high) *\/ */
489 /* begp = octp->begp; */
490 /* curp = octp->curp; */
491 /* endp = octp->endp; */
492 /* wrap = curp - begp; */
493 /* bufp = p->linbufp; */
494 /* while (curp < endp) /\* copy circbuf to linbuf *\/ */
495 /* *bufp++ = *curp++; */
496 /* for (curp=begp,len=wrap; len--; ) */
497 /* *bufp++ = *curp++; */
498 /* cosp = p->cosp; /\* get start windowed sines *\/ */
499 /* sinp = p->sinp; */
500 /* lenp = p->winlen; */
501 /* for (nfreqs=p->nfreqs; nfreqs--; ) { /\* now for each freq this oct: *\/ */
502 /* a = 0.0; */
503 /* b = 0.0; */
504 /* bufp = p->linbufp; */
505 /* for (len = *lenp++; len--; bufp++) {/\* apply windowed sine seg *\/ */
506 /* a += *bufp * *cosp++; */
507 /* b += *bufp * *sinp++; */
508 /* } */
509 /* c = a*a + b*b; /\* get magnitude squared *\/ */
510 /* if (!(p->dbout)) /\* & optionally convert *\/ */
511 /* c = sqrt(c); /\* to mag or db *\/ */
512 /* else if (p->dbout == 1) { */
513 /* if (c < .001) c = .001; */
514 /* c = 10. * log10(c); */
515 /* } */
516 /* *dftp++ = c; /\* store in out spectrum *\/ */
517 /* } */
518 /* } */
519 /* specp->ktimstamp = CS_KCNT; /\* time-stamp the output *\/ */
520 /* return OK; */
521 /* } */
522
spdspset(CSOUND * csound,SPECDISP * p)523 int32_t spdspset(CSOUND *csound, SPECDISP *p)
524 {
525 char strmsg[256];
526 /* RWD is this enough? */
527 if (UNLIKELY(p->wsig->auxch.auxp==NULL)) {
528 return csound->InitError(csound, Str("specdisp: not initialised"));
529 }
530 if (UNLIKELY((p->timcount = (int32_t)(CS_EKR * *p->iprd)) <= 0)) {
531 return csound->InitError(csound, Str("illegal iperiod"));
532 }
533 if (!(p->dwindow.windid)) {
534 SPECDAT *specp = p->wsig;
535 DOWNDAT *downp = specp->downsrcp;
536 if (downp->lofrq > FL(5.0)) {
537 snprintf(strmsg, 256,
538 Str("instr %d %s, dft (%s), %d octaves (%d - %d Hz):"),
539 (int32_t) p->h.insdshead->p1.value, p->h.optext->t.inlist->arg[0],
540 outstring[specp->dbout],
541 downp->nocts, (int32_t)downp->lofrq, (int32_t)downp->hifrq);
542 }
543 else { /* more detail if low frequency */
544 snprintf(strmsg, 256,
545 Str("instr %d %s, dft (%s), %d octaves (%3.1f - %3.1f Hz):"),
546 (int32_t) p->h.insdshead->p1.value, p->h.optext->t.inlist->arg[0],
547 outstring[specp->dbout],
548 downp->nocts, downp->lofrq, downp->hifrq);
549 }
550 csound->dispset(csound, &p->dwindow, (MYFLT*) specp->auxch.auxp,
551 (int32_t)specp->npts, strmsg, (int32_t)*p->iwtflg,
552 "specdisp");
553 }
554 p->countdown = p->timcount; /* prime the countdown */
555 return OK;
556 }
557
specdisp(CSOUND * csound,SPECDISP * p)558 int32_t specdisp(CSOUND *csound, SPECDISP *p)
559 {
560 /* RWD is this enough? */
561 if (UNLIKELY(p->wsig->auxch.auxp==NULL)) goto err1;
562 if (!(--p->countdown)) { /* on countdown */
563 csound->display(csound, &p->dwindow); /* display spect */
564 p->countdown = p->timcount; /* & reset count */
565 }
566 return OK;
567 err1:
568 return csound->PerfError(csound, &(p->h),
569 Str("specdisp: not initialised"));
570 }
571
sptrkset(CSOUND * csound,SPECPTRK * p)572 int32_t sptrkset(CSOUND *csound, SPECPTRK *p)
573 {
574 SPECDAT *inspecp = p->wsig;
575 int32_t npts, nptls, nn, lobin;
576 int32_t *dstp, ptlmax, inc;
577 MYFLT nfreqs, rolloff, *oct0p, *flop, *fhip, *fundp, *fendp, *fp;
578 MYFLT weight, weightsum, dbthresh, ampthresh;
579
580 if ((npts = inspecp->npts) != p->winpts) { /* if size has changed */
581 SPECset(csound,
582 &p->wfund, (int32_t)npts); /* realloc for wfund */
583 p->wfund.downsrcp = inspecp->downsrcp;
584 p->fundp = (MYFLT *) p->wfund.auxch.auxp;
585 p->winpts = npts;
586 }
587 if ((p->ftimcnt = (int32_t)(CS_EKR**p->ifprd)) > 0) {
588 /* if displaying wfund */
589 SPECDISP *fdp = &p->fdisplay;
590 fdp->h = p->h;
591 fdp->wsig = &p->wfund; /* pass the param pntrs */
592 fdp->iprd = p->ifprd;
593 fdp->iwtflg = p->iwtflg;
594 /* fdp->altname = "specptrk"; */
595 /* fdp->altarg = "X-corr"; */
596 p->wfund.dbout = inspecp->dbout;
597 spdspset(csound,fdp); /* & call specdisp init */
598 }
599 else p->ftimcnt = 0;
600 if (UNLIKELY((nptls = (int32_t)*p->inptls) <= 0 || nptls > MAXPTL)) {
601 return csound->InitError(csound, Str("illegal no of partials"));
602 }
603 p->nptls = nptls; /* number, whether all or odd */
604 if (*p->iodd == FL(0.0)) {
605 ptlmax = nptls;
606 inc = 1;
607 } else {
608 ptlmax = nptls * 2 - 1;
609 inc = 2;
610 }
611 dstp = p->pdist;
612 nfreqs = (MYFLT)inspecp->nfreqs;
613 for (nn = 1; nn <= ptlmax; nn += inc)
614 *dstp++ = (int32_t) ((log((double) nn) / LOGTWO) * nfreqs + 0.5);
615 if ((rolloff = *p->irolloff) == 0.0 || rolloff == 1.0 || nptls == 1) {
616 p->rolloff = 0;
617 weightsum = (MYFLT)nptls;
618 } else {
619 MYFLT *fltp = p->pmult;
620 MYFLT octdrop = (FL(1.0) - rolloff) / nfreqs;
621 weightsum = FL(0.0);
622 for (dstp = p->pdist, nn = nptls; nn--; ) {
623 weight = FL(1.0) - octdrop * *dstp++; /* rolloff * octdistance */
624 weightsum += weight;
625 *fltp++ = weight;
626 }
627 if (UNLIKELY(*--fltp < FL(0.0))) {
628 return csound->InitError(csound, Str("per oct rolloff too steep"));
629 }
630 p->rolloff = 1;
631 }
632 lobin = (int32_t)(inspecp->downsrcp->looct * nfreqs);
633 oct0p = p->fundp - lobin; /* virtual loc of oct 0 */
634
635 flop = oct0p + (int32_t)(*p->ilo * nfreqs);
636 fhip = oct0p + (int32_t)(*p->ihi * nfreqs);
637 fundp = p->fundp;
638 fendp = fundp + inspecp->npts;
639 if (flop < fundp) flop = fundp;
640 if (fhip > fendp) fhip = fendp;
641 if (UNLIKELY(flop >= fhip)) { /* chk hi-lo range valid */
642 return csound->InitError(csound, Str("illegal lo-hi values"));
643 }
644 for (fp = fundp; fp < flop; )
645 *fp++ = FL(0.0); /* clear unused lo and hi range */
646 for (fp = fhip; fp < fendp; )
647 *fp++ = FL(0.0);
648
649 csound->Warning(csound, Str("specptrk: %d freqs, %d%s ptls at "),
650 (int32_t)nfreqs, (int32_t)nptls, inc==2 ? Str(" odd") : "");
651 for (nn = 0; nn < nptls; nn++)
652 csound->Warning(csound, "\t%d", p->pdist[nn]);
653 if (p->rolloff) {
654 csound->Warning(csound, Str("\n\t\trolloff vals:"));
655 for (nn = 0; nn < nptls; nn++)
656 csound->Warning(csound, "\t%4.2f", p->pmult[nn]);
657 }
658
659 dbthresh = *p->idbthresh; /* thresholds: */
660 ampthresh = (MYFLT)exp((double)dbthresh * LOG10D20);
661 switch(inspecp->dbout) {
662 case 0: p->threshon = ampthresh; /* mag */
663 p->threshoff = ampthresh / FL(2.0);
664 break;
665 case 1: p->threshon = dbthresh; /* db */
666 p->threshoff = dbthresh - FL(6.0);
667 break;
668 case 2: p->threshon = ampthresh * ampthresh; /* mag sqrd */
669 p->threshoff = p->threshon / FL(4.0);
670 break;
671 case 3: p->threshon = (MYFLT)sqrt(ampthresh); /* root mag */
672 p->threshoff = p->threshon / FL(1.414);
673 break;
674 }
675 p->threshon *= weightsum;
676 p->threshoff *= weightsum;
677 csound->Warning(csound, Str("\n\tdbthresh %4.1f: X-corr %s "
678 "threshon %4.1f, threshoff %4.1f\n"),
679 dbthresh, outstring[inspecp->dbout],
680 p->threshon, p->threshoff);
681 p->oct0p = oct0p; /* virtual loc of oct 0 */
682 p->confact = *p->iconf;
683 p->flop = flop;
684 p->fhip = fhip;
685 p->kinterp = (*p->interp == FL(0.0)) ? 0 : 1;
686 p->playing = 0;
687 p->kvalsav = *p->istrt;
688 p->kval = p->kinc = FL(0.0);
689 p->kavl = p->kanc = FL(0.0);
690 p->jmpcount = 0;
691 return OK;
692 }
693
694 #define STARTING 1
695 #define PLAYING 2
696
specptrk(CSOUND * csound,SPECPTRK * p)697 int32_t specptrk(CSOUND *csound, SPECPTRK *p)
698 {
699 SPECDAT *inspecp = p->wsig;
700
701 if (inspecp->ktimstamp == CS_KCNT) { /* if inspectrum is new: */
702 MYFLT *inp = (MYFLT *) inspecp->auxch.auxp;
703 MYFLT *endp = inp + inspecp->npts;
704 MYFLT *inp2, sum, *fp;
705 int32_t nn, *pdist, confirms;
706 MYFLT kval, kvar, fmax, *fmaxp, absdiff, realbin;
707 MYFLT *flop, *fhip, *ilop, *ihip, a, b, c, denom, delta;
708 int32_t lobin, hibin;
709
710 if (UNLIKELY(inp==NULL)) goto err1; /* RWD fix */
711 kvar = FABS(*p->kvar);
712 kval = p->playing == PLAYING ? p->kval : p->kvalsav;
713 lobin =
714 (int32_t)((kval-kvar) * inspecp->nfreqs); /* set lim of frq interest */
715 hibin = (int32_t)((kval+kvar) * inspecp->nfreqs);
716 if ((flop = p->oct0p + lobin) < p->flop) /* as fundp bin pntrs */
717 flop = p->flop;
718 if ((fhip = p->oct0p + hibin) > p->fhip) /* within hard limits */
719 fhip = p->fhip;
720 ilop = inp + (flop - p->fundp); /* similar for input bins */
721 ihip = inp + (fhip - p->fundp);
722 if (p->ftimcnt) { /* if displaying, */
723 for (fp = p->flop; fp < flop; ) /* clr to limits */
724 *fp++ = FL(0.0);
725 for (fp = p->fhip; fp > fhip; )
726 *--fp = FL(0.0);
727 }
728 inp = ilop;
729 fp = flop;
730 if (p->rolloff) {
731 MYFLT *pmult;
732 do {
733 sum = *inp;
734 pdist = p->pdist + 1;
735 pmult = p->pmult + 1;
736 for (nn = p->nptls; --nn; ) {
737 if ((inp2 = inp + *pdist++) >= endp)
738 break;
739 sum += *inp2 * *pmult++;
740 }
741 *fp++ = sum;
742 } while (++inp < ihip);
743 }
744 else {
745 do {
746 sum = *inp;
747 pdist = p->pdist + 1;
748 for (nn = p->nptls; --nn; ) {
749 if ((inp2 = inp + *pdist++) >= endp)
750 break;
751 sum += *inp2;
752 }
753 *fp++ = sum;
754 } while (++inp < ihip);
755 }
756 fp = flop; /* now srch fbins for peak */
757 for (fmaxp = fp, fmax = *fp; ++fp<fhip; )
758 if (*fp > fmax) {
759 fmax = *fp;
760 fmaxp = fp;
761 }
762 if (!p->playing) {
763 if (fmax > p->threshon) /* not playing & threshon? */
764 p->playing = STARTING; /* prepare to turn on */
765 else goto output;
766 }
767 else {
768 if (fmax < p->threshoff) { /* playing & threshoff ? */
769 if (p->playing == PLAYING)
770 p->kvalsav = p->kval; /* save val & turn off */
771 p->kval = FL(0.0);
772 p->kavl = FL(0.0);
773 p->kinc = FL(0.0);
774 p->kanc = FL(0.0);
775 p->playing = 0;
776 goto output;
777 }
778 }
779 a = fmaxp>flop ? *(fmaxp-1) : FL(0.0); /* calc a refined bin no */
780 b = fmax;
781 c = fmaxp<fhip-1 ? *(fmaxp+1) : FL(0.0);
782 if (b < FL(2.0) * (a + c))
783 denom = b * FL(2.0) - a - c;
784 else denom = a + b + c;
785 if (denom != FL(0.0))
786 delta = FL(0.5) * (c - a) / denom;
787 else delta = FL(0.0);
788 realbin = (fmaxp - p->oct0p) + delta; /* get modified bin number */
789 kval = realbin / inspecp->nfreqs; /* & cvt to true decoct */
790
791 if (p->playing == STARTING) { /* STARTING mode: */
792 absdiff = FABS(kval - p->kvalsav);
793 confirms = (int32_t)(absdiff * p->confact); /* get interval dependency */
794 if (p->jmpcount < confirms) {
795 p->jmpcount += 1; /* if not enough confirms, */
796 goto output; /* must wait some more */
797 } else {
798 p->playing = PLAYING; /* else switch on playing */
799 p->jmpcount = 0;
800 p->kval = kval; /* but suppress interp */
801 p->kinc = FL(0.0);
802 }
803 } else { /* PLAYING mode: */
804 absdiff = FABS(kval - p->kval);
805 confirms = (int32_t)(absdiff * p->confact); /* get interval dependency */
806 if (p->jmpcount < confirms) {
807 p->jmpcount += 1; /* if not enough confirms, */
808 p->kinc = FL(0.0); /* must wait some more */
809 } else {
810 p->jmpcount = 0; /* else OK to jump interval */
811 if (p->kinterp) /* with optional interp */
812 p->kinc = (kval - p->kval) / inspecp->ktimprd;
813 else p->kval = kval;
814 }
815 }
816 fmax += delta * (c - a) / FL(4.0); /* get modified amp */
817 if (p->kinterp) /* & new kanc if interp */
818 p->kanc = (fmax - p->kavl) / inspecp->ktimprd;
819 else p->kavl = fmax;
820 }
821 output:
822 *p->koct = p->kval; /* output true decoct & amp */
823 *p->kamp = p->kavl;
824 if (p->kinterp) { /* interp if reqd */
825 p->kval += p->kinc;
826 p->kavl += p->kanc;
827 }
828 if (p->ftimcnt)
829 specdisp(csound,&p->fdisplay);
830 return OK;
831 err1:
832 return csound->PerfError(csound, &(p->h),
833 Str("specptrk: not initialised"));
834 }
835
spsumset(CSOUND * csound,SPECSUM * p)836 int32_t spsumset(CSOUND *csound, SPECSUM *p)
837 {
838 IGN(csound);
839 p->kinterp = (*p->interp == FL(0.0)) ? 0 : 1;
840 p->kinc = p->kval = FL(0.0);
841 return OK;
842 }
843
specsum(CSOUND * csound,SPECSUM * p)844 int32_t specsum(CSOUND *csound, SPECSUM *p)
845 /* sum all vals of a spectrum and put as ksig */
846 /* optionally interpolate the output */
847 {
848 SPECDAT *specp = p->wsig;
849 if (UNLIKELY(specp->auxch.auxp==NULL)) goto err1; /* RWD fix */
850 if (specp->ktimstamp == CS_KCNT) { /* if spectrum is new */
851 MYFLT *valp = (MYFLT *) specp->auxch.auxp;
852 MYFLT sum = FL(0.0);
853 int32_t n,npts = specp->npts; /* sum all the values */
854 for (n=0;n<npts;n++) {
855 sum += valp[n];
856 }
857 if (p->kinterp) /* new kinc if interp */
858 p->kinc = (sum - p->kval) / specp->ktimprd;
859 else p->kval = sum;
860 }
861 *p->ksum = p->kval; /* output current kval */
862 if (p->kinterp) /* & interp if reqd */
863 p->kval += p->kinc;
864 return OK;
865 err1:
866 return csound->PerfError(csound, &(p->h),
867 Str("specsum: not initialised"));
868 }
869
spadmset(CSOUND * csound,SPECADDM * p)870 int32_t spadmset(CSOUND *csound, SPECADDM *p)
871 {
872 SPECDAT *inspec1p = p->wsig1;
873 SPECDAT *inspec2p = p->wsig2;
874 int32_t npts;
875
876 if (UNLIKELY((npts = inspec1p->npts) != inspec2p->npts))
877 /* inspecs must agree in size */
878 return csound->InitError(csound, Str("inputs have different sizes"));
879 if (UNLIKELY(inspec1p->ktimprd != inspec2p->ktimprd))
880 /* time period */
881 return csound->InitError(csound, Str("inputs have diff. time periods"));
882 if (UNLIKELY(inspec1p->nfreqs != inspec2p->nfreqs))
883 /* frq resoltn */
884 return csound->InitError(csound,
885 Str("inputs have different freq resolution"));
886 if (UNLIKELY(inspec1p->dbout != inspec2p->dbout))
887 /* and db type */
888 return csound->InitError(csound, Str("inputs have different amptypes"));
889 if (npts != p->waddm->npts) { /* if out does not match ins */
890 SPECset(csound,
891 p->waddm, (int32_t)npts); /* reinit the out spec */
892 p->waddm->downsrcp = inspec1p->downsrcp;
893 }
894 p->waddm->ktimprd = inspec1p->ktimprd; /* pass the other specinfo */
895 p->waddm->nfreqs = inspec1p->nfreqs;
896 p->waddm->dbout = inspec1p->dbout;
897 p->waddm->ktimstamp = 0; /* mark the outspec not new */
898 return OK;
899 }
900
specaddm(CSOUND * csound,SPECADDM * p)901 int32_t specaddm(CSOUND *csound, SPECADDM *p)
902 {
903 if (UNLIKELY((p->wsig1->auxch.auxp==NULL) || /* RWD fix */
904 (p->wsig2->auxch.auxp==NULL) ||
905 (p->waddm->auxch.auxp==NULL))) goto err1;
906 if (p->wsig1->ktimstamp == CS_KCNT) { /* if inspec1 is new: */
907 MYFLT *in1p = (MYFLT *) p->wsig1->auxch.auxp;
908 MYFLT *in2p = (MYFLT *) p->wsig2->auxch.auxp;
909 MYFLT *outp = (MYFLT *) p->waddm->auxch.auxp;
910 MYFLT mul2 = p->mul2;
911 int32_t n,npts = p->wsig1->npts;
912
913 for (n=0;n<npts;n++) {
914 outp[n] = in1p[n] + in2p[n] * mul2; /* out = in1 + in2 * mul2 */
915 }
916 p->waddm->ktimstamp = CS_KCNT; /* mark the output spec as new */
917 }
918 return OK;
919 err1:
920 return csound->PerfError(csound, &(p->h),
921 Str("specaddm: not initialised"));
922 }
923
spdifset(CSOUND * csound,SPECDIFF * p)924 int32_t spdifset(CSOUND *csound, SPECDIFF *p)
925 {
926 SPECDAT *inspecp = p->wsig;
927 MYFLT *lclp;
928 MYFLT *outp;
929 int32_t npts;
930
931 if ((npts = inspecp->npts) != p->specsave.npts) { /* if inspec not matched */
932 SPECset(csound,
933 &p->specsave, (int32_t)npts); /* reinit the save spec */
934 SPECset(csound,
935 p->wdiff, (int32_t)npts); /* & the out diff spec */
936 p->wdiff->downsrcp = inspecp->downsrcp;
937 }
938 p->wdiff->ktimprd = inspecp->ktimprd; /* pass the other specinfo */
939 p->wdiff->nfreqs = inspecp->nfreqs;
940 p->wdiff->dbout = inspecp->dbout;
941 lclp = (MYFLT *) p->specsave.auxch.auxp;
942 outp = (MYFLT *) p->wdiff->auxch.auxp;
943 if (UNLIKELY(lclp==NULL || outp==NULL)) { /* RWD */
944 return csound->InitError(csound,
945 Str("specdiff: local buffers not initialised"));
946 }
947 memset(lclp, 0, npts*sizeof(MYFLT)); /* clr local & out spec bufs */
948 memset(outp, 0, npts*sizeof(MYFLT));
949 p->wdiff->ktimstamp = 0; /* mark the out spec not new */
950 return OK;
951 }
952
specdiff(CSOUND * csound,SPECDIFF * p)953 int32_t specdiff(CSOUND *csound, SPECDIFF *p)
954 {
955 SPECDAT *inspecp = p->wsig;
956
957 if (UNLIKELY((inspecp->auxch.auxp==NULL) /* RWD fix */
958 ||
959 (p->specsave.auxch.auxp==NULL)
960 ||
961 (p->wdiff->auxch.auxp==NULL))) goto err1;
962 if (inspecp->ktimstamp == CS_KCNT) { /* if inspectrum is new: */
963 MYFLT *newp = (MYFLT *) inspecp->auxch.auxp;
964 MYFLT *prvp = (MYFLT *) p->specsave.auxch.auxp;
965 MYFLT *difp = (MYFLT *) p->wdiff->auxch.auxp;
966 MYFLT newval, prvval, diff, possum = FL(0.0); /* possum not used! */
967 int32_t n,npts = inspecp->npts;
968
969 for (n=0; n<npts;n++) {
970 newval = newp[n]; /* compare new & old coefs */
971 prvval = prvp[n];
972 if ((diff = newval-prvval) > FL(0.0)) { /* if new coef > prv coef */
973 difp[n] = diff;
974 possum += diff; /* enter & accum diff */
975 }
976 else difp[n] = FL(0.0); /* else enter zero */
977 prvp[n] = newval; /* sav newval for nxt time */
978 }
979 p->wdiff->ktimstamp = CS_KCNT; /* mark the output spec as new */
980 }
981 return OK;
982 err1:
983 return csound->PerfError(csound, &(p->h),
984 Str("specdiff: not initialised"));
985 }
986
spsclset(CSOUND * csound,SPECSCAL * p)987 int32_t spsclset(CSOUND *csound, SPECSCAL *p)
988 {
989 SPECDAT *inspecp = p->wsig;
990 SPECDAT *outspecp = p->wscaled;
991 FUNC *ftp;
992 int32_t npts;
993
994 if ((npts = inspecp->npts) != outspecp->npts) { /* if size has changed, */
995 SPECset(csound,
996 outspecp, (int32_t)npts); /* realloc */
997 outspecp->downsrcp = inspecp->downsrcp;
998 csound->AuxAlloc(csound, (int32_t)npts * 2 * sizeof(MYFLT), &p->auxch);
999 }
1000 outspecp->ktimprd = inspecp->ktimprd; /* pass the source spec info */
1001 outspecp->nfreqs = inspecp->nfreqs;
1002 outspecp->dbout = inspecp->dbout;
1003 p->fscale = (MYFLT *) p->auxch.auxp; /* setup scale & thresh fn areas */
1004 if (UNLIKELY(p->fscale==NULL)) { /* RWD fix */
1005 return csound->InitError(csound,
1006 Str("specscal: local buffer not initialised"));
1007 }
1008 p->fthresh = p->fscale + npts;
1009 if (UNLIKELY((ftp=csound->FTFind(csound, p->ifscale)) == NULL)) {
1010 /* if fscale given, */
1011 return csound->InitError(csound, Str("missing fscale table"));
1012 }
1013 else {
1014 int32_t nn = npts;
1015 int32_t phs = 0;
1016 int32_t inc = (int32_t)PHMASK / npts;
1017 int32_t lobits = ftp->lobits;
1018 MYFLT *ftable = ftp->ftable;
1019 MYFLT *flp = p->fscale;
1020 for (nn=0;nn<npts;nn++) {
1021 flp[nn] = *(ftable + (phs >> lobits)); /* sample into scale area */
1022 phs += inc;
1023 }
1024 }
1025 if ((p->thresh = (int32_t)*p->ifthresh) &&
1026 (ftp=csound->FTFind(csound, p->ifthresh)) != NULL) {
1027 /* if fthresh given, */
1028 int32_t nn = npts;
1029 int32_t phs = 0;
1030 int32_t inc = (int32_t)PHMASK / npts;
1031 int32_t lobits = ftp->lobits;
1032 MYFLT *ftable = ftp->ftable;
1033 MYFLT *flp = p->fthresh;
1034 for (nn=0;nn<npts;nn++) {
1035 flp[nn] = *(ftable + (phs >> lobits)); /* sample into thresh area */
1036 phs += inc;
1037 }
1038 }
1039 else p->thresh = 0;
1040 outspecp->ktimstamp = 0; /* mark the out spec not new */
1041 return OK;
1042 }
1043
specscal(CSOUND * csound,SPECSCAL * p)1044 int32_t specscal(CSOUND *csound, SPECSCAL *p)
1045 {
1046 SPECDAT *inspecp = p->wsig;
1047 if ((inspecp->auxch.auxp==NULL) /* RWD fix */ ||
1048 (p->wscaled->auxch.auxp==NULL) ||
1049 (p->fscale==NULL)) goto err1;
1050 if (inspecp->ktimstamp == CS_KCNT) { /* if inspectrum is new: */
1051 SPECDAT *outspecp = p->wscaled;
1052 MYFLT *inp = (MYFLT *) inspecp->auxch.auxp;
1053 MYFLT *outp = (MYFLT *) outspecp->auxch.auxp;
1054 MYFLT *sclp = p->fscale;
1055 int32_t n,npts = inspecp->npts;
1056
1057 if (p->thresh) { /* if thresh requested, */
1058 MYFLT *threshp = p->fthresh;
1059 MYFLT val;
1060 for (n=0; n<npts;n++) {
1061 if ((val = inp[n] - threshp[n]) > FL(0.0)) /* for vals above thresh */
1062 outp[n] = val * sclp[n]; /* scale & write out */
1063 else outp[n] = FL(0.0); /* else output is 0. */
1064 }
1065 }
1066 else {
1067 for (n=0; n<npts;n++) {
1068 outp[n] = inp[n] * sclp[n]; /* no thresh: rescale only */
1069 }
1070 }
1071 outspecp->ktimstamp = CS_KCNT; /* mark the outspec as new */
1072 }
1073 return OK;
1074 err1:
1075 return csound->PerfError(csound, &(p->h),
1076 Str("specscal: not initialised"));
1077 }
1078
sphstset(CSOUND * csound,SPECHIST * p)1079 int32_t sphstset(CSOUND *csound, SPECHIST *p)
1080 {
1081 SPECDAT *inspecp = p->wsig;
1082 MYFLT *lclp;
1083 MYFLT *outp;
1084 int32_t npts;
1085
1086 if ((npts = inspecp->npts) != p->accumer.npts) { /* if inspec not matched */
1087 SPECset(csound,
1088 &p->accumer, (int32_t)npts); /* reinit the accum spec */
1089 SPECset(csound,
1090 p->wacout, (int32_t)npts); /* & the output spec */
1091 p->wacout->downsrcp = inspecp->downsrcp;
1092 }
1093 p->wacout->ktimprd = inspecp->ktimprd; /* pass the other specinfo */
1094 p->wacout->nfreqs = inspecp->nfreqs;
1095 p->wacout->dbout = inspecp->dbout;
1096 lclp = (MYFLT *) p->accumer.auxch.auxp;
1097 outp = (MYFLT *) p->wacout->auxch.auxp;
1098 if (UNLIKELY(lclp==NULL || outp==NULL)) { /* RWD fix */
1099 return csound->InitError(csound,
1100 Str("spechist: local buffers not initialised"));
1101 }
1102 memset(lclp,0,npts*sizeof(MYFLT)); /* clr local & out spec bufs */
1103 memset(outp,0,npts*sizeof(MYFLT));
1104 p->wacout->ktimstamp = 0; /* mark the out spec not new */
1105 return OK;
1106 }
1107
spechist(CSOUND * csound,SPECHIST * p)1108 int32_t spechist(CSOUND *csound, SPECHIST *p)
1109 {
1110 SPECDAT *inspecp = p->wsig;
1111 if (UNLIKELY((inspecp->auxch.auxp==NULL) /* RWD fix */
1112 ||
1113 (p->accumer.auxch.auxp==NULL)
1114 ||
1115 (p->wacout->auxch.auxp==NULL))) goto err1;
1116 if (inspecp->ktimstamp == CS_KCNT) { /* if inspectrum is new: */
1117 MYFLT *newp = (MYFLT *) inspecp->auxch.auxp;
1118 MYFLT *acup = (MYFLT *) p->accumer.auxch.auxp;
1119 MYFLT *outp = (MYFLT *) p->wacout->auxch.auxp;
1120 MYFLT newval;
1121 int32_t n,npts = inspecp->npts;
1122
1123 for (n=0;n<npts;n++) {
1124 newval = acup[n] + newp[n]; /* add new to old coefs */
1125 acup[n] = newval; /* sav in accumulator */
1126 outp[n] = newval; /* & copy to output */
1127 }
1128 p->wacout->ktimstamp = CS_KCNT; /* mark the output spec as new */
1129 }
1130 return OK;
1131 err1:
1132 return csound->PerfError(csound, &(p->h),
1133 Str("spechist: not initialised"));
1134 }
1135
spfilset(CSOUND * csound,SPECFILT * p)1136 int32_t spfilset(CSOUND *csound, SPECFILT *p)
1137 {
1138 SPECDAT *inspecp = p->wsig;
1139 SPECDAT *outspecp = p->wfil;
1140 FUNC *ftp;
1141 int32_t npts;
1142
1143 if ((npts = inspecp->npts) != outspecp->npts) { /* if inspec not matched */
1144 SPECset(csound,
1145 outspecp, (int32_t)npts); /* reinit the out spec */
1146 csound->AuxAlloc(csound,
1147 (size_t)npts*2* sizeof(MYFLT),
1148 &p->auxch); /* & local auxspace */
1149 p->coefs = (MYFLT *) p->auxch.auxp; /* reassign filt tbls */
1150 p->states = p->coefs + npts;
1151 }
1152 if (UNLIKELY(p->coefs==NULL || p->states==NULL)) { /* RWD fix */
1153 return csound->InitError(csound,
1154 Str("specfilt: local buffers not initialised"));
1155 }
1156 outspecp->ktimprd = inspecp->ktimprd; /* pass other spect info */
1157 outspecp->nfreqs = inspecp->nfreqs;
1158 outspecp->dbout = inspecp->dbout;
1159 outspecp->downsrcp = inspecp->downsrcp;
1160 if (UNLIKELY((ftp=csound->FTFind(csound, p->ifhtim)) == NULL)) {
1161 /* if fhtim table given, */
1162 return csound->InitError(csound, Str("missing htim ftable"));
1163 }
1164 {
1165 int32_t nn;
1166 int32_t phs = 0;
1167 int32_t inc = (int32_t)PHMASK / npts;
1168 int32_t lobits = ftp->lobits;
1169 MYFLT *ftable = ftp->ftable;
1170 MYFLT *flp = p->coefs;
1171 for (nn=0;nn<npts;nn++) {
1172 flp[nn] = *(ftable + (phs >> lobits)); /* sample into coefs area */
1173 phs += inc;
1174 }
1175 }
1176 {
1177 int32_t nn;
1178 MYFLT *flp = p->coefs;
1179 double halftim, reittim = inspecp->ktimprd * CS_ONEDKR;
1180 for (nn=0;nn<npts;nn++) {
1181 if ((halftim = flp[nn]) > 0.)
1182 flp[nn] = (MYFLT)pow(0.5, reittim/halftim);
1183 else {
1184 return csound->InitError(csound,
1185 Str("htim ftable must be all-positive"));
1186 }
1187 }
1188 }
1189 csound->Warning(csound, Str("coef range: %6.3f - %6.3f\n"),
1190 *p->coefs, *(p->coefs+npts-1));
1191 {
1192 MYFLT *flp = (MYFLT *) p->states;
1193 memset(flp,0,npts*sizeof(MYFLT)); /* clr the persist buf state mem */
1194 }
1195 outspecp->ktimstamp = 0; /* mark the output spec as not new */
1196 return OK;
1197 }
1198
specfilt(CSOUND * csound,SPECFILT * p)1199 int32_t specfilt(CSOUND *csound, SPECFILT *p)
1200 {
1201 if (p->wsig->ktimstamp == CS_KCNT) { /* if input spec is new, */
1202 SPECDAT *inspecp = p->wsig;
1203 SPECDAT *outspecp = p->wfil;
1204 MYFLT *newp = (MYFLT *) inspecp->auxch.auxp;
1205 MYFLT *outp = (MYFLT *) outspecp->auxch.auxp;
1206 MYFLT curval, *coefp = p->coefs;
1207 MYFLT *persp = p->states;
1208 int32_t n,npts = inspecp->npts;
1209
1210 if (UNLIKELY(newp==NULL || outp==NULL ||
1211 coefp==NULL || persp==NULL)) /* RWD */
1212 goto err1;
1213 for (n=0; n<npts;n++) { /* for npts of inspec: */
1214 outp[n] = curval = persp[n]; /* output current point */
1215 persp[n] = coefp[n] * curval + newp[n]; /* decay & addin newval */
1216 }
1217 outspecp->ktimstamp = CS_KCNT; /* mark output spec as new */
1218 }
1219 return OK;
1220 err1:
1221 return csound->PerfError(csound, &(p->h),
1222 Str("specfilt: not initialised"));
1223 }
1224
1225 #define S sizeof
1226
1227 static OENTRY spectra_localops[] = {
1228 { "spectrum", S(SPECTRUM),_QQ, 3, "w", "kiiiqoooo",
1229 (SUBR)spectset,(SUBR)spectrum },
1230 { "spectrum", S(SPECTRUM),_QQ, 3, "w", "aiiiqoooo",
1231 (SUBR)spectset,(SUBR)spectrum },
1232 { "specaddm", S(SPECADDM),_QQ, 3, "w", "wwp", (SUBR)spadmset, (SUBR)specaddm},
1233 { "specdiff", S(SPECDIFF),_QQ, 3, "w", "w", (SUBR)spdifset, (SUBR)specdiff},
1234 { "specscal", S(SPECSCAL),_QQ, 3, "w", "wii", (SUBR)spsclset, (SUBR)specscal},
1235 { "spechist", S(SPECHIST),_QQ, 3, "w", "w", (SUBR)sphstset, (SUBR)spechist},
1236 { "specfilt", S(SPECFILT),_QQ, 3, "w", "wi", (SUBR)spfilset, (SUBR)specfilt},
1237 { "specptrk", S(SPECPTRK),_QQ, 3, "kk", "wkiiiiiioqooo",
1238 (SUBR)sptrkset,(SUBR)specptrk},
1239 { "specsum", S(SPECSUM), _QQ, 3, "k", "wo", (SUBR)spsumset, (SUBR)specsum },
1240 { "specdisp", S(SPECDISP),_QQ, 3, "", "wio", (SUBR)spdspset, (SUBR)specdisp},
1241 { "pitch", S(PITCH), 0, 3, "kk", "aiiiiqooooojo",
1242 (SUBR)pitchset, (SUBR)pitch },
1243 { "maca", S(SUM), 0, 3, "a", "y", (SUBR)macset, (SUBR)maca },
1244 { "mac", S(SUM), 0, 3, "a", "Z", (SUBR)macset, (SUBR)mac },
1245 { "clockon", S(CLOCK), 0, 3, "", "i", (SUBR)clockset, (SUBR)clockon, NULL },
1246 { "clockoff", S(CLOCK),0, 3, "", "i", (SUBR)clockset, (SUBR)clockoff, NULL },
1247 { "readclock", S(CLKRD),0, 1, "i", "i", (SUBR)clockread, NULL, NULL },
1248 { "readscratch", S(SCRATCHPAD),0, 1, "i", "o", (SUBR)scratchread, NULL, NULL },
1249 { "writescratch", S(SCRATCHPAD),0, 1, "", "io", (SUBR)scratchwrite, NULL, NULL },
1250 { "pitchamdf",S(PITCHAMDF),0,3,"kk","aiioppoo",
1251 (SUBR)pitchamdfset, (SUBR)pitchamdf },
1252 { "hsboscil",S(HSBOSC), TR, 3, "a", "kkkiiioo",
1253
1254 (SUBR)hsboscset,(SUBR)hsboscil },
1255 { "phasorbnk", S(PHSORBNK),0,3,"a", "xkio",
1256 (SUBR)phsbnkset, (SUBR)phsorbnk },
1257 { "phasorbnk.k", S(PHSORBNK),0,3,"k", "xkio",
1258 (SUBR)phsbnkset, (SUBR)kphsorbnk, NULL},
1259 { "adsynt",S(HSBOSC), TR, 3, "a", "kkiiiio", (SUBR)adsyntset, (SUBR)adsynt },
1260 { "mpulse", S(IMPULSE), 0, 3, "a", "kko",
1261 (SUBR)impulse_set, (SUBR)impulse },
1262 { "lpf18", S(LPF18), 0, 3, "a", "axxxo", (SUBR)lpf18set, (SUBR)lpf18db },
1263 { "waveset", S(BARRI), 0, 3, "a", "ako", (SUBR)wavesetset, (SUBR)waveset},
1264 { "pinkish", S(PINKISH), 0, 3, "a", "xoooo", (SUBR)pinkset, (SUBR)pinkish },
1265 { "noise", S(VARI), 0, 3, "a", "xk", (SUBR)varicolset, (SUBR)varicol },
1266 { "transeg", S(TRANSEG),0, 3, "k", "iiim",
1267 (SUBR)trnset,(SUBR)ktrnseg, NULL},
1268 { "transeg.a", S(TRANSEG),0, 3, "a", "iiim",
1269 (SUBR)trnset,(SUBR)trnseg},
1270 { "transegb", S(TRANSEG),0, 3, "k", "iiim",
1271 (SUBR)trnset_bkpt,(SUBR)ktrnseg,(SUBR)NULL},
1272 { "transegb.a", S(TRANSEG),0, 3, "a", "iiim",
1273 (SUBR)trnset_bkpt,(SUBR)trnseg },
1274 { "transegr", S(TRANSEG),0, 3, "k", "iiim",
1275 (SUBR)trnsetr,(SUBR)ktrnsegr,(SUBR)NULL },
1276 { "transegr.a", S(TRANSEG),0, 3, "a", "iiim",
1277 (SUBR)trnsetr,(SUBR)trnsegr },
1278 { "clip", S(CLIP), 0, 3, "a", "aiiv", (SUBR)clip_set, (SUBR)clip },
1279 { "cpuprc", S(CPU_PERC),0, 1, "", "Si", (SUBR)cpuperc_S, NULL, NULL },
1280 { "maxalloc", S(CPU_PERC),0, 1, "", "Si", (SUBR)maxalloc_S, NULL, NULL },
1281 { "cpuprc", S(CPU_PERC),0, 1, "", "ii", (SUBR)cpuperc, NULL, NULL },
1282 { "maxalloc", S(CPU_PERC),0, 1, "", "ii", (SUBR)maxalloc, NULL, NULL },
1283 { "active", 0xffff },
1284 { "active.iS", S(INSTCNT),0,1, "i", "Soo", (SUBR)instcount_S, NULL, NULL },
1285 { "active.kS", S(INSTCNT),0,2, "k", "Soo", NULL, (SUBR)instcount_S, NULL },
1286 { "active.i", S(INSTCNT),0,1, "i", "ioo", (SUBR)instcount, NULL, NULL },
1287 { "active.k", S(INSTCNT),0,2, "k", "koo", NULL, (SUBR)instcount, NULL },
1288 { "p.i", S(PFUN), 0,1, "i", "i", (SUBR)pfun, NULL, NULL },
1289 { "p.k", S(PFUNK), 0,3, "k", "k",
1290 (SUBR)pfunk_init, (SUBR)pfunk, NULL },
1291 { "mute", S(MUTE), 0,1, "", "So", (SUBR)mute_inst_S },
1292 { "mute.i", S(MUTE), 0,1, "", "io", (SUBR)mute_inst },
1293 { "median", S(MEDFILT), 0, 3, "a", "akio", (SUBR)medfiltset, (SUBR)medfilt },
1294 { "mediank", S(MEDFILT), 0,3, "k", "kkio", (SUBR)medfiltset, (SUBR)kmedfilt},
1295 };
1296
1297 LINKAGE_BUILTIN(spectra_localops)
1298