src/calf/osc.h

   1 /* Calf DSP Library
   2  * Oscillators
   3  *
   4  * Copyright (C) 2001-2007 Krzysztof Foltman
   5  *
   6  * This program is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2 of the License, or (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General
  17  * Public License along with this program; if not, write to the
  18  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  19  * Boston, MA 02111-1307, USA.
  20  */
  21
  22 #ifndef CALF_OSC_H
  23 #define CALF_OSC_H
  24
  25 #include "fft.h"
  26 #include <map>
  27
  28 namespace dsp
  29 {
  30
  31 /** Very simple, non-bandlimited saw oscillator. Should not be used for anything
  32  *  else than testing/prototyping. Unless get() function is replaced with something
  33  * with "proper" oscillator code, as the frequency setting function is fine.
  34  */
  35 struct simple_oscillator
  36 {
  37     /// Phase (from 0 to 0xFFFFFFFF)
  38     uint32_t phase;
  39     /// Per-sample phase delta (phase increment), equal to 2^32*freq/sr.
  40     uint32_t phasedelta;
  41     /// Reset oscillator phase to zero.
  42     void reset()
  43     {
  44         phase = 0;
  45     }
  46     /// Set phase delta based on oscillator frequency and sample rate.
  47     void set_freq(float freq, float sr)
  48     {
  49         phasedelta = (int)(freq * 65536.0 * 256.0 * 16.0 / sr) << 4;
  50     }
  51     /// Set phase delta based on oscillator frequency and inverse of sample rate.
  52     void set_freq_odsr(float freq, double odsr)
  53     {
  54         phasedelta = (int)(freq * 65536.0 * 256.0 * 16.0 * odsr) << 4;
  55     }
  56     /// Make one phase increment
  57     inline void step()
  58     {
  59         phase += phasedelta;
  60     }
  61     /// Make one phase increment and return a value from -0.5 to 0.5
  62     inline float get()
  63     {
  64         float value = (phase >> 16 ) / 65535.0 - 0.5;
  65         step();
  66         return value;
  67     }
  68 };
  69
  70 /**
  71  * FFT-based bandlimiting helper class. Allows conversion between time and frequency domains and generating brickwall filtered
  72  * versions of a waveform given a pre-computed spectrum.
  73  * Waveform size must be a power of two, and template argument SIZE_BITS is log2 of waveform size.
  74  */
  75 template<int SIZE_BITS>
  76 struct bandlimiter
  77 {
  78     enum { SIZE = 1 << SIZE_BITS };
  79     static dsp::fft<float, SIZE_BITS> &get_fft()
  80     {
  81         static dsp::fft<float, SIZE_BITS> fft;
  82         return fft;
  83     }
  84
  85     std::complex<float> spectrum[SIZE];
  86
  87     /// Import time domain waveform and calculate spectrum from it
  88     void compute_spectrum(float input[SIZE])
  89     {
  90         dsp::fft<float, SIZE_BITS> &fft = get_fft();
  91         std::complex<float> *data = new std::complex<float>[SIZE];
  92         for (int i = 0; i < SIZE; i++)
  93             data[i] = input[i];
  94         fft.calculate(data, spectrum, false);
  95         delete []data;
  96     }
  97
  98     /// Generate the waveform from the contained spectrum.
  99     void compute_waveform(float output[SIZE])
 100     {
 101         dsp::fft<float, SIZE_BITS> &fft = get_fft();
 102         std::complex<float> *data = new std::complex<float>[SIZE];
 103         fft.calculate(spectrum, data, true);
 104         for (int i = 0; i < SIZE; i++)
 105             output[i] = data[i].real();
 106         delete []data;
 107     }
 108
 109     /// remove DC offset of the spectrum (it usually does more harm than good!)
 110     void remove_dc()
 111     {
 112         spectrum[0] = 0.f;
 113     }
 114
 115     /// Very basic bandlimiting (brickwall filter)
 116     /// might need to be improved much in future!
 117     void make_waveform(float output[SIZE], int cutoff, bool foldover = false)
 118     {
 119         dsp::fft<float, SIZE_BITS> &fft = get_fft();
 120         std::vector<std::complex<float> > new_spec, iffted;
 121         new_spec.resize(SIZE);
 122         iffted.resize(SIZE);
 123         // Copy original harmonics up to cutoff point
 124         new_spec[0] = spectrum[0];
 125         for (int i = 1; i < cutoff; i++)
 126             new_spec[i] = spectrum[i],
 127             new_spec[SIZE - i] = spectrum[SIZE - i];
 128         // Fill the rest with zeros, optionally folding over harmonics over the
 129         // cutoff point into the lower octaves while halving the amplitude.
 130         // (I think it is almost nice for bell type waveforms when the original
 131         // waveform has few widely spread harmonics)
 132         if (foldover)
 133         {
 134             std::complex<float> fatt(0.5);
 135             cutoff /= 2;
 136             if (cutoff < 2)
 137                 cutoff = 2;
 138             for (int i = SIZE / 2; i >= cutoff; i--)
 139             {
 140                 new_spec[i / 2] += new_spec[i] * fatt;
 141                 new_spec[SIZE - i / 2] += new_spec[SIZE - i] * fatt;
 142                 new_spec[i] = 0.f,
 143                 new_spec[SIZE - i] = 0.f;
 144             }
 145         }
 146         else
 147         {
 148             if (cutoff < 1)
 149                 cutoff = 1;
 150             for (int i = cutoff; i < SIZE / 2; i++)
 151                 new_spec[i] = 0.f,
 152                 new_spec[SIZE - i] = 0.f;
 153         }
 154         // convert back to time domain (IFFT) and extract only real part
 155         fft.calculate(&new_spec.front(), &iffted.front(), true);
 156         for (int i = 0; i < SIZE; i++)
 157             output[i] = iffted[i].real();
 158     }
 159 };
 160
 161 /// Set of bandlimited wavetables
 162 template<int SIZE_BITS>
 163 struct waveform_family: public std::map<uint32_t, float *>
 164 {
 165     enum { SIZE = 1 << SIZE_BITS };
 166     using std::map<uint32_t, float *>::iterator;
 167     using std::map<uint32_t, float *>::end;
 168     using std::map<uint32_t, float *>::lower_bound;
 169     float original[SIZE];
 170
 171     /// Fill the family using specified bandlimiter and original waveform. Optionally apply foldover.
 172     /// Does not produce harmonics over specified limit (limit = (SIZE / 2) / min_number_of_harmonics)
 173     void make(bandlimiter<SIZE_BITS> &bl, float input[SIZE], bool foldover = false, uint32_t limit = SIZE / 2)
 174     {
 175         memcpy(original, input, sizeof(original));
 176         bl.compute_spectrum(input);
 177         make_from_spectrum(bl, foldover);
 178     }
 179
 180     /// Fill the family using specified bandlimiter and spectrum contained within. Optionally apply foldover.
 181     /// Does not produce harmonics over specified limit (limit = (SIZE / 2) / min_number_of_harmonics)
 182     void make_from_spectrum(bandlimiter<SIZE_BITS> &bl, bool foldover = false, uint32_t limit = SIZE / 2)
 183     {
 184         bl.remove_dc();
 185
 186         uint32_t base = 1 << (32 - SIZE_BITS);
 187         uint32_t cutoff = SIZE / 2, top = SIZE / 2;
 188         float vmax = 0;
 189         for (unsigned int i = 0; i < cutoff; i++)
 190             vmax = std::max(vmax, abs(bl.spectrum[i]));
 191         float vthres = vmax / 1024.0;  // -60dB
 192         float cumul = 0.f;
 193         while(cutoff > (SIZE / limit)) {
 194             if (!foldover)
 195             {
 196                 // skip harmonics too quiet to be heard, but measure their loudness cumulatively,
 197                 // because even if a single harmonic is too quiet, a whole bunch of them may add up
 198                 // to considerable amount of space
 199                 cumul = 0.f;
 200                 while(cutoff > 1 && cumul + abs(bl.spectrum[cutoff - 1]) < vthres)
 201                 {
 202                     cumul += abs(bl.spectrum[cutoff - 1]);
 203                     cutoff--;
 204                 }
 205             }
 206             float *wf = new float[SIZE+1];
 207             bl.make_waveform(wf, cutoff, foldover);
 208             wf[SIZE] = wf[0];
 209             (*this)[base * (top / cutoff)] = wf;
 210             cutoff = (int)(0.75 * cutoff);
 211         }
 212     }
 213
 214     /// Retrieve waveform pointer suitable for specified phase_delta
 215     inline float *get_level(uint32_t phase_delta)
 216     {
 217         iterator i = upper_bound(phase_delta);
 218         if (i == end())
 219             return NULL;
 220         // printf("Level = %08x\n", i->first);
 221         return i->second;
 222     }
 223     /// Destructor, deletes the waveforms and removes them from the map.
 224     ~waveform_family()
 225     {
 226         for (iterator i = begin(); i != end(); i++)
 227             delete []i->second;
 228         clear();
 229     }
 230 };
 231
 232 #if 0
 233 // cubic interpolation
 234 static inline float cerp(float pm1, float p0, float p1, float p2, float t)
 235 {
 236     return (-t*(t-1)*(t-2) * pm1 + 3*(t+1)*(t-1)*(t-2) * p0 - 3*(t+1)*t*(t-2) * p1 + (t+1)*t*(t-1) * p2) * (1.0 / 6.0);
 237 }
 238 #endif
 239 /**
 240  * Simple table-based lerping oscillator. Uses waveform of size 2^SIZE_BITS.
 241  * Combine with waveform_family if bandlimited waveforms are needed. Because
 242  * of linear interpolation, it's usually a good idea to use large tables
 243  * (2048-4096 points), otherwise aliasing may be produced.
 244  */
 245 template<int SIZE_BITS>
 246 struct waveform_oscillator: public simple_oscillator
 247 {
 248     enum { SIZE = 1 << SIZE_BITS, MASK = SIZE - 1, SCALE = 1 << (32 - SIZE_BITS) };
 249     float *waveform;
 250     waveform_oscillator()
 251     {
 252         waveform = NULL;
 253     }
 254
 255     /// Get the value from single oscillator at current position
 256     inline float get()
 257     {
 258         uint32_t wpos = phase >> (32 - SIZE_BITS);
 259         return dsp::lerp(waveform[wpos], waveform[(wpos + 1) & MASK], (phase & (SCALE - 1)) * (1.0f / SCALE));
 260     }
 261     /// Add/substract two phase-shifted values
 262     inline float get_phaseshifted(uint32_t shift, float mix)
 263     {
 264         uint32_t wpos = phase >> (32 - SIZE_BITS);
 265         float value1 = dsp::lerp(waveform[wpos], waveform[(wpos + 1) & MASK], (phase & (SCALE - 1)) * (1.0f / SCALE));
 266         wpos = (phase + shift) >> (32 - SIZE_BITS);
 267         float value2 = dsp::lerp(waveform[wpos], waveform[(wpos + 1) & MASK], ((phase + shift) & (SCALE - 1)) * (1.0f / SCALE));
 268         return value1 + mix * value2;
 269     }
 270     /// Add/substract two phase-shifted values
 271     inline float get_phaseshifted2(uint32_t shift, int32_t gshift, float mix)
 272     {
 273         uint32_t wpos = (phase + gshift) >> (32 - SIZE_BITS);
 274         float value1 = dsp::lerp(waveform[wpos], waveform[(wpos + 1) & MASK], (phase & (SCALE - 1)) * (1.0f / SCALE));
 275         wpos = (phase + gshift + shift) >> (32 - SIZE_BITS);
 276         float value2 = dsp::lerp(waveform[wpos], waveform[(wpos + 1) & MASK], ((phase + shift) & (SCALE - 1)) * (1.0f / SCALE));
 277         return value1 + mix * value2;
 278     }
 279     /// Get the value of a hard synced osc (65536 = 1:1 ratio)
 280     inline float get_phasedist(uint32_t sync, uint32_t shift, float mix)
 281     {
 282         uint32_t phase_mod = (uint64_t(phase) * sync >> 16);
 283
 284         uint32_t wpos = phase_mod >> (32 - SIZE_BITS);
 285         float value1 = dsp::lerp(waveform[wpos], waveform[(wpos + 1) & MASK], (phase & (SCALE - 1)) * (1.0f / SCALE));
 286         wpos = (phase_mod + shift) >> (32 - SIZE_BITS);
 287         float value2 = dsp::lerp(waveform[wpos], waveform[(wpos + 1) & MASK], ((phase + shift) & (SCALE - 1)) * (1.0f / SCALE));
 288         return value1 + mix * value2;
 289     }
 290     /// One step
 291     inline void advance()
 292     {
 293         phase += phasedelta;
 294     }
 295 };
 296
 297 /**
 298  * Simple triangle LFO without any smoothing or anything of this sort.
 299  */
 300 struct triangle_lfo: public simple_oscillator
 301 {
 302     /// Previous value (not stored here, but may be used by calling code)
 303     float last;
 304
 305     triangle_lfo()
 306     {
 307         reset();
 308     }
 309     void reset()
 310     {
 311         simple_oscillator::reset();
 312         last = 0;
 313     }
 314     inline float get()
 315     {
 316         uint32_t phase2 = phase;
 317         // start at 90 degrees point of the "/\" wave (-1 to +1)
 318         phase2 += 1<<30;
 319         // if in second half, invert the wave (so it falls back into 0..0x7FFFFFFF)
 320         phase2 ^= ((int32_t)phase2)>>31;
 321
 322         float value = (phase2 >> 6) / 16777216.0 - 1.0;
 323         phase += phasedelta;
 324         return value;
 325     }
 326 };
 327
 328 /// Simple stupid inline function to normalize a waveform (by removing DC offset and ensuring max absolute value of 1).
 329 static inline void normalize_waveform(float *table, unsigned int size)
 330 {
 331     float dc = 0;
 332     for (unsigned int i = 0; i < size; i++)
 333         dc += table[i];
 334     dc /= size;
 335     for (unsigned int i = 0; i < size; i++)
 336         table[i] -= dc;
 337     float thismax = 0;
 338     for (unsigned int i = 0; i < size; i++)
 339         thismax = std::max(thismax, fabsf(table[i]));
 340     if (thismax < 0.000001f)
 341         return;
 342     double divv = 1.0 / thismax;
 343     for (unsigned int i = 0; i < size; i++)
 344         table[i] *= divv;
 345 }
 346
 347
 348
 349 };
 350
 351 #endif