RTS noise

ChangeLog

@@ -1,3 +1,12 @@
+2010-12-18 Holger Vogt
+ * rnorrexp.c, randnumb.c, 1-f-code.c, main.c,
+   fteext.h, 1-f-code.h,
+   isrcacct.c, isrcload.c, isrcpar.c,
+   vsrcacct.c, vsrcload.c, vsrcpar.c,
+   maths/misc/makefile.am,
+   visualc/vngspice.vcproj, vngspice.sln:
+   Random telegraph noise added to independent voltage and current sources
+ 
 2010-12-17 Holger Vogt
  * isrc.c, isrcacct.c, isrcload.c, isrcpar.c, isrcdefs.h:
  transient noise in independent current source

src/frontend/trannoise/1-f-code.c

@@ -20,6 +20,7 @@
 #include "fftext.h"
 #include "wallace.h"
 
+extern double exprand(double);
 
 void f_alpha(int n_pts, int n_exp, float X[], float Q_d,
 float alpha)
@@ -72,7 +73,8 @@ 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;
+            this -> NA = this -> TS = this -> NALPHA = this -> NAMP =
+               this -> RTSAM = this -> RTSCAPT = this -> RTSEMT = 0.0;
 
         if((this->NALPHA > 0.0) && (this->NAMP > 0.0)) {
 
@@ -148,7 +150,7 @@ trnoise_state_gen(struct trnoise_state *this, CKTcircuit *ckt)
 
 
 struct trnoise_state *
-trnoise_state_init(double NA, double TS, double NALPHA, double NAMP)
+trnoise_state_init(double NA, double TS, double NALPHA, double NAMP, double RTSAM, double RTSCAPT, double RTSEMT)
 {
     struct trnoise_state *this = TMALLOC(struct trnoise_state, 1);
 
@@ -156,7 +158,13 @@ trnoise_state_init(double NA, double TS, double NALPHA, double NAMP)
     this->TS = TS;
     this->NALPHA = NALPHA;
     this->NAMP = NAMP;
-
+    this->RTSAM = RTSAM;
+    this->RTSCAPT = RTSCAPT;
+    this->RTSEMT = RTSEMT;
+    if (RTSAM > 0) {
+        this->RTScapTime = exprand(RTSCAPT);
+        this->RTSemTime = this->RTScapTime + exprand(RTSEMT);
+    }
     this -> top = 0;
 
     this -> oneof = NULL;

src/include/1-f-code.h

@@ -13,14 +13,17 @@ struct trnoise_state
     double points[TRNOISE_STATE_MEM_LEN];
     size_t top;
 
-    double NA, TS, NAMP, NALPHA;
+    double NA, TS, NAMP, NALPHA, RTSAM, RTSCAPT, RTSEMT;
 
     float *oneof;
     size_t oneof_length;
+
+    double RTScapTime, RTSemTime;
+    bool RTS;
 };
 
 
-struct trnoise_state *trnoise_state_init(double NA, double TS, double NALPHA, double NAMP);
+struct trnoise_state *trnoise_state_init(double NA, double TS, double NALPHA, double NAMP, double RTSAM, double RTSCAPT, double RTSEMT);
 
 void trnoise_state_gen(struct trnoise_state *this, CKTcircuit *ckt);
 void trnoise_state_free(struct trnoise_state *this);

src/include/fteext.h

@@ -231,6 +231,9 @@ extern bool check_autostop(char *what);
 
 /* randnumb.c */
 extern void TausSeed(void);
+/* rnorrexp.c */
+extern void zigset(unsigned long jsrseed);
+
 /* resource.c */
 
 extern void ft_ckspace(void);

src/main.c

@@ -848,7 +848,8 @@ main(int argc, char **argv)
     }
     cp_program = ft_sim->simulator;
 
-    srand(getpid()); //srandom(getpid());
+    srand(getpid());
+    zigset(getpid());
     TausSeed();
 
     /* --- Process command line options --- */

src/maths/misc/Makefile.am

@@ -15,7 +15,8 @@ libmathmisc_la_SOURCES = \
 	scalb.c		\
 	norm.h		\
 	norm.c		\
-	randnumb.c
+	randnumb.c  \
+	rnorrexp.c
 
 
 EXTRA_DIST = test_accuracy.c test_erfc.c

src/maths/misc/randnumb.c

@@ -84,6 +84,7 @@ void checkseed(void)
    if (cp_getvar("rndseed", CP_NUM, &newseed)) {
       if ((newseed > 0) && (oldseed != newseed)) {
          srand(newseed); //srandom(newseed);
+         zigset(newseed);
          TausSeed();
          oldseed = newseed;
          printf("Seed value for random number generator is set to %d\n", newseed);

src/maths/misc/rnorrexp.c

@@ -0,0 +1,133 @@
+/* The ziggurat method for RNOR and REXP
+Combine the code below with the main program in which you want
+normal or exponential variates.   Then use of RNOR in any expression
+will provide a standard normal variate with mean zero, variance 1,
+while use of REXP in any expression will provide an exponential variate
+with density exp(-x),x>0.
+Before using RNOR or REXP in your main, insert a command such as
+zigset(86947731 );
+with your own choice of seed value>0, rather than 86947731.
+(If you do not invoke zigset(...) you will get all zeros for RNOR and REXP.)
+For details of the method, see Marsaglia and Tsang, "The ziggurat method
+for generating random variables", Journ. Statistical Software.
+
+http://www.jstatsoft.org/v05/i08/supp/1
+download 2010-12-18
+*/
+
+#include <math.h>
+static unsigned long jz,jsr=123456789;
+
+#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
+#define UNI (.5f + (signed) SHR3*.2328306e-9f)
+#define IUNI SHR3
+
+static long hz;
+static unsigned long iz, kn[128], ke[256];
+static float wn[128],fn[128], we[256],fe[256];
+
+#define RNOR (hz=SHR3, iz=hz&127, (fabs(hz)<kn[iz])? hz*wn[iz] : nfix())
+#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())
+
+/* prototypes */
+double exprand(double);
+float nfix(void);
+float efix(void);
+void zigset(unsigned long jsrseed);
+
+/* nfix() generates variates from the residue when rejection in RNOR occurs. */
+
+float nfix(void)
+{
+const float r = 3.442620f;     /* The start of the right tail */
+static float x, y;
+ for(;;)
+  {  x=hz*wn[iz];      /* iz==0, handles the base strip */
+     if(iz==0)
+       { do{ x=-log(UNI)*0.2904764; y=-log(UNI);}	/* .2904764 is 1/r */
+        while(y+y<x*x);
+        return (hz>0)? r+x : -r-x;
+       }
+                         /* iz>0, handle the wedges of other strips */
+      if( fn[iz]+UNI*(fn[iz-1]-fn[iz]) < exp(-.5*x*x) ) return x;
+
+     /* initiate, try to exit for(;;) for loop*/
+      hz=SHR3;
+      iz=hz&127;
+      if(fabs(hz)<kn[iz]) return (hz*wn[iz]);
+  }
+
+}
+
+/* efix() generates variates from the residue when rejection in REXP occurs. */
+float efix(void)
+{ float x;
+ for(;;)
+  {  if(iz==0) return (7.69711f-log(UNI));          /* iz==0 */
+     x=jz*we[iz]; if( fe[iz]+UNI*(fe[iz-1]-fe[iz]) < exp(-x) ) return (x);
+
+      /* initiate, try to exit for(;;) loop */
+   jz=SHR3;
+   iz=(jz&255);
+   if(jz<ke[iz]) return (jz*we[iz]);
+  }
+}
+
+/* return an exponential random number */
+double exprand(double mean)
+{
+    return (double)REXP * mean;
+}
+
+/*--------This procedure sets the seed and creates the tables------*/
+
+void zigset(unsigned long jsrseed)
+{  const double m1 = 2147483648.0, m2 = 4294967296.;
+   double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
+   double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
+   int i;
+   jsr^=jsrseed;
+
+/* Set up tables for RNOR */
+   q=vn/exp(-.5*dn*dn);
+   kn[0]=(dn/q)*m1;
+   kn[1]=0;
+
+   wn[0]=q/m1;
+   wn[127]=dn/m1;
+
+   fn[0]=1.;
+   fn[127]=exp(-.5*dn*dn);
+
+    for(i=126;i>=1;i--)
+    {dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
+     kn[i+1]=(dn/tn)*m1;
+     tn=dn;
+     fn[i]=exp(-.5*dn*dn);
+     wn[i]=dn/m1;
+    }
+
+/* Set up tables for REXP */
+    q = ve/exp(-de);
+    ke[0]=(de/q)*m2;
+    ke[1]=0;
+
+    we[0]=q/m2;
+    we[255]=de/m2;
+
+    fe[0]=1.;
+    fe[255]=exp(-de);
+
+   for(i=254;i>=1;i--)
+  {de=-log(ve/de+exp(-de));
+   ke[i+1]= (de/te)*m2;
+   te=de;
+   fe[i]=exp(-de);
+   we[i]=de/m2;
+  }
+}
+
+
+
+
+

src/spicelib/devices/isrc/isrcacct.c

@@ -15,6 +15,7 @@ Author: 1985 Thomas L. Quarles
 extern int fftInit(long M);
 extern void fftFree(void);
 extern void rffts(float *data, long M, long Rows);
+extern double exprand(double);
 
 int
 ISRCaccept(CKTcircuit *ckt, GENmodel *inModel)
@@ -202,8 +203,9 @@ INoi1 1 0  DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise
 
                         struct trnoise_state *state = here -> ISRCtrnoise_state;
                         double TS = state -> TS;
+                        double RTSAM = state ->RTSAM;
 
-                        if (TS == 0.0) // no further breakpoint if value not given
+                        if ((TS == 0.0) &&  (RTSAM == 0.0)) // no further breakpoint if value not given
                             break;
 
                         /* FIXME, dont' want this here, over to aof_get or somesuch */
@@ -228,6 +230,46 @@ INoi1 1 0  DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise
                                     return(error);
                             }
                         }
+
+
+                    if (RTSAM > 0) {
+                        double RTScapTime = state->RTScapTime;
+                        double RTSemTime = state->RTSemTime;
+                        double RTSCAPT = state->RTSCAPT;
+                        double RTSEMT = state->RTSEMT;
+
+                        if (ckt->CKTtime == 0) {
+                            if (ckt->CKTbreak) {
+                                error = CKTsetBreak(ckt, RTScapTime);
+                                if(error)
+                                    return(error);
+                            }
+                        }
+
+                        if(AlmostEqualUlps(RTScapTime, ckt->CKTtime, 3)) {
+                            if (ckt->CKTbreak) {
+                                error = CKTsetBreak(ckt, RTSemTime);
+                                if(error)
+                                    return(error);
+                            }
+                        }
+                        
+                        if(AlmostEqualUlps(RTSemTime, ckt->CKTtime, 3)) {
+                            /* new values */ 
+                            RTScapTime = here -> ISRCtrnoise_state ->RTScapTime = ckt->CKTtime + exprand(RTSCAPT);
+                            here -> ISRCtrnoise_state ->RTSemTime = RTScapTime + exprand(RTSEMT);
+
+                            if (ckt->CKTbreak) {
+                                error = CKTsetBreak(ckt, RTScapTime);
+                                if(error)
+                                    return(error);
+                            }
+
+
+                        }
+
+                    }
+
                     }
                     break;
 

src/spicelib/devices/isrc/isrcload.c

@@ -327,10 +327,13 @@ INoi1 1 0  DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise
                         struct trnoise_state *state = here -> ISRCtrnoise_state;
 
                         double TS = state -> TS;
+                        double RTSAM = state->RTSAM;
 
+                        /* no noise */
                         if(TS == 0.0) {
                             value = 0.0;
                         } else {
+                            /* 1/f and white noise */
                             size_t n1 = (size_t) floor(time / TS);
 
                             double V1 = trnoise_state_get(state, ckt, n1);
@@ -339,6 +342,14 @@ INoi1 1 0  DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise
                             value = V1 + (V2 - V1) * (time / TS - n1);
                         }
 
+                        /* RTS noise */
+                        if (RTSAM > 0) {
+                            double RTScapTime = state->RTScapTime;
+                            if (time >= RTScapTime)
+                                value += RTSAM;
+                        }
+
+                        /* DC value */
                         if(here -> ISRCdcGiven)
                             value += here->ISRCdcValue;
                     } // case

src/spicelib/devices/isrc/isrcpar.c

@@ -155,6 +155,9 @@ ISRCparam(int param, IFvalue *value, GENinstance *inst, IFvalue *select)
             double NA, TS;
             double NALPHA = 0.0;
             double NAMP   = 0.0;
+            double RTSAM   = 0.0;
+            double RTSCAPT   = 0.0;
+            double RTSEMT   = 0.0;
 
             here->ISRCfunctionType = TRNOISE;
             here->ISRCfuncTGiven = TRUE;
@@ -171,8 +174,18 @@ ISRCparam(int param, IFvalue *value, GENinstance *inst, IFvalue *select)
             if (here->ISRCfunctionOrder > 3 && NALPHA != 0.0)
                 NAMP = here->ISRCcoeffs[3];
 
+            if (here->ISRCfunctionOrder > 4)
+                RTSAM = here->ISRCcoeffs[4]; // RTS amplitude
+
+            if (here->ISRCfunctionOrder > 5 && RTSAM != 0.0)
+                RTSCAPT = here->ISRCcoeffs[5]; // RTS trap capture time
+
+            if (here->ISRCfunctionOrder > 6 && RTSAM != 0.0)
+                RTSEMT = here->ISRCcoeffs[6]; // RTS trap emission time
+
             here->ISRCtrnoise_state =
-                trnoise_state_init(NA, TS, NALPHA, NAMP);
+                trnoise_state_init(NA, TS, NALPHA, NAMP, RTSAM, RTSCAPT, RTSEMT);
+
         }
             break;	
 

src/spicelib/devices/vsrc/vsrcacct.c

@@ -15,6 +15,7 @@ Author: 1985 Thomas L. Quarles
 extern int fftInit(long M);
 extern void fftFree(void);
 extern void rffts(float *data, long M, long Rows);
+extern double exprand(double);
 
 #define SAMETIME(a,b)    (fabs((a)-(b))<= TIMETOL * PW)
 #define TIMETOL    1e-7
@@ -197,8 +198,9 @@ VNoi1 1 0  DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise
 
                     struct trnoise_state *state = here -> VSRCtrnoise_state;
                     double TS = state -> TS;
+                    double RTSAM = state ->RTSAM;
 
-                    if (TS == 0.0) // no further breakpoint if value not given
+                    if ((TS == 0.0) &&  (RTSAM == 0.0)) // no further breakpoint if value not given
                         break;
 
                     /* FIXME, dont' want this here, over to aof_get or somesuch */
@@ -223,6 +225,45 @@ VNoi1 1 0  DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise
                                 return(error);
                         }
                     }
+
+                    if (RTSAM > 0) {
+                        double RTScapTime = state->RTScapTime;
+                        double RTSemTime = state->RTSemTime;
+                        double RTSCAPT = state->RTSCAPT;
+                        double RTSEMT = state->RTSEMT;
+
+                        if (ckt->CKTtime == 0) {
+                            if (ckt->CKTbreak) {
+                                error = CKTsetBreak(ckt, RTScapTime);
+                                if(error)
+                                    return(error);
+                            }
+                        }
+
+                        if(AlmostEqualUlps(RTScapTime, ckt->CKTtime, 3)) {
+                            if (ckt->CKTbreak) {
+                                error = CKTsetBreak(ckt, RTSemTime);
+                                if(error)
+                                    return(error);
+                            }
+                        }
+                        
+                        if(AlmostEqualUlps(RTSemTime, ckt->CKTtime, 3)) {
+                            /* new values */ 
+                            RTScapTime = here -> VSRCtrnoise_state ->RTScapTime = ckt->CKTtime + exprand(RTSCAPT);
+                            here -> VSRCtrnoise_state ->RTSemTime = RTScapTime + exprand(RTSEMT);
+
+                            if (ckt->CKTbreak) {
+                                error = CKTsetBreak(ckt, RTScapTime);
+                                if(error)
+                                    return(error);
+                            }
+
+
+                        }
+
+                    }
+
                     }
                     break;
  				

src/spicelib/devices/vsrc/vsrcload.c

@@ -236,13 +236,13 @@ VSRCload(GENmodel *inModel, CKTcircuit *ckt)
                                 sin((2.0 * M_PI * FC * time) +
                                 MDI * sin(2 * M_PI * FS * time));
 #endif /* XSPICE */
-/* gtri - end - wbk - add PHASE parameters */
+    /* gtri - end - wbk - add PHASE parameters */
                     }
                     break;
                     case AM:{
     		
                         double VA, FC, MF, VO, TD;
-/* gtri - begin - wbk - add PHASE parameters */
+    /* gtri - begin - wbk - add PHASE parameters */
 #ifdef XSPICE
                         double PHASEC, PHASES;
                         double phasec;
@@ -284,7 +284,7 @@ VSRCload(GENmodel *inModel, CKTcircuit *ckt)
 #endif			
                         }
     		    
-/* gtri - end - wbk - add PHASE parameters */
+    /* gtri - end - wbk - add PHASE parameters */
                     }
                     break;
                     case PWL: {
@@ -333,16 +333,22 @@ VNoi2 2 0  DC 0 TRNOISE(10n 0.5n 0 0n) : generate gaussian distributed noise
                         rms value, time step, 0 0
 VNoi1 1 0  DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise
                         0,  time step, exponent < 2, rms value
+
+VNoi3 3 0  DC 0 TRNOISE(0 0 0 0 15m 22u 50u) : generate RTS noise
+                        0 0 0 0, amplitude, capture time, emission time
 */
                     case TRNOISE: {
 
                         struct trnoise_state *state = here -> VSRCtrnoise_state;
 
                         double TS = state -> TS;
+                        double RTSAM = state->RTSAM;
 
+                        /* no noise */
                         if(TS == 0.0) {
                             value = 0.0;
                         } else {
+                            /* 1/f and white noise */
                             size_t n1 = (size_t) floor(time / TS);
 
                             double V1 = trnoise_state_get(state, ckt, n1);
@@ -351,6 +357,14 @@ VNoi1 1 0  DC 0 TRNOISE(0n 0.5n 1 10n) : generate 1/f noise
                             value = V1 + (V2 - V1) * (time / TS - n1);
                         }
 
+                        /* RTS noise */
+                        if (RTSAM > 0) {
+                            double RTScapTime = state->RTScapTime;
+                            if (time >= RTScapTime)
+                                value += RTSAM;
+                        }
+
+                            /* DC value */
                         if(here -> VSRCdcGiven)
                             value += here->VSRCdcValue;
                     } // case

src/spicelib/devices/vsrc/vsrcpar.c

@@ -174,6 +174,9 @@ VSRCparam(int param, IFvalue *value, GENinstance *inst, IFvalue *select)
             double NA, TS;
             double NALPHA = 0.0;
             double NAMP   = 0.0;
+            double RTSAM   = 0.0;
+            double RTSCAPT   = 0.0;
+            double RTSEMT   = 0.0;
 
             here->VSRCfunctionType = TRNOISE;
             here->VSRCfuncTGiven = TRUE;
@@ -185,13 +188,22 @@ VSRCparam(int param, IFvalue *value, GENinstance *inst, IFvalue *select)
             TS = here->VSRCcoeffs[1]; // time step
 
             if (here->VSRCfunctionOrder > 2)
-                NALPHA = here->VSRCcoeffs[2];
+                NALPHA = here->VSRCcoeffs[2]; // 1/f exponent
 
             if (here->VSRCfunctionOrder > 3 && NALPHA != 0.0)
-                NAMP = here->VSRCcoeffs[3];
+                NAMP = here->VSRCcoeffs[3]; // 1/f amplitude
+
+            if (here->VSRCfunctionOrder > 4)
+                RTSAM = here->VSRCcoeffs[4]; // RTS amplitude
+
+            if (here->VSRCfunctionOrder > 5 && RTSAM != 0.0)
+                RTSCAPT = here->VSRCcoeffs[5]; // RTS trap capture time
+
+            if (here->VSRCfunctionOrder > 6 && RTSAM != 0.0)
+                RTSEMT = here->VSRCcoeffs[6]; // RTS trap emission time
 
             here->VSRCtrnoise_state =
-                trnoise_state_init(NA, TS, NALPHA, NAMP);
+                trnoise_state_init(NA, TS, NALPHA, NAMP, RTSAM, RTSCAPT, RTSEMT);
         }
             break;		
         default:

visualc/vngspice.sln

@@ -40,4 +40,9 @@ Global
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE
 	EndGlobalSection
+	GlobalSection(ExtensibilityGlobals) = postSolution
+		AMDCaProjectFile = D:\Spice_general\ng-spice-rework\visualc\CodeAnalyst\vngspice.caw
+		AMDCaPersistentStartup = vngspice
+		AMDCaPersistentConfig = Debug|Win32
+	EndGlobalSection
 EndGlobal

visualc/vngspice.vcproj

@@ -1891,11 +1891,11 @@
 				>
 			</File>
 			<File
-				RelativePath="..\src\frontend\plotting\grid.h"
+				RelativePath="..\src\include\grid.h"
 				>
 			</File>
 			<File
-				RelativePath="..\src\include\grid.h"
+				RelativePath="..\src\frontend\plotting\grid.h"
 				>
 			</File>
 			<File
@@ -2003,11 +2003,11 @@
 				>
 			</File>
 			<File
-				RelativePath="..\src\frontend\inp.h"
+				RelativePath="..\src\spicelib\parser\inp.h"
 				>
 			</File>
 			<File
-				RelativePath="..\src\spicelib\parser\inp.h"
+				RelativePath="..\src\frontend\inp.h"
 				>
 			</File>
 			<File
@@ -8130,6 +8130,10 @@
 				RelativePath="..\src\spicelib\devices\res\restemp.c"
 				>
 			</File>
+			<File
+				RelativePath="..\src\maths\misc\rnorrexp.c"
+				>
+			</File>
 			<File
 				RelativePath="..\src\frontend\runcoms.c"
 				>