1 /** 2 * Copyright: Copyright Auburn Sounds 2015 and later. 3 * License: $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0) 4 * Authors: Guillaume Piolat 5 */ 6 module dplug.dsp.wavetable; 7 8 import std.math; 9 10 //import dplug.core.nogc; 11 import dplug.core.math; 12 import dplug.core.alignedbuffer; 13 14 /// Generate a sine. 15 /// It turns out it's very stable, stable enough for table generation at least. 16 /// FUTURE: resync method 17 struct SineGenerator(T) 18 { 19 void initialize(T initPhase, T frequency, T samplerate) nothrow @nogc 20 { 21 T w = frequency * 2 * PI / samplerate; 22 _b1 = 2 * cos(w); 23 _y1 = sin(initPhase - w); 24 _y2 = sin(initPhase - 2 * w); 25 } 26 27 T next() nothrow @nogc 28 { 29 T y0 = _b1 * _y1 - _y2; 30 _y2 = _y1; 31 _y1 = y0; 32 return y0; 33 } 34 35 T _y1, _y2; 36 T _b1; 37 } 38 39 enum WaveformType 40 { 41 SINE, 42 SAWTOOTH, 43 SQUARE, 44 TRIANGLE 45 } 46 47 // wavetable with mip-maps 48 49 /// Generates anti-aliased waveform generation through 50 /// procedurally generated mipmapped tables. 51 /// FUTURE: only integer phase 52 struct Wavetable 53 { 54 nothrow: 55 @nogc: 56 57 void initialize(int largestSize, WaveformType waveform) 58 { 59 resize(largestSize); 60 generate(waveform); // regenerate tables 61 } 62 63 ~this() 64 { 65 _mipmapData.reallocBuffer(0); 66 _wholeBuffer.reallocBuffer(0); 67 } 68 69 @disable this(this); 70 71 float lookupLinear(uint phaseIntPart, float phaseFractional, int level) 72 { 73 float* mipmap0 = mipmapData(level); 74 int mask = sizeOfMipmap(level) - 1; 75 float a = mipmap0[ phaseIntPart & mask ]; 76 float b = mipmap0[ ( 1 + phaseIntPart ) & mask ]; 77 return a * (1 - phaseFractional) + phaseFractional * b; 78 } 79 80 float lookupCatmullRom(uint phaseIntPart, float phaseFractional, int level) 81 { 82 float* mipmap0 = mipmapData(level); 83 int mask = sizeOfMipmap(level) - 1; 84 float p0 = mipmap0[ (phaseIntPart - 1) & mask ]; 85 float p1 = mipmap0[ phaseIntPart & mask ]; 86 float p2 = mipmap0[ ( 1 + phaseIntPart ) & mask ]; 87 float p3 = mipmap0[ ( 2 + phaseIntPart ) & mask ]; 88 float t = phaseFractional; 89 90 return 0.5f * ((2 * p1) 91 + (-p0 + p2) * t 92 + (2 * p0 - 5 * p1 + 4 * p2 - p3) * t * t 93 + (-p0 + 3 * p1 - 3 * p2 + p3) * t * t * t); 94 } 95 96 float lookupCatmullRomMipmap(uint phaseIntPart, float phaseFractional, float phaseIncrementSamples) 97 { 98 99 float level = cast(float)log2(phaseIncrementSamples); 100 int level0 = cast(int)floor(level); 101 int level1 = level0 + 1; 102 if (level0 < 0) 103 level0 = 0; 104 if (level1 < 0) 105 level1 = 0; 106 int maxLevel = cast(int)_numTables - 1; 107 if (level0 > maxLevel) 108 level0 = maxLevel; 109 if (level1 > maxLevel) 110 level1 = maxLevel; 111 112 if (level1 == 0) 113 { 114 return lookupCatmullRom(phaseIntPart, phaseFractional, 0); 115 } 116 else 117 { 118 float fractionalLevel = level - level0; 119 120 float phaseL0 = (phaseIntPart + phaseFractional) / cast(float)(1 << level0); 121 int iPart0 = cast(int)(phaseL0); 122 float fractional0 = phaseL0 - iPart0; 123 124 float phaseL1 = (phaseIntPart + phaseFractional) / cast(float)(1 << level1); 125 int iPart1 = cast(int)(phaseL1); 126 float fractional1 = phaseL1 - iPart1; 127 128 float L0 = lookupCatmullRom(iPart0, fractional0, level0); 129 float L1 = lookupCatmullRom(iPart1, fractional1, level1); 130 return lerp(L0, L1, fractionalLevel); 131 } 132 } 133 134 // mimaps levels range from 0 to numMipmaps() - 1 135 int numMipmaps() const 136 { 137 return _numTables; 138 } 139 140 int sizeOfMipmap(int level) const 141 { 142 return _largestSize >> level; 143 } 144 145 float* mipmapData(int level) 146 { 147 return _mipmapData[level]; 148 } 149 150 private: 151 152 int _largestSize; 153 int _mast; 154 int _numTables; 155 156 float*[] _mipmapData; 157 float[] _wholeBuffer; 158 159 /// Defines the harmonic rolloff, critical to avoid aliasing around the last mipmap 160 /// This is very arbitrary and ultimately power of two mipmaps are maybe not sufficient 161 /// to have less aliasing. 162 static double rolloffHarmonic(double normalizedFrequency) pure nothrow @nogc // between 0 and 1 163 { 164 double cosF0 = cos(normalizedFrequency * PI); 165 return cosF0 * cosF0; 166 } 167 168 void resize(int largestSize) nothrow @nogc 169 { 170 assert(isPowerOfTwo(largestSize)); 171 172 _largestSize = largestSize; 173 // compute size for all mipmaps 174 int sizeNeeded = 0; 175 _numTables = 0; 176 int sizeOfTable = largestSize; 177 while (sizeOfTable > 0) 178 { 179 sizeNeeded += sizeOfTable; 180 sizeOfTable /= 2; 181 _numTables += 1; 182 } 183 184 _wholeBuffer.reallocBuffer(sizeNeeded); 185 186 // fill table pointers 187 { 188 _mipmapData.reallocBuffer(_numTables); 189 int cumulated = 0; 190 for (int level = 0; level < _numTables; ++level) 191 { 192 _mipmapData[level] = &_wholeBuffer[cumulated]; 193 cumulated += _largestSize >> level; 194 } 195 } 196 } 197 198 // fill all table with waveform 199 void generate(WaveformType waveform) nothrow @nogc 200 { 201 for (int level = 0; level < _numTables; ++level) 202 { 203 int size = sizeOfMipmap(level); 204 float* data = mipmapData(level); 205 206 for (int t = 0; t < size; ++t) 207 { 208 data[t] = 0; 209 } 210 211 int numHarmonics = size / 2; 212 213 if (size > 2) 214 { 215 for (int h = 0; h < numHarmonics; ++h) 216 { 217 double normalizedFrequency = (1 + h) / cast(double)(numHarmonics - 1); 218 double amplitude = getWaveHarmonicAmplitude(waveform, h + 1) * rolloffHarmonic(normalizedFrequency); 219 220 for (int t = 0; t < size; ++t) 221 { 222 double x = sin( cast(double)t * (2 * PI) * (h + 1) / cast(double)size ) * amplitude; 223 data[t] += cast(float)x; 224 } 225 } 226 } 227 228 for (int t = 0; t < size; ++t) 229 { 230 assert(isFinite(data[t])); 231 } 232 } 233 } 234 } 235 236 237 struct WavetableOsc 238 { 239 public: 240 void initialize(Wavetable* wavetable, float sampleRate) nothrow @nogc 241 { 242 _wavetable = wavetable; 243 _phaseIntPart = 0; 244 _phaseFractional = 0; 245 _phaseFactor = _wavetable.sizeOfMipmap(0) / sampleRate; 246 } 247 248 /// Allows dirty resync 249 void resetPhase() nothrow @nogc 250 { 251 _phaseIntPart = 0; 252 _phaseFractional = 0; 253 } 254 255 deprecated alias next = nextSample; 256 257 /// Get next sample. 258 float nextSample(float frequencyHz) nothrow @nogc 259 { 260 float phaseIncrementSamples = frequencyHz * _phaseFactor; 261 assert(phaseIncrementSamples >= 0); 262 int iPart = cast(int)(phaseIncrementSamples); 263 _phaseIntPart += iPart; 264 265 _phaseFractional += (phaseIncrementSamples - iPart); 266 if (_phaseFractional >= 1) 267 { 268 _phaseFractional -= 1.0; 269 _phaseIntPart += 1; 270 271 // in case the assert above would fail (which happen) 272 if (_phaseFractional >= 1) 273 { 274 _phaseFractional -= 1.0; 275 _phaseIntPart += 1; 276 } 277 } 278 279 assert(_phaseFractional >= 0 && _phaseFractional <= 1); 280 281 return _wavetable.lookupCatmullRomMipmap(_phaseIntPart, cast(float)_phaseFractional, phaseIncrementSamples); 282 // return _wavetable->lookupLinear(_phaseIntPart, _phaseFractional, 1); 283 } 284 285 private: 286 uint _phaseIntPart; 287 float _phaseFractional; 288 Wavetable* _wavetable; 289 float _phaseFactor; 290 } 291 292 293 private: 294 295 double getWaveHarmonicAmplitude(WaveformType waveform, int n) nothrow @nogc 296 { 297 assert(n > 0); 298 final switch(waveform) 299 { 300 case WaveformType.SINE: 301 { 302 if (n == 1) 303 return 1; 304 else 305 return 0; 306 } 307 308 case WaveformType.SAWTOOTH: 309 { 310 return 1 / cast(double)n; 311 } 312 313 case WaveformType.SQUARE: 314 { 315 if (n % 2 == 0) 316 return 0; 317 else 318 return 1 / cast(double)n; 319 } 320 321 case WaveformType.TRIANGLE: 322 { 323 if (n % 2 == 0) 324 return 0; 325 else if (n % 4 == 1) 326 return 1 / cast(double)(n*n); 327 else 328 { 329 assert(n % 4 == 3); 330 return -1 / cast(double)(n*n); 331 } 332 } 333 } 334 }