zynaddsubfx

ZynAddSubFX open source synthesizer
Log | Files | Refs | Submodules | LICENSE

OscilGen.cpp (60293B)


      1 /*
      2   ZynAddSubFX - a software synthesizer
      3 
      4   OscilGen.cpp - Waveform generator for ADnote
      5   Copyright (C) 2021 Hans Petter Selasky
      6   Copyright (C) 2002-2005 Nasca Octavian Paul
      7   Author: Nasca Octavian Paul
      8 
      9   This program is free software; you can redistribute it and/or
     10   modify it under the terms of the GNU General Public License
     11   as published by the Free Software Foundation; either version 2
     12   of the License, or (at your option) any later version.
     13 */
     14 
     15 #include "OscilGen.h"
     16 #include "../DSP/FFTwrapper.h"
     17 #include "../Synth/Resonance.h"
     18 #include "../Misc/WaveShapeSmps.h"
     19 
     20 #include <cassert>
     21 #include <cstdlib>
     22 #include <cmath>
     23 #include <cstdio>
     24 #include <cstddef>
     25 #include <complex>
     26 
     27 #include <unistd.h>
     28 
     29 #include <rtosc/ports.h>
     30 #include <rtosc/port-sugar.h>
     31 
     32 namespace zyn {
     33 
     34 
     35 #define rObject OscilGen
     36 const rtosc::Ports OscilGen::non_realtime_ports = {
     37     rSelf(OscilGen),
     38     rPaste,
     39 #undef rDefaultProps
     40 #define rDefaultProps rProp(non-realtime)
     41     //TODO ensure min/max
     42     rOption(Phmagtype, rShort("scale"),
     43             rOptions(linear,dB scale (-40),
     44                      dB scale (-60), dB scale (-80),
     45                      dB scale (-100)),
     46             rDefault(linear),
     47             "Type of magnitude for harmonics"),
     48     rOption(Pcurrentbasefunc, rShort("base"),
     49             rOptions(sine, triangle, pulse, saw, power, gauss,
     50                 diode, abssine, pulsesine, stretchsine,
     51                 chirp, absstretchsine, chebyshev, sqr,
     52                 spike, circle, powersinus), rOpt(127,use-as-base waveform),
     53             rDefault(sine),
     54             "Base Waveform for harmonics"),
     55     rParamZyn(Pbasefuncpar, rShort("shape"), rDefault(64),
     56             "Morph between possible base function shapes "
     57             "(e.g. rising sawtooth vs a falling sawtooth)"),
     58     rOption(Pbasefuncmodulation, rShort("mod"),
     59             rOptions(None, Rev, Sine, Power, Chop), rDefault(None),
     60             "Modulation applied to Base function spectra"),
     61     rParamZyn(Pbasefuncmodulationpar1, rShort("p1"), rDefault(64),
     62             "Base function modulation parameter"),
     63     rParamZyn(Pbasefuncmodulationpar2, rShort("p2"), rDefault(64),
     64             "Base function modulation parameter"),
     65     rParamZyn(Pbasefuncmodulationpar3, rShort("p3"), rDefault(32),
     66             "Base function modulation parameter"),
     67     rParamZyn(Pwaveshaping, rShort("amount"), rDefault(64),
     68             "Degree Of Waveshaping"),
     69     rOption(Pwaveshapingfunction, rShort("distort"), rDefault(Undistorted),
     70             rOptions(Undistorted,
     71                 Arctangent, Asymmetric, Pow, Sine, Quantisize,
     72                 Zigzag, Limiter, Upper Limiter, Lower Limiter,
     73                 Inverse Limiter, Clip, Asym2, Pow2, sigmoid, Tanh, Cubic, Square),
     74             "Shape of distortion to be applied"),
     75     rOption(Pfiltertype, rShort("filter"), rOptions(No Filter,
     76             lp, hp1, hp1b, bp1, bs1, lp2, hp2, bp2, bs2,
     77             cos, sin, low_shelf, s, lpsk), rDefaultId(No Filter), "Harmonic Filter"),
     78     rParamZyn(Pfilterpar1, rShort("p1"), rDefault(64), "Filter parameter"),
     79     rParamZyn(Pfilterpar2, rShort("p2"), rDefault(64), "Filter parameter"),
     80     rToggle(Pfilterbeforews, rShort("pre/post"), rDefault(false),
     81             "Filter before waveshaping spectra;"
     82             "When enabled oscilfilter(freqs); then waveshape(freqs);, "
     83             "otherwise waveshape(freqs); then oscilfilter(freqs);"),
     84     rOption(Psatype, rShort("spec. adj."), rOptions(None, Pow, ThrsD, ThrsU),
     85             rDefault(None), "Spectral Adjustment Type"),
     86     rParamZyn(Psapar, rShort("p1"), rDefault(64),
     87               "Spectral Adjustment Parameter"),
     88     rParamI(Pharmonicshift, rLinear(-64,64), rShort("shift"), rDefault(0),
     89             "Amount of shift on harmonics"),
     90     rToggle(Pharmonicshiftfirst, rShort("pre/post"), rDefault(false),
     91             "If harmonics are shifted before waveshaping/filtering"),
     92     rOption(Pmodulation, rShort("FM"), rOptions(None, Rev, Sine, Power),
     93             rDefault(None), "Frequency Modulation To Combined Spectra"),
     94     rParamZyn(Pmodulationpar1, rShort("p1"), rDefault(64),
     95               "modulation parameter"),
     96     rParamZyn(Pmodulationpar2, rShort("p2"), rDefault(64),
     97               "modulation parameter"),
     98     rParamZyn(Pmodulationpar3, rShort("p3"), rDefault(32),
     99               "modulation parameter"),
    100     rToggle(ADvsPAD, rShort("If it is used by PADSynth"),
    101             "If it is used by PADSynth (and not ADSynth)"),
    102 #undef rDefaultProps
    103 #define rDefaultProps
    104 
    105 
    106     //TODO update to rArray and test
    107     {"phase#128::c:i", rProp(parameter) rLinear(0,127) rProp(non-realtime)
    108         rDefault([64 ...])
    109         rDoc("Sets harmonic phase"),
    110         NULL, [](const char *m, rtosc::RtData &d) {
    111             const char *mm = m;
    112             while(*mm && !isdigit(*mm)) ++mm;
    113             unsigned char &phase = ((OscilGen*)d.obj)->Phphase[atoi(mm)];
    114             if(!rtosc_narguments(m))
    115                 d.reply(d.loc, "i", phase);
    116             else {
    117                 phase = rtosc_argument(m,0).i;
    118                 //XXX hack hack
    119                 char  repath[128];
    120                 strcpy(repath, d.loc);
    121                 char *edit   = strrchr(repath, '/')+1;
    122                 strcpy(edit, "prepare");
    123                 OscilGen &o = *((OscilGen*)d.obj);
    124                 FFTfreqBuffer freqs = o.fft->allocFreqBuf();
    125                 OscilGenBuffers& bfrs = o.myBuffers();
    126                 o.prepare(bfrs, freqs);
    127                 // fprintf(stderr, "sending '%p' of fft data\n", data);
    128                 d.chain(repath, "b", sizeof(fft_t*), &freqs.data);
    129                 bfrs.pendingfreqs = freqs.data;
    130                 d.broadcast(d.loc, "i", phase);
    131             }
    132         }},
    133     //TODO update to rArray and test
    134     {"magnitude#128::c:i", rProp(parameter) rLinear(0,127) rProp(non-realtime)
    135         rProp(no port checker) // buggy port
    136         rDefault([127 64 64 ...]) rDoc("Sets harmonic magnitude"),
    137         NULL, [](const char *m, rtosc::RtData &d) {
    138             //printf("I'm at '%s'\n", d.loc);
    139             const char *mm = m;
    140             while(*mm && !isdigit(*mm)) ++mm;
    141             unsigned char &mag = ((OscilGen*)d.obj)->Phmag[atoi(mm)];
    142             if(!rtosc_narguments(m))
    143                 d.reply(d.loc, "i", mag);
    144             else {
    145                 mag = rtosc_argument(m,0).i;
    146                 //printf("setting magnitude\n\n");
    147                 //XXX hack hack
    148                 char  repath[128];
    149                 strcpy(repath, d.loc);
    150                 char *edit   = strrchr(repath, '/')+1;
    151                 strcpy(edit, "prepare");
    152                 OscilGen &o = *((OscilGen*)d.obj);
    153                 FFTfreqBuffer freqs = o.fft->allocFreqBuf();
    154                 OscilGenBuffers& bfrs = o.myBuffers();
    155                 o.prepare(bfrs, freqs);
    156                 // fprintf(stderr, "sending '%p' of fft data\n", data);
    157                 d.chain(repath, "b", sizeof(fft_t*), &freqs.data);
    158                 bfrs.pendingfreqs = freqs.data;
    159                 d.broadcast(d.loc, "i", mag);
    160             }
    161         }},
    162     {"basefuncFFTfreqs::b", rProp(parameter) rProp(non-realtime)
    163         rProp(no port checker) // buggy port
    164         rBlobType(f) rDepends(magnitude, phase) rDefault([0.f 0.f ...])
    165         rDoc("FFT freqs of base function"),
    166         NULL, rBOIL_BEGIN
    167             // no. of floats to send/recv is (oscilsize/2-1) * 2,
    168             // because the first FFT bin is not used, and each bin has real+imag
    169             const int32_t bufsize = obj->synth.oscilsize-2;
    170             if(!rtosc_narguments(msg))
    171             {
    172                 float* tmpbuf = new float[bufsize];
    173                 for (int i = 0; i < bufsize/2; ++i)
    174                 {
    175                     tmpbuf[2*i]   = obj->myBuffers().basefuncFFTfreqs[i+1].real();
    176                     tmpbuf[2*i+1] = obj->myBuffers().basefuncFFTfreqs[i+1].imag();
    177                 }
    178                 data.reply(loc, "b", bufsize*sizeof(float), tmpbuf);
    179                 delete[] tmpbuf;
    180             } else {
    181                 rtosc_blob_t blob = rtosc_argument(msg, 0).b;
    182                 float* buf = (float*) blob.data;
    183                 int len = blob.len/sizeof(float);
    184                 int max = std::min(len, bufsize);
    185                 for (int i = 0; i < max/2; ++i)
    186                 {
    187                     obj->myBuffers().basefuncFFTfreqs[i+1] =
    188                         fft_t(buf[2*i], buf[2*i+1]);
    189                 }
    190                 // TODO: Simplify this code section when merging with WT branch
    191                 char  repath[128];
    192                 strcpy(repath, data.loc);
    193                 char *edit   = strrchr(repath, '/')+1;
    194                 strcpy(edit, "prepare");
    195                 FFTfreqBuffer freqs = obj->fft->allocFreqBuf();
    196                 OscilGenBuffers& bfrs = obj->myBuffers();
    197                 obj->prepare(bfrs, freqs);
    198                 data.chain(repath, "b", sizeof(fft_t*), &freqs.data);
    199                 bfrs.pendingfreqs = freqs.data;
    200                 data.broadcast(loc, "b", bufsize*sizeof(float), buf);
    201             }
    202         rBOIL_END
    203         },
    204     {"base-spectrum:", rProp(non-realtime) rDoc("Returns spectrum of base waveshape"),
    205         NULL, [](const char *, rtosc::RtData &d) {
    206             OscilGen &o = *((OscilGen*)d.obj);
    207             const unsigned n = o.synth.oscilsize / 2;
    208             float *spc = new float[n];
    209             memset(spc, 0, 4*n);
    210             ((OscilGen*)d.obj)->getspectrum(n,spc,1);
    211             d.reply(d.loc, "b", n*sizeof(float), spc);
    212             delete[] spc;
    213         }},
    214     {"base-waveform:", rProp(non-realtime) rDoc("Returns base waveshape points"),
    215         NULL, [](const char *, rtosc::RtData &d) {
    216             OscilGen &o = *((OscilGen*)d.obj);
    217             FFTsampleBuffer smps = o.fft->allocSampleBuf();
    218             const unsigned n = smps.allocSize() * sizeof(float);
    219             memset(smps.data, 0, n);
    220             ((OscilGen*)d.obj)->getcurrentbasefunction(smps);
    221             d.reply(d.loc, "b", n, smps.data);
    222             delete[] smps.data;
    223         }},
    224     {"prepare:", rProp(non-realtime) rDoc("Performs setup operation to oscillator"),
    225         NULL, [](const char *, rtosc::RtData &d) {
    226             //fprintf(stderr, "prepare: got a message from '%s'\n", m);
    227             OscilGen &o = *(OscilGen*)d.obj;
    228             FFTfreqBuffer freqs = o.fft->allocFreqBuf();
    229             OscilGenBuffers& bfrs = o.myBuffers();
    230             o.prepare(bfrs, freqs);
    231             // fprintf(stderr, "sending '%p' of fft data\n", data);
    232             d.chain(d.loc, "b", sizeof(fft_t*), &freqs.data);
    233             bfrs.pendingfreqs = freqs.data;
    234         }},
    235     {"convert2sine:", rProp(non-realtime) rDoc("Translates waveform into FS"),
    236         NULL, [](const char *, rtosc::RtData &d) {
    237             ((OscilGen*)d.obj)->convert2sine();
    238             //XXX hack hack
    239             char  repath[128];
    240             strcpy(repath, d.loc);
    241             char *edit   = strrchr(repath, '/')+1;
    242             *edit = 0;
    243             d.broadcast("/damage", "s", repath);
    244         }},
    245     {"use-as-base:", rProp(non-realtime) rDoc("Translates current waveform into base"),
    246         NULL, [](const char *, rtosc::RtData &d) {
    247             ((OscilGen*)d.obj)->useasbase();
    248             //XXX hack hack
    249             char  repath[128];
    250             strcpy(repath, d.loc);
    251             char *edit   = strrchr(repath, '/')+1;
    252             *edit = 0;
    253             d.broadcast("/damage", "s", repath);
    254         }}};
    255 
    256 const rtosc::Ports OscilGen::realtime_ports{
    257     rSelf(OscilGen),
    258     rPresetType,
    259     rParamZyn(Prand, rLinear(-64, 63), rShort("phase rnd"),
    260             rDefaultDepends(ADvsPAD), rPreset(true, 127), rPreset(false, 64),
    261             "Oscillator Phase Randomness: smaller than 0 is \""
    262             "group\", larger than 0 is for each harmonic"),
    263     rParamZyn(Pamprandpower, rShort("variance"), rDefault(64),
    264             "Variance of harmonic randomness"),
    265     rOption(Pamprandtype, rShort("distribution"), rOptions(None, Pow, Sin),
    266             rDefault(None),
    267             "Harmonic random distribution to select from"),
    268     rOption(Padaptiveharmonics, rShort("adapt")
    269             rOptions(OFF, ON, Square, 2xSub, 2xAdd, 3xSub, 3xAdd, 4xSub, 4xAdd),
    270             rDefault(OFF),
    271             "Adaptive Harmonics Mode"),
    272     rParamI(Padaptiveharmonicsbasefreq, rShort("c. freq"), rLinear(0,255),
    273             rDefault(128), "Base frequency of adaptive harmonic (30..3000Hz)"),
    274     rParamI(Padaptiveharmonicspower, rShort("amount"), rLinear(0,200),
    275             rDefault(100), "Adaptive Harmonic Strength"),
    276     rParamI(Padaptiveharmonicspar, rShort("power"), rLinear(0,100),
    277             rDefault(50), "Adaptive Harmonics Postprocessing Power"),
    278     {"waveform:", rDoc("Returns waveform points"),
    279         NULL, [](const char *, rtosc::RtData &d) {
    280             OscilGen &o = *((OscilGen*)d.obj);
    281             const unsigned n = o.synth.oscilsize;
    282             float *smps = new float[n]; // XXXRT
    283             memset(smps, 0, 4*n);
    284             //printf("%d\n", o->needPrepare());
    285             o.get(o.myBuffers(),smps,-1.0);
    286             //printf("wave: %f %f %f %f\n", smps[0], smps[1], smps[2], smps[3]);
    287             d.reply(d.loc, "b", n*sizeof(float), smps);
    288             delete[] smps;
    289         }},
    290     {"spectrum:", rDoc("Returns spectrum of waveform"),
    291         NULL, [](const char *, rtosc::RtData &d) {
    292             OscilGen &o = *((OscilGen*)d.obj);
    293             const unsigned n = o.synth.oscilsize / 2;
    294             float *spc = new float[n]; // XXXRT
    295             memset(spc, 0, 4*n);
    296             ((OscilGen*)d.obj)->getspectrum(n,spc,0);
    297             d.reply(d.loc, "b", n*sizeof(float), spc);
    298             delete[] spc;
    299         }},
    300     {"prepare:b", rProp(internal) rProp(realtime) rProp(pointer) rDoc("Sets prepared fft data"),
    301         NULL, [](const char *m, rtosc::RtData &d) {
    302             // fprintf(stderr, "prepare:b got a message from '%s'\n", m);
    303             OscilGen &o = *(OscilGen*)d.obj;
    304             OscilGenBuffers& bfrs = o.myBuffers();
    305             assert(rtosc_argument(m,0).b.len == sizeof(void*));
    306             d.reply("/free", "sb", "fft_t", sizeof(void*), &bfrs.oscilFFTfreqs.data);
    307             assert(bfrs.oscilFFTfreqs.data !=*(fft_t**)rtosc_argument(m,0).b.data);
    308             bfrs.oscilFFTfreqs.data = *(fft_t**)rtosc_argument(m,0).b.data;
    309         }},
    310 
    311 };
    312 #undef rDefaultProps
    313 
    314 const rtosc::MergePorts OscilGen::ports{
    315     &OscilGen::realtime_ports,
    316     &OscilGen::non_realtime_ports
    317 };
    318 
    319 #ifndef M_PI_2
    320 # define M_PI_2		1.57079632679489661923	/* pi/2 */
    321 #endif
    322 
    323 
    324 //operations on FFTfreqs
    325 inline void clearAll(fft_t *freqs, int oscilsize)
    326 {
    327     fft_t zero = 0;
    328     std::fill_n(freqs, oscilsize / 2, zero);
    329 }
    330 
    331 inline void clearDC(fft_t *freqs)
    332 {
    333     freqs[0] = fft_t(0.0f, 0.0f);
    334 }
    335 
    336 //return magnitude squared
    337 inline float normal(const fft_t *freqs, off_t x)
    338 {
    339     return (float)norm(freqs[x]);
    340 }
    341 
    342 //return magnitude
    343 inline float abs(const fft_t *freqs, off_t x)
    344 {
    345     return (float)abs(freqs[x]);
    346 }
    347 
    348 //return angle aka phase from a sine (not cosine wave)
    349 inline float arg(const fft_t *freqs, off_t x)
    350 {
    351     const fft_t tmp(freqs[x].imag(), freqs[x].real());
    352     return (float)arg(tmp);
    353 }
    354 
    355 /**
    356  * Take frequency spectrum and ensure values are normalized based upon
    357  * magnitude to 0<=x<=1
    358  */
    359 void normalize(fft_t *freqs, int oscilsize)
    360 {
    361     float normMax = 0.0f;
    362     for(int i = 0; i < oscilsize / 2; ++i) {
    363         //magnitude squared
    364         const float norm = normal(freqs, i);
    365         if(normMax < norm)
    366             normMax = norm;
    367     }
    368 
    369     const float max = sqrtf(normMax);
    370     if(max < 1e-8) //data is all ~zero, do not amplify noise
    371         return;
    372 
    373     for(int i = 0; i < oscilsize / 2; ++i)
    374         freqs[i] /= max;
    375 }
    376 
    377 //Full RMS normalize
    378 void rmsNormalize(fft_t *freqs, int oscilsize)
    379 {
    380     float sum = 0.0f;
    381     for(int i = 1; i < oscilsize / 2; ++i)
    382         sum += normal(freqs, i);
    383 
    384     if(sum < 0.000001f)
    385         return;  //data is all ~zero, do not amplify noise
    386 
    387     const float gain = 1.0f / sqrtf(sum);
    388 
    389     for(int i = 1; i < oscilsize / 2; ++i)
    390         freqs[i] *= gain;
    391 }
    392 
    393 #define DIFF(par) (bfrs.old ## par != P ## par)
    394 
    395 FFTfreqBuffer ctorAllocFreqs(FFTwrapper* fft, int oscilsize) {
    396     return fft ? fft->allocFreqBuf() : FFTwrapper::riskAllocFreqBufWithSize(oscilsize);
    397 }
    398 
    399 FFTsampleBuffer ctorAllocSamples(FFTwrapper* fft, int oscilsize) {
    400     return fft ? fft->allocSampleBuf() : FFTwrapper::riskAllocSampleBufWithSize(oscilsize);
    401 }
    402 
    403 OscilGenBuffers::OscilGenBuffers(zyn::OscilGenBuffersCreator c) :
    404     oscilsize(c.oscilsize),
    405     // fft_ can be nullptr in case of pasting
    406     oscilFFTfreqs(ctorAllocFreqs(c.fft, c.oscilsize)),
    407     pendingfreqs(oscilFFTfreqs.data),
    408     tmpsmps(ctorAllocSamples(c.fft, c.oscilsize)),
    409     outoscilFFTfreqs(ctorAllocFreqs(c.fft, c.oscilsize)),
    410     cachedbasefunc(ctorAllocSamples(c.fft, c.oscilsize)),
    411     cachedbasevalid(false),
    412     basefuncFFTfreqs(ctorAllocFreqs(c.fft, c.oscilsize)),
    413     scratchFreqs(ctorAllocFreqs(c.fft, c.oscilsize))
    414 {
    415     defaults();
    416 }
    417 
    418 OscilGenBuffers::~OscilGenBuffers()
    419 {
    420     delete[] tmpsmps.data;
    421     delete[] outoscilFFTfreqs.data;
    422     delete[] basefuncFFTfreqs.data;
    423     delete[] oscilFFTfreqs.data;
    424     delete[] cachedbasefunc.data;
    425     delete[] scratchFreqs.data;
    426 }
    427 
    428 zyn::OscilGenBuffersCreator OscilGen::createOscilGenBuffers() const
    429 {
    430     return OscilGenBuffersCreator(fft, synth.oscilsize);
    431 }
    432 
    433 void OscilGenBuffers::defaults()
    434 {
    435     oldbasefunc = 0;
    436     oldbasepar  = 64;
    437     oldhmagtype = 0;
    438     oldwaveshapingfunction = 0;
    439     oldwaveshaping = 64;
    440     oldbasefuncmodulation     = 0;
    441     oldharmonicshift          = 0;
    442     oldbasefuncmodulationpar1 = 0;
    443     oldbasefuncmodulationpar2 = 0;
    444     oldbasefuncmodulationpar3 = 0;
    445     oldmodulation     = 0;
    446     oldmodulationpar1 = 0;
    447     oldmodulationpar2 = 0;
    448     oldmodulationpar3 = 0;
    449 
    450     for(int i = 0; i < MAX_AD_HARMONICS; ++i) {
    451         hmag[i]    = 0.0f;
    452         hphase[i]  = 0.0f;
    453     }
    454 
    455     clearAll(oscilFFTfreqs.data, oscilsize);
    456     clearAll(basefuncFFTfreqs.data, oscilsize);
    457     oscilprepared = 0;
    458     oldfilterpars = 0;
    459     oldsapars     = 0;
    460 }
    461 
    462 OscilGen::OscilGen(const SYNTH_T &synth_, FFTwrapper *fft_, const Resonance *res_)
    463     :Presets(),
    464       m_myBuffers(OscilGenBuffersCreator(fft_, synth_.oscilsize)),
    465       fft(fft_),
    466       res(res_),
    467       synth(synth_)
    468 {
    469     if(fft_) {
    470         // FFTwrapper should operate exactly on all "oscilsize" bytes
    471         assert(fft_->fftsize() == synth_.oscilsize);
    472     } else {
    473         // this is possible if *this is a temporary paste object
    474     }
    475 
    476     setpresettype("Poscilgen");
    477 
    478     randseed = 1; // TODO: take care of random?
    479     ADvsPAD  = false;
    480 
    481     defaults();
    482 }
    483 
    484 void OscilGen::defaults()
    485 {
    486     for(int i = 0; i < MAX_AD_HARMONICS; ++i) {
    487         Phmag[i]   = 64;
    488         Phphase[i] = 64;
    489     }
    490     Phmag[0]  = 127;
    491     Phmagtype = 0;
    492     if(ADvsPAD)
    493         Prand = 127;       //max phase randomness (useful if the oscil will be imported to a ADsynth from a PADsynth
    494     else
    495         Prand = 64;  //no randomness
    496 
    497     Pcurrentbasefunc = 0;
    498     Pbasefuncpar     = 64;
    499 
    500     Pbasefuncmodulation     = 0;
    501     Pbasefuncmodulationpar1 = 64;
    502     Pbasefuncmodulationpar2 = 64;
    503     Pbasefuncmodulationpar3 = 32;
    504 
    505     Pmodulation     = 0;
    506     Pmodulationpar1 = 64;
    507     Pmodulationpar2 = 64;
    508     Pmodulationpar3 = 32;
    509 
    510     Pwaveshapingfunction = 0;
    511     Pwaveshaping    = 64;
    512     Pfiltertype     = 0;
    513     Pfilterpar1     = 64;
    514     Pfilterpar2     = 64;
    515     Pfilterbeforews = 0;
    516     Psatype = 0;
    517     Psapar  = 64;
    518 
    519     Pamprandpower = 64;
    520     Pamprandtype  = 0;
    521 
    522     Pharmonicshift      = 0;
    523     Pharmonicshiftfirst = 0;
    524 
    525     Padaptiveharmonics         = 0;
    526     Padaptiveharmonicspower    = 100;
    527     Padaptiveharmonicsbasefreq = 128;
    528     Padaptiveharmonicspar      = 50;
    529 
    530     prepare(myBuffers());
    531 }
    532 
    533 void OscilGen::convert2sine()
    534 {
    535     OscilGenBuffers& bfrs = myBuffers();
    536     float  mag[MAX_AD_HARMONICS], phase[MAX_AD_HARMONICS];
    537 
    538     {
    539         FFTwrapper *fft = new FFTwrapper(synth.oscilsize);
    540         FFTsampleBuffer oscil = fft->allocSampleBuf();
    541         get(oscil.data, -1.0f);
    542         fft->smps2freqs_noconst_input(oscil, bfrs.scratchFreqs);
    543         delete (fft);
    544     }
    545 
    546     normalize(bfrs.scratchFreqs.data, synth.oscilsize);
    547 
    548     mag[0]   = 0;
    549     phase[0] = 0;
    550     for(int i = 0; i < MAX_AD_HARMONICS; ++i) {
    551         mag[i]   = abs(bfrs.scratchFreqs.data, i + 1);
    552         phase[i] = arg(bfrs.scratchFreqs.data, i + 1);
    553     }
    554 
    555     defaults();
    556 
    557     for(int i = 0; i < MAX_AD_HARMONICS - 1; ++i) {
    558         float newmag   = mag[i];
    559         float newphase = phase[i];
    560 
    561         Phmag[i] = (int) ((newmag) * 63.0f) + 64;
    562 
    563         Phphase[i] = 64 - (int) (64.0f * newphase / PI);
    564         if(Phphase[i] > 127)
    565             Phphase[i] = 127;
    566 
    567         if(Phmag[i] == 64)
    568             Phphase[i] = 64;
    569     }
    570     prepare();
    571 }
    572 
    573 float OscilGen::userfunc(OscilGenBuffers& bfrs, float x) const
    574 {
    575     if (!fft)
    576         return 0;
    577     if (!bfrs.cachedbasevalid) {
    578         fft->freqs2smps(bfrs.basefuncFFTfreqs, bfrs.cachedbasefunc, bfrs.scratchFreqs);
    579         bfrs.cachedbasevalid = true;
    580     }
    581     return cinterpolate(bfrs.cachedbasefunc.data,
    582                         synth.oscilsize,
    583                         synth.oscilsize * (x + 1) - 1);
    584 }
    585 
    586 /*
    587  * Get the base function
    588  */
    589 void OscilGen::getbasefunction(OscilGenBuffers& bfrs, FFTsampleBuffer smps) const
    590 {
    591     float par = (Pbasefuncpar + 0.5f) / 128.0f;
    592     if(Pbasefuncpar == 64)
    593         par = 0.5f;
    594 
    595     float p1 = Pbasefuncmodulationpar1 / 127.0f,
    596           p2 = Pbasefuncmodulationpar2 / 127.0f,
    597           p3 = Pbasefuncmodulationpar3 / 127.0f;
    598 
    599     switch(Pbasefuncmodulation) {
    600         case 1:
    601             p1 = (powf(2, p1 * 5.0f) - 1.0f) / 10.0f;
    602             p3 = floorf(powf(2, p3 * 5.0f) - 1.0f);
    603             if(p3 < 0.9999f)
    604                 p3 = -1.0f;
    605             break;
    606         case 2:
    607             p1 = (powf(2, p1 * 5.0f) - 1.0f) / 10.0f;
    608             p3 = 1.0f + floorf(powf(2, p3 * 5.0f) - 1.0f);
    609             break;
    610         case 3:
    611             p1 = (powf(2, p1 * 7.0f) - 1.0f) / 10.0f;
    612             p3 = 0.01f + (powf(2, p3 * 16.0f) - 1.0f) / 10.0f;
    613             break;
    614     }
    615 
    616     base_func_t *func = getBaseFunction(Pcurrentbasefunc);
    617 
    618     for(int i = 0; i < synth.oscilsize; ++i) {
    619         float t = i * 1.0f / synth.oscilsize;
    620 
    621         switch(Pbasefuncmodulation) {
    622             case 1: //rev
    623                 t = t * p3 + sinf((t + p2) * 2.0f * PI) * p1;
    624                 break;
    625             case 2: //sine
    626                 t += sinf( (t * p3 + p2) * 2.0f * PI) * p1;
    627                 break;
    628             case 3: //power
    629                 t += powf((1.0f - cosf((t + p2) * 2.0f * PI)) * 0.5f, p3) * p1;
    630                 break;
    631             case 4: //chop
    632                 t = t * (powf(2.0, Pbasefuncmodulationpar1/32.f +
    633                               Pbasefuncmodulationpar2/2048.f)) + p3;
    634         }
    635 
    636         t = t - floorf(t);
    637 
    638         if(func)
    639             smps[i] = func(t, par);
    640         else if (Pcurrentbasefunc == 0)
    641             smps[i] = -sinf(2.0f * PI * i / synth.oscilsize);
    642         else
    643             smps[i] = userfunc(bfrs, t);
    644     }
    645 }
    646 
    647 
    648 /*
    649  * Filter the oscillator
    650  */
    651 void OscilGen::oscilfilter(fft_t *freqs) const
    652 {
    653     if(Pfiltertype == 0)
    654         return;
    655 
    656     const float par    = 1.0f - Pfilterpar1 / 128.0f;
    657     const float par2   = Pfilterpar2 / 127.0f;
    658     filter_func_t *filter = getFilter(Pfiltertype);
    659 
    660     for(int i = 1; i < synth.oscilsize / 2; ++i)
    661         freqs[i] *= filter(i, par, par2);
    662 
    663     normalize(freqs, synth.oscilsize);
    664 }
    665 
    666 
    667 /*
    668  * Change the base function
    669  */
    670 void OscilGen::changebasefunction(OscilGenBuffers& bfrs) const
    671 {
    672     if(Pcurrentbasefunc != 0) {
    673         if(Pcurrentbasefunc == 127 && !Pbasefuncmodulation) {
    674             // this would be a no-op, skip it
    675         }
    676         else {
    677             getbasefunction(bfrs, bfrs.tmpsmps);
    678             if(fft)
    679                 fft->smps2freqs_noconst_input(bfrs.tmpsmps, bfrs.basefuncFFTfreqs);
    680         }
    681         clearDC(bfrs.basefuncFFTfreqs.data);
    682     }
    683     else //in this case bfrs.basefuncFFTfreqs are not used
    684         clearAll(bfrs.basefuncFFTfreqs.data, synth.oscilsize);
    685     bfrs.oscilprepared = 0;
    686     bfrs.oldbasefunc   = Pcurrentbasefunc;
    687     bfrs.oldbasepar    = Pbasefuncpar;
    688     bfrs.oldbasefuncmodulation     = Pbasefuncmodulation;
    689     bfrs.oldbasefuncmodulationpar1 = Pbasefuncmodulationpar1;
    690     bfrs.oldbasefuncmodulationpar2 = Pbasefuncmodulationpar2;
    691     bfrs.oldbasefuncmodulationpar3 = Pbasefuncmodulationpar3;
    692 }
    693 
    694 inline void normalize(float *smps, size_t N)
    695 {
    696     //Find max
    697     float max = 0.0f;
    698     for(size_t i = 0; i < N; ++i)
    699         if(max < fabsf(smps[i]))
    700             max = fabsf(smps[i]);
    701     if(max < 0.00001f)
    702         max = 1.0f;
    703 
    704     //Normalize to +-1
    705     for(size_t i = 0; i < N; ++i)
    706         smps[i] /= max;
    707 }
    708 
    709 /*
    710  * Waveshape
    711  */
    712 void OscilGen::waveshape(OscilGenBuffers& bfrs, FFTfreqBuffer freqs) const
    713 {
    714     bfrs.oldwaveshapingfunction = Pwaveshapingfunction;
    715     bfrs.oldwaveshaping = Pwaveshaping;
    716     if(Pwaveshapingfunction == 0)
    717         return;
    718 
    719     clearDC(freqs.data);
    720     //reduce the amplitude of the freqs near the nyquist
    721     for(int i = 1; i < synth.oscilsize / 8; ++i) {
    722         float gain = i / (synth.oscilsize / 8.0f);
    723         freqs[synth.oscilsize / 2 - i] *= gain;
    724     }
    725     fft->freqs2smps_noconst_input(freqs, bfrs.tmpsmps);
    726 
    727     //Normalize
    728     normalize(bfrs.tmpsmps.data, synth.oscilsize);
    729 
    730     //Do the waveshaping
    731     waveShapeSmps(synth.oscilsize, bfrs.tmpsmps.data, Pwaveshapingfunction, Pwaveshaping);
    732 
    733     fft->smps2freqs_noconst_input(bfrs.tmpsmps, freqs); //perform FFT
    734 }
    735 
    736 
    737 /*
    738  * Do the Frequency Modulation of the Oscil
    739  */
    740 void OscilGen::modulation(OscilGenBuffers& bfrs, FFTfreqBuffer freqs) const
    741 {
    742     bfrs.oldmodulation     = Pmodulation;
    743     bfrs.oldmodulationpar1 = Pmodulationpar1;
    744     bfrs.oldmodulationpar2 = Pmodulationpar2;
    745     bfrs.oldmodulationpar3 = Pmodulationpar3;
    746     if(Pmodulation == 0)
    747         return;
    748 
    749 
    750     float modulationpar1 = Pmodulationpar1 / 127.0f,
    751           modulationpar2 = 0.5f - Pmodulationpar2 / 127.0f,
    752           modulationpar3 = Pmodulationpar3 / 127.0f;
    753 
    754     switch(Pmodulation) {
    755         case 1:
    756             modulationpar1 = (powf(2, modulationpar1 * 7.0f) - 1.0f) / 100.0f;
    757             modulationpar3 = floorf((powf(2, modulationpar3 * 5.0f) - 1.0f));
    758             if(modulationpar3 < 0.9999f)
    759                 modulationpar3 = -1.0f;
    760             break;
    761         case 2:
    762             modulationpar1 = (powf(2, modulationpar1 * 7.0f) - 1.0f) / 100.0f;
    763             modulationpar3 = 1.0f
    764                              + floorf((powf(2, modulationpar3 * 5.0f) - 1.0f));
    765             break;
    766         case 3:
    767             modulationpar1 = (powf(2, modulationpar1 * 9.0f) - 1.0f) / 100.0f;
    768             modulationpar3 = 0.01f
    769                              + (powf(2, modulationpar3 * 16.0f) - 1.0f) / 10.0f;
    770             break;
    771     }
    772 
    773     clearDC(freqs.data); //remove the DC
    774     //reduce the amplitude of the freqs near the nyquist
    775     for(int i = 1; i < synth.oscilsize / 8; ++i) {
    776         const float tmp = i / (synth.oscilsize / 8.0f);
    777         freqs[synth.oscilsize / 2 - i] *= tmp;
    778     }
    779     fft->freqs2smps_noconst_input(freqs, bfrs.tmpsmps);
    780     const int    extra_points = 2;
    781     float *in = new float[synth.oscilsize + extra_points];
    782 
    783     //Normalize
    784     normalize(bfrs.tmpsmps.data, synth.oscilsize);
    785 
    786     for(int i = 0; i < synth.oscilsize; ++i)
    787         in[i] = bfrs.tmpsmps[i];
    788     for(int i = 0; i < extra_points; ++i)
    789         in[i + synth.oscilsize] = bfrs.tmpsmps[i];
    790 
    791     //Do the modulation
    792     for(int i = 0; i < synth.oscilsize; ++i) {
    793         float t = i * 1.0f / synth.oscilsize;
    794 
    795         switch(Pmodulation) {
    796             case 1:
    797                 t = t * modulationpar3
    798                     + sinf((t + modulationpar2) * 2.0f * PI) * modulationpar1; //rev
    799                 break;
    800             case 2:
    801                 t = t
    802                     + sinf((t * modulationpar3
    803                             + modulationpar2) * 2.0f * PI) * modulationpar1; //sine
    804                 break;
    805             case 3:
    806                 t = t + powf((1.0f - cosf(
    807                                   (t + modulationpar2) * 2.0f * PI)) * 0.5f,
    808                              modulationpar3) * modulationpar1; //power
    809                 break;
    810         }
    811 
    812         t = (t - floorf(t)) * synth.oscilsize;
    813 
    814         const int   poshi = (int) t;
    815         const float poslo = t - floorf(t);
    816 
    817         bfrs.tmpsmps[i] = in[poshi] * (1.0f - poslo) + in[poshi + 1] * poslo;
    818     }
    819 
    820     delete [] in;
    821     fft->smps2freqs_noconst_input(bfrs.tmpsmps, freqs); //perform FFT
    822 }
    823 
    824 
    825 /*
    826  * Adjust the spectrum
    827  */
    828 void OscilGen::spectrumadjust(fft_t *freqs) const
    829 {
    830     if(Psatype == 0)
    831         return;
    832     float par = Psapar / 127.0f;
    833     switch(Psatype) {
    834         case 1:
    835             par = 1.0f - par * 2.0f;
    836             if(par >= 0.0f)
    837                 par = powf(5.0f, par);
    838             else
    839                 par = powf(8.0f, par);
    840             break;
    841         case 2:
    842             par = powf(10.0f, (1.0f - par) * 3.0f) * 0.001f;
    843             break;
    844         case 3:
    845             par = powf(10.0f, (1.0f - par) * 3.0f) * 0.001f;
    846             break;
    847     }
    848 
    849     normalize(freqs, synth.oscilsize);
    850 
    851     for(int i = 0; i < synth.oscilsize / 2; ++i) {
    852         float mag   = abs(freqs, i);
    853         float phase = ((float)M_PI_2) - arg(freqs, i);
    854 
    855         switch(Psatype) {
    856             case 1:
    857                 mag = powf(mag, par);
    858                 break;
    859             case 2:
    860                 if(mag < par)
    861                     mag = 0.0f;
    862                 break;
    863             case 3:
    864                 mag /= par;
    865                 if(mag > 1.0f)
    866                     mag = 1.0f;
    867                 break;
    868         }
    869         freqs[i] = FFTpolar<fftwf_real>(mag, phase);
    870     }
    871 }
    872 
    873 void OscilGen::shiftharmonics(fft_t *freqs) const
    874 {
    875     if(Pharmonicshift == 0)
    876         return;
    877 
    878     int   harmonicshift = -Pharmonicshift;
    879     fft_t h;
    880 
    881     if(harmonicshift > 0)
    882         for(int i = synth.oscilsize / 2 - 2; i >= 0; i--) {
    883             int oldh = i - harmonicshift;
    884             if(oldh < 0)
    885                 h = 0.0f;
    886             else
    887                 h = freqs[oldh + 1];
    888             freqs[i + 1] = h;
    889         }
    890     else
    891         for(int i = 0; i < synth.oscilsize / 2 - 1; ++i) {
    892         int oldh = i + ::abs(harmonicshift);
    893             if(oldh >= (synth.oscilsize / 2 - 1))
    894                 h = 0.0f;
    895             else {
    896                 h = freqs[oldh + 1];
    897                 if(abs(h) < 0.000001f)
    898                     h = 0.0f;
    899             }
    900 
    901             freqs[i + 1] = h;
    902         }
    903 
    904     clearDC(freqs);
    905 }
    906 
    907 /*
    908  * Prepare the Oscillator
    909  */
    910 void OscilGen::prepare(OscilGenBuffers& bfrs) const
    911 {
    912     prepare(bfrs, bfrs.oscilFFTfreqs);
    913 }
    914 
    915 void OscilGen::prepare(OscilGenBuffers& bfrs, FFTfreqBuffer freqs) const
    916 {
    917     if((bfrs.oldbasepar != Pbasefuncpar) || (bfrs.oldbasefunc != Pcurrentbasefunc)
    918        || DIFF(basefuncmodulation) || DIFF(basefuncmodulationpar1)
    919        || DIFF(basefuncmodulationpar2) || DIFF(basefuncmodulationpar3))
    920         changebasefunction(bfrs);
    921 
    922     for(int i = 0; i < MAX_AD_HARMONICS; ++i)
    923         bfrs.hphase[i] = (Phphase[i] - 64.0f) / 64.0f * PI / (i + 1);
    924 
    925     for(int i = 0; i < MAX_AD_HARMONICS; ++i) {
    926         const float hmagnew = 1.0f - fabsf(Phmag[i] / 64.0f - 1.0f);
    927         switch(Phmagtype) {
    928             case 1:
    929                 bfrs.hmag[i] = expf(hmagnew * logf(0.01f));
    930                 break;
    931             case 2:
    932                 bfrs.hmag[i] = expf(hmagnew * logf(0.001f));
    933                 break;
    934             case 3:
    935                 bfrs.hmag[i] = expf(hmagnew * logf(0.0001f));
    936                 break;
    937             case 4:
    938                 bfrs.hmag[i] = expf(hmagnew * logf(0.00001f));
    939                 break;
    940             default:
    941                 bfrs.hmag[i] = 1.0f - hmagnew;
    942                 break;
    943         }
    944 
    945         if(Phmag[i] < 64)
    946             bfrs.hmag[i] = -bfrs.hmag[i];
    947     }
    948 
    949     //remove the harmonics where Phmag[i]==64
    950     for(int i = 0; i < MAX_AD_HARMONICS; ++i)
    951         if(Phmag[i] == 64)
    952             bfrs.hmag[i] = 0.0f;
    953 
    954 
    955     clearAll(freqs.data, synth.oscilsize);
    956     if(Pcurrentbasefunc == 0)   //the sine case
    957         for(int i = 0; i < MAX_AD_HARMONICS - 1; ++i) {
    958             freqs[i + 1] =
    959                 std::complex<float>(-bfrs.hmag[i] * sinf(bfrs.hphase[i] * (i + 1)) / 2.0f,
    960                         bfrs.hmag[i] * cosf(bfrs.hphase[i] * (i + 1)) / 2.0f);
    961         }
    962     else
    963         for(int j = 0; j < MAX_AD_HARMONICS; ++j) {
    964             if(Phmag[j] == 64)
    965                 continue;
    966             for(int i = 1; i < synth.oscilsize / 2; ++i) {
    967                 int k = i * (j + 1);
    968                 if(k >= synth.oscilsize / 2)
    969                     break;
    970                 freqs[k] += bfrs.basefuncFFTfreqs[i] * FFTpolar<fftwf_real>(
    971                     bfrs.hmag[j],
    972                     bfrs.hphase[j] * k);
    973             }
    974         }
    975 
    976     if(Pharmonicshiftfirst != 0)
    977         shiftharmonics(freqs.data);
    978 
    979     if(Pfilterbeforews) {
    980         oscilfilter(freqs.data);
    981         waveshape(bfrs, freqs);
    982     } else {
    983         waveshape(bfrs, freqs);
    984         oscilfilter(freqs.data);
    985     }
    986 
    987     modulation(bfrs, freqs);
    988     spectrumadjust(freqs.data);
    989     if(Pharmonicshiftfirst == 0)
    990         shiftharmonics(freqs.data);
    991 
    992     clearDC(freqs.data);
    993 
    994     bfrs.oldhmagtype      = Phmagtype;
    995     bfrs.oldharmonicshift = Pharmonicshift + Pharmonicshiftfirst * 256;
    996 
    997     bfrs.oscilprepared = 1;
    998 }
    999 
   1000 fft_t operator*(float a, fft_t b)
   1001 {
   1002     return std::complex<float>((float)(a*b.real()), (float)(a*b.imag()));
   1003 }
   1004 
   1005 void OscilGen::adaptiveharmonic(fft_t *f, float freq) const
   1006 {
   1007     if(Padaptiveharmonics == 0 /*||(freq<1.0f)*/)
   1008         return;
   1009     if(freq < 1.0f)
   1010         freq = 440.0f;
   1011 
   1012     fft_t *inf = new fft_t[synth.oscilsize / 2]; // XXXRT
   1013     for(int i = 0; i < synth.oscilsize / 2; ++i)
   1014         inf[i] = f[i];
   1015     clearAll(f, synth.oscilsize);
   1016     clearDC(inf);
   1017 
   1018     float basefreq = 30.0f * powf(10.0f, Padaptiveharmonicsbasefreq / 128.0f);
   1019     float power    = (Padaptiveharmonicspower + 1.0f) / 101.0f;
   1020 
   1021     float rap = freq / basefreq;
   1022 
   1023     rap = powf(rap, power);
   1024 
   1025     bool down = false;
   1026     if(rap > 1.0f) {
   1027         rap  = 1.0f / rap;
   1028         down = true;
   1029     }
   1030 
   1031     for(int i = 0; i < synth.oscilsize / 2 - 2; ++i) {
   1032         const int   high = (int)(i * rap);
   1033         const float low  = fmodf(i * rap, 1.0f);
   1034 
   1035         if(high >= (synth.oscilsize / 2 - 2))
   1036             break;
   1037 
   1038         if(down) {
   1039             f[high] += (1.0f - low) * inf[i];
   1040             f[high + 1] += low * inf[i];
   1041         }
   1042         else {
   1043             f[i] = (1.0f - low) * inf[high] + low * inf[high + 1];
   1044         }
   1045     }
   1046     if(!down)//correct the amplitude of the first harmonic
   1047         f[0] *= rap;
   1048 
   1049     f[1] += f[0];
   1050     clearDC(f);
   1051     delete[] inf;
   1052 }
   1053 
   1054 void OscilGen::adaptiveharmonicpostprocess(fft_t *f, int size) const
   1055 {
   1056     if(Padaptiveharmonics <= 1)
   1057         return;
   1058     fft_t *inf = new fft_t[size]; // XXXRT
   1059     float  par = Padaptiveharmonicspar * 0.01f;
   1060     par = 1.0f - powf((1.0f - par), 1.5f);
   1061 
   1062     for(int i = 0; i < size; ++i) {
   1063         inf[i] = f[i] * par;
   1064         f[i]  *= (1.0f - par);
   1065     }
   1066 
   1067 
   1068     if(Padaptiveharmonics == 2) { //2n+1
   1069         for(int i = 0; i < size; ++i)
   1070             if((i % 2) == 0)
   1071                 f[i] += inf[i];  //i=0 first harmonic,etc.
   1072     }
   1073     else {  //other ways
   1074         int nh = (Padaptiveharmonics - 3) / 2 + 2;
   1075         int sub_vs_add = (Padaptiveharmonics - 3) % 2;
   1076         if(sub_vs_add == 0) {
   1077             for(int i = 0; i < size; ++i)
   1078                 if(((i + 1) % nh) == 0)
   1079                     f[i] += inf[i];
   1080         }
   1081         else
   1082             for(int i = 0; i < size / nh - 1; ++i)
   1083                 f[(i + 1) * nh - 1] += inf[i];
   1084     }
   1085 
   1086     delete [] inf;
   1087 }
   1088 
   1089 void OscilGen::newrandseed(unsigned int randseed)
   1090 {
   1091     this->randseed = randseed;
   1092 }
   1093 
   1094 bool OscilGen::needPrepare(OscilGenBuffers& bfrs) const
   1095 {
   1096     bool outdated = false;
   1097 
   1098     //Check function parameters
   1099     if((bfrs.oldbasepar != Pbasefuncpar) || (bfrs.oldbasefunc != Pcurrentbasefunc)
   1100        || DIFF(hmagtype) || DIFF(waveshaping) || DIFF(waveshapingfunction))
   1101         outdated = true;
   1102 
   1103     //Check filter parameters
   1104     if(bfrs.oldfilterpars != Pfiltertype * 256 + Pfilterpar1 + Pfilterpar2 * 65536
   1105        + Pfilterbeforews * 16777216) {
   1106         outdated      = true;
   1107         bfrs.oldfilterpars = Pfiltertype * 256 + Pfilterpar1 + Pfilterpar2 * 65536
   1108                         + Pfilterbeforews * 16777216;
   1109     }
   1110 
   1111     //Check spectrum adjustments
   1112     if(bfrs.oldsapars != Psatype * 256 + Psapar) {
   1113         outdated  = true;
   1114         bfrs.oldsapars = Psatype * 256 + Psapar;
   1115     }
   1116 
   1117     //Check function modulation
   1118     if(DIFF(basefuncmodulation) || DIFF(basefuncmodulationpar1)
   1119        || DIFF(basefuncmodulationpar2) || DIFF(basefuncmodulationpar3))
   1120         outdated = true;
   1121 
   1122     //Check overall modulation
   1123     if(DIFF(modulation) || DIFF(modulationpar1)
   1124        || DIFF(modulationpar2) || DIFF(modulationpar3))
   1125         outdated = true;
   1126 
   1127     //Check harmonic shifts
   1128     if(bfrs.oldharmonicshift != Pharmonicshift + Pharmonicshiftfirst * 256)
   1129         outdated = true;
   1130 
   1131     return outdated == true || bfrs.oscilprepared == false;
   1132 }
   1133 
   1134 /*
   1135  * Get the oscillator function
   1136  */
   1137 short int OscilGen::get(OscilGenBuffers& bfrs, float* smps, float freqHz, int resonance) const
   1138 {
   1139     if(needPrepare(bfrs))
   1140         prepare(bfrs);
   1141 
   1142     fft_t *input = freqHz > 0.0f ? bfrs.oscilFFTfreqs.data : bfrs.pendingfreqs;
   1143 
   1144     unsigned int realrnd = prng();
   1145     sprng(randseed);
   1146 
   1147     int outpos =
   1148         (int)((RND * 2.0f
   1149                - 1.0f) * synth.oscilsize_f * (Prand - 64.0f) / 64.0f);
   1150     outpos = (outpos + 2 * synth.oscilsize) % synth.oscilsize;
   1151 
   1152 
   1153     clearAll(bfrs.outoscilFFTfreqs.data, synth.oscilsize);
   1154 
   1155     int nyquist = (int)(0.5f * synth.samplerate_f / fabsf(freqHz)) + 2;
   1156     if(ADvsPAD)
   1157         nyquist = (int)(synth.oscilsize / 2);
   1158     if(nyquist > synth.oscilsize / 2)
   1159         nyquist = synth.oscilsize / 2;
   1160 
   1161     //Process harmonics
   1162     {
   1163         int realnyquist = nyquist;
   1164 
   1165         if(Padaptiveharmonics != 0)
   1166             nyquist = synth.oscilsize / 2;
   1167         for(int i = 1; i < nyquist - 1; ++i)
   1168             bfrs.outoscilFFTfreqs[i] = input[i];
   1169 
   1170         adaptiveharmonic(bfrs.outoscilFFTfreqs.data, freqHz);
   1171         adaptiveharmonicpostprocess(&bfrs.outoscilFFTfreqs[1],
   1172                                     synth.oscilsize / 2 - 1);
   1173 
   1174         nyquist = realnyquist;
   1175     }
   1176 
   1177     if(Padaptiveharmonics)   //do the antialiasing in the case of adaptive harmonics
   1178         for(int i = nyquist; i < synth.oscilsize / 2; ++i)
   1179             bfrs.outoscilFFTfreqs[i] = fft_t(0.0f, 0.0f);
   1180 
   1181     // Randomness (each harmonic), the block type is computed
   1182     // in ADnote by setting start position according to this setting
   1183     if((Prand > 64) && (freqHz >= 0.0f) && (!ADvsPAD)) {
   1184         const float rnd = PI * powf((Prand - 64.0f) / 64.0f, 2.0f);
   1185         for(int i = 1; i < nyquist - 1; ++i) //to Nyquist only for AntiAliasing
   1186             bfrs.outoscilFFTfreqs[i] *=
   1187                 FFTpolar<fftwf_real>(1.0f, (float)(rnd * i * RND));
   1188     }
   1189 
   1190     //Harmonic Amplitude Randomness
   1191     if((freqHz > 0.1f) && (!ADvsPAD)) {
   1192 
   1193         float power     = Pamprandpower / 127.0f;
   1194         float normalize = 1.0f / (1.2f - power);
   1195         switch(Pamprandtype) {
   1196             case 1:
   1197                 power = power * 2.0f - 0.5f;
   1198                 power = powf(15.0f, power);
   1199                 for(int i = 1; i < nyquist - 1; ++i)
   1200                     bfrs.outoscilFFTfreqs[i] *= powf(RND, power) * normalize;
   1201                 break;
   1202             case 2:
   1203                 power = power * 2.0f - 0.5f;
   1204                 power = powf(15.0f, power) * 2.0f;
   1205                 float rndfreq = 2 * PI * RND;
   1206                 for(int i = 1; i < nyquist - 1; ++i)
   1207                     bfrs.outoscilFFTfreqs[i] *= powf(fabsf(sinf(i * rndfreq)), power)
   1208                                            * normalize;
   1209                 break;
   1210         }
   1211     }
   1212 
   1213     if((freqHz > 0.1f) && (resonance != 0))
   1214         res->applyres(nyquist - 1, bfrs.outoscilFFTfreqs.data, freqHz);
   1215 
   1216     rmsNormalize(bfrs.outoscilFFTfreqs.data, synth.oscilsize);
   1217 
   1218     if((ADvsPAD) && (freqHz > 0.1f)) //in this case the smps will contain the freqs
   1219         for(int i = 1; i < synth.oscilsize / 2; ++i)
   1220             smps[i - 1] = abs(bfrs.outoscilFFTfreqs.data, i);
   1221     else {
   1222         fft->freqs2smps(bfrs.outoscilFFTfreqs, bfrs.tmpsmps, bfrs.scratchFreqs);
   1223         for(int i = 0; i < synth.oscilsize; ++i)
   1224             smps[i] = bfrs.tmpsmps[i] * 0.25f;            //correct the amplitude
   1225     }
   1226 
   1227     sprng(realrnd + 1);
   1228 
   1229     if(Prand < 64)
   1230         return outpos;
   1231     else
   1232         return 0;
   1233 }
   1234 
   1235 ///*
   1236 // * Get the oscillator function's harmonics
   1237 // */
   1238 //void OscilGen::getPad(float *smps, float freqHz)
   1239 //{
   1240 //    if(needPrepare())
   1241 //        prepare();
   1242 //
   1243 //    clearAll(bfrs.outoscilFFTfreqs);
   1244 //
   1245 //    const int nyquist = (synth.oscilsize / 2);
   1246 //
   1247 //    //Process harmonics
   1248 //    for(int i = 1; i < nyquist - 1; ++i)
   1249 //        bfrs.outoscilFFTfreqs[i] = bfrs.oscilFFTfreqs[i];
   1250 //
   1251 //    adaptiveharmonic(bfrs.outoscilFFTfreqs, freqHz);
   1252 //    adaptiveharmonicpostprocess(&bfrs.outoscilFFTfreqs[1], nyquist - 1);
   1253 //
   1254 //    rmsNormalize(bfrs.outoscilFFTfreqs);
   1255 //
   1256 //    for(int i = 1; i < nyquist; ++i)
   1257 //        smps[i - 1] = abs(bfrs.outoscilFFTfreqs, i);
   1258 //}
   1259 //
   1260 
   1261 /*
   1262  * Get the spectrum of the oscillator for the UI
   1263  */
   1264 void OscilGen::getspectrum(int n, float *spc, int what)
   1265 {
   1266     OscilGenBuffers& bfrs = myBuffers();
   1267 
   1268     if(n > synth.oscilsize / 2)
   1269         n = synth.oscilsize / 2;
   1270 
   1271     for(int i = 1; i < n; ++i) {
   1272         if(what == 0)
   1273             spc[i] = abs(bfrs.pendingfreqs, i);
   1274         else {
   1275             if(Pcurrentbasefunc == 0)
   1276                 spc[i] = ((i == 1) ? (1.0f) : (0.0f));
   1277             else
   1278                 spc[i] = abs(bfrs.basefuncFFTfreqs.data, i);
   1279         }
   1280     }
   1281     spc[0]=0;
   1282 
   1283     if(what == 0) {
   1284         for(int i = 0; i < n; ++i)
   1285             bfrs.outoscilFFTfreqs[i] = fft_t(spc[i], spc[i]);
   1286         fft_t zero = 0;
   1287         std::fill_n(bfrs.outoscilFFTfreqs.data + n, synth.oscilsize / 2 - n, zero);
   1288         adaptiveharmonic(bfrs.outoscilFFTfreqs.data, 0.0f);
   1289         adaptiveharmonicpostprocess(bfrs.outoscilFFTfreqs.data, n - 1);
   1290         for(int i = 0; i < n; ++i)
   1291             spc[i] = (float)bfrs.outoscilFFTfreqs[i].imag();
   1292     }
   1293 }
   1294 
   1295 
   1296 /*
   1297  * Convert the oscillator as base function
   1298  */
   1299 void OscilGen::useasbase()
   1300 {
   1301     OscilGenBuffers& bfrs = myBuffers();
   1302 
   1303     for(int i = 0; i < synth.oscilsize / 2; ++i)
   1304         bfrs.basefuncFFTfreqs[i] = bfrs.oscilFFTfreqs[i];
   1305 
   1306     bfrs.oldbasefunc = Pcurrentbasefunc = 127;
   1307     prepare(bfrs);
   1308     bfrs.cachedbasevalid = false;
   1309 }
   1310 
   1311 
   1312 /*
   1313  * Get the base function for UI
   1314  */
   1315 void OscilGen::getcurrentbasefunction(FFTsampleBuffer smps)
   1316 {
   1317     OscilGenBuffers& bfrs = myBuffers();
   1318     if(Pcurrentbasefunc != 0)
   1319         fft->freqs2smps(bfrs.basefuncFFTfreqs, smps, bfrs.scratchFreqs);
   1320     else
   1321         getbasefunction(bfrs, smps);   //the sine case
   1322 }
   1323 
   1324 #define COPY(y) this->y = o.y
   1325 void OscilGen::paste(OscilGen &o)
   1326 {
   1327     OscilGenBuffers& bfrs = myBuffers();
   1328 
   1329     //XXX Figure out a better implementation of this sensitive to RT issues...
   1330     for(int i=0; i<MAX_AD_HARMONICS; ++i) {
   1331         COPY(Phmag[i]);
   1332         COPY(Phphase[i]);
   1333     }
   1334 
   1335     COPY(Phmagtype);
   1336     COPY(Pcurrentbasefunc);
   1337     COPY(Pbasefuncpar);
   1338 
   1339     COPY(Pbasefuncmodulation);
   1340     COPY(Pbasefuncmodulationpar1);
   1341     COPY(Pbasefuncmodulationpar2);
   1342     COPY(Pbasefuncmodulationpar3);
   1343 
   1344     COPY(Pwaveshaping);
   1345     COPY(Pwaveshapingfunction);
   1346     COPY(Pfiltertype);
   1347     COPY(Pfilterpar1);
   1348     COPY(Pfilterpar2);
   1349     COPY(Pfilterbeforews);
   1350     COPY(Psatype);
   1351     COPY(Psapar);
   1352 
   1353     COPY(Pharmonicshift);
   1354     COPY(Pharmonicshiftfirst);
   1355 
   1356     COPY(Pmodulation);
   1357     COPY(Pmodulationpar1);
   1358     COPY(Pmodulationpar2);
   1359     COPY(Pmodulationpar3);
   1360 
   1361     COPY(Prand);
   1362     COPY(Pamprandpower);
   1363     COPY(Pamprandtype);
   1364     COPY(Padaptiveharmonics);
   1365     COPY(Padaptiveharmonicsbasefreq);
   1366     COPY(Padaptiveharmonicspower);
   1367     COPY(Padaptiveharmonicspar);
   1368 
   1369 
   1370     if(this->Pcurrentbasefunc)
   1371         changebasefunction(bfrs);
   1372     this->prepare(bfrs);
   1373 }
   1374 #undef COPY
   1375 
   1376 void OscilGen::add2XML(XMLwrapper& xml)
   1377 {
   1378     const OscilGenBuffers& bfrs = myBuffers();
   1379 
   1380     xml.addpar("harmonic_mag_type", Phmagtype);
   1381 
   1382     xml.addpar("base_function", Pcurrentbasefunc);
   1383     xml.addpar("base_function_par", Pbasefuncpar);
   1384     xml.addpar("base_function_modulation", Pbasefuncmodulation);
   1385     xml.addpar("base_function_modulation_par1", Pbasefuncmodulationpar1);
   1386     xml.addpar("base_function_modulation_par2", Pbasefuncmodulationpar2);
   1387     xml.addpar("base_function_modulation_par3", Pbasefuncmodulationpar3);
   1388 
   1389     xml.addpar("modulation", Pmodulation);
   1390     xml.addpar("modulation_par1", Pmodulationpar1);
   1391     xml.addpar("modulation_par2", Pmodulationpar2);
   1392     xml.addpar("modulation_par3", Pmodulationpar3);
   1393 
   1394     xml.addpar("wave_shaping", Pwaveshaping);
   1395     xml.addpar("wave_shaping_function", Pwaveshapingfunction);
   1396 
   1397     xml.addpar("filter_type", Pfiltertype);
   1398     xml.addpar("filter_par1", Pfilterpar1);
   1399     xml.addpar("filter_par2", Pfilterpar2);
   1400     xml.addpar("filter_before_wave_shaping", Pfilterbeforews);
   1401 
   1402     xml.addpar("spectrum_adjust_type", Psatype);
   1403     xml.addpar("spectrum_adjust_par", Psapar);
   1404 
   1405     xml.addpar("rand", Prand);
   1406     xml.addpar("amp_rand_type", Pamprandtype);
   1407     xml.addpar("amp_rand_power", Pamprandpower);
   1408 
   1409     xml.addpar("harmonic_shift", Pharmonicshift);
   1410     xml.addparbool("harmonic_shift_first", Pharmonicshiftfirst);
   1411 
   1412     xml.addpar("adaptive_harmonics", Padaptiveharmonics);
   1413     xml.addpar("adaptive_harmonics_base_frequency", Padaptiveharmonicsbasefreq);
   1414     xml.addpar("adaptive_harmonics_power", Padaptiveharmonicspower);
   1415     xml.addpar("adaptive_harmonics_par", Padaptiveharmonicspar);
   1416 
   1417     xml.beginbranch("HARMONICS");
   1418     for(int n = 0; n < MAX_AD_HARMONICS; ++n) {
   1419         if((Phmag[n] == 64) && (Phphase[n] == 64))
   1420             continue;
   1421         xml.beginbranch("HARMONIC", n + 1);
   1422         xml.addpar("mag", Phmag[n]);
   1423         xml.addpar("phase", Phphase[n]);
   1424         xml.endbranch();
   1425     }
   1426     xml.endbranch();
   1427 
   1428     if(Pcurrentbasefunc == 127) {
   1429         normalize(bfrs.basefuncFFTfreqs.data, synth.oscilsize);
   1430 
   1431         xml.beginbranch("BASE_FUNCTION");
   1432         for(int i = 1; i < synth.oscilsize / 2; ++i) {
   1433             float xc = (float)bfrs.basefuncFFTfreqs[i].real();
   1434             float xs = (float)bfrs.basefuncFFTfreqs[i].imag();
   1435             if((fabsf(xs) > 1e-6f) || (fabsf(xc) > 1e-6f)) {
   1436                 xml.beginbranch("BF_HARMONIC", i);
   1437                 xml.addparreal("cos", xc);
   1438                 xml.addparreal("sin", xs);
   1439                 xml.endbranch();
   1440             }
   1441         }
   1442         xml.endbranch();
   1443     }
   1444 }
   1445 
   1446 void OscilGen::getfromXML(XMLwrapper& xml)
   1447 {
   1448     OscilGenBuffers& bfrs = myBuffers();
   1449 
   1450     Phmagtype = xml.getpar127("harmonic_mag_type", Phmagtype);
   1451 
   1452     Pcurrentbasefunc = xml.getpar127("base_function", Pcurrentbasefunc);
   1453     Pbasefuncpar     = xml.getpar127("base_function_par", Pbasefuncpar);
   1454 
   1455     Pbasefuncmodulation = xml.getpar127("base_function_modulation",
   1456                                          Pbasefuncmodulation);
   1457     Pbasefuncmodulationpar1 = xml.getpar127("base_function_modulation_par1",
   1458                                              Pbasefuncmodulationpar1);
   1459     Pbasefuncmodulationpar2 = xml.getpar127("base_function_modulation_par2",
   1460                                              Pbasefuncmodulationpar2);
   1461     Pbasefuncmodulationpar3 = xml.getpar127("base_function_modulation_par3",
   1462                                              Pbasefuncmodulationpar3);
   1463 
   1464     Pmodulation     = xml.getpar127("modulation", Pmodulation);
   1465     Pmodulationpar1 = xml.getpar127("modulation_par1",
   1466                                      Pmodulationpar1);
   1467     Pmodulationpar2 = xml.getpar127("modulation_par2",
   1468                                      Pmodulationpar2);
   1469     Pmodulationpar3 = xml.getpar127("modulation_par3",
   1470                                      Pmodulationpar3);
   1471 
   1472     Pwaveshaping = xml.getpar127("wave_shaping", Pwaveshaping);
   1473     Pwaveshapingfunction = xml.getpar127("wave_shaping_function",
   1474                                           Pwaveshapingfunction);
   1475 
   1476     Pfiltertype     = xml.getpar127("filter_type", Pfiltertype);
   1477     Pfilterpar1     = xml.getpar127("filter_par1", Pfilterpar1);
   1478     Pfilterpar2     = xml.getpar127("filter_par2", Pfilterpar2);
   1479     Pfilterbeforews = xml.getpar127("filter_before_wave_shaping",
   1480                                      Pfilterbeforews);
   1481 
   1482     Psatype = xml.getpar127("spectrum_adjust_type", Psatype);
   1483     Psapar  = xml.getpar127("spectrum_adjust_par", Psapar);
   1484 
   1485     Prand = xml.getpar127("rand", Prand);
   1486     Pamprandtype  = xml.getpar127("amp_rand_type", Pamprandtype);
   1487     Pamprandpower = xml.getpar127("amp_rand_power", Pamprandpower);
   1488 
   1489     Pharmonicshift = xml.getpar("harmonic_shift",
   1490                                  Pharmonicshift,
   1491                                  -64,
   1492                                  64);
   1493     Pharmonicshiftfirst = xml.getparbool("harmonic_shift_first",
   1494                                           Pharmonicshiftfirst);
   1495 
   1496     Padaptiveharmonics = xml.getpar("adaptive_harmonics",
   1497                                      Padaptiveharmonics,
   1498                                      0,
   1499                                      127);
   1500     Padaptiveharmonicsbasefreq = xml.getpar(
   1501         "adaptive_harmonics_base_frequency",
   1502         Padaptiveharmonicsbasefreq,
   1503         0,
   1504         255);
   1505     Padaptiveharmonicspower = xml.getpar("adaptive_harmonics_power",
   1506                                           Padaptiveharmonicspower,
   1507                                           0,
   1508                                           200);
   1509     Padaptiveharmonicspar = xml.getpar("adaptive_harmonics_par",
   1510                                        Padaptiveharmonicspar,
   1511                                        0,
   1512                                        100);
   1513 
   1514 
   1515     if(xml.enterbranch("HARMONICS")) {
   1516         Phmag[0]   = 64;
   1517         Phphase[0] = 64;
   1518         for(int n = 0; n < MAX_AD_HARMONICS; ++n) {
   1519             if(xml.enterbranch("HARMONIC", n + 1) == 0)
   1520                 continue;
   1521             Phmag[n]   = xml.getpar127("mag", 64);
   1522             Phphase[n] = xml.getpar127("phase", 64);
   1523             xml.exitbranch();
   1524         }
   1525         xml.exitbranch();
   1526     }
   1527 
   1528     if(Pcurrentbasefunc != 0)
   1529         changebasefunction(bfrs);
   1530 
   1531     if(xml.enterbranch("BASE_FUNCTION")) {
   1532         for(int i = 1; i < synth.oscilsize / 2; ++i)
   1533             if(xml.enterbranch("BF_HARMONIC", i)) {
   1534                 bfrs.basefuncFFTfreqs[i] =
   1535                     std::complex<float>(xml.getparreal("cos", 0.0f),
   1536                             xml.getparreal("sin", 0.0f));
   1537                 xml.exitbranch();
   1538             }
   1539         xml.exitbranch();
   1540 
   1541         clearDC(bfrs.basefuncFFTfreqs.data);
   1542         normalize(bfrs.basefuncFFTfreqs.data, synth.oscilsize);
   1543         bfrs.cachedbasevalid = false;
   1544     }
   1545 }
   1546 
   1547 
   1548 //Define basic functions
   1549 #define FUNC(b) float basefunc_ ## b(float x, float a)
   1550 
   1551 FUNC(pulse)
   1552 {
   1553     return (fmodf(x, 1.0f) < a) ? -1.0f : 1.0f;
   1554 }
   1555 
   1556 FUNC(saw)
   1557 {
   1558     if(a < 0.00001f)
   1559         a = 0.00001f;
   1560     else
   1561     if(a > 0.99999f)
   1562         a = 0.99999f;
   1563     x = fmodf(x, 1);
   1564     if(x < a)
   1565         return x / a * 2.0f - 1.0f;
   1566     else
   1567         return (1.0f - x) / (1.0f - a) * 2.0f - 1.0f;
   1568 }
   1569 
   1570 FUNC(triangle)
   1571 {
   1572     x = fmodf(x + 0.25f, 1);
   1573     a = 1 - a;
   1574     if(a < 0.00001f)
   1575         a = 0.00001f;
   1576     if(x < 0.5f)
   1577         x = x * 4 - 1.0f;
   1578     else
   1579         x = (1.0f - x) * 4 - 1.0f;
   1580     x /= -a;
   1581     if(x < -1.0f)
   1582         x = -1.0f;
   1583     if(x > 1.0f)
   1584         x = 1.0f;
   1585     return x;
   1586 }
   1587 
   1588 FUNC(power)
   1589 {
   1590     x = fmodf(x, 1);
   1591     if(a < 0.00001f)
   1592         a = 0.00001f;
   1593     else
   1594     if(a > 0.99999f)
   1595         a = 0.99999f;
   1596     return powf(x, expf((a - 0.5f) * 10.0f)) * 2.0f - 1.0f;
   1597 }
   1598 
   1599 FUNC(gauss)
   1600 {
   1601     x = fmodf(x, 1) * 2.0f - 1.0f;
   1602     if(a < 0.00001f)
   1603         a = 0.00001f;
   1604     return expf(-x * x * (expf(a * 8) + 5.0f)) * 2.0f - 1.0f;
   1605 }
   1606 
   1607 FUNC(diode)
   1608 {
   1609     if(a < 0.00001f)
   1610         a = 0.00001f;
   1611     else
   1612     if(a > 0.99999f)
   1613         a = 0.99999f;
   1614     a = a * 2.0f - 1.0f;
   1615     x = cosf((x + 0.5f) * 2.0f * PI) - a;
   1616     if(x < 0.0f)
   1617         x = 0.0f;
   1618     return x / (1.0f - a) * 2 - 1.0f;
   1619 }
   1620 
   1621 FUNC(abssine)
   1622 {
   1623     x = fmodf(x, 1);
   1624     if(a < 0.00001f)
   1625         a = 0.00001f;
   1626     else
   1627     if(a > 0.99999f)
   1628         a = 0.99999f;
   1629     return sinf(powf(x, expf((a - 0.5f) * 5.0f)) * PI) * 2.0f - 1.0f;
   1630 }
   1631 
   1632 FUNC(pulsesine)
   1633 {
   1634     if(a < 0.00001f)
   1635         a = 0.00001f;
   1636     x = (fmodf(x, 1) - 0.5f) * expf((a - 0.5f) * logf(128));
   1637     if(x < -0.5f)
   1638         x = -0.5f;
   1639     else
   1640     if(x > 0.5f)
   1641         x = 0.5f;
   1642     x = sinf(x * PI * 2.0f);
   1643     return x;
   1644 }
   1645 
   1646 FUNC(stretchsine)
   1647 {
   1648     x = fmodf(x + 0.5f, 1) * 2.0f - 1.0f;
   1649     a = (a - 0.5f) * 4;
   1650     if(a > 0.0f)
   1651         a *= 2;
   1652     a = powf(3.0f, a);
   1653     float b = powf(fabsf(x), a);
   1654     if(x < 0)
   1655         b = -b;
   1656     return -sinf(b * PI);
   1657 }
   1658 
   1659 FUNC(chirp)
   1660 {
   1661     x = fmodf(x, 1.0f) * 2.0f * PI;
   1662     a = (a - 0.5f) * 4;
   1663     if(a < 0.0f)
   1664         a *= 2.0f;
   1665     a = powf(3.0f, a);
   1666     return sinf(x / 2.0f) * sinf(a * x * x);
   1667 }
   1668 
   1669 FUNC(absstretchsine)
   1670 {
   1671     x = fmodf(x + 0.5f, 1) * 2.0f - 1.0f;
   1672     a = (a - 0.5f) * 9;
   1673     a = powf(3.0f, a);
   1674     float b = powf(fabsf(x), a);
   1675     if(x < 0)
   1676         b = -b;
   1677     return -powf(sinf(b * PI), 2);
   1678 }
   1679 
   1680 FUNC(chebyshev)
   1681 {
   1682     a = a * a * a * 30.0f + 1.0f;
   1683     return cosf(acosf(x * 2.0f - 1.0f) * a);
   1684 }
   1685 
   1686 FUNC(sqr)
   1687 {
   1688     a = a * a * a * a * 160.0f + 0.001f;
   1689     return -atanf(sinf(x * 2.0f * PI) * a);
   1690 }
   1691 
   1692 FUNC(spike)
   1693 {
   1694     float b = a * 0.66666f; // the width of the range: if a == 0.5, b == 0.33333
   1695 
   1696     if(x < 0.5) {
   1697         if(x < (0.5 - (b / 2.0)))
   1698             return 0.0;
   1699         else {
   1700             x = (x + (b / 2) - 0.5f) * (2.f / b); // shift to zero, and expand to range from 0 to 1
   1701             return x * (2.f / b); // this is the slope: 1 / (b / 2)
   1702         }
   1703     }
   1704     else {
   1705         if(x > (0.5 + (b / 2.0)))
   1706             return 0.0;
   1707         else {
   1708             x = (x - .5f) * (2.f / b);
   1709             return (1 - x) * (2.f / b);
   1710         }
   1711     }
   1712 }
   1713 
   1714 FUNC(circle)
   1715 {
   1716     // a is parameter: 0 -> 0.5 -> 1 // O.5 = circle
   1717     float b, y;
   1718 
   1719     b = 2 - (a * 2); // b goes from 2 to 0
   1720     x = x * 4;
   1721 
   1722     if(x < 2) {
   1723         x = x - 1; // x goes from -1 to 1
   1724         if((x < -b) || (x > b))
   1725             y = 0;
   1726         else
   1727             y = sqrtf(1 - (powf(x, 2) / powf(b, 2)));  // normally * a^2, but a stays 1
   1728     }
   1729     else {
   1730         x = x - 3; // x goes from -1 to 1 as well
   1731         if((x < -b) || (x > b))
   1732             y = 0;
   1733         else
   1734             y = -sqrtf(1 - (powf(x, 2) / powf(b, 2)));
   1735     }
   1736     return y;
   1737 }
   1738 
   1739 static float
   1740 power_cosinus_32(float _x, double _power)
   1741 {
   1742     uint32_t x = (_x - floorf(_x)) * (1ULL << 32);
   1743     double retval;
   1744     uint8_t num;
   1745 
   1746     /* Handle special cases, if any */
   1747     switch (x) {
   1748     case 0xFFFFFFFFU:
   1749     case 0x00000000U:
   1750         return (1.0f);
   1751     case 0x3FFFFFFFU:
   1752     case 0x40000000U:
   1753     case 0xBFFFFFFFU:
   1754     case 0xC0000000U:
   1755         return (0.0f);
   1756     case 0x7FFFFFFFU:
   1757     case 0x80000000U:
   1758         return (-1.0f);
   1759     }
   1760 
   1761     /* Apply "grey" encoding */
   1762     for (uint32_t mask = 1U << 31; mask != 1; mask /= 2) {
   1763         if (x & mask)
   1764             x ^= (mask - 1);
   1765     }
   1766 
   1767     /* Find first set bit */
   1768     for (num = 0; num != 30; num++) {
   1769         if (x & (1U << num)) {
   1770             num++;
   1771             break;
   1772         }
   1773     }
   1774 
   1775     /* Initialize return value */
   1776     retval = 0.0;
   1777 
   1778     /* Compute the rest of the power series */
   1779     for (; num != 30; num++) {
   1780         if (x & (1U << num))
   1781             retval = pow((1.0 - retval) / 2.0, _power);
   1782         else
   1783             retval = pow((1.0 + retval) / 2.0, _power);
   1784     }
   1785 
   1786     /* Check if halfway */
   1787     if (x & (1ULL << 30))
   1788         retval = -retval;
   1789 
   1790     return (retval);
   1791 }
   1792 
   1793 //
   1794 // power argument magic values:
   1795 //     0.0: Converges to a square wave
   1796 //     0.5: Sinus wave
   1797 //     1.0: Triangle wave
   1798 // x: phase value [0..1>
   1799 //
   1800 static float
   1801 power_sinus_32(float _x, double _power)
   1802 {
   1803     return (power_cosinus_32(_x + 0.75f, _power));
   1804 }
   1805 
   1806 FUNC(powersinus)
   1807 {
   1808     return (power_sinus_32(x, 2.0 * a));
   1809 }
   1810 
   1811 base_func_t *getBaseFunction(unsigned char func)
   1812 {
   1813     static base_func_t * const functions[] = {
   1814         basefunc_triangle,
   1815         basefunc_pulse,
   1816         basefunc_saw,
   1817         basefunc_power,
   1818         basefunc_gauss,
   1819         basefunc_diode,
   1820         basefunc_abssine,
   1821         basefunc_pulsesine,
   1822         basefunc_stretchsine,
   1823         basefunc_chirp,
   1824         basefunc_absstretchsine,
   1825         basefunc_chebyshev,
   1826         basefunc_sqr,
   1827         basefunc_spike,
   1828         basefunc_circle,
   1829         basefunc_powersinus,
   1830     };
   1831 
   1832     if(!func)
   1833         return NULL;
   1834 
   1835     if(func == 127) //should be the custom wave
   1836         return NULL;
   1837 
   1838     func--;
   1839     assert(func < (sizeof(functions)/sizeof(functions[0])));
   1840     return functions[func];
   1841 }
   1842 
   1843 //And filters
   1844 
   1845 #define FILTER(x) float osc_ ## x(unsigned int i, float par, float par2)
   1846 FILTER(lp)
   1847 {
   1848     float gain = powf(1.0f - par * par * par * 0.99f, i);
   1849     float tmp  = par2 * par2 * par2 * par2 * 0.5f + 0.0001f;
   1850     if(gain < tmp)
   1851         gain = powf(gain, 10.0f) / powf(tmp, 9.0f);
   1852     return gain;
   1853 }
   1854 
   1855 FILTER(hp1)
   1856 {
   1857     float gain = 1.0f - powf(1.0f - par * par, i + 1);
   1858     return powf(gain, par2 * 2.0f + 0.1f);
   1859 }
   1860 
   1861 FILTER(hp1b)
   1862 {
   1863     if(par < 0.2f)
   1864         par = par * 0.25f + 0.15f;
   1865     float gain = 1.0f - powf(1.0f - par * par * 0.999f + 0.001f,
   1866                              i * 0.05f * i + 1.0f);
   1867     float tmp = powf(5.0f, par2 * 2.0f);
   1868     return powf(gain, tmp);
   1869 }
   1870 
   1871 FILTER(bp1)
   1872 {
   1873     float gain = i + 1 - powf(2, (1.0f - par) * 7.5f);
   1874     gain = 1.0f / (1.0f + gain * gain / (i + 1.0f));
   1875     float tmp = powf(5.0f, par2 * 2.0f);
   1876     gain = powf(gain, tmp);
   1877     if(gain < 1e-5)
   1878         gain = (float)1e-5;
   1879     return gain;
   1880 }
   1881 
   1882 FILTER(bs1)
   1883 {
   1884     float gain = i + 1 - powf(2, (1.0f - par) * 7.5f);
   1885     gain = powf(atanf(gain / (i / 10.0f + 1)) / 1.57f, 6);
   1886     return powf(gain, par2 * par2 * 3.9f + 0.1f);
   1887 }
   1888 
   1889 FILTER(lp2)
   1890 {
   1891     return (i + 1 >
   1892             powf(2, (1.0f - par) * 10) ? 0.0f : 1.0f) * par2 + (1.0f - par2);
   1893 }
   1894 
   1895 FILTER(hp2)
   1896 {
   1897     if(par == 1)
   1898         return 1.0f;
   1899     return (i + 1 >
   1900             powf(2, (1.0f - par) * 7) ? 1.0f : 0.0f) * par2 + (1.0f - par2);
   1901 }
   1902 
   1903 FILTER(bp2)
   1904 {
   1905     return (fabsf(powf(2,
   1906                       (1.0f
   1907                        - par)
   1908                       * 7)
   1909                  - i) > i / 2 + 1 ? 0.0f : 1.0f) * par2 + (1.0f - par2);
   1910 }
   1911 
   1912 FILTER(bs2)
   1913 {
   1914     return (fabsf(powf(2,
   1915                       (1.0f
   1916                        - par)
   1917                       * 7)
   1918                  - i) < i / 2 + 1 ? 0.0f : 1.0f) * par2 + (1.0f - par2);
   1919 }
   1920 
   1921 bool floatEq(float a, float b)
   1922 {
   1923     const float fudge = .01f;
   1924     return a + fudge > b && a - fudge < b;
   1925 }
   1926 
   1927 FILTER(cos)
   1928 {
   1929     float tmp = powf(5.0f, par2 * 2.0f - 1.0f);
   1930     tmp = powf(i / 32.0f, tmp) * 32.0f;
   1931     if(floatEq(par2 * 127.0f, 64.0f))
   1932         tmp = i;
   1933     float gain = cosf(par * par * PI / 2.0f * tmp);
   1934     gain *= gain;
   1935     return gain;
   1936 }
   1937 
   1938 FILTER(sin)
   1939 {
   1940     float tmp = powf(5.0f, par2 * 2.0f - 1.0f);
   1941     tmp = powf(i / 32.0f, tmp) * 32.0f;
   1942     if(floatEq(par2 * 127.0f, 64.0f))
   1943         tmp = i;
   1944     float gain = sinf(par * par * PI / 2.0f * tmp);
   1945     gain *= gain;
   1946     return gain;
   1947 }
   1948 
   1949 FILTER(low_shelf)
   1950 {
   1951     float p2 = 1.0f - par + 0.2f;
   1952     float x  = i / (64.0f * p2 * p2);
   1953     if(x < 0.0f)
   1954         x = 0.0f;
   1955     else
   1956     if(x > 1.0f)
   1957         x = 1.0f;
   1958     float tmp = powf(1.0f - par2, 2.0f);
   1959     return cosf(x * PI) * (1.0f - tmp) + 1.01f + tmp;
   1960 }
   1961 
   1962 FILTER(s)
   1963 {
   1964     unsigned int tmp = (int) (powf(2.0f, (1.0f - par) * 7.2f));
   1965     float gain = 1.0f;
   1966     if(i == tmp)
   1967         gain = powf(2.0f, par2 * par2 * 8.0f);
   1968     return gain;
   1969 }
   1970 
   1971 FILTER(lpsk)
   1972 {
   1973     float tmp2PIf = 2.0f * PI * (1.05f-par)*64.0f;
   1974     std::complex<float> s  (0.0f,2.0f*PI*i);
   1975     float vOut = tmp2PIf * tmp2PIf;
   1976     std::complex<float> vIn = s*s + tmp2PIf*s/((par2)+(2.0f*par*par2)+0.5f) + tmp2PIf*tmp2PIf;
   1977     return std::abs((vOut*vOut*vOut) / (vIn*vIn*vIn));
   1978 
   1979 }
   1980 #undef FILTER
   1981 
   1982 filter_func_t *getFilter(unsigned char func)
   1983 {
   1984     static filter_func_t * const functions[] = {
   1985         osc_lp,
   1986         osc_hp1,
   1987         osc_hp1b,
   1988         osc_bp1,
   1989         osc_bs1,
   1990         osc_lp2,
   1991         osc_hp2,
   1992         osc_bp2,
   1993         osc_bs2,
   1994         osc_cos,
   1995         osc_sin,
   1996         osc_low_shelf,
   1997         osc_s,
   1998         osc_lpsk
   1999     };
   2000 
   2001     if(!func)
   2002         return NULL;
   2003 
   2004     func--;
   2005     assert(func < (sizeof(functions)/sizeof(functions[0])));
   2006     return functions[func];
   2007 }
   2008 
   2009 }