From 86f00e38fc176ee147cff7775d90ee7cfb2e34a6 Mon Sep 17 00:00:00 2001 From: Victor Bombi Date: Tue, 23 Jun 2015 19:57:04 +0200 Subject: [PATCH 1/3] dwglib changes --- source/DWGUGens/DWGBowed.cpp | 72 +++++++++++++++++++++++++++++++--- source/DWGUGens/DWGPlucked.cpp | 6 +-- source/DWGUGens/dwglib/DWG.cpp | 26 ++++++++++++ source/DWGUGens/dwglib/DWG.hpp | 29 +++++++++++++- 4 files changed, 124 insertions(+), 9 deletions(-) diff --git a/source/DWGUGens/DWGBowed.cpp b/source/DWGUGens/DWGBowed.cpp index f51ef5640a..6af1d86157 100644 --- a/source/DWGUGens/DWGBowed.cpp +++ b/source/DWGUGens/DWGBowed.cpp @@ -50,22 +50,23 @@ DWGBowedSimple::DWGBowedSimple(Unit* unit){ } void DWGBowedSimple::Release(float trig,float *out,int NumSamples){ - if(this->m_trig <=0 && trig > 0){ + if(this->m_trig <=0.0f && trig > 0.0f){ this->m_trig = trig; } - if((this->m_trig >0 && trig <= 0)){ + if((this->m_trig >0.0f && trig <= 0.0f)){ int relcount = this->relcount; float rellevel = this->rellevel; float rellevelstep = this->rellevelstep; for(int i=0; i 0){ + //if(relcount > 0){ rellevel -= rellevelstep; + rellevel = std::max(0.0f, rellevel); relcount--; - } + //} out[i] *=rellevel; } - if(relcount <=0) + if(relcount <=0 && rellevel == 0.0f) DoneAction(2,this); this->relcount = relcount; @@ -73,6 +74,13 @@ void DWGBowedSimple::Release(float trig,float *out,int NumSamples){ } } //////////////////////////////////////////////////// +struct DWGBowedStk : public DWGBowedSimple +{ + DWGBowedStk(Unit* unit); +}; +SCWrapClass(DWGBowedStk); +DWGBowedStk::DWGBowedStk(Unit* unit):DWGBowedSimple(unit){ SETCALC(DWGBowedStk_next);} +//////////////////////////////////////////////////// struct DWGBowed : public DWGBowedSimple { @@ -169,6 +177,16 @@ float Bow(float vb,float fb,float vsr_plus_vsl){ float ref = BowTable(vdeltap,fb); return ref * vdeltap; } + +float BowStk(float vb,float fb,float vsr_plus_vsl){ + float vdeltap = vb - (vsr_plus_vsl); + float slope = ( 5.0 - (4.0 * fb) ); + float sample = vdeltap + 0.001; //offset + sample *= slope; + sample = fabs( sample ) + 0.75; + sample = pow( sample, -4.0 ); + return sample; +} ///////////////////////////////////////// float hyperbolicbow(float deltav,float fb){ float mus = 0.8; @@ -477,6 +495,49 @@ void DWGBowedSimple_next(DWGBowedSimple *unit, int inNumSamples) unit->Release(trig,out,inNumSamples); } +void DWGBowedStk_next(DWGBowedStk *unit, int inNumSamples) +{ + + float *out = OUT(0); + float freq = ZIN0(0); + float bowvelocity = ZIN0(1); + float bowforce = ZIN0(2); + float trig = ZIN0(3); + float pos = ZIN0(4); + + float c1 = ZIN0(6); + float c3 = std::max(ZIN0(7),(float)1e-9); + + + unit->Loss.setcoeffs(freq,c1,c3); + float lossdelay = unit->Loss.groupdelay(freq,SAMPLERATE); + float deltot = SAMPLERATE/freq; + float del1 = (deltot - lossdelay)*0.5 - 1; + + + float PMAS,PMAS2; + float PMENOS; + for (int i=0; i < inNumSamples; ++i) + { + float vel = unit->DWGF[0].get(pos*del1) + unit->DWGF[1].get(del1*(1-pos)); + vel = BowStk(bowvelocity,bowforce,vel); + //vel = unit->Bow2(bowvelocity,bowforce,vel); + unit->DWGF[0].add(vel,pos*del1); + unit->DWGF[1].add(vel,del1*(1-pos)); + + PMAS = unit->DWGF[0].delay(del1); + PMAS2 = unit->Loss.pushConvol(PMAS); + PMENOS = unit->DWGF[1].delay(del1); + + unit->DWGF[1].push(-PMAS2); + unit->DWGF[0].push(-PMENOS); + + out[i] = PMAS;// + PMAS2; + + } + unit->Release(trig,out,inNumSamples); +} + ///////////////////////////////////////// PluginLoad(DWGBowed) { @@ -484,5 +545,6 @@ PluginLoad(DWGBowed) DefineDtorUnit(DWGBowed); DefineDtorUnit(DWGBowedTor); DefineDtorUnit(DWGBowedSimple); + DefineDtorUnit(DWGBowedStk); DefineDtorUnit(DWGSoundBoard); } diff --git a/source/DWGUGens/DWGPlucked.cpp b/source/DWGUGens/DWGPlucked.cpp index b177f8f882..31136dc2d8 100644 --- a/source/DWGUGens/DWGPlucked.cpp +++ b/source/DWGUGens/DWGPlucked.cpp @@ -133,10 +133,10 @@ void DWGPluckedStiff_next(DWGPluckedStiff *unit, int inNumSamples) float B = ZIN0(8)/100000; unit->disper.setcoeffs(freq,B); - float disperdelay = unit->disper.groupdelay(SAMPLERATE); - + float disperdelay = unit->disper.phasedelay(SAMPLERATE); + //Print("diff %f\n",disperdelay - unit->disper.groupdelay(SAMPLERATE)); unit->Loss.setcoeffs(freq,c1,c3); - float lossdelay = unit->Loss.groupdelay(freq,SAMPLERATE); + float lossdelay = unit->Loss.phasedelay(freq,SAMPLERATE); float deltot = SAMPLERATE/freq; float del1 = (deltot - lossdelay - disperdelay)*0.5 - 1; diff --git a/source/DWGUGens/dwglib/DWG.cpp b/source/DWGUGens/dwglib/DWG.cpp index b0ac32a964..b910ba2a55 100644 --- a/source/DWGUGens/dwglib/DWG.cpp +++ b/source/DWGUGens/dwglib/DWG.cpp @@ -100,6 +100,32 @@ float groupdelay(float f,float *B,int sizeB,float *A,int sizeA,float FS){ float c[2] ; complex_divide(Num,AxB,c); return c[0]; } +//// +void evalpoly(float omega,float* B,int sizeB,float* A,int sizeA,float H[2]){ + float HB[2]; + float HA[2]; + evalpolyB(omega,B,sizeB,HB); + evalpolyA(omega,A,sizeA,HA); + complex_divide(HB,HA,H); +} +float PhaseIIR(float omega,float* B,int sizeB,float* A,int sizeA){ + float C[2]; + evalpoly(omega,B,sizeB,A,sizeA,C); + return atan2(C[1],C[0]); +} +float PhaseDelayDerive(float omega,float* B,int sizeB,float* A,int sizeA,float delta){ + float omega1 = omega - delta; + float omega2 = omega + delta; + float p1 = PhaseIIR(omega1,B,sizeB,A,sizeA); + float p2 = PhaseIIR(omega2,B,sizeB,A,sizeA); + return (-omega1*p2 + omega2*p1)/(2*delta*omega1*omega2); +} +float PhaseDelay(float f,float* B,int sizeB,float* A,int sizeA,float FS){ + float grpdel = ::groupdelay(f,B,sizeB,A,sizeA,FS); + float omega = 2.0*M_PI*f/FS; + return grpdel - omega*PhaseDelayDerive(omega,B,sizeB,A,sizeA); +} +/// long Nchoose(long n, long k) { long divisor = 1; long multiplier = n; diff --git a/source/DWGUGens/dwglib/DWG.hpp b/source/DWGUGens/dwglib/DWG.hpp index c526f609ef..93ccaeab8d 100644 --- a/source/DWGUGens/dwglib/DWG.hpp +++ b/source/DWGUGens/dwglib/DWG.hpp @@ -74,6 +74,8 @@ inline bool approximatelyEqual(float a, float b, float epsilon = 1e-7f) float absb = fabs(b); return fabs(a - b) <= ( (absa < absb ? absb : absa) * epsilon); } +float PhaseDelay(float f,float* B,int sizeB,float* A,int sizeA,float FS); +float PhaseDelayDerive(float omega,float* B,int sizeB,float* A,int sizeA,float delta=0.0005); float groupdelay(float f,float *B,int sizeB,float *A,int sizeA,float FS); long Nchoose(long n, long k); float ValimakiDispersion(float B, float f, int M); @@ -364,6 +366,16 @@ class LTITv } return grdel; } + float phasedelay(float f,float FS){ + //return ::PhaseDelay(f,KernelB,kernel_sizeB,KernelA,kernel_sizeA,FS); + //if(dirty_phdel){ + float grpdel = groupdelay(f,FS); + float omega = 2.0*M_PI*f/FS; + float phdel = grpdel - omega*::PhaseDelayDerive(omega,KernelB,kernel_sizeB,KernelA,kernel_sizeA); + //dirty_phdel = false; + //} + return phdel; + } }; //////////////specialization template<> @@ -410,6 +422,16 @@ class LTITv<1,1> } return grdel; } + float phasedelay(float f,float FS){ + //return ::PhaseDelay(f,&KernelB,1,&KernelA,1,FS); + //if(dirty_phdel){ + float grpdel = groupdelay(f,FS); + float omega = 2.0*M_PI*f/FS; + float phdel = grpdel - omega*::PhaseDelayDerive(omega,&KernelB,1,&KernelA,1); + //dirty_phdel = false; + //} + return phdel; + } }; ////////////////////////////////////////////////////////// template @@ -779,7 +801,7 @@ struct ThirianDispersion{ return; float D = ValimakiDispersion(B,freq,M); //if(D <=1) - // printf("D es %g\n",D); + //Print("D es %g\n",D); for(int i=0; ifreq = freq; @@ -791,6 +813,11 @@ struct ThirianDispersion{ return 0; return M*dispersion[0].groupdelay(freq,FS); } + float phasedelay(float FS){ + if(B==0) + return 0; + return M*dispersion[0].phasedelay(freq,FS); + } float filter(float a){ if(B==0) return a; From 197628e507fbad9deecaeddd41c0d049774ff261 Mon Sep 17 00:00:00 2001 From: Victor Bombi Date: Thu, 9 Jul 2015 13:52:56 +0200 Subject: [PATCH 2/3] dwglib: add KLjunctions and NJunctions and TUBE --- source/DWGUGens/dwglib/DWG.cpp | 12 ++ source/DWGUGens/dwglib/DWG.hpp | 275 +++++++++++++++++++++++++++++++++ 2 files changed, 287 insertions(+) diff --git a/source/DWGUGens/dwglib/DWG.cpp b/source/DWGUGens/dwglib/DWG.cpp index b910ba2a55..b9a285aab4 100644 --- a/source/DWGUGens/dwglib/DWG.cpp +++ b/source/DWGUGens/dwglib/DWG.cpp @@ -24,6 +24,18 @@ void volcar(const char *_Message, const char *_File, unsigned _Line) { Print("assertion %s ,%s:%d\n",_Message,_File,_Line); } +void kill_denormals(float &val) +{ + static const float anti_denormal = 1e-18; + val += anti_denormal; + val -= anti_denormal; +} +void kill_denormals(double &val) +{ + static const double anti_denormal = 1e-18; + val += anti_denormal; + val -= anti_denormal; +} float ValimakiDispersion(float B, float f, int M) { float C1,C2,k1,k2,k3; if(M==4) { diff --git a/source/DWGUGens/dwglib/DWG.hpp b/source/DWGUGens/dwglib/DWG.hpp index 93ccaeab8d..30a9582aff 100644 --- a/source/DWGUGens/dwglib/DWG.hpp +++ b/source/DWGUGens/dwglib/DWG.hpp @@ -79,6 +79,8 @@ float PhaseDelayDerive(float omega,float* B,int sizeB,float* A,int sizeA,float d float groupdelay(float f,float *B,int sizeB,float *A,int sizeA,float FS); long Nchoose(long n, long k); float ValimakiDispersion(float B, float f, int M); +void kill_denormals(float &val); +void kill_denormals(double &val); /////////////////////////// //Circular BufferT template @@ -224,7 +226,49 @@ class LagrangeT : public CircularBuffer2POWSizedT float delay(){return delay(actdelay);} }; +//1st order lagrange fractional delay +template +class Lagrange1T : public CircularBuffer2POWSizedT +{ + public: + float lastdelay; + float h[2]; + int ptL; + float actdelay; + Lagrange1T() + { + actdelay = 0; + lastdelay = 0; + ptL = CalcCoeffs(0); + //Print("LAGRANGE construc\n"); + } + int CalcCoeffs(float delay) + { + int intd =(int) delay; + float Dm1 = delay - (float)intd; + h[0] = (1.0 - Dm1); + h[1] = Dm1; + return intd ; + } + float delay(float pos) + { + assertv(size >= pos); + if (pos != lastdelay){ + ptL = CalcCoeffs(pos); + lastdelay = pos; + } + //return this->Buffer[(this->pointer - (int)pos) & this->mask]; + float sum = 0; + for(int i=0; i < 2; i++){ + sum += this->Buffer[(this->pointer + ptL + i) & this->mask]*h[i]; + } + //DUMPONNAN(sum); + return sum; + } + void setdelay(float del){actdelay = del;} + float delay(){return delay(actdelay);} +}; template class ConvolverT @@ -841,4 +885,235 @@ struct DCBlocker:public LTITv<2,1>{ this->dirty_grdel = true; } }; + +struct KL_Junction{ + Unit * unit; + float k; + float inR; + float inL; + float outR; + float outL; + KL_Junction(Unit * unit):inR(0),inL(0),outR(0),outL(0),k(0){this->unit = unit;}; + void* operator new(size_t sz, Unit* unit) { + //Print("KL_Junction new\n"); + void * ptr = NULL; + ptr = RTAlloc(unit->mWorld,sz); + /* + if(ptr==NULL) + Print("Null pointer\n"); + else + Print("not Null pointer\n");*/ + return ptr; + } + void operator delete(void* pObject) { + //Print("KL_Junction delete\n"); + KL_Junction * obj = (KL_Junction*)pObject; + RTFree(obj->unit->mWorld, pObject); + } + void set_areas(float a1,float a2){ + float num = a1 - a2; + if (num == 0.0) + k = 0.0; + else + k = num/(a1 + a2); + /* + if(a1 < 1e-18) + lossR = 0.0; + else + lossR = 1.0 - std::min(lossF/sqrt(a1),1.0); + if(a2 < 1e-18) + lossL = 0.0; + else + lossL = 1.0 - std::min(lossF/sqrt(a2),1.0); + */ + } + void go(){ + //inR *= lossR; + //inL *= lossL; + outR = inR*(1.0 + k) - inL*k; + outL = inL*(1.0 - k) + inR*k; + kill_denormals(outR); + kill_denormals(outL); + } +}; + +struct KL_JunctionU{ + Unit * unit; + float k; + float inR; + float inL; + float outR; + float outL; + KL_JunctionU(Unit * unit):inR(0),inL(0),outR(0),outL(0),k(0){this->unit = unit;}; + void* operator new(size_t sz, Unit* unit) { + //Print("KL_Junction new\n"); + void * ptr = NULL; + ptr = RTAlloc(unit->mWorld,sz); + /* + if(ptr==NULL) + Print("Null pointer\n"); + else + Print("not Null pointer\n");*/ + return ptr; + } + void operator delete(void* pObject) { + //Print("KL_Junction delete\n"); + KL_Junction * obj = (KL_Junction*)pObject; + RTFree(obj->unit->mWorld, pObject); + } + void set_areas(float a1,float a2){ + float num = a1 - a2; + if (num == 0.0) + k = 0.0; + else + k = num/(a1 + a2); + /* + if(a1 < 1e-18) + lossR = 0.0; + else + lossR = 1.0 - std::min(lossF/sqrt(a1),1.0); + if(a2 < 1e-18) + lossL = 0.0; + else + lossL = 1.0 - std::min(lossF/sqrt(a2),1.0); + */ + } + void go(){ + //inR *= lossR; + //inL *= lossL; + outR = inR*(1.0 - k) - inL*k; + outL = inL*(1.0 + k) + inR*k; + kill_denormals(outR); + kill_denormals(outL); + } +}; +template +class NJunction{ + public: + float in[N]; + float out[N]; + float alpha[N]; + /* + void* operator new(size_t sz, Unit* unit) { + return RTAlloc(unit->mWorld,sz); + } + void operator delete(void* pObject) { + //Print("KL_Junction delete\n"); + NJunction * obj = (NJunction*)pObject; + RTFree(obj->unit->mWorld, pObject); + }*/ + NJunction(){ + Print("NJunction const\n"); + for(int i=0;i +class NJunctionU{ + public: + float in[N]; + float out[N]; + float alpha[N]; + /* + void* operator new(size_t sz, Unit* unit) { + return RTAlloc(unit->mWorld,sz); + } + void operator delete(void* pObject) { + //Print("KL_Junction delete\n"); + NJunction * obj = (NJunction*)pObject; + RTFree(obj->unit->mWorld, pObject); + }*/ + NJunctionU(){ + Print("NJunctionU const\n"); + for(int i=0;i Rline; + Lagrange1T<1024> Lline; + float inR; + float inL; + float outR; + float outL; + float del; + float loss; + void* operator new(size_t sz, Unit* unit) { + //Print("TUBE new\n"); + void * ptr = NULL; + ptr = RTAlloc(unit->mWorld,sz); + + return ptr; + } + void operator delete(void* pObject) { + TUBE * obj = (TUBE*)pObject; + RTFree(obj->unit->mWorld, pObject); + } + TUBE(Unit * unit):inR(0),inL(0),outR(0),outL(0),del(0.0),loss(1){ + //Print("TUBE constructor\n"); + this->unit = unit; + } + void set_area(float a1,float lossF) + { + if(a1 < 1e-18) + loss = 0.0; + else + loss = 1.0 - std::min(lossF/sqrt(a1),1.0); + + } + void go(){ + Rline.push(inR); + Lline.push(inL); + outR = loss*Rline.delay(del); + outL = loss*Lline.delay(del); + } +}; #endif \ No newline at end of file From 3cd0ed4bc50f944813d7be697b53a62b1a75195c Mon Sep 17 00:00:00 2001 From: Victor Bombi Date: Thu, 15 Oct 2015 18:24:21 +0200 Subject: [PATCH 3/3] dwglib: implement multitap delay --- source/DWGUGens/dwglib/DWG.hpp | 67 ++++++++++++++++++++++++++++++++-- 1 file changed, 64 insertions(+), 3 deletions(-) diff --git a/source/DWGUGens/dwglib/DWG.hpp b/source/DWGUGens/dwglib/DWG.hpp index 30a9582aff..d0e65baf24 100644 --- a/source/DWGUGens/dwglib/DWG.hpp +++ b/source/DWGUGens/dwglib/DWG.hpp @@ -225,6 +225,61 @@ class LagrangeT : public CircularBuffer2POWSizedT void setdelay(float del){actdelay = del;} float delay(){return delay(actdelay);} }; +//3rd order lagrange fractional multitap delay +template +class LagrangeMtapT : public CircularBuffer2POWSizedT +{ + + public: + float lastdelay[taps]; + float h[taps][4]; + int ptL[taps]; + //float actdelay; + LagrangeMtapT() + { + //actdelay = 0; + for(int i=0;i= pos); + if (pos != lastdelay[tap]){ + ptL[tap] = CalcCoeffs(pos,tap); + lastdelay[tap] = pos; + } + //return this->Buffer[(this->pointer - (int)pos) & this->mask]; + float sum = 0; + for(int i=0; i < 4; i++){ + sum += this->Buffer[(this->pointer + ptL[tap] + i) & this->mask]*h[tap][i]; + } + //DUMPONNAN(sum); + return sum; + } + //void setdelay(float del){actdelay = del;} + //float delay(){return delay(actdelay);} +}; //1st order lagrange fractional delay template @@ -309,6 +364,10 @@ class LTIT void push(float a){ cbuf.push(a); } + float filter(float a){ + push(a); + return Convol(); + } float pushConvol(float a){ push(a); return Convol(); @@ -888,7 +947,7 @@ struct DCBlocker:public LTITv<2,1>{ struct KL_Junction{ Unit * unit; - float k; + float k,one_plus_k,one_minus_k; float inR; float inL; float outR; @@ -916,6 +975,8 @@ struct KL_Junction{ k = 0.0; else k = num/(a1 + a2); + one_plus_k = 1.0 + k; + one_minus_k = 1.0 - k; /* if(a1 < 1e-18) lossR = 0.0; @@ -930,8 +991,8 @@ struct KL_Junction{ void go(){ //inR *= lossR; //inL *= lossL; - outR = inR*(1.0 + k) - inL*k; - outL = inL*(1.0 - k) + inR*k; + outR = inR*(one_plus_k) - inL*k; + outL = inL*(one_minus_k) + inR*k; kill_denormals(outR); kill_denormals(outL); }