WaveShapeSmps.cpp (11305B)
1 /* 2 ZynAddSubFX - a software synthesizer 3 4 WaveShapeSmps.cpp - Sample Distortion 5 Copyright (C) 2002-2005 Nasca Octavian Paul 6 Author: Nasca Octavian Paul 7 8 This program is free software; you can redistribute it and/or 9 modify it under the terms of the GNU General Public License 10 as published by the Free Software Foundation; either version 2 11 of the License, or (at your option) any later version. 12 */ 13 14 #include "WaveShapeSmps.h" 15 #include <cmath> 16 17 namespace zyn { 18 19 float polyblampres(float smp, float ws, float dMax) 20 { 21 // Formula from: Esqueda, Välimäki, Bilbao (2015): ALIASING REDUCTION IN SOFT-CLIPPING ALGORITHMS 22 // http://dafx16.vutbr.cz/dafxpapers/18-DAFx-16_paper_33-PN.pdf pg 123, table 1 23 // Four-point polyBLAMP residual: 24 // [−2T, T] d^5/120 25 // [−T, 0] −d^5/40 + d^4/24 + d^3/12 + d^2/12 + d/24 + 1/120 26 // [0, T] d^5/40 − d^4/12 + d^2/3 − d/2 + 7/30 27 // [T, 2T] −d^5/120 + d^4/24 − d^3/12 + d^2/12 − d/24 + 1/120 28 if (dMax == 0) return 0.0; 29 float dist = fabsf(smp) - ws; 30 float res, d1, d2, d3, d4, d5; 31 if(fabsf(dist) < dMax) { 32 if(dist < -dMax/2.0f) { 33 d1 = (dist + dMax)/dMax*2.0f; // [-dMax ... -dMax/2] -> [0 ... 1] 34 res = d1 * d1 * d1 * d1 * d1 / 120.0f; 35 } 36 else if(dist < 0.0f) { 37 d1 = (dist + dMax/2.0f)/dMax*2.0f; // [-dMax/2 ... 0] -> [0 ... 1] 38 d2 = d1*d1; 39 d3 = d2*d1; 40 d4 = d3*d1; 41 d5 = d4*d1; 42 res = (-d5/40.0f) + (d4/24.0f) + (d3/12.0f) + (d2/12.0f) + (d1/24.0f) + (1.0f/120.0f); 43 } 44 else if(dist < dMax/2.0f) { 45 d1 = dist/dMax*2.0f; //[0 ... dMax/2] -> [0 ... 1] 46 d2 = d1*d1; 47 d4 = d2*d2; 48 d5 = d4*d1; 49 res = (d5/40.0f) - (d4/12.0f) + (d2/3.0f) - (d1/2.0f) + (7.0f/30.0f); 50 } 51 else { 52 d1 = (dist - dMax/2.0f)/dMax*2.0f; //[dMax/2 ... dMax] -> [0 ... 1] 53 d2 = d1*d1; 54 d3 = d2*d1; 55 d4 = d3*d1; 56 d5 = d4*d1; 57 res = (-d5/120.0f) + (d4/24.0f) - (d3/12.0f) + (d2/12.0f) - (d1/24.0f) + (1.0f/120.0f); 58 } 59 } 60 else 61 res = 0.0f; 62 63 return res*dMax/2.0f; 64 } 65 66 void waveShapeSmps(int n, 67 float *smps, 68 unsigned char type, 69 unsigned char drive, 70 unsigned char offset, 71 unsigned char funcpar) 72 { 73 int i; 74 float ws = drive / 127.0f; 75 float par = funcpar / 127.0f; 76 float offs = (offset - 64.0f) / 64.0f; 77 float tmpv; 78 79 switch(type) { 80 case 1: 81 ws = powf(10, ws * ws * 3.0f) - 1.0f + 0.001f; //Arctangent 82 for(i = 0; i < n; ++i) { 83 smps[i] += offs; 84 smps[i] = atanf(smps[i] * ws) / atanf(ws); 85 smps[i] -= offs; 86 } 87 break; 88 case 2: 89 ws = ws * ws * 32.0f + 0.0001f; //Asymmetric 90 if(ws < 1.0f) 91 tmpv = sinf(ws) + 0.1f; 92 else 93 tmpv = 1.1f; 94 for(i = 0; i < n; ++i) 95 smps[i] = sinf(smps[i] * (0.1f + ws - ws * smps[i])) / tmpv; 96 ; 97 break; 98 case 3: 99 ws = ws * ws * ws * 20.0f + 0.0001f; //Pow 100 for(i = 0; i < n; ++i) { 101 smps[i] *= ws; 102 if(fabsf(smps[i]) < 1.0f) { 103 smps[i] = (smps[i] - smps[i] * smps[i] * smps[i]) * 3.0f; 104 if(ws < 1.0f) 105 smps[i] /= ws; 106 } 107 else 108 smps[i] = 0.0f; 109 } 110 break; 111 case 4: 112 ws = ws * ws * ws * 32.0f + 0.0001f; //Sine 113 if(ws < 1.57f) 114 tmpv = sinf(ws); 115 else 116 tmpv = 1.0f; 117 for(i = 0; i < n; ++i) 118 smps[i] = sinf(smps[i] * ws) / tmpv; 119 break; 120 case 5: 121 ws = ws * ws + 0.000001f; //Quantisize 122 for(i = 0; i < n; ++i) 123 smps[i] = floor(smps[i] / ws + 0.5f) * ws; 124 break; 125 case 6: 126 ws = ws * ws * ws * 32 + 0.0001f; //Zigzag 127 if(ws < 1.0f) 128 tmpv = sinf(ws); 129 else 130 tmpv = 1.0f; 131 for(i = 0; i < n; ++i) 132 smps[i] = asinf(sinf(smps[i] * ws)) / tmpv; 133 break; 134 case 7: 135 ws = powf(2.0f, -ws * ws * 8.0f); //Limiter 136 par = par/4; 137 if (par > ws - 0.01) par = ws - 0.01; 138 for(i = 0; i < n; ++i) { 139 // add the offset: x = smps[i] + offs 140 smps[i] += offs; 141 float res = polyblampres(smps[i], ws, par); 142 // now apply the polyblamped limiter: y = f(x) 143 if (smps[i]>=0) 144 smps[i] = ( smps[i] > ws ? ws-res : smps[i]-res ); 145 else 146 smps[i] = ( smps[i] < -ws ? -ws+res : smps[i]+res ); 147 // and subtract the polyblamp-limited offset again: smps[i] = y - f(offs) 148 res = polyblampres(offs, ws, par); 149 if (offs>=0) 150 smps[i] -= ( offs >= ws ? ws-res : offs-res ); 151 else 152 smps[i] -= ( offs <= -ws ? -ws+res : offs+res ); 153 // divide through the drive factor: prevents limited signals to get low 154 smps[i] /= ws; 155 156 } 157 break; 158 case 8: 159 ws = powf(2.0f, -ws * ws * 8.0f); //Upper Limiter 160 for(i = 0; i < n; ++i) { 161 float tmp = smps[i]; 162 if(tmp > ws) 163 smps[i] = ws; 164 smps[i] *= 2.0f; 165 } 166 break; 167 case 9: 168 ws = powf(2.0f, -ws * ws * 8.0f); //Lower Limiter 169 for(i = 0; i < n; ++i) { 170 float tmp = smps[i]; 171 if(tmp < -ws) 172 smps[i] = -ws; 173 smps[i] *= 2.0f; 174 } 175 break; 176 case 10: 177 ws = (powf(2.0f, ws * 6.0f) - 1.0f) / powf(2.0f, 6.0f); //Inverse Limiter 178 if (par > ws - 0.01) par = ws - 0.01; 179 for(i = 0; i < n; ++i) { 180 smps[i] += offs; 181 float res = polyblampres(smps[i], ws, par); 182 if (smps[i]>=0) 183 smps[i] = ( smps[i] > ws ? smps[i]-ws+res : res ); 184 else 185 smps[i] = ( smps[i] < -ws ? smps[i]+ws-res : -res ); 186 smps[i] -= offs; 187 } 188 break; 189 case 11: 190 ws = powf(5, ws * ws * 1.0f) - 1.0f; //Clip 191 for(i = 0; i < n; ++i) 192 smps[i] = smps[i] 193 * (ws + 0.5f) * 0.9999f - floor( 194 0.5f + smps[i] * (ws + 0.5f) * 0.9999f); 195 break; 196 case 12: 197 ws = ws * ws * ws * 30 + 0.001f; //Asym2 198 if(ws < 0.3f) 199 tmpv = ws; 200 else 201 tmpv = 1.0f; 202 for(i = 0; i < n; ++i) { 203 float tmp = smps[i] * ws; 204 if((tmp > -2.0f) && (tmp < 1.0f)) 205 smps[i] = tmp * (1.0f - tmp) * (tmp + 2.0f) / tmpv; 206 else 207 smps[i] = 0.0f; 208 } 209 break; 210 case 13: 211 ws = ws * ws * ws * 32.0f + 0.0001f; //Pow2 212 if(ws < 1.0f) 213 tmpv = ws * (1 + ws) / 2.0f; 214 else 215 tmpv = 1.0f; 216 for(i = 0; i < n; ++i) { 217 float tmp = smps[i] * ws; 218 if((tmp > -1.0f) && (tmp < 1.618034f)) 219 smps[i] = tmp * (1.0f - tmp) / tmpv; 220 else 221 if(tmp > 0.0f) 222 smps[i] = -1.0f; 223 else 224 smps[i] = -2.0f; 225 } 226 break; 227 case 14: 228 ws = powf(ws, 5.0f) * 80.0f + 0.0001f; //sigmoid 229 if(ws > 10.0f) 230 tmpv = 0.5f; 231 else 232 tmpv = 0.5f - 1.0f / (expf(ws) + 1.0f); 233 for(i = 0; i < n; ++i) { 234 smps[i] += offs; //add offset 235 // calculate sigmoid function 236 float tmp = smps[i] * ws; 237 if(tmp < -10.0f) 238 tmp = -10.0f; 239 else 240 if(tmp > 10.0f) 241 tmp = 10.0f; 242 tmp = 0.5f - 1.0f / (expf(tmp) + 1.0f); 243 // calculate the same for offset value 244 float tmpo = offs * ws; 245 if(tmpo < -10.0f) 246 tmpo = -10.0f; 247 else 248 if(tmpo > 10.0f) 249 tmpo = 10.0f; 250 tmpo = 0.5f - 1.0f / (expf(tmpo) + 1.0f); 251 252 smps[i] = tmp / tmpv; 253 smps[i] -= tmpo / tmpv; // subtract offset 254 } 255 break; 256 case 15: // tanh soft limiter 257 // f(x) = x / ((1+|x|^n)^(1/n)) // tanh approximation for n=2.5 258 // Formula from: Yeh, Abel, Smith (2007): SIMPLIFIED, PHYSICALLY-INFORMED MODELS OF DISTORTION AND OVERDRIVE GUITAR EFFECTS PEDALS 259 par = (20.0f) * par * par + (0.1f) * par + 1.0f; //Pfunpar=32 -> n=2.5 260 ws = ws * ws * 35.0f + 1.0f; 261 for(i = 0; i < n; ++i) { 262 smps[i] *= ws;// multiply signal to drive it in the saturation of the function 263 smps[i] += offs; // add dc offset 264 smps[i] = smps[i] / powf(1+powf(fabsf(smps[i]), par), 1/par); 265 smps[i] -= offs / powf(1+powf(fabsf(offs), par), 1/par); 266 } 267 break; 268 case 16: //cubic distortion 269 // f(x) = 1.5 * (x-(x^3/3)) 270 // Formula from: https://ccrma.stanford.edu/~jos/pasp/Soft_Clipping.html 271 // modified with factor 1.5 to go through [1,1] and [-1,-1] 272 ws = ws * ws * ws * 20.0f + 0.168f; // plain cubic at drive=44 273 for(i = 0; i < n; ++i) { 274 smps[i] *= ws; // multiply signal to drive it in the saturation of the function 275 smps[i] += offs; // add dc offset 276 if(fabsf(smps[i]) < 1.0f) 277 smps[i] = 1.5 * (smps[i] - (smps[i]*smps[i]*smps[i] / 3.0) ); 278 else 279 smps[i] = (smps[i] > 0 ? 1.0f : -1.0f); 280 //subtract offset with distortion function applied 281 smps[i] -= 1.5 * (offs - (offs*offs*offs / 3.0)); 282 } 283 break; 284 case 17: //square distortion 285 // f(x) = x*(2-abs(x)) 286 // Formula of cubic changed to square but still going through [1,1] and [-1,-1] 287 ws = ws * ws * ws * 20.0f + 0.168f; // plain square at drive=44 288 for(i = 0; i < n; ++i) { 289 smps[i] *= ws; // multiply signal to drive it in the saturation of the function 290 smps[i] += offs; // add dc offset 291 if(fabsf(smps[i]) < 1.0f) 292 smps[i] = smps[i]*(2-fabsf(smps[i])); 293 else 294 smps[i] = (smps[i] > 0 ? 1.0f : -1.0f); 295 //subtract offset with distortion function applied 296 smps[i] -= offs*(2-fabsf(offs)); 297 } 298 break; 299 } 300 } 301 302 }