server/plugins/ML_SpecStats.cpp

   1 /*
   2
   3 Spectral statistics UGens for SuperCollider, by Dan Stowell.
   4 Copyright (c) Dan Stowell 2006-2007.
   5 Now part of SuperCollider 3, (c) James McCartney.
   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 "SC_PlugIn.h"
  24 #include "SCComplex.h"
  25 #include "FFT_UGens.h"
  26
  27 #include "ML.h"
  28
  29 //////////////////////////////////////////////////////////////////////////////////////////////////
  30
  31 /*
  32 struct FFTAnalyser_Unit : Unit
  33 {
  34         float outval;
  35
  36         // Not always used: multipliers which convert from bin indices to freq vals, and vice versa.
  37         // See also the macros for deriving these.
  38         float m_bintofreq; // , m_freqtobin;
  39 };
  40
  41 struct FFTAnalyser_OutOfPlace : FFTAnalyser_Unit
  42 {
  43         int m_numbins;
  44         float *m_tempbuf;
  45 };
  46
  47 struct SpecFlatness : FFTAnalyser_Unit
  48 {
  49 };
  50
  51 struct SpecPcile : FFTAnalyser_OutOfPlace
  52 {
  53         bool m_interpolate;
  54 };
  55
  56 struct SpecCentroid : FFTAnalyser_Unit
  57 {
  58 };
  59 */
  60
  61 //////////////////////////////////////////////////////////////////////////////////////////////////
  62
  63
  64 // for operation on one buffer
  65 // just like PV_GET_BUF except it outputs unit->outval rather than -1 when FFT not triggered
  66 #define FFTAnalyser_GET_BUF \
  67         float fbufnum = ZIN0(0); \
  68         if (fbufnum < 0.f) { ZOUT0(0) = unit->outval; return; } \
  69         ZOUT0(0) = fbufnum; \
  70         uint32 ibufnum = (uint32)fbufnum; \
  71         World *world = unit->mWorld; \
  72         SndBuf *buf; \
  73         if (!(ibufnum < world->mNumSndBufs)) { \
  74                 int localBufNum = ibufnum - world->mNumSndBufs; \
  75                 Graph *parent = unit->mParent; \
  76                 if(!(localBufNum > parent->localBufNum)) { \
  77                         buf = parent->mLocalSndBufs + localBufNum; \
  78                 } else { \
  79                         buf = world->mSndBufs; \
  80                 } \
  81         } else { \
  82                 buf = world->mSndBufs + ibufnum; \
  83         } \
  84         LOCK_SNDBUF(buf); \
  85         int numbins = buf->samples - 2 >> 1;
  86
  87 // Copied from FFT_UGens.cpp
  88 #define GET_BINTOFREQ \
  89         if(unit->m_bintofreq==0.f){ \
  90                 unit->m_bintofreq = world->mFullRate.mSampleRate / buf->samples; \
  91         } \
  92         float bintofreq = unit->m_bintofreq;
  93 /*
  94 #define GET_FREQTOBIN \
  95         if(unit->m_freqtobin==0.f){ \
  96                 unit->m_freqtobin = buf->samples / world->mFullRate.mSampleRate; \
  97         } \
  98         float freqtobin = unit->m_freqtobin;
  99 */
 100
 101 //////////////////////////////////////////////////////////////////////////////////////////////////
 102 /*
 103 extern "C"
 104 {
 105 //      void load(InterfaceTable *inTable);
 106
 107         void SpecFlatness_Ctor(SpecFlatness *unit);
 108         void SpecFlatness_next(SpecFlatness *unit, int inNumSamples);
 109
 110         void SpecPcile_Ctor(SpecPcile *unit);
 111         void SpecPcile_next(SpecPcile *unit, int inNumSamples);
 112         void SpecPcile_Dtor(SpecPcile *unit);
 113
 114         void SpecCentroid_Ctor(SpecCentroid *unit);
 115         void SpecCentroid_next(SpecCentroid *unit, int inNumSamples);
 116
 117 }
 118 */
 119 /*
 120 SCPolarBuf* ToPolarApx(SndBuf *buf)
 121 {
 122         if (buf->coord == coord_Complex) {
 123                 SCComplexBuf* p = (SCComplexBuf*)buf->data;
 124                 int numbins = buf->samples - 2 >> 1;
 125                 for (int i=0; i<numbins; ++i) {
 126                         p->bin[i].ToPolarApxInPlace();
 127                 }
 128                 buf->coord = coord_Polar;
 129         }
 130
 131         return (SCPolarBuf*)buf->data;
 132 }
 133
 134 SCComplexBuf* ToComplexApx(SndBuf *buf)
 135 {
 136         if (buf->coord == coord_Polar) {
 137                 SCPolarBuf* p = (SCPolarBuf*)buf->data;
 138                 int numbins = buf->samples - 2 >> 1;
 139                 for (int i=0; i<numbins; ++i) {
 140                         p->bin[i].ToComplexApxInPlace();
 141                 }
 142                 buf->coord = coord_Complex;
 143         }
 144         return (SCComplexBuf*)buf->data;
 145 }
 146
 147 InterfaceTable *ft;
 148
 149 void init_SCComplex(InterfaceTable *inTable);
 150 */
 151 //////////////////////////////////////////////////////////////////////////////////////////////////
 152
 153 void SpecFlatness_Ctor(SpecFlatness *unit)
 154 {
 155         SETCALC(SpecFlatness_next);
 156         ZOUT0(0) = unit->outval = 0.;
 157         unit->m_oneovern = 0.;
 158 }
 159
 160 void SpecFlatness_next(SpecFlatness *unit, int inNumSamples)
 161 {
 162         FFTAnalyser_GET_BUF
 163         if(unit->m_oneovern == 0.)
 164                 unit->m_oneovern = 1./(numbins + 2);
 165
 166         SCComplexBuf *p = ToComplexApx(buf);
 167
 168         // Spectral Flatness Measure is geometric mean divided by arithmetic mean.
 169         //
 170         // In order to calculate geom mean without hitting the precision limit,
 171         //  we use the trick of converting to log, taking the average, then converting back from log.
 172         double geommean = std::log(sc_abs(p->dc)) + std::log(sc_abs(p->nyq));
 173         double mean     = sc_abs(p->dc)      + sc_abs(p->nyq);
 174
 175         for (int i=0; i<numbins; ++i) {
 176                 float rabs = (p->bin[i].real);
 177                 float iabs = (p->bin[i].imag);
 178                 float amp = std::sqrt((rabs*rabs) + (iabs*iabs));
 179                 if(amp != 0.f) { // zeroes lead to NaNs
 180                         geommean += std::log(amp);
 181                         mean += amp;
 182                 }
 183         }
 184
 185         double oneovern = unit->m_oneovern;
 186         geommean = exp(geommean * oneovern); // Average and then convert back to linear
 187         mean *= oneovern;
 188
 189         // Store the val for output in future calls
 190         unit->outval = geommean / mean;
 191
 192         ZOUT0(0) = unit->outval;
 193 }
 194
 195 ////////////////////////////////////////////////////////////////////////////////////
 196
 197 void SpecPcile_Ctor(SpecPcile *unit)
 198 {
 199         SETCALC(SpecPcile_next);
 200
 201         unit->m_interpolate = ZIN0(2) > 0.f;
 202
 203         ZOUT0(0) = unit->outval = 0.;
 204         unit->m_tempbuf = 0;
 205 }
 206
 207 void SpecPcile_next(SpecPcile *unit, int inNumSamples)
 208 {
 209         FFTAnalyser_GET_BUF
 210
 211         // Used to be MAKE_TEMP_BUF but we can handle it more cleanly in this specific case:
 212         if (!unit->m_tempbuf) {
 213                 unit->m_tempbuf = (float*)RTAlloc(unit->mWorld, numbins * sizeof(float));
 214                 unit->m_numbins = numbins;
 215                 unit->m_halfnyq_over_numbinsp2 = ((float)unit->mWorld->mSampleRate) * 0.5f / (float)(numbins+2);
 216         } else if (numbins != unit->m_numbins) return;
 217
 218         // Percentile value as a fraction. eg: 0.5 == 50-percentile (median).
 219         float fraction = ZIN0(1);
 220         bool interpolate = unit->m_interpolate;
 221
 222         // The magnitudes in *p will be converted to cumulative sum values and stored in *q temporarily
 223         SCComplexBuf *p = ToComplexApx(buf);
 224         float *q = (float*)unit->m_tempbuf;
 225
 226         float cumul = sc_abs(p->dc);
 227
 228         for (int i=0; i<numbins; ++i) {
 229                 float real = p->bin[i].real;
 230                 float imag = p->bin[i].imag;
 231                 cumul += std::sqrt(real*real + imag*imag);
 232
 233                 // A convenient place to store the mag values...
 234                 q[i] = cumul;
 235         }
 236
 237         cumul += sc_abs(p->nyq);
 238
 239         float target = cumul * fraction; // The target cumul value, stored somewhere in q
 240
 241         float bestposition = 0; // May be linear-interpolated between bins, but not implemented yet
 242                    // NB If nothing beats the target (e.g. if fraction is -1) zero Hz is returned
 243         float binpos;
 244         for(int i=0; i<numbins; ++i) {
 245                 //Print("Testing %g, at position %i", q->bin[i].real, i);
 246                 if(!(q[i] < target)){ // this is a ">=" comparison, done more efficiently as "!(<)"
 247                         if(interpolate && i!=0) {
 248                                 binpos = ((float)i) + 1.f - (q[i] - target) / (q[i] - q[i-1]);
 249                         } else {
 250                                 binpos = ((float)i) + 1.f;
 251                         }
 252                         bestposition = binpos * unit->m_halfnyq_over_numbinsp2;
 253                         //Print("Target %g beaten by %g (at position %i), equating to freq %g\n",
 254                         //                              target, p->bin[i].real, i, bestposition);
 255                         break;
 256                 }
 257         }
 258
 259         // Store the val for output in future calls
 260         unit->outval = bestposition;
 261
 262         ZOUT0(0) = unit->outval;
 263 }
 264
 265
 266 void SpecPcile_Dtor(SpecPcile *unit)
 267 {
 268         if(unit->m_tempbuf) RTFree(unit->mWorld, unit->m_tempbuf);
 269 }
 270
 271 ////////////////////////////////////////////////////////////////////////////////////////////////////////
 272
 273 void SpecCentroid_Ctor(SpecCentroid *unit)
 274 {
 275         SETCALC(SpecCentroid_next);
 276         ZOUT0(0) = unit->outval = 0.;
 277
 278         unit->m_bintofreq = 0.f;
 279 }
 280
 281 void SpecCentroid_next(SpecCentroid *unit, int inNumSamples)
 282 {
 283         FFTAnalyser_GET_BUF
 284
 285         SCPolarBuf *p = ToPolarApx(buf);
 286
 287         GET_BINTOFREQ
 288
 289
 290         double num   = sc_abs(p->nyq) * (numbins+1);
 291         double denom = sc_abs(p->nyq);
 292
 293         for (int i=0; i<numbins; ++i) {
 294                 num   += sc_abs(p->bin[i].mag) * (i+1);
 295                 denom += sc_abs(p->bin[i].mag);
 296         }
 297
 298         ZOUT0(0) = unit->outval = denom == 0.0 ? 0.f : (float) (bintofreq * num/denom);
 299 }
 300
 301 //////////////////////////////////////////////////////////////////////////////////////////////////
 302 /*
 303 void load(InterfaceTable *inTable)
 304 {
 305         ft= inTable;
 306
 307         //(*ft->fDefineUnit)("SpecFlatness", sizeof(FFTAnalyser_Unit), (UnitCtorFunc)&SpecFlatness_Ctor, 0, 0);
 308         //(*ft->fDefineUnit)("SpecPcile", sizeof(SpecPcile_Unit), (UnitCtorFunc)&SpecPcile_Ctor, (UnitDtorFunc)&SpecPcile_Dtor, 0);
 309
 310         DefineSimpleUnit(SpecFlatness);
 311         DefineDtorUnit(SpecPcile);
 312         DefineSimpleUnit(SpecCentroid);
 313 }
 314 */