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 SpecFlatness_Ctor(SpecFlatness *unit);
 106         void SpecFlatness_next(SpecFlatness *unit, int inNumSamples);
 107
 108         void SpecPcile_Ctor(SpecPcile *unit);
 109         void SpecPcile_next(SpecPcile *unit, int inNumSamples);
 110         void SpecPcile_Dtor(SpecPcile *unit);
 111
 112         void SpecCentroid_Ctor(SpecCentroid *unit);
 113         void SpecCentroid_next(SpecCentroid *unit, int inNumSamples);
 114
 115 }
 116 */
 117 /*
 118 SCPolarBuf* ToPolarApx(SndBuf *buf)
 119 {
 120         if (buf->coord == coord_Complex) {
 121                 SCComplexBuf* p = (SCComplexBuf*)buf->data;
 122                 int numbins = buf->samples - 2 >> 1;
 123                 for (int i=0; i<numbins; ++i) {
 124                         p->bin[i].ToPolarApxInPlace();
 125                 }
 126                 buf->coord = coord_Polar;
 127         }
 128
 129         return (SCPolarBuf*)buf->data;
 130 }
 131
 132 SCComplexBuf* ToComplexApx(SndBuf *buf)
 133 {
 134         if (buf->coord == coord_Polar) {
 135                 SCPolarBuf* p = (SCPolarBuf*)buf->data;
 136                 int numbins = buf->samples - 2 >> 1;
 137                 for (int i=0; i<numbins; ++i) {
 138                         p->bin[i].ToComplexApxInPlace();
 139                 }
 140                 buf->coord = coord_Complex;
 141         }
 142         return (SCComplexBuf*)buf->data;
 143 }
 144
 145 InterfaceTable *ft;
 146
 147 void init_SCComplex(InterfaceTable *inTable);
 148 */
 149 //////////////////////////////////////////////////////////////////////////////////////////////////
 150
 151 void SpecFlatness_Ctor(SpecFlatness *unit)
 152 {
 153         SETCALC(SpecFlatness_next);
 154         ZOUT0(0) = unit->outval = 0.;
 155         unit->m_oneovern = 0.;
 156 }
 157
 158 void SpecFlatness_next(SpecFlatness *unit, int inNumSamples)
 159 {
 160         FFTAnalyser_GET_BUF
 161         if(unit->m_oneovern == 0.)
 162                 unit->m_oneovern = 1./(numbins + 2);
 163
 164         SCComplexBuf *p = ToComplexApx(buf);
 165
 166         // Spectral Flatness Measure is geometric mean divided by arithmetic mean.
 167         //
 168         // In order to calculate geom mean without hitting the precision limit,
 169         //  we use the trick of converting to log, taking the average, then converting back from log.
 170         double geommean = std::log(sc_abs(p->dc)) + std::log(sc_abs(p->nyq));
 171         double mean     = sc_abs(p->dc)      + sc_abs(p->nyq);
 172
 173         for (int i=0; i<numbins; ++i) {
 174                 float rabs = (p->bin[i].real);
 175                 float iabs = (p->bin[i].imag);
 176                 float amp = std::sqrt((rabs*rabs) + (iabs*iabs));
 177                 if(amp != 0.f) { // zeroes lead to NaNs
 178                         geommean += std::log(amp);
 179                         mean += amp;
 180                 }
 181         }
 182
 183         double oneovern = unit->m_oneovern;
 184         geommean = exp(geommean * oneovern); // Average and then convert back to linear
 185         mean *= oneovern;
 186
 187         // Store the val for output in future calls
 188         unit->outval = geommean / mean;
 189
 190         ZOUT0(0) = unit->outval;
 191 }
 192
 193 ////////////////////////////////////////////////////////////////////////////////////
 194
 195 void SpecPcile_Ctor(SpecPcile *unit)
 196 {
 197         SETCALC(SpecPcile_next);
 198
 199         unit->m_interpolate = ZIN0(2) > 0.f;
 200
 201         ZOUT0(0) = unit->outval = 0.;
 202         unit->m_tempbuf = 0;
 203 }
 204
 205 void SpecPcile_next(SpecPcile *unit, int inNumSamples)
 206 {
 207         FFTAnalyser_GET_BUF
 208
 209         // Used to be MAKE_TEMP_BUF but we can handle it more cleanly in this specific case:
 210         if (!unit->m_tempbuf) {
 211                 unit->m_tempbuf = (float*)RTAlloc(unit->mWorld, numbins * sizeof(float));
 212                 unit->m_numbins = numbins;
 213                 unit->m_halfnyq_over_numbinsp2 = ((float)unit->mWorld->mSampleRate) * 0.5f / (float)(numbins+2);
 214         } else if (numbins != unit->m_numbins) return;
 215
 216         // Percentile value as a fraction. eg: 0.5 == 50-percentile (median).
 217         float fraction = ZIN0(1);
 218         bool interpolate = unit->m_interpolate;
 219
 220         // The magnitudes in *p will be converted to cumulative sum values and stored in *q temporarily
 221         SCComplexBuf *p = ToComplexApx(buf);
 222         float *q = (float*)unit->m_tempbuf;
 223
 224         float cumul = sc_abs(p->dc);
 225
 226         for (int i=0; i<numbins; ++i) {
 227                 float real = p->bin[i].real;
 228                 float imag = p->bin[i].imag;
 229                 cumul += std::sqrt(real*real + imag*imag);
 230
 231                 // A convenient place to store the mag values...
 232                 q[i] = cumul;
 233         }
 234
 235         cumul += sc_abs(p->nyq);
 236
 237         float target = cumul * fraction; // The target cumul value, stored somewhere in q
 238
 239         float bestposition = 0; // May be linear-interpolated between bins, but not implemented yet
 240                    // NB If nothing beats the target (e.g. if fraction is -1) zero Hz is returned
 241         float binpos;
 242         for(int i=0; i<numbins; ++i) {
 243                 //Print("Testing %g, at position %i", q->bin[i].real, i);
 244                 if(!(q[i] < target)){ // this is a ">=" comparison, done more efficiently as "!(<)"
 245                         if(interpolate && i!=0) {
 246                                 binpos = ((float)i) + 1.f - (q[i] - target) / (q[i] - q[i-1]);
 247                         } else {
 248                                 binpos = ((float)i) + 1.f;
 249                         }
 250                         bestposition = binpos * unit->m_halfnyq_over_numbinsp2;
 251                         //Print("Target %g beaten by %g (at position %i), equating to freq %g\n",
 252                         //                              target, p->bin[i].real, i, bestposition);
 253                         break;
 254                 }
 255         }
 256
 257         // Store the val for output in future calls
 258         unit->outval = bestposition;
 259
 260         ZOUT0(0) = unit->outval;
 261 }
 262
 263
 264 void SpecPcile_Dtor(SpecPcile *unit)
 265 {
 266         if(unit->m_tempbuf) RTFree(unit->mWorld, unit->m_tempbuf);
 267 }
 268
 269 ////////////////////////////////////////////////////////////////////////////////////////////////////////
 270
 271 void SpecCentroid_Ctor(SpecCentroid *unit)
 272 {
 273         SETCALC(SpecCentroid_next);
 274         ZOUT0(0) = unit->outval = 0.;
 275
 276         unit->m_bintofreq = 0.f;
 277 }
 278
 279 void SpecCentroid_next(SpecCentroid *unit, int inNumSamples)
 280 {
 281         FFTAnalyser_GET_BUF
 282
 283         SCPolarBuf *p = ToPolarApx(buf);
 284
 285         GET_BINTOFREQ
 286
 287
 288         double num   = sc_abs(p->nyq) * (numbins+1);
 289         double denom = sc_abs(p->nyq);
 290
 291         for (int i=0; i<numbins; ++i) {
 292                 num   += sc_abs(p->bin[i].mag) * (i+1);
 293                 denom += sc_abs(p->bin[i].mag);
 294         }
 295
 296         ZOUT0(0) = unit->outval = denom == 0.0 ? 0.f : (float) (bintofreq * num/denom);
 297 }
 298
 299 //////////////////////////////////////////////////////////////////////////////////////////////////
 300 /*
 301 void load(InterfaceTable *inTable)
 302 {
 303         ft= inTable;
 304
 305         //(*ft->fDefineUnit)("SpecFlatness", sizeof(FFTAnalyser_Unit), (UnitCtorFunc)&SpecFlatness_Ctor, 0, 0);
 306         //(*ft->fDefineUnit)("SpecPcile", sizeof(SpecPcile_Unit), (UnitCtorFunc)&SpecPcile_Ctor, (UnitDtorFunc)&SpecPcile_Dtor, 0);
 307
 308         DefineSimpleUnit(SpecFlatness);
 309         DefineDtorUnit(SpecPcile);
 310         DefineSimpleUnit(SpecCentroid);
 311 }
 312 */