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 audiosize = unit->m_audiosize * sizeof(float);
 156
 157         int hopsize = (int)(sc_max(sc_min(ZIN0(2), 1.f), 0.f) * unit->m_audiosize);
 158         if (hopsize < unit->mWorld->mFullRate.mBufLength) {
 159                 Print("FFT_Ctor: hopsize smaller than SC's block size (%i) - automatically corrected.\n", hopsize, unit->mWorld->mFullRate.mBufLength);
 160                 hopsize = unit->mWorld->mFullRate.mBufLength;
 161         } else if (((int)(hopsize / unit->mWorld->mFullRate.mBufLength)) * unit->mWorld->mFullRate.mBufLength
 162                                 != hopsize) {
 163                 Print("FFT_Ctor: hopsize (%i) not an exact multiple of SC's block size (%i) - automatically corrected.\n", hopsize, unit->mWorld->mFullRate.mBufLength);
 164                 hopsize = ((int)(hopsize / unit->mWorld->mFullRate.mBufLength)) * unit->mWorld->mFullRate.mBufLength;
 165         }
 166         unit->m_hopsize = hopsize;
 167         unit->m_shuntsize = unit->m_audiosize - hopsize;
 168
 169         unit->m_inbuf = (float*)RTAlloc(unit->mWorld, audiosize);
 170
 171         SCWorld_Allocator alloc(ft, unit->mWorld);
 172         unit->m_scfft = scfft_create(unit->m_fullbufsize, unit->m_audiosize, (SCFFT_WindowFunction)unit->m_wintype, unit->m_inbuf,
 173                                                                  unit->m_fftsndbuf->data, kForward, alloc);
 174
 175         memset(unit->m_inbuf, 0, audiosize);
 176
 177         //Print("FFT_Ctor: hopsize %i, shuntsize %i, bufsize %i, wintype %i, \n",
 178         //      unit->m_hopsize, unit->m_shuntsize, unit->m_bufsize, unit->m_wintype);
 179
 180         if (INRATE(1) == calc_FullRate) {
 181                 unit->m_numSamples = unit->mWorld->mFullRate.mBufLength;
 182         } else {
 183                 unit->m_numSamples = 1;
 184         }
 185
 186         SETCALC(FFT_next);
 187 }
 188
 189 void FFT_Dtor(FFT *unit)
 190 {
 191         SCWorld_Allocator alloc(ft, unit->mWorld);
 192         if(unit->m_scfft)
 193                 scfft_destroy(unit->m_scfft, alloc);
 194
 195         if(unit->m_inbuf)
 196                 RTFree(unit->mWorld, unit->m_inbuf);
 197 }
 198
 199 // Ordinary ClearUnitOutputs outputs zero, potentially telling the IFFT (+ PV UGens) to act on buffer zero, so let's skip that:
 200 void FFT_ClearUnitOutputs(FFT *unit, int wrongNumSamples)
 201 {
 202         ZOUT0(0) = -1;
 203 }
 204
 205 void FFT_next(FFT *unit, int wrongNumSamples)
 206 {
 207         float *in = IN(1);
 208         float *out = unit->m_inbuf + unit->m_pos + unit->m_shuntsize;
 209
 210 //      int numSamples = unit->mWorld->mFullRate.mBufLength;
 211         int numSamples = unit->m_numSamples;
 212
 213         // copy input
 214         memcpy(out, in, numSamples * sizeof(float));
 215
 216         unit->m_pos += numSamples;
 217
 218         bool gate = ZIN0(4) > 0.f; // Buffer shunting continues, but no FFTing
 219
 220         if (unit->m_pos != unit->m_hopsize || !unit->m_fftsndbuf->data || unit->m_fftsndbuf->samples != unit->m_fullbufsize) {
 221                 if(unit->m_pos == unit->m_hopsize)
 222                         unit->m_pos = 0;
 223                 ZOUT0(0) = -1.f;
 224         } else {
 225
 226                 unit->m_pos = 0;
 227                 if(gate){
 228                         scfft_dofft(unit->m_scfft);
 229                         unit->m_fftsndbuf->coord = coord_Complex;
 230                         ZOUT0(0) = unit->m_fftbufnum;
 231                 } else {
 232                         ZOUT0(0) = -1;
 233                 }
 234                 // Shunt input buf down
 235                 memcpy(unit->m_inbuf, unit->m_inbuf + unit->m_hopsize, unit->m_shuntsize * sizeof(float));
 236         }
 237 }
 238
 239 /////////////////////////////////////////////////////////////////////////////////////////////
 240
 241 void IFFT_Ctor(IFFT* unit){
 242         unit->m_wintype = (int)ZIN0(1); // wintype may be used by the base ctor
 243         if(!FFTBase_Ctor(unit, 2)){
 244                 SETCALC(*ClearUnitOutputs);
 245                 // These zeroes are to prevent the dtor freeing things that don't exist:
 246                 unit->m_olabuf = 0;
 247                 return;
 248         }
 249
 250         // This will hold the transformed and progressively overlap-added data ready for outputting.
 251         unit->m_olabuf = (float*)RTAlloc(unit->mWorld, unit->m_audiosize * sizeof(float));
 252         memset(unit->m_olabuf, 0, unit->m_audiosize * sizeof(float));
 253
 254         SCWorld_Allocator alloc(ft, unit->mWorld);
 255         unit->m_scfft = scfft_create(unit->m_fullbufsize, unit->m_audiosize, (SCFFT_WindowFunction)unit->m_wintype, unit->m_fftsndbuf->data,
 256                                                                  unit->m_fftsndbuf->data, kBackward, alloc);
 257
 258         // "pos" will be reset to zero when each frame comes in. Until then, the following ensures silent output at first:
 259         unit->m_pos = 0; //unit->m_audiosize;
 260
 261         if (unit->mCalcRate == calc_FullRate) {
 262                 unit->m_numSamples = unit->mWorld->mFullRate.mBufLength;
 263         } else {
 264                 unit->m_numSamples = 1;
 265         }
 266
 267         SETCALC(IFFT_next);
 268 }
 269
 270 void IFFT_Dtor(IFFT* unit)
 271 {
 272         if(unit->m_olabuf)
 273                 RTFree(unit->mWorld, unit->m_olabuf);
 274
 275         SCWorld_Allocator alloc(ft, unit->mWorld);
 276         if(unit->m_scfft)
 277                 scfft_destroy(unit->m_scfft, alloc);
 278 }
 279
 280 void IFFT_next(IFFT *unit, int wrongNumSamples)
 281 {
 282         float *out = OUT(0); // NB not ZOUT0
 283
 284         // Load state from struct into local scope
 285         int pos     = unit->m_pos;
 286         int audiosize = unit->m_audiosize;
 287 //      int numSamples = unit->mWorld->mFullRate.mBufLength;
 288         int numSamples = unit->m_numSamples;
 289         float *olabuf = unit->m_olabuf;
 290         float fbufnum = ZIN0(0);
 291
 292         // Only run the IFFT if we're receiving a new block of input data - otherwise just output data already received
 293         if (fbufnum >= 0.f){
 294                 // Ensure it's in cartesian format, not polar
 295                 ToComplexApx(unit->m_fftsndbuf);
 296
 297                 float* fftbuf = unit->m_fftsndbuf->data;
 298
 299                 scfft_doifft(unit->m_scfft);
 300
 301                 // Then shunt the "old" time-domain output down by one hop
 302                 int hopsamps = pos;
 303                 int shuntsamps = audiosize - hopsamps;
 304                 if(hopsamps != audiosize)  // There's only copying to be done if the position isn't all the way to the end of the buffer
 305                         memcpy(olabuf, olabuf+hopsamps, shuntsamps * sizeof(float));
 306
 307                 // Then mix the "new" time-domain data in - adding at first, then just setting (copying) where the "old" is supposed to be zero.
 308                 #if defined(__APPLE__) && !defined(SC_IPHONE)
 309                         vDSP_vadd(olabuf, 1, fftbuf, 1, olabuf, 1, shuntsamps);
 310                 #else
 311                         // NB we re-use the "pos" variable temporarily here for write rather than read
 312                         for(pos = 0; pos < shuntsamps; ++pos){
 313                                 olabuf[pos] += fftbuf[pos];
 314                         }
 315                 #endif
 316                 memcpy(olabuf + shuntsamps, fftbuf + shuntsamps, (hopsamps) * sizeof(float));
 317
 318                 // Move the pointer back to zero, which is where playback will next begin
 319                 pos = 0;
 320
 321         } // End of has-the-chain-fired
 322
 323         // Now we can output some stuff, as long as there is still data waiting to be output.
 324         // If there is NOT data waiting to be output, we output zero. (Either irregular/negative-overlap
 325         //     FFT firing, or FFT has given up, or at very start of execution.)
 326         if(pos >= audiosize)
 327                 ClearUnitOutputs(unit, numSamples);
 328         else {
 329                 memcpy(out, olabuf + pos, numSamples * sizeof(float));
 330                 pos += numSamples;
 331         }
 332         unit->m_pos = pos;
 333 }
 334
 335 /////////////////////////////////////////////////////////////////////////////////////////////
 336
 337 void FFTTrigger_Ctor(FFTTrigger *unit)
 338 {
 339         World *world = unit->mWorld;
 340
 341 /*
 342         uint32 bufnum = (uint32)IN0(0);
 343         Print("FFTTrigger_Ctor: bufnum is %i\n", bufnum);
 344         if (bufnum >= world->mNumSndBufs) bufnum = 0;
 345         SndBuf *buf = world->mSndBufs + bufnum;
 346 */
 347
 348
 349         uint32 bufnum = (uint32)IN0(0);
 350         //Print("FFTTrigger_Ctor: bufnum is %i\n", bufnum);
 351         SndBuf *buf;
 352         if (bufnum >= world->mNumSndBufs) {
 353                 int localBufNum = bufnum - world->mNumSndBufs;
 354                 Graph *parent = unit->mParent;
 355                 if(localBufNum <= parent->localMaxBufNum) {
 356                         buf = parent->mLocalSndBufs + localBufNum;
 357                 } else {
 358                         bufnum = 0;
 359                         buf = world->mSndBufs + bufnum;
 360                 }
 361         } else {
 362                 buf = world->mSndBufs + bufnum;
 363         }
 364         LOCK_SNDBUF(buf);
 365
 366
 367         unit->m_fftsndbuf = buf;
 368         unit->m_fftbufnum = bufnum;
 369         unit->m_fullbufsize = buf->samples;
 370
 371         int numSamples = unit->mWorld->mFullRate.mBufLength;
 372         float dataHopSize = IN0(1);
 373         unit->m_numPeriods = unit->m_periodsRemain = (int)(((float)unit->m_fullbufsize * dataHopSize) / numSamples) - 1;
 374
 375         buf->coord = (IN0(2) == 1.f) ? coord_Polar : coord_Complex;
 376
 377         OUT0(0) = IN0(0);
 378         SETCALC(FFTTrigger_next);
 379 }
 380
 381 void FFTTrigger_next(FFTTrigger *unit, int inNumSamples)
 382 {
 383         if (unit->m_periodsRemain > 0) {
 384                 ZOUT0(0) = -1.f;
 385                 unit->m_periodsRemain--;
 386         } else {
 387                 ZOUT0(0) = unit->m_fftbufnum;
 388                 unit->m_pos = 0;
 389                 unit->m_periodsRemain = unit->m_numPeriods;
 390         }
 391
 392 }
 393
 394 void initFFT(InterfaceTable *inTable)
 395 {
 396         ft = inTable;
 397
 398         DefineDtorUnit(FFT);
 399         DefineDtorUnit(IFFT);
 400         DefineSimpleUnit(FFTTrigger);
 401 }
 402