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., 59 Temple Place, Suite 330,
  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     inline float get()
  57     {
  58         float value = (phase >> 16 ) / 65535.0 - 0.5;
  59         phase += phasedelta;
  60         return value;
  61     }
  62 };
  63
  64 /**
  65  * FFT-based bandlimiting helper class. Allows conversion between time and frequency domains and generating brickwall filtered
  66  * versions of a waveform given a pre-computed spectrum.
  67  * Waveform size must be a power of two, and template argument SIZE_BITS is log2 of waveform size.
  68  */
  69 template<int SIZE_BITS>
  70 struct bandlimiter
  71 {
  72     enum { SIZE = 1 << SIZE_BITS };
  73     static dsp::fft<float, SIZE_BITS> &get_fft()
  74     {
  75         static dsp::fft<float, SIZE_BITS> fft;
  76         return fft;
  77     }
  78
  79     std::complex<float> spectrum[SIZE];
  80
  81     /// Import time domain waveform and calculate spectrum from it
  82     void compute_spectrum(float input[SIZE])
  83     {
  84         dsp::fft<float, SIZE_BITS> &fft = get_fft();
  85         std::complex<float> *data = new std::complex<float>[SIZE];
  86         for (int i = 0; i < SIZE; i++)
  87             data[i] = input[i];
  88         fft.calculate(data, spectrum, false);
  89         delete []data;
  90     }
  91
  92     /// Generate the waveform from the contained spectrum.
  93     void compute_waveform(float output[SIZE])
  94     {
  95         dsp::fft<float, SIZE_BITS> &fft = get_fft();
  96         std::complex<float> *data = new std::complex<float>[SIZE];
  97         fft.calculate(spectrum, data, true);
  98         for (int i = 0; i < SIZE; i++)
  99             output[i] = data[i].real();
 100         delete []data;
 101     }
 102
 103     /// remove DC offset of the spectrum (it usually does more harm than good!)
 104     void remove_dc()
 105     {
 106         spectrum[0] = 0.f;
 107     }
 108
 109     /// Very basic bandlimiting (brickwall filter)
 110     /// might need to be improved much in future!
 111     void make_waveform(float output[SIZE], int cutoff, bool foldover = false)
 112     {
 113         dsp::fft<float, SIZE_BITS> &fft = get_fft();
 114         std::vector<std::complex<float> > new_spec, iffted;
 115         new_spec.resize(SIZE);
 116         iffted.resize(SIZE);
 117         // Copy original harmonics up to cutoff point
 118         new_spec[0] = spectrum[0];
 119         for (int i = 1; i < cutoff; i++)
 120             new_spec[i] = spectrum[i],
 121             new_spec[SIZE - i] = spectrum[SIZE - i];
 122         // Fill the rest with zeros, optionally folding over harmonics over the
 123         // cutoff point into the lower octaves while halving the amplitude.
 124         // (I think it is almost nice for bell type waveforms when the original
 125         // waveform has few widely spread harmonics)
 126         if (foldover)
 127         {
 128             std::complex<float> fatt(0.5);
 129             cutoff /= 2;
 130             if (cutoff < 2)
 131                 cutoff = 2;
 132             for (int i = SIZE / 2; i >= cutoff; i--)
 133             {
 134                 new_spec[i / 2] += new_spec[i] * fatt;
 135                 new_spec[SIZE - i / 2] += new_spec[SIZE - i] * fatt;
 136                 new_spec[i] = 0.f,
 137                 new_spec[SIZE - i] = 0.f;
 138             }
 139         }
 140         else
 141         {
 142             if (cutoff < 1)
 143                 cutoff = 1;
 144             for (int i = cutoff; i < SIZE / 2; i++)
 145                 new_spec[i] = 0.f,
 146                 new_spec[SIZE - i] = 0.f;
 147         }
 148         // convert back to time domain (IFFT) and extract only real part
 149         fft.calculate(new_spec.data(), iffted.data(), true);
 150         for (int i = 0; i < SIZE; i++)
 151             output[i] = iffted[i].real();
 152     }
 153 };
 154
 155 /// Set of bandlimited wavetables
 156 template<int SIZE_BITS>
 157 struct waveform_family: public std::map<uint32_t, float *>
 158 {
 159     enum { SIZE = 1 << SIZE_BITS };
 160     using std::map<uint32_t, float *>::iterator;
 161     using std::map<uint32_t, float *>::end;
 162     using std::map<uint32_t, float *>::lower_bound;
 163     float original[SIZE];
 164
 165     /// Fill the family using specified bandlimiter and original waveform. Optionally apply foldover.
 166     /// Does not produce harmonics over specified limit (limit = (SIZE / 2) / min_number_of_harmonics)
 167     void make(bandlimiter<SIZE_BITS> &bl, float input[SIZE], bool foldover = false, uint32_t limit = SIZE / 2)
 168     {
 169         memcpy(original, input, sizeof(original));
 170         bl.compute_spectrum(input);
 171         make_from_spectrum(bl, foldover);
 172     }
 173
 174     /// Fill the family using specified bandlimiter and spectrum contained within. Optionally apply foldover.
 175     /// Does not produce harmonics over specified limit (limit = (SIZE / 2) / min_number_of_harmonics)
 176     void make_from_spectrum(bandlimiter<SIZE_BITS> &bl, bool foldover = false, uint32_t limit = SIZE / 2)
 177     {
 178         bl.remove_dc();
 179
 180         uint32_t base = 1 << (32 - SIZE_BITS);
 181         uint32_t cutoff = SIZE / 2, top = SIZE / 2;
 182         float vmax = 0;
 183         for (unsigned int i = 0; i < cutoff; i++)
 184             vmax = std::max(vmax, abs(bl.spectrum[i]));
 185         float vthres = vmax / 1024.0;  // -60dB
 186         while(cutoff > (SIZE / limit)) {
 187             if (!foldover)
 188             {
 189                 while(cutoff > 1 && abs(bl.spectrum[cutoff - 1]) < vthres)
 190                     cutoff--;
 191             }
 192             float *wf = new float[SIZE+1];
 193             bl.make_waveform(wf, cutoff, foldover);
 194             wf[SIZE] = wf[0];
 195             (*this)[base * (top / cutoff)] = wf;
 196             cutoff = (int)(0.75 * cutoff);
 197         }
 198     }
 199
 200     /// Retrieve waveform pointer suitable for specified phase_delta
 201     inline float *get_level(uint32_t phase_delta)
 202     {
 203         iterator i = upper_bound(phase_delta);
 204         if (i == end())
 205             return NULL;
 206         // printf("Level = %08x\n", i->first);
 207         return i->second;
 208     }
 209     /// Destructor, deletes the waveforms and removes them from the map.
 210     ~waveform_family()
 211     {
 212         for (iterator i = begin(); i != end(); i++)
 213             delete []i->second;
 214         clear();
 215     }
 216 };
 217
 218 /**
 219  * Simple table-based lerping oscillator. Uses waveform of size 2^SIZE_BITS.
 220  * Combine with waveform_family if bandlimited waveforms are needed. Because
 221  * of linear interpolation, it's usually a good idea to use large tables
 222  * (2048-4096 points), otherwise aliasing may be produced.
 223  */
 224 template<int SIZE_BITS>
 225 struct waveform_oscillator: public simple_oscillator
 226 {
 227     enum { SIZE = 1 << SIZE_BITS, MASK = SIZE - 1, SCALE = 1 << (32 - SIZE_BITS) };
 228     float *waveform;
 229     inline float get()
 230     {
 231         uint32_t wpos = phase >> (32 - SIZE_BITS);
 232         float value = dsp::lerp(waveform[wpos], waveform[(wpos + 1) & MASK], (phase & (SCALE - 1)) * (1.0f / SCALE));
 233         phase += phasedelta;
 234         return value;
 235     }
 236 };
 237
 238 /// Simple stupid inline function to normalize a waveform (by removing DC offset and ensuring max absolute value of 1).
 239 static inline void normalize_waveform(float *table, unsigned int size)
 240 {
 241     float dc = 0;
 242     for (unsigned int i = 0; i < size; i++)
 243         dc += table[i];
 244     dc /= size;
 245     for (unsigned int i = 0; i < size; i++)
 246         table[i] -= dc;
 247     float thismax = 0;
 248     for (unsigned int i = 0; i < size; i++)
 249         thismax = std::max(thismax, fabsf(table[i]));
 250     if (thismax < 0.000001f)
 251         return;
 252     double divv = 1.0 / thismax;
 253     for (unsigned int i = 0; i < size; i++)
 254         table[i] *= divv;
 255 }
 256
 257
 258
 259 };
 260
 261 #endif