1 /*
2 ugens9.c:
3
4 Copyright (C) 1996 Greg Sullivan, 2004 ma++ ingalls
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 "stdopcod.h" /* UGENS9.C */
25 #include <math.h>
26 #include "convolve.h"
27 #include "ugens9.h"
28 #include "soundio.h"
29 #include <inttypes.h>
30
cvset_(CSOUND * csound,CONVOLVE * p,int32_t stringname)31 static int32_t cvset_(CSOUND *csound, CONVOLVE *p, int32_t stringname)
32 {
33 char cvfilnam[MAXNAME];
34 MEMFIL *mfp;
35 MYFLT *fltp;
36 CVSTRUCT *cvh;
37 int32_t siz;
38 int32 Hlenpadded = 1, obufsiz, Hlen;
39 uint32_t nchanls;
40 uint32_t nsmps = CS_KSMPS;
41
42 if (UNLIKELY(csound->oparms->odebug))
43 csound->Message(csound, CONVOLVE_VERSION_STRING);
44
45 if (stringname==0){
46 if (csound->ISSTRCOD(*p->ifilno))
47 strNcpy(cvfilnam,get_arg_string(csound, *p->ifilno), MAXNAME-1);
48 else csound->strarg2name(csound, cvfilnam,p->ifilno, "convolve.",0);
49 }
50 else strNcpy(cvfilnam, ((STRINGDAT *)p->ifilno)->data, MAXNAME-1);
51
52 if ((mfp = p->mfp) == NULL || strcmp(mfp->filename, cvfilnam) != 0) {
53 /* if file not already readin */
54 if (UNLIKELY((mfp = csound->ldmemfile2withCB(csound, cvfilnam,
55 CSFTYPE_CVANAL,NULL))
56 == NULL)) {
57 return csound->InitError(csound,
58 Str("CONVOLVE cannot load %s"), cvfilnam);
59 }
60 }
61 cvh = (CVSTRUCT *)mfp->beginp;
62 if (UNLIKELY(cvh->magic != CVMAGIC)) {
63 return csound->InitError(csound,
64 Str("%s not a CONVOLVE file (magic %"PRIi32")"),
65 cvfilnam, cvh->magic);
66 }
67
68 nchanls = (cvh->channel == ALLCHNLS ? cvh->src_chnls : 1);
69
70 if (*p->channel == FL(0.0)) {
71 if (LIKELY(p->OUTOCOUNT == nchanls))
72 p->nchanls = nchanls;
73 else {
74 return csound->InitError(csound,
75 Str("CONVOLVE: output channels not equal "
76 "to number of channels in source"));
77 }
78 }
79 else {
80 if (*p->channel <= nchanls) {
81 if (UNLIKELY(p->OUTOCOUNT != 1)) {
82 return csound->InitError(csound,
83 Str("CONVOLVE: output channels not equal "
84 "to number of channels in source"));
85 }
86 else
87 p->nchanls = 1;
88 }
89 else {
90 return csound->InitError(csound,
91 Str("CONVOLVE: channel number greater than "
92 "number of channels in source"));
93 }
94 }
95 Hlen = p->Hlen = cvh->Hlen;
96 while (Hlenpadded < 2*Hlen-1)
97 Hlenpadded <<= 1;
98 p->Hlenpadded = Hlenpadded;
99 p->H = (MYFLT *) ((char *)cvh+cvh->headBsize);
100 if ((p->nchanls == 1) && (*p->channel > 0))
101 p->H += (Hlenpadded + 2) * (int32_t)(*p->channel - 1);
102
103 if (UNLIKELY(cvh->samplingRate != CS_ESR)) {
104 /* & chk the data */
105 csound->Warning(csound, Str("%s's srate = %8.0f, orch's srate = %8.0f"),
106 cvfilnam, cvh->samplingRate, CS_ESR);
107 }
108 if (UNLIKELY(cvh->dataFormat != CVMYFLT)) {
109 return csound->InitError(csound,
110 Str("unsupported CONVOLVE data "
111 "format %"PRIi32" in %s"),
112 cvh->dataFormat, cvfilnam);
113 }
114
115 /* Determine size of circular output buffer */
116 if (Hlen >= (int32)nsmps)
117 obufsiz = (int32) CEIL((MYFLT) Hlen / nsmps) * nsmps;
118 else
119 obufsiz = (int32) CEIL(CS_KSMPS / (MYFLT) Hlen) * Hlen;
120
121 siz = ((Hlenpadded + 2) + p->nchanls * ((Hlen - 1) + obufsiz)
122 + (p->nchanls > 1 ? (Hlenpadded + 2) : 0));
123 if (p->auxch.auxp == NULL || p->auxch.size < siz*sizeof(MYFLT)) {
124 /* if no buffers yet, alloc now */
125 csound->AuxAlloc(csound, (size_t) siz*sizeof(MYFLT), &p->auxch);
126 fltp = (MYFLT *) p->auxch.auxp;
127 p->fftbuf = fltp; fltp += (Hlenpadded + 2); /* and insert addresses */
128 p->olap = fltp; fltp += p->nchanls*(Hlen - 1);
129 p->outbuf = fltp; fltp += p->nchanls*obufsiz;
130 p->X = fltp;
131 }
132 else {
133 fltp = (MYFLT *) p->auxch.auxp;
134 memset(fltp, 0, sizeof(MYFLT)*siz);
135 /* for(i=0; i < siz; i++) fltp[i] = FL(0.0); */
136 }
137 p->obufsiz = obufsiz;
138 p->outcnt = obufsiz;
139 p->incount = 0;
140 p->obufend = p->outbuf + obufsiz - 1;
141 p->outhead = p->outail = p->outbuf;
142 p->fwdsetup = csound->RealFFT2Setup(csound, Hlenpadded, FFT_FWD);
143 p->invsetup = csound->RealFFT2Setup(csound, Hlenpadded, FFT_INV);
144 return OK;
145 }
146
cvset(CSOUND * csound,CONVOLVE * p)147 static int32_t cvset(CSOUND *csound, CONVOLVE *p){
148 return cvset_(csound,p,0);
149
150 }
151
cvset_S(CSOUND * csound,CONVOLVE * p)152 static int32_t cvset_S(CSOUND *csound, CONVOLVE *p){
153 return cvset_(csound,p,1);
154
155 }
156 /* Write from a circular buffer into a linear output buffer without
157 clearing data
158 UPDATES SOURCE & DESTINATION POINTERS TO REFLECT NEW POSITIONS */
writeFromCircBuf(MYFLT ** sce,MYFLT ** dst,MYFLT * sceStart,MYFLT * sceEnd,int32 numToDo)159 static void writeFromCircBuf(
160 MYFLT **sce,
161 MYFLT **dst, /* Circular source and linear destination */
162 MYFLT *sceStart,
163 MYFLT *sceEnd, /* Address of start & end of source buffer */
164 int32 numToDo) /* How many points to write (<= circBufSize) */
165 {
166 MYFLT *srcindex = *sce;
167 MYFLT *dstindex = *dst;
168 int32 breakPoint; /* how many points to add before having to wrap */
169
170 breakPoint = sceEnd - srcindex + 1;
171 if (numToDo >= breakPoint) { /* we will do 2 in 1st loop, rest in 2nd. */
172 numToDo -= breakPoint;
173 for (; breakPoint > 0; --breakPoint) {
174 *dstindex++ = *srcindex++;
175 }
176 srcindex = sceStart;
177 }
178 for (; numToDo > 0; --numToDo) {
179 *dstindex++ = *srcindex++;
180 }
181 *sce = srcindex;
182 *dst = dstindex;
183 return;
184 }
185
convolve(CSOUND * csound,CONVOLVE * p)186 static int32_t convolve(CSOUND *csound, CONVOLVE *p)
187 {
188 int32_t nsmpso=CS_KSMPS,nsmpsi=CS_KSMPS,outcnt_sav;
189 int32_t nchm1 = p->nchanls - 1,chn;
190 int32 i,j;
191 MYFLT *ar[4];
192 MYFLT *ai = p->ain;
193 MYFLT *fftbufind;
194 int32 outcnt = p->outcnt;
195 int32 incount=p->incount;
196 int32 Hlen = p->Hlen;
197 int32 Hlenm1 = Hlen - 1;
198 int32 obufsiz = p->obufsiz;
199 MYFLT *outhead = NULL;
200 MYFLT *outail = p->outail;
201 MYFLT *olap;
202 MYFLT *X;
203 int32 Hlenpadded = p->Hlenpadded;
204 MYFLT scaleFac;
205 uint32_t offset = p->h.insdshead->ksmps_offset;
206 uint32_t early = p->h.insdshead->ksmps_no_end;
207 uint32_t nn,nsmpso_sav;
208
209 scaleFac = csound->GetInverseRealFFTScale(csound, (int32_t) Hlenpadded);
210 ar[0] = p->ar1;
211 ar[1] = p->ar2;
212 ar[2] = p->ar3;
213 ar[3] = p->ar4;
214
215 if (UNLIKELY(p->auxch.auxp==NULL)) goto err1;
216 /* First dump as much pre-existing audio in output buffer as possible */
217 if (outcnt > 0) {
218 if (outcnt <= (int32_t)CS_KSMPS)
219 i = outcnt;
220 else
221 i = CS_KSMPS;
222 nsmpso -= i; outcnt -= i;
223 for (chn = nchm1;chn >= 0;chn--) {
224 outhead = p->outhead + chn*obufsiz;
225 writeFromCircBuf(&outhead,&ar[chn],p->outbuf+chn*obufsiz,
226 p->obufend+chn*obufsiz,i);
227 }
228 p->outhead = outhead;
229 }
230 while (nsmpsi > 0) {
231 /* Read input audio and place into work buffer. */
232
233 fftbufind = p->fftbuf + incount;
234 if ((incount + nsmpsi) <= Hlen)
235 i = nsmpsi;
236 else
237 i = Hlen - incount;
238 nsmpsi -= i;
239 incount += i;
240 nsmpso_sav = CS_KSMPS-early;
241 for (nn=0; i>0; nn++, i--) {
242 if (nn<offset|| nn>(uint32_t)nsmpso_sav)
243 *fftbufind++ = FL(0.0);
244 else
245 *fftbufind++ = scaleFac * ai[nn];
246 }
247 if (incount == Hlen) {
248 /* We have enough audio for a convolution. */
249 incount = 0;
250 /* FFT the input (to create X) */
251 /*csound->Message(csound, "CONVOLVE: ABOUT TO FFT\n"); */
252 csound->RealFFT2(csound, p->fwdsetup, p->fftbuf);
253 p->fftbuf[Hlenpadded] = p->fftbuf[1];
254 p->fftbuf[1] = p->fftbuf[Hlenpadded + 1L] = FL(0.0);
255 /* save the result if multi-channel */
256 if (nchm1) {
257 fftbufind = p->fftbuf;
258 X = p->X;
259 for (i = Hlenpadded + 2;i > 0;i--)
260 *X++ = *fftbufind++;
261 }
262 nsmpso_sav = nsmpso;
263 outcnt_sav = outcnt;
264 for (chn = nchm1;chn >= 0;chn--) {
265 outhead = p->outhead + chn*obufsiz;
266 outail = p->outail + chn*obufsiz;
267 olap = p->olap + chn*Hlenm1;
268 if (chn < nchm1) {
269 fftbufind = p->fftbuf;
270 X = p->X;
271 for (i = Hlenpadded + 2;i> 0;i--)
272 *fftbufind++ = *X++;
273 }
274 /*csound->Message(csound, "CONVOLVE: ABOUT TO MULTIPLY\n"); */
275 /* Multiply H * X, point for point */
276
277 {
278 MYFLT *a, *b, re, im;
279 int32_t i;
280 a = (MYFLT*) p->H + (int32_t) (chn * (Hlenpadded + 2));
281 b = (MYFLT*) p->fftbuf;
282 for (i = 0; i <= (int32_t) Hlenpadded; i += 2) {
283 re = a[i + 0] * b[i + 0] - a[i + 1] * b[i + 1];
284 im = a[i + 0] * b[i + 1] + a[i + 1] * b[i + 0];
285 b[i + 0] = re;
286 b[i + 1] = im;
287 }
288 }
289
290 /*csound->Message(csound, "CONVOLVE: ABOUT TO IFFT\n"); */
291 /* Perform inverse FFT on X */
292
293 p->fftbuf[1] = p->fftbuf[Hlenpadded];
294 p->fftbuf[Hlenpadded] = p->fftbuf[Hlenpadded + 1L] = FL(0.0);
295 csound->RealFFT2(csound, p->invsetup, p->fftbuf);
296
297 /* Take the first Hlen output samples and output them to
298 either the real audio output buffer or the local circular
299 buffer */
300 nsmpso = nsmpso_sav;
301 outcnt = outcnt_sav;
302 fftbufind = p->fftbuf;
303 if ( (nsmpso > 0)&&(outcnt == 0) ) {
304 /* csound->Message(csound, "Outputting to audio buffer proper\n");*/
305 /* space left in output buffer, and nothing currently in circular
306 buffer, so write as much as possible to output buffer first */
307 if (nsmpso >= Hlenm1) {
308 nsmpso -= Hlenm1;
309 for (i=Hlenm1;i > 0;--i)
310 *ar[chn]++ = *fftbufind++ + *olap++;
311 if (nsmpso > 0) {
312 *ar[chn]++ = *fftbufind++;
313 --nsmpso;
314 }
315 }
316 else {
317 for (;nsmpso > 0;--nsmpso)
318 *ar[chn]++ = *fftbufind++ + *olap++;
319 }
320 }
321 /* Any remaining output must go into circular buffer */
322 /*csound->Message(csound, "Outputting to circ. buffer\n");*/
323 i = Hlen - (fftbufind - p->fftbuf);
324 outcnt += i;
325 i--; /* do first Hlen -1 samples with overlap */
326 j = p->obufend+chn*obufsiz - outail + 1;
327 if (i >= j) {
328 i -= j;
329 for (;j > 0;--j)
330 *outail++ = *fftbufind++ + *olap++;
331 outail = p->outbuf+chn*obufsiz;
332 }
333 for (;i > 0;--i)
334 *outail++ = *fftbufind++ + *olap++;
335 /* just need to do sample at Hlen now */
336 *outail++ = *fftbufind++;
337 if (outail > p->obufend+chn*obufsiz)
338 outail = p->outbuf+chn*obufsiz;
339
340 /* Pad the rest to zero, and save first remaining (Hlen - 1) to overlap
341 buffer */
342 olap = p->olap+chn*Hlenm1;
343 for (i = Hlenm1;i > 0;--i) {
344 *olap++ = *fftbufind;
345 *fftbufind++ = FL(0.0);
346 }
347 olap = p->olap+chn*Hlenm1;
348 /* Now pad the rest to zero as well. In theory, this shouldn't be
349 necessary, however it's conceivable that rounding errors may
350 creep in, and these cells won't be exactly zero. So, let's
351 make absolutely sure */
352 for (i = Hlenpadded - (Hlen+Hlenm1);i > 0;--i)
353 *fftbufind++ = FL(0.0);
354 } /* end main for loop */
355 p->outhead = outhead;
356 p->outail = outail;
357 }
358 } /* end while */
359
360 /* update state in p */
361 p->incount = incount;
362 p->outcnt = outcnt;
363 p->outhead = outhead;
364 p->outail = outail;
365 return OK;
366 err1:
367 return csound->PerfError(csound, &(p->h),
368 Str("convolve: not initialised"));
369 }
370
371 /* partitioned (low latency) overlap-save convolution.
372 we break up the IR into separate blocks, then perform
373 an FFT on each partition. For this reason we ONLY accept
374 soundfiles as input, and do all of the traditional 'cvanal'
375 processing at i-time. it would be nice to eventually
376 have cvanal create a partitioned format, which in turn would
377 allow this opcode to accept .con files.
378 -ma++ april 2004 */
379
pconvset_(CSOUND * csound,PCONVOLVE * p,int32_t stringname)380 static int32_t pconvset_(CSOUND *csound, PCONVOLVE *p, int32_t stringname)
381 {
382 int32_t channel = (*(p->channel) <= 0 ? ALLCHNLS : (int32_t) *(p->channel));
383 SNDFILE *infd;
384 SOUNDIN IRfile;
385 MYFLT *inbuf, *fp1,*fp2;
386 int32 i, j, read_in, part;
387 MYFLT *IRblock;
388 MYFLT ainput_dur, scaleFac;
389 MYFLT partitionSize;
390
391 /* IV - 2005-04-06: fixed bug: was uninitialised */
392 memset(&IRfile, 0, sizeof(SOUNDIN));
393 /* open impulse response soundfile [code derived from SAsndgetset()] */
394 IRfile.skiptime = FL(0.0);
395
396 if (stringname==0){
397 if (csound->ISSTRCOD(*p->ifilno))
398 strNcpy(IRfile.sfname,get_arg_string(csound, *p->ifilno), 511);
399 else csound->strarg2name(csound, IRfile.sfname, p->ifilno, "soundin.",0);
400 }
401 else strNcpy(IRfile.sfname, ((STRINGDAT *)p->ifilno)->data, 511);
402
403 IRfile.sr = 0;
404 if (UNLIKELY(channel < 1 || ((channel > 4) && (channel != ALLCHNLS)))) {
405 return csound->InitError(csound, Str("channel request %d illegal"), channel);
406 }
407 IRfile.channel = channel;
408 IRfile.analonly = 1;
409 if (UNLIKELY((infd = csound->sndgetset(csound, &IRfile)) == NULL)) {
410 return csound->InitError(csound, Str("pconvolve: error while impulse file"));
411 }
412
413 if (UNLIKELY(IRfile.framesrem < 0)) {
414 csound->Warning(csound, Str("undetermined file length, "
415 "will attempt requested duration"));
416 ainput_dur = FL(0.0); /* This is probably wrong -- JPff */
417 }
418 else {
419 IRfile.getframes = IRfile.framesrem;
420 if (UNLIKELY(IRfile.sr==0)) return csound->InitError(csound, Str("SR zero"));
421 ainput_dur = (MYFLT) IRfile.getframes / IRfile.sr;
422 }
423
424 csound->Warning(csound, Str("analyzing %ld sample frames (%3.1f secs)\n"),
425 (long) IRfile.getframes, ainput_dur);
426
427 p->nchanls = (channel != ALLCHNLS ? 1 : IRfile.nchanls);
428 if (UNLIKELY(p->nchanls != (int32_t)p->OUTOCOUNT)) {
429 return csound->InitError(csound, Str("PCONVOLVE: number of output channels "
430 "not equal to input channels"));
431 }
432
433 if (UNLIKELY(IRfile.sr != CS_ESR)) {
434 /* ## RWD suggests performing sr conversion here! */
435 csound->Warning(csound, Str("IR srate != orch's srate"));
436 }
437
438 /* make sure the partition size is nonzero and a power of 2 */
439 if (*p->partitionSize <= 0)
440 partitionSize = csound->oparms->outbufsamps / csound->GetNchnls(csound);
441 else
442 partitionSize = *p->partitionSize;
443
444 p->Hlen = 1;
445 while (p->Hlen < partitionSize)
446 p->Hlen <<= 1;
447
448 p->Hlenpadded = 2*p->Hlen;
449
450 /* determine the number of partitions */
451 p->numPartitions = CEIL((MYFLT)(IRfile.getframes) / (MYFLT)p->Hlen);
452
453 /* set up FFT tables */
454 inbuf = (MYFLT *) csound->Malloc(csound,
455 p->Hlen * p->nchanls * sizeof(MYFLT));
456 csound->AuxAlloc(csound, p->numPartitions * (p->Hlenpadded + 2) *
457 sizeof(MYFLT) * p->nchanls, &p->H);
458 IRblock = (MYFLT *)p->H.auxp;
459 p->fwdsetup = csound->RealFFT2Setup(csound,p->Hlenpadded, FFT_FWD);
460 p->invsetup = csound->RealFFT2Setup(csound,p->Hlenpadded, FFT_INV);
461 /* form each partition and take its FFT */
462 for (part = 0; part < p->numPartitions; part++) {
463 /* get the block of input samples and normalize -- soundin code
464 handles finding the right channel */
465 if (UNLIKELY((read_in = csound->getsndin(csound, infd, inbuf,
466 p->Hlen*p->nchanls, &IRfile)) <= 0))
467 return csound->InitError(csound,
468 Str("PCONVOLVE: less sound than expected!"));
469
470 /* take FFT of each channel */
471 scaleFac = csound->dbfs_to_float
472 * csound->GetInverseRealFFTScale(csound, (int32_t) p->Hlenpadded);
473 for (i = 0; i < p->nchanls; i++) {
474 fp1 = inbuf + i;
475 fp2 = IRblock;
476 for (j = 0; j < read_in/p->nchanls; j++) {
477 *fp2++ = *fp1 * scaleFac;
478 fp1 += p->nchanls;
479 }
480
481 csound->RealFFT2(csound, p->fwdsetup, IRblock);
482 IRblock[p->Hlenpadded] = IRblock[1];
483 IRblock[1] = IRblock[p->Hlenpadded + 1L] = FL(0.0);
484 IRblock += (p->Hlenpadded + 2);
485 }
486 }
487
488 csound->Free(csound, inbuf);
489 csound->FileClose(csound, IRfile.fd);
490
491 /* allocate the buffer saving recent input samples */
492 csound->AuxAlloc(csound, p->Hlen * sizeof(MYFLT), &p->savedInput);
493 p->inCount = 0;
494
495 /* allocate the convolution work buffer */
496 csound->AuxAlloc(csound, (p->Hlenpadded+2) * sizeof(MYFLT), &p->workBuf);
497 p->workWrite = (MYFLT *)p->workBuf.auxp + p->Hlen;
498
499 /* allocate the buffer holding recent past convolutions */
500 csound->AuxAlloc(csound, (p->Hlenpadded+2) * p->numPartitions
501 * p->nchanls * sizeof(MYFLT), &p->convBuf);
502 p->curPart = 0;
503
504 /* allocate circular output sample buffer */
505 p->outBufSiz = sizeof(MYFLT) * p->nchanls *
506 (p->Hlen >= (int32_t)CS_KSMPS ? p->Hlenpadded : 2*(int32_t)CS_KSMPS);
507 csound->AuxAlloc(csound, p->outBufSiz, &p->output);
508 p->outRead = (MYFLT *)p->output.auxp;
509
510 /* if ksmps < hlen, we have to pad initial output to prevent a possible
511 empty ksmps pass after a few initial generated buffers. There is
512 probably an equation to figure this out to reduce the delay, but
513 I can't seem to figure it out */
514 if (p->Hlen > (int32_t)CS_KSMPS) {
515 p->outCount = p->Hlen + CS_KSMPS;
516 p->outWrite = p->outRead + (p->nchanls * p->outCount);
517 }
518 else {
519 p->outCount = 0;
520 p->outWrite = p->outRead;
521 }
522 return OK;
523 }
524
pconvset(CSOUND * csound,PCONVOLVE * p)525 static int32_t pconvset(CSOUND *csound, PCONVOLVE *p){
526 return pconvset_(csound,p,0);
527 }
528
pconvset_S(CSOUND * csound,PCONVOLVE * p)529 static int32_t pconvset_S(CSOUND *csound, PCONVOLVE *p){
530 return pconvset_(csound,p,1);
531 }
532
pconvolve(CSOUND * csound,PCONVOLVE * p)533 static int32_t pconvolve(CSOUND *csound, PCONVOLVE *p)
534 {
535 uint32_t nn, nsmps = CS_KSMPS;
536 uint32_t offset = p->h.insdshead->ksmps_offset;
537 uint32_t early = nsmps - p->h.insdshead->ksmps_no_end;
538 MYFLT *ai = p->ain;
539 MYFLT *buf;
540 MYFLT *input = (MYFLT*) p->savedInput.auxp, *workWrite = p->workWrite;
541 MYFLT *a1 = p->ar1, *a2 = p->ar2, *a3 = p->ar3, *a4 = p->ar4;
542 int32 i, j, count = p->inCount;
543 int32 hlenpaddedplus2 = p->Hlenpadded+2;
544
545 for (nn=0; nn<nsmps; nn++) {
546 /* Read input audio and place into buffer. */
547 input[count++] = *workWrite++ = (nn<offset||nn>early? FL(0.0) : ai[nn]);
548
549 /* We have enough audio for a convolution. */
550 if (count == p->Hlen) {
551 MYFLT *dest = (MYFLT*) p->convBuf.auxp
552 + p->curPart * (p->Hlenpadded + 2) * p->nchanls;
553 MYFLT *h = (MYFLT*) p->H.auxp;
554 MYFLT *workBuf = (MYFLT*) p->workBuf.auxp;
555
556 /* FFT the input (to create X) */
557 *workWrite = FL(0.0); /* zero out nyquist bin from last fft result
558 - maybe is ignored for input(?) but just in case.. */
559 csound->RealFFT2(csound, p->fwdsetup, workBuf);
560 workBuf[p->Hlenpadded] = workBuf[1];
561 workBuf[1] = workBuf[p->Hlenpadded + 1L] = FL(0.0);
562
563 /* for every IR partition convolve and add to previous convolves */
564 for (i = 0; i < p->numPartitions*p->nchanls; i++) {
565 MYFLT *src = workBuf;
566 int32_t n;
567 for (n = 0; n <= (int32_t) p->Hlenpadded; n += 2) {
568 dest[n + 0] += (h[n + 0] * src[n + 0]) - (h[n + 1] * src[n + 1]);
569 dest[n + 1] += (h[n + 1] * src[n + 0]) + (h[n + 0] * src[n + 1]);
570 }
571 h += n; dest += n;
572 if (UNLIKELY(dest == (MYFLT*)p->convBuf.endp))
573 dest = (MYFLT*)p->convBuf.auxp;
574 }
575
576 /* Perform inverse FFT of the ondeck partion block */
577 buf = (MYFLT*) p->convBuf.auxp
578 + p->curPart * p->nchanls * hlenpaddedplus2;
579 for (i = 0; i < p->nchanls; i++) {
580 MYFLT *bufp;
581 bufp = buf + i * hlenpaddedplus2;
582 bufp[1] = bufp[p->Hlenpadded];
583 bufp[p->Hlenpadded] = bufp[p->Hlenpadded + 1L] = FL(0.0);
584 csound->RealFFT2(csound, p->invsetup, bufp);
585 }
586 /* We only take only the last Hlen output samples so we first zero out
587 the first half for next time, then we copy the rest to output buffer
588 */
589 for (j = 0; j < p->nchanls; j++) {
590 MYFLT *outp = p->outWrite + j;
591 for (i = 0; i < p->Hlen; i++)
592 *buf++ = FL(0.0);
593 for (i = 0; i < p->Hlen; i++) {
594 *outp = *buf;
595 *buf++ = FL(0.0);
596 outp += p->nchanls;
597 if (outp >= (MYFLT *)p->output.endp)
598 outp = (MYFLT *)p->output.auxp + j;
599 }
600 buf += 2;
601 }
602 p->outWrite += p->Hlen*p->nchanls;
603 if (p->outWrite >= (MYFLT *)p->output.endp)
604 p->outWrite -= p->outBufSiz/sizeof(MYFLT);
605 p->outCount += p->Hlen;
606 if (++p->curPart == p->numPartitions)
607 /* advance to the next partition */
608 p->curPart = 0;
609 /* copy the saved input into the work buffer for next time around */
610 memcpy(p->workBuf.auxp, input, p->Hlen * sizeof(MYFLT));
611 count = 0;
612 workWrite = (MYFLT *)p->workBuf.auxp + p->Hlen;
613 }
614 } /* end while */
615
616 /* copy to output if we have enough samples [we always should
617 except the first Hlen samples] */
618 if (p->outCount >= (int32_t)CS_KSMPS) {
619 uint32_t n;
620 p->outCount -= CS_KSMPS;
621 for (n=0; n < CS_KSMPS; n++) {
622 switch (p->nchanls) {
623 case 1:
624 *a1++ = *p->outRead++;
625 break;
626 case 2:
627 *a1++ = *p->outRead++;
628 *a2++ = *p->outRead++;
629 break;
630 case 3:
631 *a1++ = *p->outRead++;
632 *a2++ = *p->outRead++;
633 *a3++ = *p->outRead++;
634 break;
635 case 4:
636 *a1++ = *p->outRead++;
637 *a2++ = *p->outRead++;
638 *a3++ = *p->outRead++;
639 *a4++ = *p->outRead++;
640 break;
641 }
642 if (p->outRead == p->output.endp)
643 p->outRead = p->output.auxp;
644 }
645 }
646
647 /* update struct */
648 p->inCount = count;
649 p->workWrite = workWrite;
650 return OK;
651 }
652
653 static OENTRY localops[] =
654 {
655 { "convolve", sizeof(CONVOLVE), 0, 3, "mmmm", "aSo",
656 (SUBR) cvset_S, (SUBR) convolve },
657 { "convle", sizeof(CONVOLVE), 0, 3, "mmmm", "aSo",
658 (SUBR) cvset_S, (SUBR) convolve },
659 { "pconvolve",sizeof(PCONVOLVE), 0, 3, "mmmm", "aSoo",
660 (SUBR) pconvset_S, (SUBR) pconvolve },
661 { "convolve.i", sizeof(CONVOLVE), 0, 3, "mmmm", "aio",
662 (SUBR) cvset, (SUBR) convolve },
663 { "convle.i", sizeof(CONVOLVE), 0, 3, "mmmm", "aio",
664 (SUBR) cvset, (SUBR) convolve },
665 { "pconvolve.i",sizeof(PCONVOLVE), 0, 3, "mmmm", "aioo",
666 (SUBR) pconvset, (SUBR) pconvolve }
667 };
668
ugens9_init_(CSOUND * csound)669 int32_t ugens9_init_(CSOUND *csound)
670 {
671 return csound->AppendOpcodes(csound, &(localops[0]),
672 (int32_t) (sizeof(localops) / sizeof(OENTRY)));
673 }
674