server/plugins/FFT_UGens.cpp

   1 /*
   2         Improved FFT and IFFT UGens for SuperCollider 3
   3         Copyright (c) 2007-2008 Dan Stowell, incorporating code from
   4         SuperCollider 3 Copyright (c) 2002 James McCartney.
   5         All rights reserved.
   6
   7     This program is free software; you can redistribute it and/or modify
   8     it under the terms of the GNU General Public License as published by
   9     the Free Software Foundation; either version 2 of the License, or
  10     (at your option) any later version.
  11
  12     This program is distributed in the hope that it will be useful,
  13     but WITHOUT ANY WARRANTY; without even the implied warranty of
  14     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15     GNU General Public License for more details.
  16
  17     You should have received a copy of the GNU General Public License
  18     along with this program; if not, write to the Free Software
  19     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA
  20 */
  21
  22
  23 #include "FFT_UGens.h"
  24
  25 // We include vDSP even if not using for FFT, since we want to use some vectorised add/mul tricks
  26 #if defined(__APPLE__) && !defined(SC_IPHONE)
  27 #include "vecLib/vDSP.h"
  28 #endif
  29
  30 struct FFTBase : public Unit
  31 {
  32         SndBuf *m_fftsndbuf;
  33         float *m_fftbuf;
  34
  35         int m_pos, m_fullbufsize, m_audiosize; // "fullbufsize" includes any zero-padding, "audiosize" does not.
  36         int m_log2n_full, m_log2n_audio;
  37
  38         uint32 m_fftbufnum;
  39
  40         scfft* m_scfft;
  41
  42         int m_hopsize, m_shuntsize; // These add up to m_audiosize
  43         int m_wintype;
  44
  45         int m_numSamples;
  46 };
  47
  48 struct FFT : public FFTBase
  49 {
  50         float *m_inbuf;
  51 };
  52
  53 struct IFFT : public FFTBase
  54 {
  55         float *m_olabuf;
  56         int m_numSamples;
  57 };
  58
  59 struct FFTTrigger : public FFTBase
  60 {
  61         int m_numPeriods, m_periodsRemain, m_polar;
  62 };
  63
  64 //////////////////////////////////////////////////////////////////////////////////////////////////
  65
  66 extern "C"
  67 {
  68         void FFT_Ctor(FFT* unit);
  69         void FFT_ClearUnitOutputs(FFT *unit, int wrongNumSamples);
  70         void FFT_next(FFT *unit, int inNumSamples);
  71         void FFT_Dtor(FFT* unit);
  72
  73         void IFFT_Ctor(IFFT* unit);
  74         void IFFT_next(IFFT *unit, int inNumSamples);
  75         void IFFT_Dtor(IFFT* unit);
  76
  77         void FFTTrigger_Ctor(FFTTrigger* unit);
  78         void FFTTrigger_next(FFTTrigger* unit, int inNumSamples);
  79 }
  80
  81 //////////////////////////////////////////////////////////////////////////////////////////////////
  82
  83 static int FFTBase_Ctor(FFTBase *unit, int frmsizinput)
  84 {
  85         World *world = unit->mWorld;
  86
  87         uint32 bufnum = (uint32)ZIN0(0);
  88         SndBuf *buf;
  89         if (bufnum >= world->mNumSndBufs) {
  90                 int localBufNum = bufnum - world->mNumSndBufs;
  91                 Graph *parent = unit->mParent;
  92                 if(localBufNum <= parent->localMaxBufNum) {
  93                         buf = parent->mLocalSndBufs + localBufNum;
  94                 } else {
  95                         if(unit->mWorld->mVerbosity > -1){ Print("FFTBase_Ctor error: invalid buffer number: %i.\n", bufnum); }
  96                         return 0;
  97                 }
  98         } else {
  99                 buf = world->mSndBufs + bufnum;
 100         }
 101
 102         if (!buf->data) {
 103                 if(unit->mWorld->mVerbosity > -1){ Print("FFTBase_Ctor error: Buffer %i not initialised.\n", bufnum); }
 104                 return 0;
 105         }
 106
 107         unit->m_fftsndbuf = buf;
 108         unit->m_fftbufnum = bufnum;
 109         unit->m_fullbufsize = buf->samples;
 110         int framesize = (int)ZIN0(frmsizinput);
 111         if(framesize < 1)
 112                 unit->m_audiosize = buf->samples;
 113         else
 114                 unit->m_audiosize = sc_min(buf->samples, framesize);
 115
 116         unit->m_log2n_full  = LOG2CEIL(unit->m_fullbufsize);
 117         unit->m_log2n_audio = LOG2CEIL(unit->m_audiosize);
 118
 119
 120         // Although FFTW allows non-power-of-two buffers (vDSP doesn't), this would complicate the windowing, so we don't allow it.
 121         if (!ISPOWEROFTWO(unit->m_fullbufsize)) {
 122                 Print("FFTBase_Ctor error: buffer size (%i) not a power of two.\n", unit->m_fullbufsize);
 123                 return 0;
 124         }
 125         else if (!ISPOWEROFTWO(unit->m_audiosize)) {
 126                 Print("FFTBase_Ctor error: audio frame size (%i) not a power of two.\n", unit->m_audiosize);
 127                 return 0;
 128         }
 129         else if (unit->m_audiosize < SC_FFT_MINSIZE ||
 130                         (((int)(unit->m_audiosize / unit->mWorld->mFullRate.mBufLength))
 131                                         * unit->mWorld->mFullRate.mBufLength != unit->m_audiosize)) {
 132                 Print("FFTBase_Ctor error: audio frame size (%i) not a multiple of the block size (%i).\n", unit->m_audiosize, unit->mWorld->mFullRate.mBufLength);
 133                 return 0;
 134         }
 135
 136         unit->m_pos = 0;
 137
 138         ZOUT0(0) = ZIN0(0);
 139
 140         return 1;
 141 }
 142
 143 //////////////////////////////////////////////////////////////////////////////////////////////////
 144
 145 void FFT_Ctor(FFT *unit)
 146 {
 147         unit->m_wintype = (int)ZIN0(3); // wintype may be used by the base ctor
 148         if(!FFTBase_Ctor(unit, 5)){
 149                 SETCALC(FFT_ClearUnitOutputs);
 150                 // These zeroes are to prevent the dtor freeing things that don't exist:
 151                 unit->m_inbuf = 0;
 152                 unit->m_scfft = 0;
 153                 return;
 154         }
 155         int fullbufsize = unit->m_fullbufsize * sizeof(float);
 156         int audiosize = unit->m_audiosize * sizeof(float);
 157
 158         int hopsize = (int)(sc_max(sc_min(ZIN0(2), 1.f), 0.f) * unit->m_audiosize);
 159         if (hopsize < unit->mWorld->mFullRate.mBufLength) {
 160                 Print("FFT_Ctor: hopsize smaller than SC's block size (%i) - automatically corrected.\n", hopsize, unit->mWorld->mFullRate.mBufLength);
 161                 hopsize = unit->mWorld->mFullRate.mBufLength;
 162         } else if (((int)(hopsize / unit->mWorld->mFullRate.mBufLength)) * unit->mWorld->mFullRate.mBufLength
 163                                 != hopsize) {
 164                 Print("FFT_Ctor: hopsize (%i) not an exact multiple of SC's block size (%i) - automatically corrected.\n", hopsize, unit->mWorld->mFullRate.mBufLength);
 165                 hopsize = ((int)(hopsize / unit->mWorld->mFullRate.mBufLength)) * unit->mWorld->mFullRate.mBufLength;
 166         }
 167         unit->m_hopsize = hopsize;
 168         unit->m_shuntsize = unit->m_audiosize - hopsize;
 169
 170         unit->m_inbuf = (float*)RTAlloc(unit->mWorld, audiosize);
 171
 172         SCWorld_Allocator alloc(ft, unit->mWorld);
 173         unit->m_scfft = scfft_create(unit->m_fullbufsize, unit->m_audiosize, (SCFFT_WindowFunction)unit->m_wintype, unit->m_inbuf,
 174                                                                  unit->m_fftsndbuf->data, kForward, alloc);
 175
 176         memset(unit->m_inbuf, 0, audiosize);
 177
 178         //Print("FFT_Ctor: hopsize %i, shuntsize %i, bufsize %i, wintype %i, \n",
 179         //      unit->m_hopsize, unit->m_shuntsize, unit->m_bufsize, unit->m_wintype);
 180
 181         if (INRATE(1) == calc_FullRate) {
 182                 unit->m_numSamples = unit->mWorld->mFullRate.mBufLength;
 183         } else {
 184                 unit->m_numSamples = 1;
 185         }
 186
 187         SETCALC(FFT_next);
 188 }
 189
 190 void FFT_Dtor(FFT *unit)
 191 {
 192         SCWorld_Allocator alloc(ft, unit->mWorld);
 193         if(unit->m_scfft)
 194                 scfft_destroy(unit->m_scfft, alloc);
 195
 196         if(unit->m_inbuf)
 197                 RTFree(unit->mWorld, unit->m_inbuf);
 198 }
 199
 200 // Ordinary ClearUnitOutputs outputs zero, potentially telling the IFFT (+ PV UGens) to act on buffer zero, so let's skip that:
 201 void FFT_ClearUnitOutputs(FFT *unit, int wrongNumSamples)
 202 {
 203         ZOUT0(0) = -1;
 204 }
 205
 206 void FFT_next(FFT *unit, int wrongNumSamples)
 207 {
 208         float *in = IN(1);
 209         float *out = unit->m_inbuf + unit->m_pos + unit->m_shuntsize;
 210
 211 //      int numSamples = unit->mWorld->mFullRate.mBufLength;
 212         int numSamples = unit->m_numSamples;
 213
 214         // copy input
 215         memcpy(out, in, numSamples * sizeof(float));
 216
 217         unit->m_pos += numSamples;
 218
 219         bool gate = ZIN0(4) > 0.f; // Buffer shunting continues, but no FFTing
 220
 221         if (unit->m_pos != unit->m_hopsize || !unit->m_fftsndbuf->data || unit->m_fftsndbuf->samples != unit->m_fullbufsize) {
 222                 if(unit->m_pos == unit->m_hopsize)
 223                         unit->m_pos = 0;
 224                 ZOUT0(0) = -1.f;
 225         } else {
 226
 227                 unit->m_pos = 0;
 228                 if(gate){
 229                         scfft_dofft(unit->m_scfft);
 230                         unit->m_fftsndbuf->coord = coord_Complex;
 231                         ZOUT0(0) = unit->m_fftbufnum;
 232                 } else {
 233                         ZOUT0(0) = -1;
 234                 }
 235                 // Shunt input buf down
 236                 memcpy(unit->m_inbuf, unit->m_inbuf + unit->m_hopsize, unit->m_shuntsize * sizeof(float));
 237         }
 238 }
 239
 240 /////////////////////////////////////////////////////////////////////////////////////////////
 241
 242 void IFFT_Ctor(IFFT* unit){
 243         unit->m_wintype = (int)ZIN0(1); // wintype may be used by the base ctor
 244         if(!FFTBase_Ctor(unit, 2)){
 245                 SETCALC(*ClearUnitOutputs);
 246                 // These zeroes are to prevent the dtor freeing things that don't exist:
 247                 unit->m_olabuf = 0;
 248                 return;
 249         }
 250
 251         // This will hold the transformed and progressively overlap-added data ready for outputting.
 252         unit->m_olabuf = (float*)RTAlloc(unit->mWorld, unit->m_audiosize * sizeof(float));
 253         memset(unit->m_olabuf, 0, unit->m_audiosize * sizeof(float));
 254
 255         SCWorld_Allocator alloc(ft, unit->mWorld);
 256         unit->m_scfft = scfft_create(unit->m_fullbufsize, unit->m_audiosize, (SCFFT_WindowFunction)unit->m_wintype, unit->m_fftsndbuf->data,
 257                                                                  unit->m_fftsndbuf->data, kBackward, alloc);
 258
 259         // "pos" will be reset to zero when each frame comes in. Until then, the following ensures silent output at first:
 260         unit->m_pos = 0; //unit->m_audiosize;
 261
 262         if (unit->mCalcRate == calc_FullRate) {
 263                 unit->m_numSamples = unit->mWorld->mFullRate.mBufLength;
 264         } else {
 265                 unit->m_numSamples = 1;
 266         }
 267
 268         SETCALC(IFFT_next);
 269 }
 270
 271 void IFFT_Dtor(IFFT* unit)
 272 {
 273         if(unit->m_olabuf)
 274                 RTFree(unit->mWorld, unit->m_olabuf);
 275
 276         SCWorld_Allocator alloc(ft, unit->mWorld);
 277         if(unit->m_scfft)
 278                 scfft_destroy(unit->m_scfft, alloc);
 279 }
 280
 281 void IFFT_next(IFFT *unit, int wrongNumSamples)
 282 {
 283         float *out = OUT(0); // NB not ZOUT0
 284
 285         // Load state from struct into local scope
 286         int pos     = unit->m_pos;
 287         int fullbufsize  = unit->m_fullbufsize;
 288         int audiosize = unit->m_audiosize;
 289 //      int numSamples = unit->mWorld->mFullRate.mBufLength;
 290         int numSamples = unit->m_numSamples;
 291         float *olabuf = unit->m_olabuf;
 292         float fbufnum = ZIN0(0);
 293
 294         // Only run the IFFT if we're receiving a new block of input data - otherwise just output data already received
 295         if (fbufnum >= 0.f){
 296                 // Ensure it's in cartesian format, not polar
 297                 ToComplexApx(unit->m_fftsndbuf);
 298
 299                 float* fftbuf = unit->m_fftsndbuf->data;
 300
 301                 scfft_doifft(unit->m_scfft);
 302
 303                 // Then shunt the "old" time-domain output down by one hop
 304                 int hopsamps = pos;
 305                 int shuntsamps = audiosize - hopsamps;
 306                 if(hopsamps != audiosize)  // There's only copying to be done if the position isn't all the way to the end of the buffer
 307                         memcpy(olabuf, olabuf+hopsamps, shuntsamps * sizeof(float));
 308
 309                 // Then mix the "new" time-domain data in - adding at first, then just setting (copying) where the "old" is supposed to be zero.
 310                 #if defined(__APPLE__) && !defined(SC_IPHONE)
 311                         vDSP_vadd(olabuf, 1, fftbuf, 1, olabuf, 1, shuntsamps);
 312                 #else
 313                         // NB we re-use the "pos" variable temporarily here for write rather than read
 314                         for(pos = 0; pos < shuntsamps; ++pos){
 315                                 olabuf[pos] += fftbuf[pos];
 316                         }
 317                 #endif
 318                 memcpy(olabuf + shuntsamps, fftbuf + shuntsamps, (hopsamps) * sizeof(float));
 319
 320                 // Move the pointer back to zero, which is where playback will next begin
 321                 pos = 0;
 322
 323         } // End of has-the-chain-fired
 324
 325         // Now we can output some stuff, as long as there is still data waiting to be output.
 326         // If there is NOT data waiting to be output, we output zero. (Either irregular/negative-overlap
 327         //     FFT firing, or FFT has given up, or at very start of execution.)
 328         if(pos >= audiosize)
 329                 ClearUnitOutputs(unit, numSamples);
 330         else {
 331                 memcpy(out, olabuf + pos, numSamples * sizeof(float));
 332                 pos += numSamples;
 333         }
 334         unit->m_pos = pos;
 335 }
 336
 337 /////////////////////////////////////////////////////////////////////////////////////////////
 338
 339 void FFTTrigger_Ctor(FFTTrigger *unit)
 340 {
 341         World *world = unit->mWorld;
 342
 343 /*
 344         uint32 bufnum = (uint32)IN0(0);
 345         Print("FFTTrigger_Ctor: bufnum is %i\n", bufnum);
 346         if (bufnum >= world->mNumSndBufs) bufnum = 0;
 347         SndBuf *buf = world->mSndBufs + bufnum;
 348 */
 349
 350
 351         uint32 bufnum = (uint32)IN0(0);
 352         //Print("FFTTrigger_Ctor: bufnum is %i\n", bufnum);
 353         SndBuf *buf;
 354         if (bufnum >= world->mNumSndBufs) {
 355                 int localBufNum = bufnum - world->mNumSndBufs;
 356                 Graph *parent = unit->mParent;
 357                 if(localBufNum <= parent->localMaxBufNum) {
 358                         buf = parent->mLocalSndBufs + localBufNum;
 359                 } else {
 360                         bufnum = 0;
 361                         buf = world->mSndBufs + bufnum;
 362                 }
 363         } else {
 364                 buf = world->mSndBufs + bufnum;
 365         }
 366         LOCK_SNDBUF(buf);
 367
 368
 369         unit->m_fftsndbuf = buf;
 370         unit->m_fftbufnum = bufnum;
 371         unit->m_fullbufsize = buf->samples;
 372
 373         int numSamples = unit->mWorld->mFullRate.mBufLength;
 374         float dataHopSize = IN0(1);
 375         int initPolar = unit->m_polar = (int)IN0(2);
 376         unit->m_numPeriods = unit->m_periodsRemain = (int)(((float)unit->m_fullbufsize * dataHopSize) / numSamples) - 1;
 377
 378         buf->coord = (IN0(2) == 1.f) ? coord_Polar : coord_Complex;
 379
 380         OUT0(0) = IN0(0);
 381         SETCALC(FFTTrigger_next);
 382 }
 383
 384 void FFTTrigger_next(FFTTrigger *unit, int inNumSamples)
 385 {
 386         if (unit->m_periodsRemain > 0) {
 387                 ZOUT0(0) = -1.f;
 388                 unit->m_periodsRemain--;
 389         } else {
 390                 ZOUT0(0) = unit->m_fftbufnum;
 391                 unit->m_pos = 0;
 392                 unit->m_periodsRemain = unit->m_numPeriods;
 393         }
 394
 395 }
 396
 397 void initFFT(InterfaceTable *inTable)
 398 {
 399         ft = inTable;
 400
 401         DefineDtorUnit(FFT);
 402         DefineDtorUnit(IFFT);
 403         DefineSimpleUnit(FFTTrigger);
 404 }
 405