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