1 /*
2
3 FFT analysis and phase vocoder UGens for SuperCollider, by Dan Stowell.
4 (c) Dan Stowell 2006-2010.
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19
20 */
21
22 #include "SC_fftlib.h"
23 #include "SC_PlugIn.h"
24 #include "FFT_UGens.h"
25 #include <stdio.h>
26
27 // Used by PV_MagLog
28 #define SMALLEST_NUM_FOR_LOG 2e-42
29
30 #define PI 3.1415926535898f
31 #define MPI -3.1415926535898f
32 #define TWOPI 6.28318530717952646f
33 #define THREEPI 9.4247779607694f
34
35 /* Rewrap phase into +-pi domain: essentially mod(phase+pi,-2pi)+pi */
36 #define PHASE_REWRAP(phase) ((phase) + TWOPI * (1.f + floorf(-((phase)+PI)/TWOPI)))
37
38 //////////////////////////////////////////////////////////////////////////////////////////////////
39
40
41 struct FFTAnalyser_Unit : Unit
42 {
43 float outval;
44
45 // Not always used: multipliers which convert from bin indices to freq vals, and vice versa.
46 // See also the macros for deriving these.
47 float m_bintofreq, m_freqtobin;
48 };
49
50 struct FFTAnalyser_OutOfPlace : FFTAnalyser_Unit
51 {
52 int m_numbins;
53 float *m_tempbuf;
54 };
55
56 struct PV_MagSubtract : Unit {
57 };
58
59 struct FFTFlux_Unit : FFTAnalyser_OutOfPlace
60 {
61 float m_yesternorm;
62 float m_yesterdc;
63 float m_yesternyq;
64 bool m_normalise;
65 };
66
67 struct FFTPower : FFTAnalyser_Unit
68 {
69 float m_normfactor;
70 bool m_square;
71 };
72
73 struct FFTSubbandPower : FFTAnalyser_Unit
74 {
75 float m_normfactor;
76 bool m_square;
77
78 int m_numbands;
79 int *m_cutoffs; // Will hold bin indices corresponding to frequencies
80 float *m_outvals;
81 bool m_cutoff_inited;
82
83 int m_scalemode;
84 };
85
86 struct FFTPhaseDev : FFTAnalyser_OutOfPlace
87 {
88 bool m_weight; // Whether or not to weight the phase deviations according to the magnitudes
89 };
90
91 struct FFTComplexDev : FFTAnalyser_OutOfPlace
92 {
93 bool m_rectify; // Whether or not to ignore bins whose power is decreasing
94 };
95
96 struct FFTMKL : FFTAnalyser_OutOfPlace
97 {
98 };
99
100 struct PV_Whiten : Unit {
101 float m_relaxcoef, m_floor, m_smear;
102 int m_bindownsample;
103 };
104
105
106 struct FFTSubbandFlatness : FFTAnalyser_Unit
107 {
108 int m_numbands;
109 int *m_cutoffs; // Will hold bin indices corresponding to frequencies
110 float *m_outvals;
111 bool m_cutoff_inited;
112 };
113 struct FFTCrest : FFTAnalyser_Unit
114 {
115 int m_frombin; // Will hold bin index
116 int m_tobinp1; // Will hold bin index
117 bool m_cutoffneedsinit;
118 };
119 struct FFTSpread : FFTAnalyser_Unit
120 {
121 };
122 struct FFTSlope : FFTAnalyser_Unit
123 {
124 };
125 struct FFTPeak : FFTAnalyser_Unit
126 {
127 float outval2, maxfreq, minfreq;
128 int maxbin, minbin;
129 };
130 struct PV_MagSmooth : Unit {
131 float *m_memory;
132 };
133 struct PV_ExtractRepeat : Unit {
134 float* m_logmags;
135 int m_cursor;
136 float m_fbufnum;
137 SndBuf* m_buf;
138 };
139
140
141 // for operation on one buffer
142 // just like PV_GET_BUF except it outputs unit->outval rather than -1 when FFT not triggered
143 #define FFTAnalyser_GET_BUF \
144 float fbufnum = ZIN0(0); \
145 if (fbufnum < 0.f) { ZOUT0(0) = unit->outval; return; } \
146 ZOUT0(0) = fbufnum; \
147 uint32 ibufnum = (uint32)fbufnum; \
148 World *world = unit->mWorld; \
149 SndBuf *buf; \
150 if (ibufnum >= world->mNumSndBufs) { \
151 int localBufNum = ibufnum - world->mNumSndBufs; \
152 Graph *parent = unit->mParent; \
153 if(localBufNum <= parent->localBufNum) { \
154 buf = parent->mLocalSndBufs + localBufNum; \
155 } else { \
156 buf = world->mSndBufs; \
157 } \
158 } else { \
159 buf = world->mSndBufs + ibufnum; \
160 } \
161 int numbins = (buf->samples - 2) >> 1;
162
163 // Same as above; but with "output2" as well as "output"
164 #define FFTAnalyser_GET_BUF_TWOOUTS \
165 float fbufnum = ZIN0(0); \
166 if (fbufnum < 0.f) { ZOUT0(0) = unit->outval; ZOUT0(1) = unit->outval2; return; } \
167 ZOUT0(0) = fbufnum; \
168 uint32 ibufnum = (uint32)fbufnum; \
169 World *world = unit->mWorld; \
170 SndBuf *buf; \
171 if (ibufnum >= world->mNumSndBufs) { \
172 int localBufNum = ibufnum - world->mNumSndBufs; \
173 Graph *parent = unit->mParent; \
174 if(localBufNum <= parent->localBufNum) { \
175 buf = parent->mLocalSndBufs + localBufNum; \
176 } else { \
177 buf = world->mSndBufs; \
178 } \
179 } else { \
180 buf = world->mSndBufs + ibufnum; \
181 } \
182 int numbins = (buf->samples - 2) >> 1;
183
184 // As above; but for operation on two input buffers
185 #define FFTAnalyser_GET_BUF2 \
186 float fbufnum1 = ZIN0(0); \
187 float fbufnum2 = ZIN0(1); \
188 if (fbufnum1 < 0.f || fbufnum2 < 0.f) { ZOUT0(0) = unit->outval; return; } \
189 uint32 ibufnum1 = (int)fbufnum1; \
190 uint32 ibufnum2 = (int)fbufnum2; \
191 World *world = unit->mWorld; \
192 SndBuf *buf1; \
193 if (ibufnum1 >= world->mNumSndBufs) { \
194 int localBufNum = ibufnum1 - world->mNumSndBufs; \
195 Graph *parent = unit->mParent; \
196 if(localBufNum <= parent->localBufNum) { \
197 buf1 = parent->mLocalSndBufs + localBufNum; \
198 } else { \
199 buf1 = world->mSndBufs; \
200 } \
201 } else { \
202 buf1 = world->mSndBufs + ibufnum1; \
203 } \
204 SndBuf *buf2; \
205 if (ibufnum2 >= world->mNumSndBufs) { \
206 int localBufNum = ibufnum2 - world->mNumSndBufs; \
207 Graph *parent = unit->mParent; \
208 if(localBufNum <= parent->localBufNum) { \
209 buf2 = parent->mLocalSndBufs + localBufNum; \
210 } else { \
211 buf2 = world->mSndBufs; \
212 } \
213 } else { \
214 buf2 = world->mSndBufs + ibufnum2; \
215 } \
216 if (buf1->samples != buf2->samples) return; \
217 int numbins = (buf1->samples - 2) >> 1;
218
219 // Copied from FFT_UGens.cpp
220 #define MAKE_TEMP_BUF \
221 if (!unit->m_tempbuf) { \
222 unit->m_tempbuf = (float*)RTAlloc(unit->mWorld, buf->samples * sizeof(float)); \
223 unit->m_numbins = numbins; \
224 } else if (numbins != unit->m_numbins) return;
225
226 #define GET_BINTOFREQ \
227 if(unit->m_bintofreq==0.f){ \
228 unit->m_bintofreq = world->mFullRate.mSampleRate / buf->samples; \
229 } \
230 float bintofreq = unit->m_bintofreq;
231 #define GET_FREQTOBIN \
232 if(unit->m_freqtobin==0.f){ \
233 unit->m_freqtobin = buf->samples / world->mFullRate.mSampleRate; \
234 } \
235 float freqtobin = unit->m_freqtobin;
236
237 //////////////////////////////////////////////////////////////////////////////////////////////////
238
239 extern "C"
240 {
241 void FFTPower_Ctor(FFTPower *unit);
242 void FFTPower_next(FFTPower *unit, int inNumSamples);
243
244 void FFTFlux_Ctor(FFTFlux_Unit *unit);
245 void FFTFlux_next(FFTFlux_Unit *unit, int inNumSamples);
246 void FFTFlux_Dtor(FFTFlux_Unit *unit);
247 void FFTFluxPos_Ctor(FFTFlux_Unit *unit);
248 void FFTFluxPos_next(FFTFlux_Unit *unit, int inNumSamples);
249 void FFTFluxPos_Dtor(FFTFlux_Unit *unit);
250
251 void FFTDiffMags_Ctor(FFTAnalyser_Unit *unit);
252 void FFTDiffMags_next(FFTAnalyser_Unit *unit, int inNumSamples);
253
254 // void PV_DiffAndToDC_Ctor(PV_Unit *unit);
255 // void PV_DiffAndToDC_next(PV_Unit *unit, int inNumSamples);
256
257 void PV_MagSubtract_Ctor(PV_Unit *unit);
258 void PV_MagSubtract_next(PV_Unit *unit, int inNumSamples);
259
260 void FFTSubbandPower_Ctor(FFTSubbandPower *unit);
261 void FFTSubbandPower_next(FFTSubbandPower *unit, int inNumSamples);
262 void FFTSubbandPower_Dtor(FFTSubbandPower *unit);
263
264 void PV_MagLog_Ctor(PV_Unit *unit);
265 void PV_MagLog_next(PV_Unit *unit, int inNumSamples);
266
267 void PV_MagExp_Ctor(PV_Unit *unit);
268 void PV_MagExp_next(PV_Unit *unit, int inNumSamples);
269
270 void FFTPhaseDev_Ctor(FFTPhaseDev *unit);
271 void FFTPhaseDev_Dtor(FFTPhaseDev *unit);
272 void FFTPhaseDev_next(FFTPhaseDev *unit, int inNumSamples);
273
274 void FFTComplexDev_Ctor(FFTComplexDev *unit);
275 void FFTComplexDev_Dtor(FFTComplexDev *unit);
276 void FFTComplexDev_next(FFTComplexDev *unit, int inNumSamples);
277
278 void FFTMKL_Ctor(FFTMKL *unit);
279 void FFTMKL_Dtor(FFTMKL *unit);
280 void FFTMKL_next(FFTMKL *unit, int inNumSamples);
281
282 void PV_Whiten_Ctor(PV_Whiten *unit);
283 void PV_Whiten_next(PV_Whiten *unit, int inNumSamples);
284
285 void FFTSubbandFlatness_Ctor(FFTSubbandFlatness *unit);
286 void FFTSubbandFlatness_next(FFTSubbandFlatness *unit, int inNumSamples);
287 void FFTSubbandFlatness_Dtor(FFTSubbandFlatness *unit);
288
289 void FFTCrest_Ctor(FFTCrest *unit);
290 void FFTCrest_next(FFTCrest *unit, int inNumSamples);
291
292 void FFTSpread_Ctor(FFTSpread *unit);
293 void FFTSpread_next(FFTSpread *unit, int inNumSamples);
294
295 void FFTSlope_Ctor(FFTSlope *unit);
296 void FFTSlope_next(FFTSlope *unit, int inNumSamples);
297
298 void FFTPeak_Ctor(FFTPeak *unit);
299 void FFTPeak_next(FFTPeak *unit, int inNumSamples);
300
301 void PV_MagSmooth_Ctor(PV_MagSmooth *unit);
302 void PV_MagSmooth_next(PV_MagSmooth *unit, int inNumSamples);
303 void PV_MagSmooth_Dtor(PV_MagSmooth *unit);
304
305 void PV_MagMulAdd_Ctor(PV_Unit *unit);
306 void PV_MagMulAdd_next(PV_Unit *unit, int inNumSamples);
307
308 void PV_ExtractRepeat_Ctor(PV_ExtractRepeat *unit);
309 void PV_ExtractRepeat_next(PV_ExtractRepeat *unit, int inNumSamples);
310 void PV_ExtractRepeat_Dtor(PV_ExtractRepeat *unit);
311
312 }
313
314 InterfaceTable *ft;
315
316 void init_SCComplex(InterfaceTable *inTable);
317
318 //////////////////////////////////////////////////////////////////////////////////////////////////
319
FFTPower_Ctor(FFTPower * unit)320 void FFTPower_Ctor(FFTPower *unit)
321 {
322 SETCALC(FFTPower_next);
323 ZOUT0(0) = unit->outval = 0.;
324
325 unit->m_square = ZIN0(1) > 0.f;
326 unit->m_normfactor = 0.f;
327 }
328
FFTPower_next(FFTPower * unit,int inNumSamples)329 void FFTPower_next(FFTPower *unit, int inNumSamples)
330 {
331 FFTAnalyser_GET_BUF
332
333 float normfactor = unit->m_normfactor;
334 bool square = unit->m_square;
335 if(normfactor == 0.f){
336 if(square)
337 unit->m_normfactor = normfactor = 1.f / powf(numbins + 2.f, 1.5f);
338 else
339 unit->m_normfactor = normfactor = 1.f / (numbins + 2.f);
340 }
341
342
343 SCComplexBuf *p = ToComplexApx(buf);
344 // SCPolarBuf *p = ToPolarApx(buf);
345
346 float total;
347 if(square){
348 total = sc_abs(p->dc) * sc_abs(p->dc) + sc_abs(p->nyq) * sc_abs(p->nyq);
349
350 for (int i=0; i<numbins; ++i) {
351 float rabs = (p->bin[i].real);
352 float iabs = (p->bin[i].imag);
353 total += (rabs*rabs) + (iabs*iabs);
354 }
355 }else{
356 total = sc_abs(p->dc) + sc_abs(p->nyq);
357
358 for (int i=0; i<numbins; ++i) {
359 float rabs = (p->bin[i].real);
360 float iabs = (p->bin[i].imag);
361 total += sqrt((rabs*rabs) + (iabs*iabs));
362 }
363 // for (int i=0; i<numbins; ++i) {
364 // total += sc_abs(p->bin[i].mag);
365 // }
366 }
367
368 // Store the val for output in future calls
369 unit->outval = total * normfactor;
370
371 ZOUT0(0) = unit->outval;
372 }
373
374 //////////////////////////////////////////////////////////////////////////////////////////////////
375
FFTSubbandPower_Ctor(FFTSubbandPower * unit)376 void FFTSubbandPower_Ctor(FFTSubbandPower *unit)
377 {
378 SETCALC(FFTSubbandPower_next);
379 ZOUT0(0) = unit->outval = 0.;
380
381 unit->m_square = ZIN0(2) > 0.f;
382 unit->m_normfactor = 0.f;
383
384 // ZIN0(1) tells us how many cutoffs we're looking for
385 int numcutoffs = (int)ZIN0(1);
386 int numbands = numcutoffs+1;
387
388 unit->m_scalemode = (int)ZIN0(3);
389
390 float * outvals = (float*)RTAlloc(unit->mWorld, numbands * sizeof(float));
391 for(int i=0; i<numbands; i++) {
392 outvals[i] = 0.f;
393 }
394 unit->m_outvals = outvals;
395
396 unit->m_cutoffs = (int*)RTAlloc(unit->mWorld, numcutoffs * sizeof(int));
397 unit->m_cutoff_inited = false;
398
399 unit->m_numbands = numbands;
400 }
401
FFTSubbandPower_next(FFTSubbandPower * unit,int inNumSamples)402 void FFTSubbandPower_next(FFTSubbandPower *unit, int inNumSamples)
403 {
404 int numbands = unit->m_numbands;
405 int numcutoffs = numbands - 1;
406
407 // Multi-output equiv of FFTAnalyser_GET_BUF
408 float fbufnum = ZIN0(0);
409 if (fbufnum < 0.f) {
410 for(int i=0; i<numbands; i++){
411 ZOUT0(i) = unit->m_outvals[i];
412 }
413 return;
414 }
415 uint32 ibufnum = (uint32)fbufnum;
416 World *world = unit->mWorld;
417 SndBuf *buf;
418 if (ibufnum >= world->mNumSndBufs) {
419 int localBufNum = ibufnum - world->mNumSndBufs;
420 Graph *parent = unit->mParent;
421 if(localBufNum <= parent->localBufNum) {
422 buf = parent->mLocalSndBufs + localBufNum;
423 } else {
424 buf = world->mSndBufs;
425 }
426 } else {
427 buf = world->mSndBufs + ibufnum;
428 }
429 int numbins = (buf->samples - 2) >> 1;
430 // End: Multi-output equiv of FFTAnalyser_GET_BUF
431
432 int scalemode = unit->m_scalemode;
433
434 float normfactor = unit->m_normfactor;
435 bool square = unit->m_square;
436 if(normfactor == 0.f){
437 if(square)
438 unit->m_normfactor = normfactor = 1.f / powf(numbins + 2.f, 1.5f);
439 else
440 unit->m_normfactor = normfactor = 1.f / (numbins + 2.f);
441 }
442
443 // Now we create the integer lookup list, if it doesn't already exist
444 int * cutoffs = unit->m_cutoffs;
445 if(!unit->m_cutoff_inited){
446
447 float srate = world->mFullRate.mSampleRate;
448 for(int i=0; i < numcutoffs; ++i) {
449 cutoffs[i] = (int)(buf->samples * ZIN0(4 + i) / srate);
450 //Print("Allocated bin cutoff #%d, at bin %d\n", i, cutoffs[i]);
451 }
452
453 unit->m_cutoff_inited = true;
454 }
455
456 SCComplexBuf *p = ToComplexApx(buf);
457
458 // Now we can actually calculate the bandwise subtotals
459 float total = sc_abs(p->dc);
460 if(square){
461 total *= total; // square the DC val
462 }
463 int binaddcount = 1; // Counts how many bins contributed to the current band (1 because of the DC value)
464 int curband = 0;
465 float * outvals = unit->m_outvals;
466 float magsq;
467 for (int i=0; i<numbins; ++i) {
468 if((curband != numbands) && (i >= cutoffs[curband])){
469 if(scalemode==1){
470 outvals[curband] = total * normfactor;
471 }else{
472 if(square)
473 outvals[curband] = total / powf((float)binaddcount, 1.5f);
474 else
475 outvals[curband] = total / binaddcount;
476 }
477 //Print("Finished off band %i while in bin %i\n", curband, i);
478 ++curband;
479 total = 0.f;
480 binaddcount = 0;
481 }
482
483 float rabs = (p->bin[i].real);
484 float iabs = (p->bin[i].imag);
485 magsq = ((rabs*rabs) + (iabs*iabs));
486 if(square)
487 total += magsq;
488 else
489 total += std::sqrt(magsq);
490 ++binaddcount;
491 }
492 // Remember to output the very last (highest) band
493 if(square)
494 total += p->nyq * p->nyq;
495 else
496 total += sc_abs(p->nyq);
497 // Pop the last one onto the end of the lovely list
498 if(scalemode==1){
499 outvals[curband] = total * normfactor;
500 }else{
501 if(square)
502 outvals[curband] = total / powf((float)binaddcount + 1.f, 1.5f); // Plus one because of the nyq value
503 else
504 outvals[curband] = total / (binaddcount + 1); // Plus one because of the nyq value
505 }
506
507 // Now we can output the vals
508 for(int i=0; i<numbands; i++) {
509 ZOUT0(i) = outvals[i];
510 }
511 }
512
FFTSubbandPower_Dtor(FFTSubbandPower * unit)513 void FFTSubbandPower_Dtor(FFTSubbandPower *unit)
514 {
515 if(unit->m_cutoffs) RTFree(unit->mWorld, unit->m_cutoffs);
516 if(unit->m_outvals) RTFree(unit->mWorld, unit->m_outvals);
517 }
518
519 ////////////////////////////////////////////////////////////////////////////////////
520
FFTFlux_Ctor(FFTFlux_Unit * unit)521 void FFTFlux_Ctor(FFTFlux_Unit *unit)
522 {
523 SETCALC(FFTFlux_next);
524
525 unit->m_tempbuf = 0;
526
527 unit->m_yesternorm = 1.0f;
528 unit->m_yesterdc = 0.0f;
529 unit->m_yesternyq = 0.0f;
530
531 unit->m_normalise = ZIN0(1) > 0.f; // Whether we want to normalise or not
532
533 unit->outval = 0.f;
534 ZOUT0(0) = 0.f;
535 }
536
FFTFlux_next(FFTFlux_Unit * unit,int inNumSamples)537 void FFTFlux_next(FFTFlux_Unit *unit, int inNumSamples)
538 {
539 FFTAnalyser_GET_BUF
540
541 // Modified form of MAKE TEMP BUF, here used for storing last frame's magnitudes:
542 if (!unit->m_tempbuf) {
543 unit->m_tempbuf = (float*)RTAlloc(unit->mWorld, numbins * sizeof(float));
544 unit->m_numbins = numbins;
545
546 // Must also ensure the yester is zero'ed
547 // because it will be compared before being filled in FFTFlux calculation
548 memset(unit->m_tempbuf, 0, numbins * sizeof(float));
549 } else if (numbins != unit->m_numbins) return;
550
551 SCPolarBuf *p = ToPolarApx(buf); // Current frame
552 float* yestermags = unit->m_tempbuf; // This is an array storing the yester magnitudes
553
554 float yesternorm = unit->m_yesternorm; // Should have been calculated on prev cycle
555
556 float currnorm;
557 if(unit->m_normalise){
558 // First iteration is to find the sum of magnitudes (to find the normalisation factor):
559 currnorm = (p->dc * p->dc) + (p->nyq * p->nyq);
560 for (int i=0; i<numbins; ++i) {
561 currnorm += p->bin[i].mag * p->bin[i].mag;
562 }
563 // The normalisation factor is 1-over-sum
564 if(currnorm != 0.0f) {
565 currnorm = 1.0f / currnorm;
566 }
567 } else {
568 currnorm = 1.f;
569 }
570
571 // This iteration is the meat of the algorithm. Compare current (normed) bins against prev.
572 float onebindiff = sc_abs(p->dc * currnorm) - sc_abs(unit->m_yesterdc * yesternorm);
573 float fluxsquared = (onebindiff * onebindiff);
574 onebindiff = sc_abs(p->nyq * currnorm) - sc_abs(unit->m_yesternyq * yesternorm);
575 fluxsquared += (onebindiff * onebindiff);
576 // Now the bins
577 for (int i=0; i<numbins; ++i) {
578 // Sum the squared difference of normalised mags onto the cumulative value
579 onebindiff = (p->bin[i].mag * currnorm) - (yestermags[i] * yesternorm);
580 fluxsquared += (onebindiff * onebindiff);
581 // Overwrite yestermag values with current values, so they're available next time
582 yestermags[i] = p->bin[i].mag;
583 }
584
585 // Store the just-calc'ed norm as yesternorm
586 unit->m_yesternorm = currnorm;
587 unit->m_yesterdc = p->dc;
588 unit->m_yesternyq = p->nyq;
589
590 // Store the val for output in future calls
591 unit->outval = sqrt(fluxsquared);
592
593 ZOUT0(0) = unit->outval;
594
595 }
FFTFlux_Dtor(FFTFlux_Unit * unit)596 void FFTFlux_Dtor(FFTFlux_Unit *unit)
597 {
598 RTFree(unit->mWorld, unit->m_tempbuf);
599 }
600
601
602 ////////////////////////////////////////////////////////////////////////////////////
603
FFTFluxPos_Ctor(FFTFlux_Unit * unit)604 void FFTFluxPos_Ctor(FFTFlux_Unit *unit)
605 {
606 SETCALC(FFTFluxPos_next);
607
608 unit->m_tempbuf = 0;
609
610 unit->m_yesternorm = 1.0f;
611 unit->m_yesterdc = 0.0f;
612 unit->m_yesternyq = 0.0f;
613
614 unit->outval = 0.f;
615 ZOUT0(0) = 0.f;
616 }
617
FFTFluxPos_next(FFTFlux_Unit * unit,int inNumSamples)618 void FFTFluxPos_next(FFTFlux_Unit *unit, int inNumSamples)
619 {
620 FFTAnalyser_GET_BUF
621
622 // Modified form of MAKE TEMP BUF, here used for storing last frame's magnitudes:
623 if (!unit->m_tempbuf) {
624 unit->m_tempbuf = (float*)RTAlloc(unit->mWorld, numbins * sizeof(float));
625 unit->m_numbins = numbins;
626
627 // Must also ensure the yester is zero'ed
628 // because it will be compared before being filled in FFTFlux calculation
629 memset(unit->m_tempbuf, 0, numbins * sizeof(float));
630 } else if (numbins != unit->m_numbins) return;
631
632 SCPolarBuf *p = ToPolarApx(buf); // Current frame
633 float* yestermags = unit->m_tempbuf; // This is an array storing the yester magnitudes
634
635 float yesternorm = unit->m_yesternorm; // Should have been calculated on prev cycle
636
637 float currnorm;
638 if(unit->m_normalise){
639 // First iteration is to find the sum of magnitudes (to find the normalisation factor):
640 currnorm = (p->dc * p->dc) + (p->nyq * p->nyq);
641 for (int i=0; i<numbins; ++i) {
642 currnorm += p->bin[i].mag * p->bin[i].mag;
643 }
644 // The normalisation factor is 1-over-sum
645 if(currnorm != 0.0f) {
646 currnorm = 1.0f / currnorm;
647 }
648 } else {
649 currnorm = 1.f;
650 }
651
652
653 // This iteration is the meat of the algorithm. Compare current (normed) bins against prev.
654 float onebindiff = sc_abs(p->dc * currnorm) - sc_abs(unit->m_yesterdc * yesternorm);
655 float fluxsquared = 0.f;
656 if(onebindiff > 0.f)
657 fluxsquared += (onebindiff * onebindiff);
658 onebindiff = sc_abs(p->nyq * currnorm) - sc_abs(unit->m_yesternyq * yesternorm);
659 if(onebindiff > 0.f)
660 fluxsquared += (onebindiff * onebindiff);
661 // Now the bins
662 for (int i=0; i<numbins; ++i) {
663 // Sum the squared difference of normalised mags onto the cumulative value
664 onebindiff = (p->bin[i].mag * currnorm) - (yestermags[i] * yesternorm);
665
666 ////////// THIS IS WHERE FFTFluxPos DIFFERS FROM FFTFlux - THE SIMPLE ADDITION OF AN "IF": //////////
667 if(onebindiff > 0.f) // The IF only applies to the next line - formatting is a bit weird to keep in line with the other func
668
669 fluxsquared += (onebindiff * onebindiff);
670 // Overwrite yestermag values with current values, so they're available next time
671 yestermags[i] = p->bin[i].mag;
672 }
673
674 // Store the just-calc'ed norm as yesternorm
675 unit->m_yesternorm = currnorm;
676 unit->m_yesterdc = p->dc;
677 unit->m_yesternyq = p->nyq;
678
679 // Store the val for output in future calls
680 unit->outval = sqrt(fluxsquared);
681
682 ZOUT0(0) = unit->outval;
683
684 }
FFTFluxPos_Dtor(FFTFlux_Unit * unit)685 void FFTFluxPos_Dtor(FFTFlux_Unit *unit)
686 {
687 RTFree(unit->mWorld, unit->m_tempbuf);
688 }
689
690 ////////////////////////////////////////////////////////////////////////////////////
691
FFTDiffMags_next(FFTAnalyser_Unit * unit,int inNumSamples)692 void FFTDiffMags_next(FFTAnalyser_Unit *unit, int inNumSamples)
693 {
694
695 FFTAnalyser_GET_BUF2
696
697 SCComplexBuf *p = ToComplexApx(buf1);
698 SCComplexBuf *q = ToComplexApx(buf2);
699
700 // First the DC and nyquist.
701 float diffsum = sc_abs(p->dc - q->dc) + sc_abs(p->nyq - q->nyq);
702
703 for (int i=0; i<numbins; ++i) {
704 float rdiff = p->bin[i].real - q->bin[i].real;
705 float idiff = p->bin[i].imag - q->bin[i].imag;
706
707 diffsum += sqrt(rdiff*rdiff + idiff*idiff);
708 }
709
710 // Store the val for output in future calls
711 unit->outval = diffsum / (numbins + 2);
712
713 //Print("FFTDiffMags_next: output is %g\n", unit->outval);
714
715 ZOUT0(0) = unit->outval;
716 }
717
FFTDiffMags_Ctor(FFTAnalyser_Unit * unit)718 void FFTDiffMags_Ctor(FFTAnalyser_Unit *unit)
719 {
720 SETCALC(FFTDiffMags_next);
721 ZOUT0(0) = ZIN0(0);
722 }
723
724 ////////////////////////////////////////////////////////////////////////////////////
725
PV_MagSubtract_next(PV_Unit * unit,int inNumSamples)726 void PV_MagSubtract_next(PV_Unit *unit, int inNumSamples)
727 {
728 PV_GET_BUF2
729
730 SCPolarBuf *p = ToPolarApx(buf1);
731 SCPolarBuf *q = ToPolarApx(buf2);
732
733 bool zerolimit = ZIN0(2) > 0.f;
734
735 if(zerolimit){
736 // First the DC and nyquist
737 if(p->dc > q->dc)
738 p->dc = p->dc - q->dc;
739 else
740 p->dc = 0.f;
741
742 if(p->nyq > q->nyq)
743 p->nyq = p->nyq - q->nyq;
744 else
745 p->nyq = 0.f;
746
747 for (int i=0; i<numbins; ++i) {
748 if(p->bin[i].mag > q->bin[i].mag)
749 p->bin[i].mag -= q->bin[i].mag;
750 else
751 p->bin[i].mag = 0.f;
752 }
753 }else{
754 // First the DC and nyquist
755 p->dc = p->dc - q->dc;
756 p->nyq = p->nyq - q->nyq;
757 for (int i=0; i<numbins; ++i) {
758 p->bin[i].mag -= q->bin[i].mag;
759 }
760 }
761 }
762
PV_MagSubtract_Ctor(PV_Unit * unit)763 void PV_MagSubtract_Ctor(PV_Unit *unit)
764 {
765 SETCALC(PV_MagSubtract_next);
766 ZOUT0(0) = ZIN0(0);
767 }
768
769 ////////////////////////////////////////////////////////////////////////////////////////////////////////
770
PV_MagLog_Ctor(PV_Unit * unit)771 void PV_MagLog_Ctor(PV_Unit *unit)
772 {
773 SETCALC(PV_MagLog_next);
774 ZOUT0(0) = ZIN0(0);
775 }
776
PV_MagLog_next(PV_Unit * unit,int inNumSamples)777 void PV_MagLog_next(PV_Unit *unit, int inNumSamples)
778 {
779 PV_GET_BUF
780
781 SCPolarBuf *p = ToPolarApx(buf);
782
783 for (int i=0; i<numbins; ++i) {
784 float mag = p->bin[i].mag;
785 p->bin[i].mag = log(mag > SMALLEST_NUM_FOR_LOG ? mag : SMALLEST_NUM_FOR_LOG);
786 }
787 }
788
789 ////
790
PV_MagExp_Ctor(PV_Unit * unit)791 void PV_MagExp_Ctor(PV_Unit *unit)
792 {
793 SETCALC(PV_MagExp_next);
794 ZOUT0(0) = ZIN0(0);
795 }
796
PV_MagExp_next(PV_Unit * unit,int inNumSamples)797 void PV_MagExp_next(PV_Unit *unit, int inNumSamples)
798 {
799 PV_GET_BUF
800
801 SCPolarBuf *p = ToPolarApx(buf);
802
803 for (int i=0; i<numbins; ++i) {
804 float mag = p->bin[i].mag;
805 p->bin[i].mag = exp(mag);
806 }
807 }
808
809 ////////////////////////////////////////////////////////////////////////////////////////////////////////
810
FFTPhaseDev_next(FFTPhaseDev * unit,int inNumSamples)811 void FFTPhaseDev_next(FFTPhaseDev *unit, int inNumSamples)
812 {
813 FFTAnalyser_GET_BUF
814
815 // Get the current frame, as polar. NB We don't care about DC or nyquist in this UGen.
816 SCPolarBuf *p = ToPolarApx(buf);
817 int tbpointer;
818
819 float powthresh = ZIN0(2);
820
821 // MAKE_TEMP_BUF but modified:
822 if (!unit->m_tempbuf) {
823 unit->m_tempbuf = (float*)RTAlloc(unit->mWorld, numbins * 2 * sizeof(float));
824 memset(unit->m_tempbuf, 0, numbins * 2 * sizeof(float)); // Ensure it's zeroed
825 // Ensure the initial values don't cause some weird jump in the output - set them to vals which will produce deviation of zero
826 tbpointer = 0;
827 for (int i=0; i<numbins; ++i) {
828 unit->m_tempbuf[tbpointer++] = p->bin[i].phase;
829 unit->m_tempbuf[tbpointer++] = 0.f;
830 }
831 unit->m_numbins = numbins;
832 } else if (numbins != unit->m_numbins) return;
833
834 // Retrieve state
835 bool useweighting = unit->m_weight;
836 float *storedvals = unit->m_tempbuf;
837
838 // Note: temp buf is stored in this format: phase[0],d_phase[0],phase[1],d_phase[1], ...
839
840 // Print("\nbin[10] phase: %g\nbin[10] yesterphase: %g\nbin[10] yesterdiff: %g\n",
841 // /*PHASE_REWRAP(*/p->bin[10].phase/*)*/, storedvals[20], storedvals[21]);
842
843 //Print("\npowthresh is %g", powthresh);
844
845 // Iterate through, calculating the deviation from expected value.
846 double totdev = 0.0;
847 tbpointer = 0;
848 float deviation;
849 for (int i=0; i<numbins; ++i) {
850 // Thresholding as Brossier did - discard bin's phase deviation if bin's power is minimal
851 if(p->bin[i].mag > powthresh) {
852
853 // Deviation is the *second difference* of the phase, which is calc'ed as curval - yesterval - yesterfirstdiff
854 deviation = p->bin[i].phase - storedvals[tbpointer] - storedvals[tbpointer+1];
855 tbpointer += 2;
856 // Wrap onto +-PI range
857 deviation = PHASE_REWRAP(deviation);
858
859 if(useweighting){
860 totdev += fabs(deviation * p->bin[i].mag);
861 } else {
862 totdev += fabs(deviation);
863 }
864 }
865 }
866
867 // totdev will be the output, but first we need to fill tempbuf with today's values, ready for tomorrow.
868 tbpointer = 0;
869 float diff;
870 for (int i=0; i<numbins; ++i) {
871 diff = p->bin[i].phase - storedvals[tbpointer]; // Retrieving yesterphase from buf
872 storedvals[tbpointer++] = p->bin[i].phase; // Storing phase
873 // Wrap onto +-PI range
874 diff = PHASE_REWRAP(diff);
875
876 storedvals[tbpointer++] = diff; // Storing first diff to buf
877
878 }
879
880 // Store the val for output in future calls
881 unit->outval = (float)totdev;
882
883 ZOUT0(0) = unit->outval;
884 }
885
FFTPhaseDev_Ctor(FFTPhaseDev * unit)886 void FFTPhaseDev_Ctor(FFTPhaseDev *unit)
887 {
888 SETCALC(FFTPhaseDev_next);
889
890 unit->m_weight = (ZIN0(1) > 0.f) ? true : false;
891
892 ZOUT0(0) = unit->outval = 0.;
893 unit->m_tempbuf = 0;
894 }
895
FFTPhaseDev_Dtor(FFTPhaseDev * unit)896 void FFTPhaseDev_Dtor(FFTPhaseDev *unit)
897 {
898 RTFree(unit->mWorld, unit->m_tempbuf);
899 }
900
901 ////////////////////////////////////////////////////////////////////////////////////////////////////////
902
FFTComplexDev_next(FFTComplexDev * unit,int inNumSamples)903 void FFTComplexDev_next(FFTComplexDev *unit, int inNumSamples)
904 {
905 FFTAnalyser_GET_BUF
906
907 // Get the current frame, as polar. NB We don't care about DC or nyquist in this UGen.
908 SCPolarBuf *p = ToPolarApx(buf);
909 int tbpointer;
910
911 float powthresh = ZIN0(2);
912
913 // MAKE_TEMP_BUF but modified:
914 if (!unit->m_tempbuf) {
915 unit->m_tempbuf = (float*)RTAlloc(unit->mWorld, numbins * 3 * sizeof(float));
916 memset(unit->m_tempbuf, 0, numbins * 3 * sizeof(float)); // Ensure it's zeroed
917 // Ensure the initial values don't cause some weird jump in the output - set them to vals which will produce deviation of zero
918 tbpointer = 0;
919 for (int i=0; i<numbins; ++i) {
920 unit->m_tempbuf[tbpointer++] = p->bin[i].phase;
921 unit->m_tempbuf[tbpointer++] = 0.f;
922 }
923 unit->m_numbins = numbins;
924 } else if (numbins != unit->m_numbins) return;
925
926 // Retrieve state
927 float *storedvals = unit->m_tempbuf;
928 bool rectify = unit->m_rectify;
929
930 // Note: temp buf is stored in this format: mag[0],phase[0],d_phase[0],mag[1],phase[1],d_phase[1], ...
931
932
933 //Print("\npowthresh is %g", powthresh);
934
935 // Iterate through, calculating the deviation from expected value.
936 double totdev = 0.0;
937 tbpointer = 0;
938 float curmag, predmag, predphase, yesterphase, yesterphasediff;
939 float deviation;
940 for (int i=0; i<numbins; ++i) {
941 curmag = p->bin[i].mag;
942
943 // Predict mag as yestermag
944 predmag = storedvals[tbpointer++];
945 yesterphase = storedvals[tbpointer++];
946 yesterphasediff = storedvals[tbpointer++];
947
948 // Thresholding as Brossier did - discard bin's deviation if bin's power is minimal
949 if(curmag > powthresh) {
950 // If rectifying, ignore decreasing bins
951 if((!rectify) || (curmag >= predmag)){
952
953 // Predict phase as yesterval + yesterfirstdiff
954 predphase = yesterphase + yesterphasediff;
955
956 // Deviation is Euclidean distance between predicted and actual.
957 // In polar coords: sqrt(r1^2 + r2^2 - r1r2 cos (theta1 - theta2))
958 deviation = sqrt(predmag * predmag + curmag * curmag
959 - predmag * predmag * cos(PHASE_REWRAP(predphase - p->bin[i].phase))
960 );
961
962 totdev += deviation;
963 }
964 }
965 }
966
967 // totdev will be the output, but first we need to fill tempbuf with today's values, ready for tomorrow.
968 tbpointer = 0;
969 float diff;
970 for (int i=0; i<numbins; ++i) {
971 storedvals[tbpointer++] = p->bin[i].mag; // Storing mag
972 diff = p->bin[i].phase - storedvals[tbpointer]; // Retrieving yesterphase from buf
973 storedvals[tbpointer++] = p->bin[i].phase; // Storing phase
974 // Wrap onto +-PI range
975 diff = PHASE_REWRAP(diff);
976
977 storedvals[tbpointer++] = diff; // Storing first diff to buf
978
979 }
980
981 // Store the val for output in future calls
982 unit->outval = (float)totdev;
983
984 ZOUT0(0) = unit->outval;
985 }
986
FFTComplexDev_Ctor(FFTComplexDev * unit)987 void FFTComplexDev_Ctor(FFTComplexDev *unit)
988 {
989 SETCALC(FFTComplexDev_next);
990
991 unit->m_rectify = (ZIN0(1) > 0.f) ? true : false;
992
993 ZOUT0(0) = unit->outval = 0.;
994 unit->m_tempbuf = 0;
995 }
996
FFTComplexDev_Dtor(FFTComplexDev * unit)997 void FFTComplexDev_Dtor(FFTComplexDev *unit)
998 {
999 RTFree(unit->mWorld, unit->m_tempbuf);
1000 }
1001
1002
1003 ////////////////////////////////////////////////////////////////////////////////////////////////////////
1004
FFTMKL_next(FFTMKL * unit,int inNumSamples)1005 void FFTMKL_next(FFTMKL *unit, int inNumSamples)
1006 {
1007 FFTAnalyser_GET_BUF
1008
1009 // Get the current frame, as polar. NB We don't care about DC or nyquist in this UGen.
1010 SCPolarBuf *p = ToPolarApx(buf);
1011 int tbpointer;
1012 float eta = ZIN0(1);
1013
1014 // MAKE_TEMP_BUF but modified:
1015 if (!unit->m_tempbuf) {
1016 unit->m_tempbuf = (float*)RTAlloc(unit->mWorld, numbins * 1 * sizeof(float));
1017 memset(unit->m_tempbuf, 0, numbins * 1 * sizeof(float)); // Ensure it's zeroed
1018 // Ensure the initial values don't cause some weird jump in the output - set them to vals which will produce deviation of zero
1019 tbpointer = 0;
1020 for (int i=0; i<numbins; ++i) {
1021 unit->m_tempbuf[tbpointer++] = p->bin[i].mag;
1022 }
1023 unit->m_numbins = numbins;
1024 } else if (numbins != unit->m_numbins) return;
1025
1026 // Retrieve state
1027 float *storedvals = unit->m_tempbuf;
1028
1029 // Note: for this UGen, temp buf is just mag[0],mag[1],... - we ain't interested in phase etc
1030
1031 // Iterate through, calculating the Modified Kullback-Liebler distance
1032 double totdev = 0.0;
1033 tbpointer = 0;
1034 float curmag, yestermag;
1035 float deviation;
1036 for (int i=0; i<numbins; ++i) {
1037 curmag = p->bin[i].mag;
1038 yestermag = storedvals[tbpointer];
1039
1040 // Here's the main implementation of Brossier's MKL eq'n (eqn 2.9 from his thesis):
1041 deviation = sc_abs(curmag) / (sc_abs(yestermag) + eta);
1042 totdev += log(1.f + deviation);
1043
1044 // Store the mag as yestermag
1045 storedvals[tbpointer++] = curmag;
1046 }
1047
1048 // Store the val for output in future calls
1049 unit->outval = (float)totdev;
1050
1051 ZOUT0(0) = unit->outval;
1052 }
1053
FFTMKL_Ctor(FFTMKL * unit)1054 void FFTMKL_Ctor(FFTMKL *unit)
1055 {
1056 SETCALC(FFTMKL_next);
1057
1058 ZOUT0(0) = unit->outval = 0.;
1059 unit->m_tempbuf = 0;
1060 }
1061
FFTMKL_Dtor(FFTMKL * unit)1062 void FFTMKL_Dtor(FFTMKL *unit)
1063 {
1064 RTFree(unit->mWorld, unit->m_tempbuf);
1065 }
1066
1067
1068 ////////////////////////////////////////////////////////////////////////////////
1069
PV_Whiten_Ctor(PV_Whiten * unit)1070 void PV_Whiten_Ctor(PV_Whiten *unit){
1071
1072 SETCALC(PV_Whiten_next);
1073
1074 ZOUT0(0) = ZIN0(0);
1075 }
1076
1077
PV_Whiten_next(PV_Whiten * unit,int inNumSamples)1078 void PV_Whiten_next(PV_Whiten *unit, int inNumSamples){
1079
1080 float fbufnum1 = ZIN0(0);
1081 float fbufnum2 = ZIN0(1);
1082 if (fbufnum1 < 0.f || fbufnum2 < 0.f) { ZOUT0(0) = -1.f; return; }
1083 // Print("\nfbufnum1: %g; fbufnum2: %g", fbufnum1, fbufnum2);
1084 uint32 ibufnum1 = (int)fbufnum1;
1085 uint32 ibufnum2 = (int)fbufnum2;
1086 World *world = unit->mWorld;
1087 SndBuf *buf1;
1088 if (ibufnum1 >= world->mNumSndBufs) {
1089 int localBufNum = ibufnum1 - world->mNumSndBufs;
1090 Graph *parent = unit->mParent;
1091 if(localBufNum <= parent->localBufNum) {
1092 buf1 = parent->mLocalSndBufs + localBufNum;
1093 } else {
1094 buf1 = world->mSndBufs;
1095 }
1096 } else {
1097 buf1 = world->mSndBufs + ibufnum1;
1098 }
1099 SndBuf *buf2;
1100 if (ibufnum2 >= world->mNumSndBufs) {
1101 int localBufNum = ibufnum2 - world->mNumSndBufs;
1102 Graph *parent = unit->mParent;
1103 if(localBufNum <= parent->localBufNum) {
1104 buf2 = parent->mLocalSndBufs + localBufNum;
1105 } else {
1106 buf2 = world->mSndBufs;
1107 }
1108 } else {
1109 buf2 = world->mSndBufs + ibufnum2;
1110 }
1111 int numbins = (buf1->samples - 2) >> 1;
1112 // Print("\nibufnum1: %d; ibufnum2: %d", ibufnum1, ibufnum2);
1113 // if (buf1->samples != buf2->samples) return;
1114
1115 // Print("\nnumbins: %d", numbins);
1116
1117 // memcpy(buf2->data, buf1->data, buf1->samples * sizeof(float));
1118
1119 SCPolarBuf *indata = ToPolarApx(buf1);
1120
1121 // This buffer stores numbins+2 amplitude tracks, in "logical" order (DC, bin1, ... nyquist), not in the order produced by the FFT
1122 float *pkdata = buf2->data;
1123
1124 // Update the parameters
1125 float relax = ZIN0(2);
1126 float relaxcoef = (relax == 0.0f) ? 0.0f : exp(log1/(relax * SAMPLERATE));
1127 float floor = ZIN0(3);
1128 float smear = ZIN0(4);
1129 // unsigned int bindownsample = (int)ZIN0(5);
1130
1131 float val,oldval;
1132 ////////////////////// Now for each bin, update the record of the peak value /////////////////////
1133
1134 val = fabs(indata->dc); // Grab current magnitude
1135 oldval = pkdata[0];
1136 // If it beats the amplitude stored then that's our new amplitude; otherwise our new amplitude is a decayed version of the old one
1137 if(val < oldval) {
1138 val = val + (oldval - val) * relaxcoef;
1139 }
1140 pkdata[0] = val; // Store the "amplitude trace" back
1141
1142 val = fabs(indata->nyq);
1143 oldval = pkdata[numbins+1];
1144 if(val < oldval) {
1145 val = val + (oldval - val) * relaxcoef;
1146 }
1147 pkdata[numbins+1] = val;
1148
1149 // Print("-----------Peaks-------\n");
1150 for(int i=0; i<numbins; ++i){
1151 val = fabs(indata->bin[i].mag);
1152 oldval = pkdata[i+1];
1153 if(val < oldval) {
1154 val = val + (oldval - val) * relaxcoef;
1155 }
1156 pkdata[i+1] = val;
1157 //Print("%g, ", val);
1158 }
1159 // Print("\n");
1160
1161 // Perform smearing now
1162 if(smear != 0.f){
1163 float oldval, newval;
1164 oldval = pkdata[0];
1165 // What we want is to keep the largest of curval, prevval*smear, nextval*smear.
1166 // We do this in two steps, by keeping the biggest of prevval and nextval, then keeping the largest of (biggest*smear, curval)
1167 for(int i=1; i<=numbins; i++){
1168 oldval = sc_max(oldval, pkdata[i+1]);
1169 newval = sc_max(oldval * smear, pkdata[i]);
1170
1171 oldval = pkdata[i]; // For next iter
1172 pkdata[i] = newval;
1173 }
1174 }
1175
1176 //////////////////////////// Now for each bin, rescale the current magnitude ////////////////////////////
1177 indata->dc /= sc_max(floor, pkdata[0]);
1178 indata->nyq /= sc_max(floor, pkdata[numbins+1]);
1179 for(int i=0; i<numbins; ++i){
1180 indata->bin[i].mag /= sc_max(floor, pkdata[i+1]);
1181 }
1182
1183 ZOUT0(0) = fbufnum1;
1184
1185 }
1186
1187 //////////////////////////////////////////////////////////////////////////////////////////////////
1188
FFTSubbandFlatness_next(FFTSubbandFlatness * unit,int inNumSamples)1189 void FFTSubbandFlatness_next(FFTSubbandFlatness *unit, int inNumSamples)
1190 {
1191 int numbands = unit->m_numbands;
1192 int numcutoffs = numbands - 1;
1193
1194 // Multi-output equiv of FFTAnalyser_GET_BUF
1195 float fbufnum = ZIN0(0);
1196 if (fbufnum < 0.f) {
1197 for(int i=0; i<numbands; i++){
1198 ZOUT0(i) = unit->m_outvals[i];
1199 }
1200 return;
1201 }
1202 uint32 ibufnum = (uint32)fbufnum;
1203 World *world = unit->mWorld;
1204 SndBuf *buf;
1205 if (ibufnum >= world->mNumSndBufs) {
1206 int localBufNum = ibufnum - world->mNumSndBufs;
1207 Graph *parent = unit->mParent;
1208 if(localBufNum <= parent->localBufNum) {
1209 buf = parent->mLocalSndBufs + localBufNum;
1210 } else {
1211 buf = world->mSndBufs;
1212 }
1213 } else {
1214 buf = world->mSndBufs + ibufnum;
1215 }
1216 int numbins = (buf->samples - 2) >> 1;
1217 // End: Multi-output equiv of FFTAnalyser_GET_BUF
1218
1219 // Now we create the integer lookup list, if it doesn't already exist
1220 int * cutoffs = unit->m_cutoffs;
1221 if(!unit->m_cutoff_inited){
1222
1223 float srate = world->mFullRate.mSampleRate;
1224 for(int i=0; i < numcutoffs; i++) {
1225 cutoffs[i] = (int)(buf->samples * ZIN0(2 + i) / srate);
1226 //Print("Allocated bin cutoff #%d, at bin %d\n", i, cutoffs[i]);
1227 }
1228
1229 unit->m_cutoff_inited = true;
1230 }
1231
1232 SCPolarBuf *p = ToPolarApx(buf);
1233
1234 // Now we can actually calculate the bandwise stuff
1235 int binaddcount = 0; // Counts how many bins contributed to the current band
1236 int curband = 0;
1237 float * outvals = unit->m_outvals;
1238
1239 double geommean = 0.0, arithmean = 0.0;
1240
1241 for (int i=0; i<numbins; ++i) {
1242 if(i == cutoffs[curband]){
1243 // Finish off the mean calculations
1244 geommean = exp(geommean / binaddcount);
1245 arithmean /= binaddcount;
1246 outvals[curband] = (float)(geommean / arithmean);
1247 curband++;
1248 geommean = arithmean = 0.0;
1249 binaddcount = 0;
1250 }
1251
1252 float mag = (p->bin[i].mag);
1253 geommean += log(mag);
1254 arithmean += mag;
1255
1256 binaddcount++;
1257 }
1258
1259 // Remember to output the very last (highest) band
1260 // Do the nyquist
1261 geommean += log(sc_abs(p->nyq));
1262 arithmean += sc_abs(p->nyq);
1263 binaddcount++;
1264 // Finish off the mean calculations
1265 geommean = exp(geommean / binaddcount);
1266 arithmean /= binaddcount;
1267 outvals[curband] = (float)(geommean / arithmean);
1268
1269 // Now we can output the vals
1270 for(int i=0; i<numbands; i++) {
1271 ZOUT0(i) = outvals[i];
1272 }
1273 }
1274
FFTSubbandFlatness_Ctor(FFTSubbandFlatness * unit)1275 void FFTSubbandFlatness_Ctor(FFTSubbandFlatness *unit)
1276 {
1277 SETCALC(FFTSubbandFlatness_next);
1278
1279 // ZIN0(1) tells us how many cutoffs we're looking for
1280 int numcutoffs = (int)ZIN0(1);
1281 int numbands = numcutoffs+1;
1282
1283 float * outvals = (float*)RTAlloc(unit->mWorld, numbands * sizeof(float));
1284 for(int i=0; i<numbands; i++) {
1285 outvals[i] = 0.f;
1286 }
1287 unit->m_outvals = outvals;
1288
1289 unit->m_cutoffs = (int*)RTAlloc(unit->mWorld,
1290 numcutoffs * sizeof(int)
1291 );
1292
1293 unit->m_cutoff_inited = false;
1294
1295 unit->m_numbands = numbands;
1296 ZOUT0(0) = unit->outval = 0.;
1297 }
1298
FFTSubbandFlatness_Dtor(FFTSubbandFlatness * unit)1299 void FFTSubbandFlatness_Dtor(FFTSubbandFlatness *unit)
1300 {
1301 RTFree(unit->mWorld, unit->m_cutoffs);
1302 RTFree(unit->mWorld, unit->m_outvals);
1303 }
1304
1305 //////////////////////////////////////////////////////////////////////////////////////////////////
1306
FFTCrest_Ctor(FFTCrest * unit)1307 void FFTCrest_Ctor(FFTCrest *unit)
1308 {
1309 SETCALC(FFTCrest_next);
1310
1311 unit->m_cutoffneedsinit = true;
1312 unit->m_freqtobin = 0.f;
1313
1314 ZOUT0(0) = unit->outval = 1.f;
1315 }
1316
FFTCrest_next(FFTCrest * unit,int inNumSamples)1317 void FFTCrest_next(FFTCrest *unit, int inNumSamples)
1318 {
1319 float freqlo = IN0(1);
1320 float freqhi = IN0(2);
1321
1322 FFTAnalyser_GET_BUF
1323
1324 //SCPolarBuf *p = ToPolarApx(buf); // Seems buggy...?
1325 SCComplexBuf *p = ToComplexApx(buf);
1326
1327 GET_FREQTOBIN
1328
1329 if(unit->m_cutoffneedsinit){
1330 // Get desired range, convert to bin index
1331 unit->m_frombin = (int)(freqtobin * freqlo);
1332 unit->m_tobinp1 = (int)(freqtobin * freqhi);
1333 if(unit->m_frombin < 0)
1334 unit->m_frombin = 0;
1335 if(unit->m_tobinp1 > numbins)
1336 unit->m_tobinp1 = numbins;
1337
1338 unit->m_cutoffneedsinit = false;
1339 }
1340 int frombin = unit->m_frombin;
1341 int tobinp1 = unit->m_tobinp1;
1342
1343 float total = 0.f, scf, sqrmag, peak=0.f;
1344 for (int i=frombin; i<tobinp1; ++i) {
1345 //sqrmag = p->bin[i].mag * p->bin[i].mag;
1346 sqrmag = (p->bin[i].real * p->bin[i].real) + (p->bin[i].imag * p->bin[i].imag);
1347 // (1) Check if it's the peak
1348 if(sqrmag >= peak){
1349 peak = sqrmag;
1350 }
1351 // (2) Add to subtotal
1352 total = total + sqrmag;
1353 }
1354
1355 // SCF defined as peak val divided by mean val; in other words, peak * count / total
1356 if(total == 0.f)
1357 scf = 1.f; // If total==0, peak==0, so algo output is indeterminate; but 1 indicates a perfectly flat spectrum, so we use that
1358 else
1359 scf = peak * ((float)(tobinp1 - frombin - 1)) / total;
1360
1361 // Store the val for output in future calls
1362 unit->outval = scf;
1363
1364 ZOUT0(0) = unit->outval;
1365 }
1366
1367 ////////////////////////////////////////////////////////////////////////////////////////////////////////
1368
FFTSpread_Ctor(FFTSpread * unit)1369 void FFTSpread_Ctor(FFTSpread *unit)
1370 {
1371 SETCALC(FFTSpread_next);
1372 ZOUT0(0) = unit->outval = 0.;
1373
1374 unit->m_bintofreq = 0.f;
1375 }
1376
FFTSpread_next(FFTSpread * unit,int inNumSamples)1377 void FFTSpread_next(FFTSpread *unit, int inNumSamples)
1378 {
1379 FFTAnalyser_GET_BUF
1380
1381 SCPolarBuf *p = ToPolarApx(buf);
1382
1383 GET_BINTOFREQ
1384
1385 float centroid = ZIN0(1);
1386
1387 float distance = ((numbins + 1) * bintofreq) - centroid;
1388 float mag = sc_abs(p->nyq);
1389 double num = mag * distance * distance;
1390 double denom = mag;
1391 for (int i=0; i<numbins; ++i) {
1392 distance = ((i+1) * bintofreq) - centroid;
1393 mag = sc_abs(p->bin[i].mag);
1394 num += mag * distance * distance;
1395 denom += mag;
1396 }
1397
1398 ZOUT0(0) = unit->outval = (denom==0.0 ? 0.f : (float) num/denom);
1399
1400 }
1401
1402 ////////////////////////////////////////////////////////////////////////////////////////////////////////
1403
FFTSlope_Ctor(FFTSlope * unit)1404 void FFTSlope_Ctor(FFTSlope *unit)
1405 {
1406 SETCALC(FFTSlope_next);
1407 ZOUT0(0) = unit->outval = 0.;
1408
1409 unit->m_bintofreq = 0.f;
1410 }
1411
FFTSlope_next(FFTSlope * unit,int inNumSamples)1412 void FFTSlope_next(FFTSlope *unit, int inNumSamples)
1413 {
1414 FFTAnalyser_GET_BUF
1415
1416 SCPolarBuf *p = ToPolarApx(buf);
1417
1418 GET_BINTOFREQ
1419
1420 // This is a straightforward linear regression slope on the magnitudes.
1421
1422 // These vars accumulate as we iter the bins. We start by putting in their values from the DC & Nyquist oddities.
1423 double sumx = (numbins+1) * bintofreq;
1424 double sumx2 = sumx * sumx;
1425 double sumxy = sumx * sc_abs(p->nyq);
1426 double sumy = sc_abs(p->dc) + sc_abs(p->nyq);
1427 double mag, freq;
1428
1429 for(int i=0; i<numbins; ++i){
1430 mag = p->bin[i].mag;
1431 freq = (i+1) * bintofreq;
1432
1433 sumxy += (freq * mag);
1434 sumx += freq;
1435 sumy += mag;
1436 sumx2 += (freq*freq);
1437 };
1438
1439 float slope = (float)((numbins * sumxy - sumx * sumy)
1440 / (numbins * sumx2 - sumx * sumx));
1441
1442 ZOUT0(0) = unit->outval = slope;
1443
1444 }
1445
1446 ////////////////////////////////////////////////////////////////////////////////////////////////////////
1447
FFTPeak_Ctor(FFTPeak * unit)1448 void FFTPeak_Ctor(FFTPeak *unit)
1449 {
1450 SETCALC(FFTPeak_next);
1451 ZOUT0(0) = unit->outval = 0.;
1452
1453 unit->m_bintofreq = 0.f;
1454 unit->m_freqtobin = 0.f;
1455
1456 unit->minbin = -99; // flag for the _next func to initialise it
1457 unit->minfreq = ZIN0(1);
1458 if(unit->minfreq < 0.f) unit->minfreq = 0.f;
1459 unit->maxfreq = ZIN0(2);
1460 if(unit->maxfreq < 0.f) unit->maxfreq = 0.f;
1461 }
1462
FFTPeak_next(FFTPeak * unit,int inNumSamples)1463 void FFTPeak_next(FFTPeak *unit, int inNumSamples)
1464 {
1465 FFTAnalyser_GET_BUF_TWOOUTS
1466
1467 SCPolarBuf *p = ToPolarApx(buf);
1468
1469 GET_BINTOFREQ
1470 GET_FREQTOBIN
1471
1472 int minbin = unit->minbin;
1473 int maxbin = unit->maxbin;
1474 if(minbin == -99){
1475 // Initialise bin ranges from frequency ranges
1476 minbin = unit->minbin = ((int)(unit->minfreq * freqtobin)) - 1;
1477 if(minbin >= numbins)
1478 minbin = unit->minbin = numbins - 1;
1479 maxbin = unit->maxbin = ((int)(unit->maxfreq * freqtobin)) - 1;
1480 if(maxbin > numbins)
1481 maxbin = unit->maxbin = numbins;
1482 }
1483
1484 // Start off assuming DC is the best...
1485 int peakbin=-1; // "-1" meaning DC here, since the DC is not included in the main list
1486 float peakmag;
1487 if(minbin == -1){
1488 peakmag = sc_abs(p->dc);
1489 minbin = 0;
1490 }else{
1491 peakmag = -9999.f; // Will ensure DC gets beaten if not in desired range
1492 }
1493
1494 // ...then check all the others. We neglect nyquist for efficiency purposes. Sorry.
1495 float mag;
1496 for(int i = minbin; i < maxbin; ++i){
1497 mag = p->bin[i].mag;
1498
1499 if(peakmag < mag){
1500 peakmag = mag;
1501 peakbin = i;
1502 }
1503 };
1504
1505 ZOUT0(0) = unit->outval = (peakbin+1) * bintofreq;
1506 ZOUT0(1) = unit->outval2 = peakmag;
1507
1508 }
1509
1510 ////////////////////////////////////////////////////////////////////////////////
1511
PV_MagSmooth_Ctor(PV_MagSmooth * unit)1512 void PV_MagSmooth_Ctor(PV_MagSmooth *unit)
1513 {
1514 SETCALC(PV_MagSmooth_next);
1515 ZOUT0(0) = ZIN0(0);
1516 unit->m_memory = NULL;
1517 }
1518
PV_MagSmooth_next(PV_MagSmooth * unit,int inNumSamples)1519 void PV_MagSmooth_next(PV_MagSmooth *unit, int inNumSamples)
1520 {
1521 PV_GET_BUF
1522
1523 SCPolarBuf *p = ToPolarApx(buf);
1524
1525 float* memory = unit->m_memory;
1526 if(memory==NULL){
1527 memory = unit->m_memory = (float*)RTAlloc(unit->mWorld, (numbins+2) * sizeof(float));
1528 // Now copy the first frame into the memory
1529 for (int i=0; i<numbins; ++i) {
1530 memory[i] = p->bin[i].mag;
1531 }
1532 memory[numbins] = p->dc;
1533 memory[numbins+1] = p->nyq;
1534 }
1535
1536 float factor = ZIN0(1);
1537 float onemfactor = 1.f - factor;
1538
1539 // The actual smoothing calculation:
1540 for (int i=0; i<numbins; ++i) {
1541 memory[i] = p->bin[i].mag = (memory[i ] * factor) + (p->bin[i].mag * onemfactor);
1542 }
1543 memory[numbins] = p->dc = (memory[numbins ] * factor) + (p->dc * onemfactor);
1544 memory[numbins+1] = p->nyq = (memory[numbins+1] * factor) + (p->nyq * onemfactor);
1545 }
1546
PV_MagSmooth_Dtor(PV_MagSmooth * unit)1547 void PV_MagSmooth_Dtor(PV_MagSmooth *unit)
1548 {
1549 if(unit->m_memory!=NULL){
1550 RTFree(unit->mWorld, unit->m_memory);
1551 }
1552 }
1553
1554 //////////////////////////////////////////////////////////////////////////////////////////////////
1555
PV_MagMulAdd_Ctor(PV_Unit * unit)1556 void PV_MagMulAdd_Ctor(PV_Unit *unit)
1557 {
1558 SETCALC(PV_MagMulAdd_next);
1559 ZOUT0(0) = ZIN0(0);
1560 }
1561
PV_MagMulAdd_next(PV_Unit * unit,int inNumSamples)1562 void PV_MagMulAdd_next(PV_Unit *unit, int inNumSamples)
1563 {
1564 PV_GET_BUF
1565
1566 SCPolarBuf *p = ToPolarApx(buf);
1567
1568 float m = ZIN0(1);
1569 float a = ZIN0(2);
1570
1571 p->dc = p->dc * m + a;
1572 p->nyq = p->nyq * m + a;
1573 for (int i=0; i<numbins; ++i) {
1574 p->bin[i].mag = p->bin[i].mag * m + a;
1575 }
1576 }
1577
1578 ////////////////////////////////////////////////////////////////////////////////////////////////////////
1579
1580 #define GET_BUF_AT1 \
1581 float fbufnum1 = ZIN0(1); \
1582 if (fbufnum1 < 0.f) { fbufnum1 = 0.f; } \
1583 if (fbufnum1 != unit->m_fbufnum) { \
1584 uint32 bufnum = (int)fbufnum1; \
1585 World *world = unit->mWorld; \
1586 if (bufnum >= world->mNumSndBufs) { \
1587 int localBufNum = bufnum - world->mNumSndBufs; \
1588 Graph *parent = unit->mParent; \
1589 if(localBufNum <= parent->localBufNum) { \
1590 unit->m_buf = parent->mLocalSndBufs + localBufNum; \
1591 } else { \
1592 bufnum = 0; \
1593 unit->m_buf = world->mSndBufs + bufnum; \
1594 } \
1595 } else { \
1596 unit->m_buf = world->mSndBufs + bufnum; \
1597 } \
1598 unit->m_fbufnum = fbufnum1; \
1599 } \
1600 SndBuf* buf1 = unit->m_buf; \
1601 LOCK_SNDBUF(buf1); \
1602 float *bufData __attribute__((__unused__)) = buf1->data; \
1603 uint32 bufChannels __attribute__((__unused__)) = buf1->channels; \
1604 uint32 bufSamples __attribute__((__unused__)) = buf1->samples; \
1605 uint32 bufFrames = buf1->frames; \
1606 int mask __attribute__((__unused__)) = buf1->mask; \
1607 int guardFrame __attribute__((__unused__)) = bufFrames - 2;
1608
PV_ExtractRepeat_Ctor(PV_ExtractRepeat * unit)1609 void PV_ExtractRepeat_Ctor(PV_ExtractRepeat *unit)
1610 {
1611 unit->m_logmags = NULL;
1612 unit->m_cursor = 0;
1613 unit->m_fbufnum = -1;
1614 SETCALC(PV_ExtractRepeat_next);
1615 ZOUT0(0) = ZIN0(0);
1616 }
1617
PV_ExtractRepeat_next(PV_ExtractRepeat * unit,int inNumSamples)1618 void PV_ExtractRepeat_next(PV_ExtractRepeat *unit, int inNumSamples)
1619 {
1620 PV_GET_BUF
1621 SCPolarBuf *p = ToPolarApx(buf);
1622 GET_BUF_AT1
1623
1624 if((numbins+2) != bufChannels){
1625 printf("PV_ExtractRepeat error: fft magnitude size != bufChannels, %i > %i\n", (numbins+2), bufChannels);
1626 return;
1627 }
1628
1629 float loopdursecs = ZIN0(2);
1630 float memorytimesecs = ZIN0(3);
1631 bool which = ZIN0(4) > 0.f;
1632 float ffthopsize = ZIN0(5);
1633 float thresh = ZIN0(6);
1634
1635 int loopdurframes = loopdursecs * FULLRATE / ((numbins+numbins+2) * ffthopsize);
1636 //printf("PV_ExtractRepeat info: loopdurframes is %i\n", bufFrames);
1637 if(loopdurframes > bufFrames){
1638 printf("PV_ExtractRepeat warning: loopdurframes > bufFrames, %i > %i\n", loopdurframes, bufFrames);
1639 loopdurframes = bufFrames;
1640 }
1641
1642 float* logmags = unit->m_logmags;
1643 if(logmags==NULL){
1644 logmags = unit->m_logmags = (float*)RTAlloc(unit->mWorld, (numbins+2) * sizeof(float));
1645 // Here we zero the buffer, not the logmags, but we do it here as first-run initialisation.
1646 memset(bufData, 0, bufChannels * bufFrames * sizeof(float));
1647 }
1648
1649 float mag;
1650 for (int i=0; i<numbins; ++i) {
1651 mag = p->bin[i].mag;
1652 logmags[i] = log(mag > SMALLEST_NUM_FOR_LOG ? mag : SMALLEST_NUM_FOR_LOG);
1653 }
1654 // Note, we store the dc and nyquist ON THE END
1655 mag = sc_abs(p->dc);
1656 logmags[numbins ] = log(mag > SMALLEST_NUM_FOR_LOG ? mag : SMALLEST_NUM_FOR_LOG);
1657 mag = sc_abs(p->nyq);
1658 logmags[numbins+1] = log(mag > SMALLEST_NUM_FOR_LOG ? mag : SMALLEST_NUM_FOR_LOG);
1659
1660 // Advance the recursive-buffer-cursor by one. Let's call the 'current' frame in the recursive buffer R.
1661 int cursor = unit->m_cursor + 1;
1662 if(cursor >= loopdurframes)
1663 cursor = 0;
1664 unit->m_cursor = cursor;
1665 float* curframe = bufData + (cursor * bufChannels);
1666
1667 // find (X-R) (the original authors would have taken abs of this too, but that seems harmful to me).
1668 // find which bins are (X-R)<t where t=1. Zero them out (in S).
1669 // (OR do the reverse; user-settable flag, applied here using xors.)
1670 for (int i=0; i<numbins; ++i)
1671 if((logmags[i] - curframe[i] < thresh) ^ which)
1672 p->bin[i].mag = 0;
1673 if((logmags[numbins ] - curframe[numbins ] < thresh) ^ which)
1674 p->dc = 0;
1675 if((logmags[numbins+1] - curframe[numbins+1] < thresh) ^ which)
1676 p->nyq = 0;
1677
1678 // Finally update the recursive buffer by writing X onto it crossfade-style.
1679 float writestrength = (memorytimesecs == 0.f) ? 0.f :
1680 exp(log001 / (memorytimesecs * FULLRATE / ((numbins+numbins+2) * ffthopsize)));
1681 float antiwritestrength = 1.f - writestrength;
1682 //printf("PV_ExtractRepeat info: memorytimesecs %g gives write strength %g\n", memorytimesecs, writestrength);
1683 for (int i=0; i<numbins+2; ++i) {
1684 curframe[i] = (curframe[i] * antiwritestrength) + (logmags[i] * writestrength);
1685 }
1686 }
1687
PV_ExtractRepeat_Dtor(PV_ExtractRepeat * unit)1688 void PV_ExtractRepeat_Dtor(PV_ExtractRepeat *unit)
1689 {
1690 if(unit->m_logmags) RTFree(unit->mWorld, unit->m_logmags);
1691 }
1692
1693 ////////////////////////////////////////////////////////////////////////////////////////////////////////
1694
PluginLoad(MCLDFFT)1695 PluginLoad(MCLDFFT)
1696 {
1697 ft= inTable;
1698
1699 init_SCComplex(inTable);
1700
1701 (*ft->fDefineUnit)("FFTPower", sizeof(FFTPower), (UnitCtorFunc)&FFTPower_Ctor, 0, 0);
1702 (*ft->fDefineUnit)("FFTFlux", sizeof(FFTFlux_Unit), (UnitCtorFunc)&FFTFlux_Ctor, (UnitDtorFunc)&FFTFlux_Dtor, 0);
1703 (*ft->fDefineUnit)("FFTFluxPos", sizeof(FFTFlux_Unit), (UnitCtorFunc)&FFTFluxPos_Ctor, (UnitDtorFunc)&FFTFluxPos_Dtor, 0);
1704 (*ft->fDefineUnit)("FFTDiffMags", sizeof(FFTAnalyser_Unit), (UnitCtorFunc)&FFTDiffMags_Ctor, 0, 0);
1705 DefineSimpleUnit(PV_MagSubtract);
1706 (*ft->fDefineUnit)("PV_MagLog", sizeof(PV_Unit), (UnitCtorFunc)&PV_MagLog_Ctor, 0, 0);
1707 (*ft->fDefineUnit)("PV_MagExp", sizeof(PV_Unit), (UnitCtorFunc)&PV_MagExp_Ctor, 0, 0);
1708 DefineDtorUnit(FFTSubbandPower);
1709
1710 DefineDtorUnit(FFTPhaseDev);
1711 DefineDtorUnit(FFTComplexDev);
1712 DefineDtorUnit(FFTMKL);
1713
1714 DefineSimpleUnit(PV_Whiten);
1715
1716 DefineSimpleUnit(FFTCrest);
1717 DefineSimpleUnit(FFTSpread);
1718 DefineSimpleUnit(FFTSlope);
1719
1720 DefineDtorUnit(FFTSubbandFlatness);
1721
1722 DefineSimpleUnit(FFTPeak);
1723
1724 DefineDtorUnit(PV_MagSmooth);
1725
1726 (*ft->fDefineUnit)("PV_MagMulAdd", sizeof(PV_Unit), (UnitCtorFunc)&PV_MagMulAdd_Ctor, 0, 0);
1727
1728 DefineDtorUnit(PV_ExtractRepeat);
1729 }
1730