diff --git a/ChangeLog b/ChangeLog index a3ce7b6cc..2be745e26 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,7 +1,19 @@ +2010-12-12 Robert Larice + * src/frontend/miscvars.c , + * src/frontend/trannoise/1-f-code.c , + * src/include/1-f-code.h , + * src/include/fftext.h , + * src/spicelib/devices/vsrc/vsrcacct.c , + * src/spicelib/devices/vsrc/vsrcdefs.h , + * src/spicelib/devices/vsrc/vsrcload.c , + * src/spicelib/devices/vsrc/vsrcpar.c : + rewrite TRNOISE, with the intention to separate the noise sequenze computation + from its use in the VSRC device. + 2010-12-12 Holger Vogt * vsrcacct.c: Patch von R. Larice für sichere Breakpoints bei TRNOISE - + 2010-12-11 Robert Larice * src/include/bool.h : Change bool from `unsigned char' to `int' diff --git a/src/frontend/miscvars.c b/src/frontend/miscvars.c index 0c57f4010..f603772fc 100644 --- a/src/frontend/miscvars.c +++ b/src/frontend/miscvars.c @@ -94,6 +94,7 @@ char *ft_setkwords[] = { "noprintscale", "nosort", "nosubckt", + "notrnoise", "numdgt", "opts", "pivrel", diff --git a/src/frontend/trannoise/1-f-code.c b/src/frontend/trannoise/1-f-code.c index e321a901f..ce0257f90 100644 --- a/src/frontend/trannoise/1-f-code.c +++ b/src/frontend/trannoise/1-f-code.c @@ -12,8 +12,10 @@ #include #include #include // var. argumente -#include "1-f-code.h" #include "ngspice.h" +#include "cpextern.h" +#include "cktdefs.h" +#include "1-f-code.h" #include "fftext.h" #include "wallace.h" @@ -60,3 +62,104 @@ float alpha) fprintf(stdout,"%d (2e%d) one over f values created\n", n_pts, n_exp); } + + +/*-----------------------------------------------------------------------------*/ + +void +trnoise_state_gen(struct trnoise_state *this, CKTcircuit *ckt) +{ + if(this->top == 0) { + + if(cp_getvar("notrnoise", CP_BOOL, NULL)) + this -> NA = this -> TS = this -> NALPHA = this -> NAMP = 0.0; + + if((this->NALPHA > 0.0) && (this->NAMP > 0.0)) { + + // add 10 steps for start up sequence + size_t nosteps = (size_t) (ckt->CKTfinalTime / this->TS) + 10; + + size_t newsteps = 1; + long int newexp = 0; + // generate number of steps as power of 2 + while(newsteps < nosteps) { + newsteps <<= 1; + newexp++; + } + + this->oneof = TMALLOC(float, newsteps); + this->oneof_length = newsteps; + + f_alpha((int) newsteps, newexp, + this -> oneof, + (float) this -> NAMP, + (float) this -> NALPHA); + } + + trnoise_state_push(this, 0.0); /* first is deterministic */ + return; + } + + + // make use of two random variables per call to rgauss() + { + double ra1, ra2; + double NA = this -> NA; + + if(NA != 0.0) { + +#ifdef FastRand + // use FastNorm3 + ra1 = NA * FastNorm; + ra2 = NA * FastNorm; +#elif defined (WaGauss) + // use WallaceHV + ra1 = NA * GaussWa; + ra2 = NA * GaussWa; +#else + rgauss(&ra1, &ra2); + ra1 *= NA; + ra2 *= NA; +#endif + + } else { + + ra1 = 0.0; + ra2 = 0.0; + + } + + if(this -> oneof) { + + if(this->top + 1 >= this->oneof_length) { + fprintf(stderr,"ouch, noise data exhausted\n"); + exit(1); + } + + ra1 += this->oneof[this->top] - this->oneof[0]; + ra2 += this->oneof[this->top + 1] - this->oneof[0]; + } + + trnoise_state_push(this, ra1); + trnoise_state_push(this, ra2); + } + +} + + +struct trnoise_state * +trnoise_state_init(double NA, double TS, double NALPHA, double NAMP) +{ + struct trnoise_state *this = TMALLOC(struct trnoise_state, 1); + + this->NA = NA; + this->TS = TS; + this->NALPHA = NALPHA; + this->NAMP = NAMP; + + this -> top = 0; + + this -> oneof = NULL; + + return this; +} diff --git a/src/include/1-f-code.h b/src/include/1-f-code.h index c7de952d5..04a30dfff 100644 --- a/src/include/1-f-code.h +++ b/src/include/1-f-code.h @@ -5,3 +5,45 @@ void f_alpha(int n_pts, int n_exp, float X[], float Q_d, float alpha); void rvfft(float X[], unsigned long int n); + + +#define TRNOISE_STATE_MEM_LEN 4 +struct trnoise_state +{ + double points[TRNOISE_STATE_MEM_LEN]; + size_t top; + + double NA, TS, NAMP, NALPHA; + + float *oneof; + size_t oneof_length; +}; + + +struct trnoise_state *trnoise_state_init(double NA, double TS, double NALPHA, double NAMP); + +void trnoise_state_gen(struct trnoise_state *this, CKTcircuit *ckt); +void trnoise_state_free(struct trnoise_state *this); + + +static inline void +trnoise_state_push(struct trnoise_state *this, double val) +{ + this->points[this->top++ % TRNOISE_STATE_MEM_LEN] = val; +} + + +static inline double +trnoise_state_get(struct trnoise_state *this, CKTcircuit *ckt, size_t index) +{ + while(index >= this->top) + trnoise_state_gen(this, ckt); + + if(index + TRNOISE_STATE_MEM_LEN < this->top) { + fprintf(stderr, "ouch, trying to fetch from the past %d %d\n", + index, this->top); + exit(1); + } + + return this->points[index % TRNOISE_STATE_MEM_LEN]; +} diff --git a/src/include/fftext.h b/src/include/fftext.h index 25dfd7386..a70432826 100644 --- a/src/include/fftext.h +++ b/src/include/fftext.h @@ -16,7 +16,7 @@ int fftInit(long M); /* OUTPUTS */ /* private cosine and bit reversed tables */ -void fftFree(); +void fftFree(void); // release storage for all private cosine and bit reversed tables void ffts(float *data, long M, long Rows); diff --git a/src/spicelib/devices/vsrc/vsrcacct.c b/src/spicelib/devices/vsrc/vsrcacct.c index 18897d286..f3e47057b 100644 --- a/src/spicelib/devices/vsrc/vsrcacct.c +++ b/src/spicelib/devices/vsrc/vsrcacct.c @@ -10,6 +10,7 @@ Author: 1985 Thomas L. Quarles #include "sperror.h" #include "suffix.h" #include "missing_math.h" +#include "1-f-code.h" extern int fftInit(long M); extern void fftFree(void); @@ -193,21 +194,18 @@ VNoi1 1 0 DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise 0, time step, exponent < 2, rms value */ case TRNOISE: { - double NA, NT, NALPHA, NAMP, TS; - NA = here->VSRCcoeffs[0]; // input is rms value - NT = here->VSRCcoeffs[1]; // time step + struct trnoise_state *state = here -> VSRCtrnoise_state; + double TS = state -> TS; - NALPHA = here->VSRCfunctionOrder > 2 - ? here->VSRCcoeffs[2] : 0.0; - NAMP = here->VSRCfunctionOrder > 3 - && here->VSRCcoeffs[3] != 0.0 - && here->VSRCcoeffs[2] != 0.0 - ? here->VSRCcoeffs[3] : 0.0; - - if ((NT == 0.) || ((NA == 0.) && (NAMP == 0.))) // no further breakpoint if value not given + if (TS == 0.0) // no further breakpoint if value not given break; - TS = NT; + + /* FIXME, dont' want this here, over to aof_get or somesuch */ + if (ckt->CKTtime == 0.0) { + printf("VSRC: free fft tables\n"); + fftFree(); + } if(ckt->CKTbreak) { diff --git a/src/spicelib/devices/vsrc/vsrcdefs.h b/src/spicelib/devices/vsrc/vsrcdefs.h index c6ea14abb..eaa7d00c4 100644 --- a/src/spicelib/devices/vsrc/vsrcdefs.h +++ b/src/spicelib/devices/vsrc/vsrcdefs.h @@ -11,6 +11,8 @@ Author: 1985 Thomas L. Quarles #include "gendefs.h" #include "complex.h" +struct trnoise_state; + /* * structures to describe independent voltage sources */ @@ -49,14 +51,7 @@ typedef struct sVSRCinstance { double VSRCdF1phase; /* distortion f1 phase */ double VSRCdF2phase; /* distortion f2 phase */ - /*transient noise*/ - double VSRCprevTime; /*last time a new random value was issued*/ - double VSRCprevVal; /*last value issued at prevTime*/ - double VSRCnewVal; /*new value issued at prevTime*/ - double VSRCsecRand; /*second random value not yet used*/ - float *VSRConeof; /*pointer to array of 1 over f noise values */ - long int VSRCncount; /* counter to retrieve noise values */ - /*end of noise*/ + struct trnoise_state *VSRCtrnoise_state; /* transient noise */ double VSRCr; /* pwl repeat */ double VSRCrdelay; /* pwl delay period */ diff --git a/src/spicelib/devices/vsrc/vsrcload.c b/src/spicelib/devices/vsrc/vsrcload.c index ee1455992..20f725693 100644 --- a/src/spicelib/devices/vsrc/vsrcload.c +++ b/src/spicelib/devices/vsrc/vsrcload.c @@ -11,14 +11,6 @@ $Id$ #include "trandefs.h" #include "sperror.h" #include "suffix.h" -#undef WaGauss -#ifdef FastRand -#include "FastNorm3.h" -#elif defined (WaGauss) -#include "wallace.h" -#else -extern void rgauss(double* py1, double* py2); -#endif #include "1-f-code.h" #ifdef XSPICE_EXP @@ -343,204 +335,24 @@ VNoi1 1 0 DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise 0, time step, exponent < 2, rms value */ case TRNOISE: { - /* Generate voltage point every TS with amplitude NA * ra, - where ra is drawn from a random number generator with - gaussian distribution with mean 0 and standard deviation 1 - */ -//#define PRVAL -// typedef int bool; - - double newval=0.0, lastval=0.0, lasttime=0.0; - double NA, NT, TS; - double V1, V2, basetime = 0.; - double scalef, ra1, ra2; - float NALPHA, NAMP; - - long int nosteps, newsteps = 1, newexp = 0; - - bool aof = FALSE; - - NA = here->VSRCcoeffs[0]; // input is rms value - NT = here->VSRCcoeffs[1]; // time step + struct trnoise_state *state = here -> VSRCtrnoise_state; - scalef = NA; -// scalef = NA*1.32; + double TS = state -> TS; - NALPHA = here->VSRCfunctionOrder > 2 - ? (float)here->VSRCcoeffs[2] : 0.0f; - NAMP = here->VSRCfunctionOrder > 3 - && here->VSRCcoeffs[3] != 0.0 - && here->VSRCcoeffs[2] != 0.0 - ? (float)here->VSRCcoeffs[3] : 0.0f; + if(TS == 0.0) { + value = 0.0; + } else { + size_t n1 = (size_t) floor(time / TS); - if ((NT == 0.) || ((NA == 0.) && (NAMP == 0.))) { - value = here->VSRCdcValue; - goto noiDone; - } - else - TS = NT; /* time step for noise */ - - if ((NALPHA > 0.0) && (NAMP > 0.0)) aof = TRUE; - - lasttime = here->VSRCprevTime; - lastval = here->VSRCprevVal; - newval = here->VSRCnewVal; - /* set all data: DC, white, 1of */ - if (time <= 0 /*ckt->CKTstep*/) { - /* data are already set */ - if ((here->VSRCprevVal != 0) || (here->VSRCnewVal != 0)) { - value = here->VSRCprevVal; - goto noiDone; - } - lasttime = 0.0; - here->VSRCsecRand = 2.; /* > 1, invalid number out of the random number range */ - /* get two random samples */ -#ifdef FastRand - // use FastNorm3 - here->VSRCprevVal = scalef * GaussWa; - here->VSRCnewVal = scalef * GaussWa; -#elif defined (WaGauss) - // use WallaceHV - here->VSRCprevVal = scalef * GaussWa; - here->VSRCnewVal = scalef * GaussWa; -#else - // make use of two random variables per call to rgauss() - rgauss(&ra1, &ra2); - here->VSRCprevVal = scalef * ra1; - // choose to set start value to 0 - here->VSRCprevVal = 0; - here->VSRCnewVal = scalef * ra2; -#endif - /* generate 1 over f noise at time 0 */ - if (aof) { - if (here->VSRCncount==0) { - // add 10 steps for start up sequence - nosteps = (long)((ckt->CKTfinalTime)/TS) + 10; - // generate number of steps as power of 2 - while(newsteps < nosteps) { - newsteps <<= 1; - newexp++; - } - here->VSRConeof = TMALLOC(float, newsteps); //(float *)tmalloc(sizeof(float) * newsteps); -#ifdef PRVAL - printf("ALPHA: %f, GAIN: %e\n", NALPHA, NAMP); -#endif - f_alpha(newsteps, newexp, here->VSRConeof, NAMP, NALPHA); -#ifdef PRVAL - printf("Noi1: %e, Noi2: %e\n", here->VSRConeof[10], here->VSRConeof[100]); -#endif - here->VSRCprevVal += here->VSRConeof[here->VSRCncount]; - here->VSRCncount++; - here->VSRCnewVal += here->VSRConeof[here->VSRCncount]; - here->VSRCncount++; - value = newval; - // add DC - here->VSRCprevVal += here->VSRCdcValue; - here->VSRCnewVal += here->VSRCdcValue; - value = here->VSRCprevVal; -#ifdef PRVAL - printf("start1, time: %e, outp: %e, rnd: %e\n", time, newval, testval); -#endif - } else { // here->VSRCncount > 0 - // add DC - here->VSRCprevVal += here->VSRCdcValue; - here->VSRCnewVal += here->VSRCdcValue; - value = here->VSRCprevVal; -#ifdef PRVAL - printf("start2, time: %e, outp: %e, rnd: %e\n", time, here->VSRCprevVal, testval); -#endif - } -#ifdef PRVAL - printf("time 0 value: %e for %s\n", here->VSRCprevVal, here->VSRCname); -#endif - goto loadDone; - } //aof - // add DC - here->VSRCprevVal += here->VSRCdcValue; - here->VSRCnewVal += here->VSRCdcValue; - value = here->VSRCprevVal; - here->VSRCprevTime = 0.; - goto loadDone; - } // time < 0 + double V1 = trnoise_state_get(state, ckt, n1); + double V2 = trnoise_state_get(state, ckt, n1+1); - V1 = here->VSRCprevVal; - V2 = here->VSRCnewVal; - if (here->VSRCprevTime == ckt->CKTtime) { - value = here->VSRCprevVal; - goto noiDone; + value = V1 + (V2 - V1) * (time / TS - n1); } - if (time > 0 && time < TS) { - value = V1 + (V2 - V1) * (time) / TS; - } - else if (time >= TS) { - /* repeating signal - figure out where we are in period */ - /* numerical correction to avoid basetime less than - next step, e.g. 4.99999999999999995 delivers a floor - of 4 instead of 5 */ - basetime = TS * floor(time*1.000000000001/TS); - time -= basetime; - -#define NSAMETIME(a,b) (fabs((a)-(b))<= NTIMETOL * TS) -#define NTIMETOL 1e-7 - - if NSAMETIME(time,0.) { - - /* get new random number */ -#ifdef FastRand - // use FastNorm3 - newval = scalef * FastNorm; -#elif defined (WaGauss) - // use WallaceHV - newval = scalef * GaussWa; -#else - // make use of two random variables per call to rgauss() - if (here->VSRCsecRand == 2.0) { - rgauss(&ra1, &ra2); - newval = scalef * ra1; - here->VSRCsecRand = scalef * ra2; - } - else { - newval = here->VSRCsecRand; - here->VSRCsecRand = 2.0; - } -#endif - V1 = here->VSRCprevVal = here->VSRCnewVal; - V2 = newval; // scale factor t.b.d. - if(here->VSRCdcGiven) V2 += here->VSRCdcValue; - if (aof) { - V2 += here->VSRConeof[here->VSRCncount]; -#ifdef PRVAL - printf("aof: %d\n", here->VSRCncount); -#endif - } - here->VSRCncount++; - value = V1; - here->VSRCnewVal = V2; - } else if (time > 0 && time < TS) { - value = V1 + (V2 - V1) * (time) / TS; -#ifdef PRVAL - printf("if1, time: %e, outp: %e, rnd: %e\n", ckt->CKTtime, - V1 + (V2 - V1) * (time) / TS, V2); -#endif - } else { /* time > TS should be never reached */ - value = V1 + (V2 - V1) * (time-TS) / TS; -#ifdef PRVAL - printf("if2, time: %e, outp: %e, rnd: %e\n", ckt->CKTtime, - V1 + (V2 - V1) * (time-TS) / TS, V2); -#endif - } - here->VSRCprevTime = ckt->CKTtime; - } -noiDone: - if (time >=ckt->CKTfinalTime) { - /* free the 1of memory */ - if (here->VSRConeof) tfree(here->VSRConeof); - /* reset the 1of counter */ - here->VSRCncount = 0; - } - goto loadDone; + if(here -> VSRCdcGiven) + value += here->VSRCdcValue; } // case break; } // switch diff --git a/src/spicelib/devices/vsrc/vsrcpar.c b/src/spicelib/devices/vsrc/vsrcpar.c index ea37d2db1..7a4f9bd22 100644 --- a/src/spicelib/devices/vsrc/vsrcpar.c +++ b/src/spicelib/devices/vsrc/vsrcpar.c @@ -11,6 +11,7 @@ Modified: 2000 AlansFixes #include "ifsim.h" #include "sperror.h" #include "suffix.h" +#include "1-f-code.h" /* ARGSUSED */ @@ -169,12 +170,29 @@ VSRCparam(int param, IFvalue *value, GENinstance *inst, IFvalue *select) return(E_BADPARM); } break; - case VSRC_TRNOISE: + case VSRC_TRNOISE: { + double NA, TS; + double NALPHA = 0.0; + double NAMP = 0.0; + here->VSRCfunctionType = TRNOISE; here->VSRCfuncTGiven = TRUE; here->VSRCcoeffs = value->v.vec.rVec; here->VSRCfunctionOrder = value->v.numValue; here->VSRCcoeffsGiven = TRUE; + + NA = here->VSRCcoeffs[0]; // input is rms value + TS = here->VSRCcoeffs[1]; // time step + + if (here->VSRCfunctionOrder > 2) + NALPHA = here->VSRCcoeffs[2]; + + if (here->VSRCfunctionOrder > 3 && NALPHA != 0.0) + NAMP = here->VSRCcoeffs[3]; + + here->VSRCtrnoise_state = + trnoise_state_init(NA, TS, NALPHA, NAMP); + } break; default: return(E_BADPARM);