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