1 /*
2 mxfft.c:
3
4 Copyright (C) 2002 Trevor Wishart, Keith Henderson
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 /* This program converted from the FORTRAN routines by Singleton in
25 * Section 1.4 of "Programs for Digital Signal Processing", IEEE Press, 1979.
26 * Conversion by Trevor Wishart and Keith Henderson, York Univ.
27 */
28 /*
29 static char *rcsid = "$Id$";
30 */
31 /*
32 * $Log$
33 * Revision 1.15 2008-12-04 14:55:04 jpff
34 * Branch prediction
35 *
36 * Revision 1.14 2005/08/13 14:44:37 istvanv
37 * Minor code changes
38 *
39 * Revision 1.13 2005/08/12 19:01:23 istvanv
40 * Renamed ENVIRON to CSOUND
41 *
42 * Revision 1.12 2005/08/10 09:57:07 istvanv
43 * Use CSOUND* type for Csound instance pointers instead of void*
44 *
45 * Revision 1.11 2005/07/15 10:13:28 istvanv
46 * Removed cs.h
47 *
48 * Revision 1.10 2005/06/05 16:36:24 istvanv
49 * Minor code improvements
50 *
51 * Revision 1.9 2005/06/05 13:07:08 istvanv
52 * Added mxfft.c functions to API
53 *
54 * Revision 1.8 2005/05/28 13:00:02 istvanv
55 * Minor code changes (tabs etc.)
56 *
57 * Revision 1.7 2005/04/17 17:56:25 jpff
58 * Warnings
59 *
60 * Revision 1.6 2005/02/18 16:21:14 istvanv
61 * added csound pointer to csound->Malloc, auxalloc, and other functions
62 *
63 * Revision 1.5 2005/01/27 19:22:50 istvanv
64 * Merged changes from 4.24.1, including new localization system,
65 * timers, and allow use of underscore character in opcode names
66 *
67 * Revision 1.4 2004/09/27 05:52:31 jpff
68 * Minor coding
69 *
70 * Revision 1.3 2004/06/07 11:33:09 jpff
71 * line endings
72 *
73 * Revision 1.2 2004/05/31 15:53:06 jpff
74 * Removing warnings mainly
75 *
76 * Revision 1.1.1.1 2003/06/19 11:11:53 jpff
77 * Initial upload of sources
78 *
79 * Revision 1.2 2003/05/21 11:31:06 jpff
80 * Added Copyright notices
81 *
82 * Revision 3.4 1994/10/31 17:37:28 martin
83 * Starting with rcs
84 *
85 */
86 #include "csoundCore.h"
87 #include <math.h>
88 #include <assert.h>
89
90 static void fft_(CSOUND *,MYFLT *, MYFLT *, int32_t, int32_t, int32_t, int32_t);
91 static void fftmx(MYFLT *, MYFLT *, int32_t, int32_t, int32_t, int32_t, int32_t,
92 int32_t*, MYFLT *, MYFLT *, MYFLT *, MYFLT *, int32_t *,
93 int32_t[]);
94 static void reals_(CSOUND *,MYFLT *, MYFLT *, int32_t, int32_t);
95
96 /*
97 *-----------------------------------------------------------------------
98 * subroutine: fft
99 * multivariate complex fourier transform, computed in place
100 * using mixed-radix fast fourier transform algorithm.
101 *-----------------------------------------------------------------------
102 *
103 * this is the call from C:
104 * fft_(anal,banal,&one,&N2,&one,&mtwo);
105 * CHANGED TO:-
106 * fft_(csound,anal,banal,one,N2,one,mtwo);
107 */
108
fft_(CSOUND * csound,MYFLT * a,MYFLT * b,int32_t nseg,int32_t n,int32_t nspn,int32_t isn)109 static void fft_(CSOUND *csound, MYFLT *a, MYFLT *b,
110 int32_t nseg, int32_t n, int32_t nspn, int32_t isn)
111 /* *a, pointer to array 'anal' */
112 /* *b; pointer to array 'banal' */
113 {
114 int32_t nfac[32]; /* These are one bigger than needed */
115 /* because wish to use Fortran array */
116 /* index which runs 1 to n, not 0 to n */
117
118 int32_t m = 0,
119 k,
120 kt,
121 jj,
122 j,
123 nf,
124 ntot,
125 maxf, maxp=-1;
126
127 /* work space pointers */
128 void *buf;
129 MYFLT *at, *ck, *bt, *sk;
130 int32_t *np;
131
132 /* reduce the pointers to input arrays - by doing this, FFT uses FORTRAN
133 indexing but retains compatibility with C arrays */
134 a--; b--;
135
136 /*
137 * determine the factors of n
138 */
139 k = nf = abs(n);
140 if (nf==1)
141 return;
142
143 nspn = abs(nf*nspn);
144 ntot = abs(nspn*nseg);
145
146 for (m=0; !(k%16); nfac[++m]=4,k/=16);
147 assert(m<16);
148 for (j=3,jj=9; jj<=k; j+=2,jj=j*j)
149 for (; !(k%jj); nfac[++m]=j,k/=jj);
150
151 if (k<=4) {
152 kt = m;
153 nfac[m+1] = k;
154 if (k != 1)
155 m++;
156 }
157 else {
158 if (k%4==0) {
159 nfac[++m]=2;
160 k/=4;
161 }
162
163 kt = m;
164 maxp = (kt+kt+2 > k-1 ? kt+kt+2 : k-1);
165 for (j=2; j<=k; j=1+((j+1)/2)*2)
166 if (k%j==0) {
167 nfac[++m]=j;
168 k/=j;
169 }
170 }
171 if (m <= kt+1)
172 maxp = m + kt + 1;
173 if (UNLIKELY(m+kt > 15)) {
174 csound->Warning(csound, Str("\nerror - fft parameter n has "
175 "more than 15 factors : %d"), n);
176 return;
177 }
178 if (kt!=0) {
179 j = kt;
180 while (j)
181 nfac[++m]=nfac[j--];
182 }
183 maxf = nfac[m-kt];
184 if (kt > 0 && maxf <nfac[kt])
185 maxf = nfac[kt];
186
187 /* allocate workspace - assume no errors! */
188 buf = csound->Calloc(csound, sizeof(MYFLT) * 4 * maxf + sizeof(int32_t) * maxp);
189 at = (MYFLT*) buf;
190 ck = (MYFLT*) at + (int32_t) maxf;
191 bt = (MYFLT*) ck + (int32_t) maxf;
192 sk = (MYFLT*) bt + (int32_t) maxf;
193 np = (int32_t*) ((void*) ((MYFLT*) sk + (int32_t) maxf));
194
195 /* decrement pointers to allow FORTRAN type usage in fftmx */
196 at--; bt--; ck--; sk--; np--;
197
198 /* call fft driver */
199 fftmx(a, b, ntot, nf, nspn, isn, m, &kt, at, ck, bt, sk, np, nfac);
200
201 /* release working storage before returning - assume no problems */
202 csound->Free(csound,buf);
203 }
204
205 /*
206 *-----------------------------------------------------------------------
207 * subroutine: fftmx
208 * called by subroutine 'fft' to compute mixed-radix fourier transform
209 *-----------------------------------------------------------------------
210 */
211
fftmx(MYFLT * a,MYFLT * b,int32_t ntot,int32_t n,int32_t nspan,int32_t isn,int32_t m,int32_t * kt,MYFLT * at,MYFLT * ck,MYFLT * bt,MYFLT * sk,int32_t * np,int32_t nfac[])212 static void fftmx(MYFLT *a, MYFLT *b,
213 int32_t ntot, int32_t n, int32_t nspan, int32_t isn, int32_t m,
214 int32_t *kt, MYFLT *at, MYFLT *ck, MYFLT *bt, MYFLT *sk,
215 int32_t *np, int32_t nfac[])
216 {
217 int32_t i,inc,
218 j,jc,jf, jj,
219 k, k1, k2, k3=0, k4,
220 kk,klim,ks,kspan, kspnn,
221 lim,
222 maxf,mm,
223 nn,nt;
224 double aa, aj, ajm, ajp, ak, akm, akp,
225 bb, bj, bjm, bjp, bk, bkm, bkp,
226 c1, c2=0, c3=0, c72, cd,
227 dr,
228 rad,
229 sd, s1, s2=0, s3=0, s72, s120;
230
231 double xx; /****** ADDED APRIL 1991 *********/
232 inc=abs(isn);
233 nt = inc*ntot;
234 ks = inc*nspan;
235 /******************* REPLACED MARCH 29: ***********************
236 rad = atan(1.0);
237 **************************************************************/
238 rad = 0.785398163397448278900;
239 /******************* REPLACED MARCH 29: ***********************
240 s72 = rad/0.625;
241 c72 = cos(s72);
242 s72 = sin(s72);
243 **************************************************************/
244 c72 = 0.309016994374947451270;
245 s72 = 0.951056516295153531190;
246 /******************* REPLACED MARCH 29: ***********************
247 s120 = sqrt(0.75);
248 **************************************************************/
249 s120 = 0.866025403784438707600;
250
251 /* scale by 1/n for isn > 0 ( reverse transform ) */
252
253 if (isn < 0) {
254 s72 = -s72;
255 s120 = -s120;
256 rad = -rad;}
257 else {
258 ak = 1.0/(double)n;
259 for (j=1; j<=nt;j += inc) {
260 a[j] *= (MYFLT)ak;
261 b[j] *= (MYFLT)ak;
262 }
263 }
264 kspan = ks;
265 nn = nt - inc;
266 jc = ks/n;
267
268 /* sin, cos values are re-initialised each lim steps */
269
270 lim = 32;
271 klim = lim * jc;
272 i = 0;
273 jf = 0;
274 maxf = m - (*kt);
275 maxf = nfac[maxf];
276 if ((*kt) > 0 && maxf < nfac[*kt])
277 maxf = nfac[*kt];
278
279 /*
280 * compute Fourier transform
281 */
282
283 lbl40:
284 dr = (8.0 * (double)jc)/((double)kspan);
285 /*************************** APRIL 1991 POW & POW2 not WORKING.. REPLACE *******
286 cd = 2.0 * (pow2 ( sin(0.5 * dr * rad)) );
287 *******************************************************************************/
288 xx = sin(0.5 * dr * rad);
289 cd = 2.0 * xx * xx;
290 sd = sin(dr * rad);
291 kk = 1;
292 if (nfac[++i]!=2) goto lbl110;
293 /*
294 * transform for factor of 2 (including rotation factor)
295 */
296 kspan /= 2;
297 k1 = kspan + 2;
298 do {
299 do {
300 k2 = kk + kspan;
301 ak = a[k2];
302 bk = b[k2];
303 a[k2] = (a[kk]) - (MYFLT)ak;
304 b[k2] = (b[kk]) - (MYFLT)bk;
305 a[kk] = (a[kk]) + (MYFLT)ak;
306 b[kk] = (b[kk]) + (MYFLT)bk;
307 kk = k2 + kspan;
308 } while (kk <= nn);
309 kk -= nn;
310 } while (kk <= jc);
311 if (kk > kspan) goto lbl350;
312 lbl60:
313 c1 = 1.0 - cd;
314 s1 = sd;
315 mm = (k1/2 < klim ? k1/2 :klim);
316 goto lbl80;
317 lbl70:
318 ak = c1 - ((cd*c1)+(sd*s1));
319 s1 = ((sd*c1)-(cd*s1)) + s1;
320 c1 = ak;
321 lbl80:
322 do {
323 do {
324 k2 = kk + kspan;
325 ak = a[kk] - a[k2];
326 bk = b[kk] - b[k2];
327 a[kk] = a[kk] + a[k2];
328 b[kk] = b[kk] + b[k2];
329 a[k2] = (MYFLT)((c1 * ak) - (s1 * bk));
330 b[k2] = (MYFLT)((s1 * ak) + (c1 * bk));
331 kk = k2 + kspan;
332 } while (kk < nt);
333 k2 = kk - nt;
334 c1 = -c1;
335 kk = k1 - k2;
336 } while (kk > k2);
337 kk += jc;
338 if (kk <= mm) goto lbl70;
339 if (kk < k2) goto lbl90;
340 k1 += (inc + inc);
341 kk = ((k1-kspan)/2) + jc;
342 if (kk <= (jc+jc)) goto lbl60;
343 goto lbl40;
344 lbl90:
345 s1 = ((double)((kk-1)/jc)) * dr * rad;
346 c1 = cos(s1);
347 s1 = sin(s1);
348 mm = (k1/2 < mm+klim ? k1/2 : mm+klim);
349 goto lbl80;
350 /*
351 * transform for factor of 3 (optional code)
352 */
353
354 lbl100:
355 k1 = kk + kspan;
356 k2 = k1 + kspan;
357 ak = a[kk];
358 bk = b[kk];
359 aj = a[k1] + a[k2];
360 bj = b[k1] + b[k2];
361 a[kk] = (MYFLT)(ak + aj);
362 b[kk] = (MYFLT)(bk + bj);
363 ak += (-0.5 * aj);
364 bk += (-0.5 * bj);
365 aj = (a[k1] - a[k2]) * s120;
366 bj = (b[k1] - b[k2]) * s120;
367 a[k1] = (MYFLT)(ak - bj);
368 b[k1] = (MYFLT)(bk + aj);
369 a[k2] = (MYFLT)(ak + bj);
370 b[k2] = (MYFLT)(bk - aj);
371 kk = k2 + kspan;
372 if (kk < nn) goto lbl100;
373 kk -= nn;
374 if (kk <= kspan) goto lbl100;
375 goto lbl290;
376
377 /*
378 * transform for factor of 4
379 */
380
381 lbl110:
382 if (nfac[i] != 4) goto lbl230;
383 kspnn = kspan;
384 kspan = kspan/4;
385 lbl120:
386 c1 = 1.0;
387 s1 = 0;
388 mm = (kspan < klim ? kspan : klim);
389 goto lbl150;
390 lbl130:
391 c2 = c1 - ((cd*c1)+(sd*s1));
392 s1 = ((sd*c1)-(cd*s1)) + s1;
393 /*
394 * the following three statements compensate for truncation
395 * error. if rounded arithmetic is used, substitute
396 * c1=c2
397 *
398 * c1 = (0.5/(pow2(c2)+pow2(s1))) + 0.5;
399 * s1 = c1*s1;
400 * c1 = c1*c2;
401 */
402 c1 = c2;
403 lbl140:
404 c2 = (c1 * c1) - (s1 * s1);
405 s2 = c1 * s1 * 2.0;
406 c3 = (c2 * c1) - (s2 * s1);
407 s3 = (c2 * s1) + (s2 * c1);
408 lbl150:
409 k1 = kk + kspan;
410 k2 = k1 + kspan;
411 k3 = k2 + kspan;
412 akp = a[kk] + a[k2];
413 akm = a[kk] - a[k2];
414 ajp = a[k1] + a[k3];
415 ajm = a[k1] - a[k3];
416 a[kk] = (MYFLT)(akp + ajp);
417 ajp = akp - ajp;
418 bkp = b[kk] + b[k2];
419 bkm = b[kk] - b[k2];
420 bjp = b[k1] + b[k3];
421 bjm = b[k1] - b[k3];
422 b[kk] = (MYFLT)(bkp + bjp);
423 bjp = bkp - bjp;
424 if (isn < 0) goto lbl180;
425 akp = akm - bjm;
426 akm = akm + bjm;
427 bkp = bkm + ajm;
428 bkm = bkm - ajm;
429 if (s1 == 0.0) goto lbl190;
430 lbl160:
431 a[k1] = (MYFLT)((akp*c1) - (bkp*s1));
432 b[k1] = (MYFLT)((akp*s1) + (bkp*c1));
433 a[k2] = (MYFLT)((ajp*c2) - (bjp*s2));
434 b[k2] = (MYFLT)((ajp*s2) + (bjp*c2));
435 a[k3] = (MYFLT)((akm*c3) - (bkm*s3));
436 b[k3] = (MYFLT)((akm*s3) + (bkm*c3));
437 kk = k3 + kspan;
438 if (kk <= nt) goto lbl150;
439 lbl170:
440 kk -= (nt - jc);
441 if (kk <= mm) goto lbl130;
442 if (kk < kspan) goto lbl200;
443 kk -= (kspan - inc);
444 if (kk <= jc) goto lbl120;
445 if (kspan==jc) goto lbl350;
446 goto lbl40;
447 lbl180:
448 akp = akm + bjm;
449 akm = akm - bjm;
450 bkp = bkm - ajm;
451 bkm = bkm + ajm;
452 if (s1 != 0.0) goto lbl160;
453 lbl190:
454 a[k1] = (MYFLT)akp;
455 b[k1] = (MYFLT)bkp;
456 a[k2] = (MYFLT)ajp;
457 b[k2] = (MYFLT)bjp;
458 a[k3] = (MYFLT)akm;
459 b[k3] = (MYFLT)bkm;
460 kk = k3 + kspan;
461 if (kk <= nt) goto lbl150;
462 goto lbl170;
463 lbl200:
464 s1 = ((double)((kk-1)/jc)) * dr * rad;
465 c1 = cos(s1);
466 s1 = sin(s1);
467 mm = (kspan < mm+klim ? kspan : mm+klim);
468 goto lbl140;
469
470 /*
471 * transform for factor of 5 (optional code)
472 */
473
474 lbl210:
475 c2 = (c72*c72) - (s72*s72);
476 s2 = 2.0 * c72 * s72;
477 lbl220:
478 k1 = kk + kspan;
479 k2 = k1 + kspan;
480 k3 = k2 + kspan;
481 k4 = k3 + kspan;
482 akp = a[k1] + a[k4];
483 akm = a[k1] - a[k4];
484 bkp = b[k1] + b[k4];
485 bkm = b[k1] - b[k4];
486 ajp = a[k2] + a[k3];
487 ajm = a[k2] - a[k3];
488 bjp = b[k2] + b[k3];
489 bjm = b[k2] - b[k3];
490 aa = a[kk];
491 bb = b[kk];
492 a[kk] = (MYFLT)(aa + akp + ajp);
493 b[kk] = (MYFLT)(bb + bkp + bjp);
494 ak = (akp*c72) + (ajp*c2) + aa;
495 bk = (bkp*c72) + (bjp*c2) + bb;
496 aj = (akm*s72) + (ajm*s2);
497 bj = (bkm*s72) + (bjm*s2);
498 a[k1] = (MYFLT)(ak - bj);
499 a[k4] = (MYFLT)(ak + bj);
500 b[k1] = (MYFLT)(bk + aj);
501 b[k4] = (MYFLT)(bk - aj);
502 ak = (akp*c2) + (ajp*c72) + aa;
503 bk = (bkp*c2) + (bjp*c72) + bb;
504 aj = (akm*s2) - (ajm*s72);
505 bj = (bkm*s2) - (bjm*s72);
506 a[k2] = (MYFLT)(ak - bj);
507 a[k3] = (MYFLT)(ak + bj);
508 b[k2] = (MYFLT)(bk + aj);
509 b[k3] = (MYFLT)(bk - aj);
510 kk = k4 + kspan;
511 if (kk < nn) goto lbl220;
512 kk -= nn;
513 if (kk <= kspan) goto lbl220;
514 goto lbl290;
515
516 /*
517 * transform for odd factors
518 */
519
520 lbl230:
521 k = nfac[i];
522 kspnn = kspan;
523 kspan /= k;
524 if (k==3) goto lbl100;
525 if (k==5) goto lbl210;
526 if (k==jf) goto lbl250;
527 jf = k;
528 s1 = rad/(((double)(k))/8.0);
529 c1 = cos(s1);
530 s1 = sin(s1);
531 ck[jf] = FL(1.0);
532 sk[jf] = FL(0.0);
533 for (j=1; j<k ; j++) {
534 ck[j] = (MYFLT)((ck[k])*c1 + (sk[k])*s1);
535 sk[j] = (MYFLT)((ck[k])*s1 - (sk[k])*c1);
536 k--;
537 ck[k] = ck[j];
538 sk[k] = -(sk[j]);
539 }
540 lbl250:
541 k1 = kk;
542 k2 = kk + kspnn;
543 aa = a[kk];
544 bb = b[kk];
545 ak = aa;
546 bk = bb;
547 j = 1;
548 k1 += kspan;
549 do {
550 k2 -= kspan;
551 j++;
552 at[j] = a[k1] + a[k2];
553 ak = at[j] + ak;
554 bt[j] = b[k1] + b[k2];
555 bk = bt[j] + bk;
556 j++;
557 at[j] = a[k1] - a[k2];
558 bt[j] = b[k1] - b[k2];
559 k1 += kspan;
560 } while (k1 < k2);
561 a[kk] = (MYFLT)ak;
562 b[kk] = (MYFLT)bk;
563 k1 = kk;
564 k2 = kk + kspnn;
565 j = 1;
566 lbl270:
567 k1 += kspan;
568 k2 -= kspan;
569 jj = j;
570 ak = aa;
571 bk = bb;
572 aj = 0.0;
573 bj = 0.0;
574 k = 1;
575 do {
576 k++;
577 ak = (at[k] * ck[jj]) + ak;
578 bk = (bt[k] * ck[jj]) + bk;
579 k++;
580 aj = (at[k] * sk[jj]) + aj;
581 bj = (bt[k] * sk[jj]) + bj;
582 jj += j;
583 if (jj > jf)
584 jj -= jf;
585 } while (k < jf);
586 k = jf - j;
587 a[k1] = (MYFLT)(ak - bj);
588 b[k1] = (MYFLT)(bk + aj);
589 a[k2] = (MYFLT)(ak + bj);
590 b[k2] = (MYFLT)(bk - aj);
591 j++;
592 if (j < k) goto lbl270;
593 kk += kspnn;
594 if (kk <= nn) goto lbl250;
595 kk -= nn;
596 if (kk<=kspan) goto lbl250;
597
598 /*
599 * multiply by rotation factor (except for factors of 2 and 4)
600 */
601
602 lbl290:
603 if (i==m) goto lbl350;
604 kk = jc + 1;
605 lbl300:
606 c2 = 1.0 - cd;
607 s1 = sd;
608 mm = (kspan < klim ? kspan : klim);
609 goto lbl320;
610 lbl310:
611 c2 = c1 - ((cd*c1) + (sd*s1));
612 s1 = s1 + ((sd*c1) - (cd*s1));
613 lbl320:
614 c1 = c2;
615 s2 = s1;
616 kk += kspan;
617 lbl330:
618 ak = a[kk];
619 a[kk] = (MYFLT)((c2*ak) - (s2 * b[kk]));
620 b[kk] = (MYFLT)((s2*ak) + (c2 * b[kk]));
621 kk += kspnn;
622 if (kk <= nt) goto lbl330;
623 ak = s1*s2;
624 s2 = (s1*c2) + (c1*s2);
625 c2 = (c1*c2) - ak;
626 kk -= (nt - kspan);
627 if (kk <= kspnn) goto lbl330;
628 kk -= (kspnn - jc);
629 if (kk <= mm) goto lbl310;
630 if (kk < kspan) goto lbl340;
631 kk -= (kspan - jc - inc);
632 if (kk <= (jc+jc)) goto lbl300;
633 goto lbl40;
634 lbl340:
635 s1 = ((double)((kk-1)/jc)) * dr * rad;
636 c2 = cos(s1);
637 s1 = sin(s1);
638 mm = (kspan < mm+klim ? kspan :mm+klim);
639 goto lbl320;
640
641 /*
642 * permute the results to normal order---done in two stages
643 * permutation for square factors of n
644 */
645
646 lbl350:
647 np[1] = ks;
648 if (!(*kt)) goto lbl440;
649 k = *kt + *kt + 1;
650 if (m < k)
651 k--;
652 np[k+1] = jc;
653 for (j=1; j < k; j++,k--) {
654 np[j+1] = np[j] / nfac[j];
655 np[k] = np[k+1] * nfac[j];
656 }
657 k3 = np[k+1];
658 kspan = np[2];
659 kk = jc + 1;
660 k2 = kspan + 1;
661 j = 1;
662 if (n != ntot) goto lbl400;
663 /*
664 * permutation for single-variate transform (optional code)
665 */
666 lbl370:
667 do {
668 ak = a[kk];
669 a[kk] = a[k2];
670 a[k2] = (MYFLT)ak;
671 bk = b[kk];
672 b[kk] = b[k2];
673 b[k2] = (MYFLT)bk;
674 kk += inc;
675 k2 += kspan;
676 } while (k2 < ks);
677 lbl380:
678 do {
679 k2 -= np[j++];
680 k2 += np[j+1];
681 } while (k2 > np[j]);
682 j = 1;
683 lbl390:
684 if (kk < k2) {
685 goto lbl370;
686 }
687 kk += inc;
688 k2 += kspan;
689 if (k2 < ks) goto lbl390;
690 if (kk < ks) goto lbl380;
691 jc = k3;
692 goto lbl440;
693 /*
694 * permutation for multivariate transform
695 */
696 lbl400:
697 do {
698 do {
699 k = kk + jc;
700 do {
701 ak = a[kk];
702 a[kk] = a[k2];
703 a[k2] = (MYFLT)ak;
704 bk = b[kk];
705 b[kk] = b[k2];
706 b[k2] = (MYFLT)bk;
707 kk += inc;
708 k2 += inc;
709 } while (kk < k);
710 kk += (ks - jc);
711 k2 += (ks - jc);
712 } while (kk < nt);
713 k2 -= (nt - kspan);
714 kk -= (nt - jc);
715 } while (k2 < ks);
716 lbl420:
717 do {
718 k2 -= np[j++];
719 k2 += np[j+1];
720 } while (k2 > np[j]);
721 j = 1;
722 lbl430:
723 if (kk < k2) goto lbl400;
724 kk += jc;
725 k2 += kspan;
726 if (k2 < ks) goto lbl430;
727 if (kk < ks) goto lbl420;
728 jc = k3;
729 lbl440:
730 if ((2*(*kt))+1 >= m)
731 return;
732
733 kspnn = *(np + *(kt) + 1);
734 j = m - *kt;
735 nfac[j+1] = 1;
736 lbl450:
737 nfac[j] = nfac[j] * nfac[j+1];
738 j--;
739 if (j != *kt) goto lbl450;
740 *kt = *(kt) + 1;
741 nn = nfac[*kt] - 1;
742 jj = 0;
743 j = 0;
744 goto lbl480;
745 lbl460:
746 jj -= k2;
747 k2 = kk;
748 kk = nfac[++k];
749 lbl470:
750 jj += kk;
751 if (jj >= k2) goto lbl460;
752 np[j] = jj;
753 lbl480:
754 k2 = nfac[*kt];
755 k = *kt + 1;
756 kk = nfac[k];
757 j++;
758 if (j <= nn) goto lbl470;
759 /* Determine permutation cycles of length greater than 1 */
760 j = 0;
761 goto lbl500;
762 lbl490:
763 k = kk;
764 kk = np[k];
765 np[k] = -kk;
766 if (kk != j) goto lbl490;
767 k3 = kk;
768 lbl500:
769 kk = np[++j];
770 if (kk < 0) goto lbl500;
771 if (kk != j) goto lbl490;
772 np[j] = -j;
773 if (j != nn) goto lbl500;
774 maxf *= inc;
775 /* Perform reordering following permutation cycles */
776 goto lbl570;
777 lbl510:
778 j--;
779 if (np[j] < 0) goto lbl510;
780 jj = jc;
781 lbl520:
782 kspan = jj;
783 if (jj > maxf)
784 kspan = maxf;
785 jj -= kspan;
786 k = np[j];
787 kk = (jc*k) + i + jj;
788 k1 = kk + kspan;
789 k2 = 0;
790 lbl530:
791 k2++;
792 at[k2] = a[k1];
793 bt[k2] = b[k1];
794 k1 -= inc;
795 if (k1 != kk) goto lbl530;
796 lbl540:
797 k1 = kk + kspan;
798 k2 = k1 - (jc * (k + np[k]));
799 k = -(np[k]);
800 lbl550:
801 a[k1] = a[k2];
802 b[k1] = b[k2];
803 k1 -= inc;
804 k2 -= inc;
805 if (k1 != kk) goto lbl550;
806 kk = k2;
807 if (k != j) goto lbl540;
808 k1 = kk + kspan;
809 k2 = 0;
810 lbl560:
811 k2++;
812 a[k1] = at[k2];
813 b[k1] = bt[k2];
814 k1 -= inc;
815 if (k1 != kk) goto lbl560;
816 if (jj) goto lbl520;
817 if (j != 1) goto lbl510;
818 lbl570:
819 j = k3 + 1;
820 nt -= kspnn;
821 i = nt - inc + 1;
822 if (nt >= 0) goto lbl510;
823 return;
824 }
825
826 /*
827 *-----------------------------------------------------------------------
828 * subroutine:
829 * reals
830 * used with 'fft' to compute fourier transform or inverse for real data
831 *-----------------------------------------------------------------------
832 * this is the call from C:
833 *
834 * reals_(anal,banal,N2,mtwo);
835 * which has been changed from CARL call
836 * reals_(csound,anal,banal,&N2,&mtwo);
837 */
838
reals_(CSOUND * csound,MYFLT * a,MYFLT * b,int32_t n,int32_t isn)839 static void reals_(CSOUND *csound, MYFLT *a, MYFLT *b, int32_t n, int32_t isn)
840
841 /* *a, a refers to an array of floats 'anal' */
842 /* *b; b refers to an array of floats 'banal' */
843 /* See IEEE book for a long comment here on usage */
844
845 {
846 IGN(csound);
847 int32_t inc,
848 j,
849 k,
850 lim,
851 mm,ml,
852 nf,nk,nh;
853
854 double aa,ab,
855 ba,bb,
856 cd,cn,
857 dr,
858 em,
859 rad,re,
860 sd,sn;
861 double xx; /******* ADDED APRIL 1991 ******/
862 /* adjust input array pointers (called from C) */
863 a--; b--;
864 inc = abs(isn);
865 nf = abs(n);
866 nk = (nf*inc) + 2;
867 nh = nk/2;
868 /*****************************
869 rad = atan((double)1.0);
870 ******************************/
871 rad = 0.785398163397448278900;
872 dr = -4.0/(double)(nf);
873 /********************************** POW2 REMOVED APRIL 1991 *****************
874 cd = 2.0 * (pow2(sin((double)0.5 * dr * rad)));
875 *****************************************************************************/
876 xx = sin((double)0.5 * dr * rad);
877 cd = 2.0 * xx * xx;
878 sd = sin(dr * rad);
879 /*
880 * sin,cos values are re-initialised each lim steps
881 */
882 lim = 32;
883 mm = lim;
884 ml = 0;
885 sn = 0.0;
886 if (isn<0) {
887 cn = 1.0;
888 a[nk-1] = a[1];
889 b[nk-1] = b[1]; }
890 else {
891 cn = -1.0;
892 sd = -sd;
893 }
894 for (j=1;j<=nh;j+=inc) {
895 k = nk - j;
896 aa = a[j] + a[k];
897 ab = a[j] - a[k];
898 ba = b[j] + b[k];
899 bb = b[j] - b[k];
900 re = (cn*ba) + (sn*ab);
901 em = (sn*ba) - (cn*ab);
902 b[k] = (MYFLT)((em-bb)*0.5);
903 b[j] = (MYFLT)((em+bb)*0.5);
904 a[k] = (MYFLT)((aa-re)*0.5);
905 a[j] = (MYFLT)((aa+re)*0.5);
906 ml++;
907 if (ml!=mm) {
908 aa = cn - ((cd*cn)+(sd*sn));
909 sn = ((sd*cn) - (cd*sn)) + sn;
910 cn = aa;}
911 else {
912 mm +=lim;
913 sn = ((MYFLT)ml) * dr * rad;
914 cn = cos(sn);
915 if (isn>0)
916 cn = -cn;
917 sn = sin(sn);
918 }
919 }
920 return;
921 }
922
923 /**
924 * Compute in-place real FFT, allowing non power of two FFT sizes.
925 *
926 * buf: array of FFTsize + 2 MYFLT values; output is in interleaved
927 * real/imaginary format (note: the real part of the Nyquist
928 * frequency is stored in buf[FFTsize], and not in buf[1]).
929 * FFTsize: FFT length in samples; not required to be an integer power of two,
930 * but should be even and not have too many factors.
931 */
csoundRealFFTnp2(CSOUND * csound,MYFLT * buf,int32_t FFTsize)932 void csoundRealFFTnp2(CSOUND *csound, MYFLT *buf, int32_t FFTsize)
933 {
934 if (!(FFTsize & (FFTsize - 1))) {
935 /* if FFT size is power of two: */
936 csound->RealFFT(csound, buf, FFTsize);
937 buf[FFTsize] = buf[1];
938 }
939 else {
940 if (UNLIKELY(FFTsize < 2 || (FFTsize & 1))) {
941 csound->Warning(csound,
942 Str("csoundRealFFTnp2(): invalid FFT size, %d"), FFTsize);
943 return;
944 }
945 buf[FFTsize] = buf[FFTsize + 1] = FL(0.0);
946 fft_(csound, buf, &(buf[1]), 1, (FFTsize >> 1), 1, -2);
947 reals_(csound, buf, &(buf[1]), (FFTsize >> 1), -2);
948 }
949 buf[1] = buf[FFTsize + 1] = FL(0.0);
950 }
951
952 /**
953 * Compute in-place inverse real FFT, allowing non power of two FFT sizes.
954 * The output does not need to be scaled.
955 *
956 * buf: array of FFTsize + 2 MYFLT values, in interleaved real/imaginary
957 * format (note: the real part of the Nyquist frequency is stored
958 * in buf[FFTsize], and not in buf[1]).
959 * FFTsize: FFT length in samples; not required to be an integer power of two,
960 * but should be even and not have too many factors.
961 */
csoundInverseRealFFTnp2(CSOUND * csound,MYFLT * buf,int32_t FFTsize)962 void csoundInverseRealFFTnp2(CSOUND *csound, MYFLT *buf, int32_t FFTsize)
963 {
964 if (UNLIKELY(FFTsize < 2 || (FFTsize & 1))){
965 csound->Warning(csound, Str("csoundInverseRealFFTnp2(): invalid FFT size"));
966 return;
967 }
968 buf[1] = buf[FFTsize + 1] = FL(0.0);
969 reals_(csound, buf, &(buf[1]), (FFTsize >> 1), 2);
970 fft_(csound, buf, &(buf[1]), 1, (FFTsize >> 1), 1, 2);
971 buf[FFTsize] = buf[FFTsize + 1] = FL(0.0);
972 }
973
csoundInverseComplexFFTnp2(CSOUND * csound,MYFLT * buf,int32_t FFTsize)974 void csoundInverseComplexFFTnp2(CSOUND *csound, MYFLT *buf, int32_t FFTsize)
975 {
976 if (UNLIKELY(FFTsize < 2 || (FFTsize & 1))){
977 csound->Warning(csound, Str("csoundInverseRealFFTnp2(): invalid FFT size"));
978 return;
979 }
980 fft_(csound, buf, buf, 1, FFTsize, 1, 2);
981 }
982
csoundComplexFFTnp2(CSOUND * csound,MYFLT * buf,int32_t FFTsize)983 void csoundComplexFFTnp2(CSOUND *csound, MYFLT *buf, int32_t FFTsize)
984 {
985 if (UNLIKELY(FFTsize < 2 || (FFTsize & 1))) {
986 csound->Warning(csound, Str("csoundRealFFTnp2(): invalid FFT size"));
987 return;
988 }
989 fft_(csound, buf, buf, 1, FFTsize, 1, -2);
990 }
991