From f543aa8c0698a9b66c6fcb858b025d9647054178 Mon Sep 17 00:00:00 2001
From: pnenzi <pnenzi>
Date: Thu, 15 Jan 2009 21:08:09 +0000
Subject: [PATCH] Improvements in vector derivative computation and new
 functions to compute group delay and moving average. From espice (A. Roldan).

---
 ChangeLog                 |  17 +++
 src/frontend/cpitf.c      |   4 +-
 src/frontend/evaluate.c   |   4 +-
 src/frontend/parse.c      |   2 +
 src/frontend/postcoms.c   |  62 ++++++++--
 src/frontend/spiceif.c    |   2 +-
 src/frontend/vectors.c    |   3 +-
 src/include/fteext.h      |   3 +
 src/maths/cmaths/cmath2.c |  45 ++++++++
 src/maths/cmaths/cmath4.c | 234 +++++++++++++++++++++++++++++---------
 src/maths/cmaths/cmath4.h |   2 +
 11 files changed, 312 insertions(+), 66 deletions(-)

diff --git a/ChangeLog b/ChangeLog
index 33564f21f..be2a447d6 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,20 @@
+2009-01-15 Paolo Nenzi
+ * src/frontend/vectors.c:
+   57: Fixed  I(vx), before the if I(*) (upper case) was not recognized as the
+       function to plot the current of vx. A. Roldan - Espice
+ * src/frontend/postcoms.c:
+   7:  Fixed plot number after "destroy all" command. A. Roldan - Espice 
+ * src/include/fteext.h
+ * src/frontend/evaluate.c, src/frontend/cpitf.c, src/frontend/parse.c 
+ * src/maths/cmaths/cmath2.c, src/cmaths/cmath4.c, src/cmaths/cmath4.h:
+   16: New function to compute the group delay has been implemented.
+       Group delay is defined as -(dphase/dfrequency) and can be printed or
+       plotted by writing vg(x), where x is a complex vector. A. Roldan - Espice
+   15: Fixed existing problems in this function due to the complex nature of 
+       the frequency vector. To get the data from frequency[i], the real part
+       must be accessed. A. Roldan - Espice
+   14: New function to compute the moving average. A. Roldan - Espice
+
 2009-01-15 Paolo Nenzi
  * src/frontend/{spiceif.c, spiceif.h, vectors.c}, src/include/fteext.h,
  * src/main.c: 
diff --git a/src/frontend/cpitf.c b/src/frontend/cpitf.c
index b880b009f..23bd459ce 100644
--- a/src/frontend/cpitf.c
+++ b/src/frontend/cpitf.c
@@ -55,7 +55,9 @@ ft_cpinit(void)
         "vi(x,y)",  "im(v(x) - v(y))",
         "vm(x)",    "mag(v(x))",
         "vm(x,y)",  "mag(v(x) - v(y))",
-        "vp(x)",    "ph(v(x))",
+        "vg(x)",    "group_delay(v(x))", //A.Rroldan 10/06/05 group delay new function
+        "gd(x)",    "group_delay(v(x))", //A.Rroldan 10/06/05 group delay new function        
+	    "vp(x)",    "ph(v(x))",
         "vp(x,y)",  "ph(v(x) - v(y))",
         "vr(x)",    "re(v(x))",
         "vr(x,y)",  "re(v(x) - v(y))"
diff --git a/src/frontend/evaluate.c b/src/frontend/evaluate.c
index 339ab8f23..a4d0bee53 100644
--- a/src/frontend/evaluate.c
+++ b/src/frontend/evaluate.c
@@ -707,8 +707,8 @@ apply_func(struct func *func, struct pnode *arg)
         }
         (void) signal(SIGILL, (SIGNAL_FUNCTION) sig_matherr);
 
-        if (eq(func->fu_name, "interpolate")
-            || eq(func->fu_name, "deriv"))       /* Ack */
+        /* Modified for passing necessary parameters to the derive function - A.Roldan */
+	if (eq(func->fu_name, "interpolate") || eq(func->fu_name, "deriv") || eq(func->fu_name, "group_delay"))       /* Ack */
 	{
 	    void *(*f)(void *data, short int type, int length,
                        int *newlength, short int *newtype, ...)=func->fu_func;
diff --git a/src/frontend/parse.c b/src/frontend/parse.c
index a249f0e1c..69f806185 100644
--- a/src/frontend/parse.c
+++ b/src/frontend/parse.c
@@ -698,6 +698,8 @@ struct func ft_funcs[] = {
         { "rnd",    cx_rnd } ,
         { "pos",    cx_pos } ,
         { "mean",   cx_mean } ,
+        { "avg",   cx_avg } ,     //A.Rroldan 03/06/05 incremental average  new function
+        { "group_delay",  cx_group_delay } , //A.Rroldan 10/06/05 group delay new function
         { "vector", cx_vector } ,
         { "unitvec",    cx_unitvec } ,
         { "length", cx_length } ,
diff --git a/src/frontend/postcoms.c b/src/frontend/postcoms.c
index 4a8ac6fc5..9cde88e82 100644
--- a/src/frontend/postcoms.c
+++ b/src/frontend/postcoms.c
@@ -260,7 +260,22 @@ nextpage:
             if (isreal(v))
                 (void) sprintf(buf2, "%-16.15s", v->v_name);
             else
-                (void) sprintf(buf2, "%-32.31s", v->v_name);
+		{
+         /* The frequency vector is complex but often with imaginary part = 0, 
+          * this prevents to print two columns.
+          */
+		 if(eq(v->v_name, "frequency"))
+		 {
+		  if(imagpart(&v->v_compdata[0])==0.0)
+		  {
+                   (void) sprintf(buf2, "%-16.15s", v->v_name);
+		  }
+		  else
+                   (void) sprintf(buf2, "%-32.31s", v->v_name);
+		 }
+		 else
+                  (void) sprintf(buf2, "%-32.31s", v->v_name);
+		}
             (void) strcat(buf, buf2);   
         }
         lineno = 3;
@@ -290,14 +305,34 @@ loop:
                     else
                         out_send("\t\t\t\t");
                 } else {
-                    if (isreal(v)) {
-                        sprintf(out_pbuf, "%e\t", 
-                        	v->v_realdata[j]);
+                    if (isreal(v)) 
+		    {
+                        printnum(numbuf,  v->v_realdata[j]);
+                        //sprintf(out_pbuf, "%e\t",v->v_realdata[j]);
+                        (void) sprintf(out_pbuf, "%s\t",numbuf);
 			out_send(out_pbuf);
-                    } else {
-                        sprintf(out_pbuf, "%e,\t%e\t",
-                        	realpart(&v->v_compdata[j]),
-                        	imagpart(&v->v_compdata[j]));
+                    }
+		    else
+		    {
+            /* In case of a single frequency and have a real part avoids print imaginary part equals 0. */
+			if(eq(v->v_name, "frequency"))
+			{
+			 if(imagpart(&v->v_compdata[j])==0.0)
+			 {
+                          printnum(numbuf,  realpart(&v->v_compdata[j]));
+                          (void) sprintf(out_pbuf, "%s\t",numbuf);
+			 }
+			}
+			else
+			{
+                          printnum(numbuf,  realpart(&v->v_compdata[j]));
+                          printnum(numbuf2, imagpart(&v->v_compdata[j]));
+                          (void) sprintf(out_pbuf, "%s,\t%s\t",numbuf,numbuf2);
+/*                       sprintf(out_pbuf, "%e,\t%e\t",
+ *                        	realpart(&v->v_compdata[j]),
+ *                        	imagpart(&v->v_compdata[j]));
+ */
+			}
 			out_send(out_pbuf);
 		    }
                 }
@@ -582,7 +617,13 @@ com_destroy(wordlist *wl)
         for (pl = plot_list; pl; pl = npl) {
             npl = pl->pl_next;
             if (!eq(pl->pl_typename, "const"))
-                killplot(pl);
+	    {
+             killplot(pl);
+	    }
+	    else
+	    {
+	     plot_num=1;
+	    } 
         }
     } else {
         while (wl) {
@@ -590,7 +631,10 @@ com_destroy(wordlist *wl)
                 if (eq(pl->pl_typename, wl->wl_word))
                     break;
             if (pl)
+	    {
                 killplot(pl);
+		plot_num--;
+	    }
             else
                 fprintf(cp_err, "Error: no such plot %s\n",
                         wl->wl_word);
diff --git a/src/frontend/spiceif.c b/src/frontend/spiceif.c
index 19f1e205e..3b9437057 100644
--- a/src/frontend/spiceif.c
+++ b/src/frontend/spiceif.c
@@ -675,7 +675,7 @@ spif_getparam_special(void *ckt,char **name,char *param,int ind,int do_model)
 		 * 2. Copy the type of variable of IFparm into a variable (thus parm-to-var)
 		 * vv->va_type = opt->dataType
 		 * The long description of the parameter:
-		 * IFparm MOS_SGTmPTable[] = { /* model parameters */
+		 * IFparm MOS_SGTmPTable[] = { // model parameters //
 		 * OP("type",   MOS_SGT_MOD_TYPE,  IF_STRING, "N-channel or P-channel MOS") 
          * goes into tv->va_name to put braces around the parameter of the model
 		 * tv->va_name += device->modelParms[i].keyword;
diff --git a/src/frontend/vectors.c b/src/frontend/vectors.c
index 6cffa70b6..6f4b0b628 100644
--- a/src/frontend/vectors.c
+++ b/src/frontend/vectors.c
@@ -275,7 +275,8 @@ vec_fromplot(char *word, struct plot *plot)
         /* ( */ ((s =strrchr(buf, ')')) != NULL) &&
             (*(s + 1) == '\0')) {
         *s = '\0';
-        if (prefix("i(", /* ) */ word)) {
+	if (prefix("i(", /* ) */ word) || prefix("I(", /* ) */ word))
+	{
             /* Spice dependency... */
             (void) sprintf(buf2, "%s#branch", buf);
             (void) strcpy(buf, buf2);
diff --git a/src/include/fteext.h b/src/include/fteext.h
index 49d5dfde6..228b0f2fc 100644
--- a/src/include/fteext.h
+++ b/src/include/fteext.h
@@ -100,6 +100,7 @@ extern void *cx_norm(void *, short int , int , int *, short int *, ...);
 extern void *cx_uminus(void *, short int , int , int *, short int *, ...);
 extern void *cx_rnd(void *, short int , int , int *, short int *, ...);
 extern void *cx_mean(void *, short int , int , int *, short int *, ...);
+extern void *cx_avg(void *, short int , int , int *, short int *, ...);
 extern void *cx_length(void *, short int , int , int *, short int *, ...);
 extern void *cx_vector(void *, short int , int , int *, short int *, ...);
 extern void *cx_unitvec(void *, short int , int , int *, short int *, ...);
@@ -133,6 +134,8 @@ extern void *cx_or(void *, void *, short int , short int , int, ...);
 extern void *cx_not(void *, short int , int , int *, short int * , ...);
 extern void *cx_interpolate(void *, short int , int , int *, short int *, ...); /* struct plot *, struct plot *, int ); */
 extern void *cx_deriv(void *, short int , int , int *, short int *, ...); /*struct plot *, struct plot *, int );*/
+extern void *cx_group_delay(void *, short int , int , int *, short int *, ...); /*struct plot *, struct plot *, int );*/
+
 
 /* cmdtab.c */
 
diff --git a/src/maths/cmaths/cmath2.c b/src/maths/cmaths/cmath2.c
index 67cd6b838..cc1b0fc6b 100644
--- a/src/maths/cmaths/cmath2.c
+++ b/src/maths/cmaths/cmath2.c
@@ -244,6 +244,51 @@ cx_rnd(void *data, short int type, int length, int *newlength, short int *newtyp
   }
 }
 
+/* Compute the avg of a vector.
+   Created by A.M.Roldan 2005-05-21  */
+   
+void 
+*cx_avg(void *data, short int type, int length, int *newlength, short int *newtype, ...)
+{
+    complex *c;
+    double *d, sum_real = 0.0,sum_imag = 0.0;
+    complex *cc = (complex *) data;
+    double *dd = (double *) data;
+    int i;
+
+    if (type == VF_REAL) {
+        d = alloc_d(length);
+        *newtype = VF_REAL;
+    } else {
+        c = alloc_c(length);
+        *newtype = VF_COMPLEX;
+    }
+    *newlength = length;
+
+    if (type == VF_COMPLEX)
+    {
+       for (i = 0; i < length; i++)
+       {
+        sum_real= sum_real + realpart(&cc[i]);
+        realpart(&c[i]) = sum_real / (double)(i+1);
+
+        sum_imag = sum_imag + imagpart(&cc[i]);
+        imagpart(&c[i]) = sum_imag / (double)(i+1);
+       }
+       return ((char *) c);
+    }
+    else
+    {      //VF_REAL
+       for (i = 0; i < length; i++)
+       {
+        sum_real= sum_real + dd[i];
+        d[i] = sum_real / (double)(i+1);
+       }
+       return ((char *) d);
+    }
+}
+
+
 /* Compute the mean of a vector. */
 
 void *
diff --git a/src/maths/cmaths/cmath4.c b/src/maths/cmaths/cmath4.c
index c9501f887..ac6cf2567 100644
--- a/src/maths/cmaths/cmath4.c
+++ b/src/maths/cmaths/cmath4.c
@@ -28,7 +28,10 @@ Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group
 #include "cmath.h"
 #include "cmath4.h"
 
+#include "sim.h" /* To get SV_TIME */
+
 #include "../../frontend/variable.h" /* for VT_NUM in cx_interpolate */
+extern bool cx_degrees;
 
 void *
 cx_and(void *data1, void *data2, short int datatype1, short int datatype2, int length)
@@ -262,67 +265,75 @@ cx_deriv(void *data, short int type, int length, int *newlength, short int *newt
 	scale = alloc_d(length);	/* XXX */
 	if (pl->pl_scale->v_type == VF_COMPLEX)
 	    /* Not ideal */
+	  for (i = 0; i < length; i++)
+		   scale[i] = realpart(&pl->pl_scale->v_compdata[i]);
+	 else
 	    for (i = 0; i < length; i++)
-		scale[i] = realpart(&pl->pl_scale->v_compdata[i]);
-	else
-	    for (i = 0; i < length; i++)
-		scale[i] = pl->pl_scale->v_realdata[i];
-	for (base = 0; base < length; base += grouping) {
+		   scale[i] = pl->pl_scale->v_realdata[i];
+
+	 for (base = 0; base < length; base += grouping)
+    {
 	    k = 0;
-	    for (i = degree; i < grouping; i += 1) {
+	    for (i = degree; i < grouping; i += 1)
+       {
 
-		/* real */
-		for (j = 0; j < n; j++)
+		  /* real */
+		  for (j = 0; j < n; j++)
 		    spare[j] = c_indata[j + i + base].cx_real;
-		if (!ft_polyfit(scale + i + base - degree,
+		  if (!ft_polyfit(scale + i + base - degree,
 		    spare, r_coefs, degree, scratch))
-		{
+		   {
 		    fprintf(stderr, "ft_polyfit @ %d failed\n", i);
-		}
-		ft_polyderiv(r_coefs, degree);
+		   }
+		  ft_polyderiv(r_coefs, degree);
 
-		/* for loop gets the beginning part */
-		for (j = k; j <= i - degree / 2; j++) {
+		  /* for loop gets the beginning part */
+		  for (j = k; j <= i + degree / 2; j++)
+        {
 		    x = scale[j + base];
 		    c_outdata[j + base].cx_real =
 			ft_peval(x, r_coefs, degree - 1);
-		}
+		  }
 
-		/* imag */
-		for (j = 0; j < n; j++)
+		  /* imag */
+		  for (j = 0; j < n; j++)
 		    spare[j] = c_indata[j + i + base].cx_imag;
-		if (!ft_polyfit(scale + i - degree + base,
+		  if (!ft_polyfit(scale + i - degree + base,
 		    spare, i_coefs, degree, scratch))
-		{
+		  {
 		    fprintf(stderr, "ft_polyfit @ %d failed\n", i);
-		}
-		ft_polyderiv(i_coefs, degree);
+		  }
+		  ft_polyderiv(i_coefs, degree);
 
-		/* for loop gets the beginning part */
-		for (j = k; j <= i - degree / 2; j++) {
+		  /* for loop gets the beginning part */
+        for (j = k; j <= i - degree / 2; j++)
+        {
 		    x = scale[j + base];
 		    c_outdata[j + base].cx_imag =
-			ft_peval(x, i_coefs, degree - 1);
-		}
-		k = j;
-	    }
+		    ft_peval(x, i_coefs, degree - 1);
+		  }
+		 k = j;
+	   }
 
 	    /* get the tail */
-	    for (j = k; j < length; j++) {
-		x = scale[j + base];
-		/* real */
-		c_outdata[j + base].cx_real = ft_peval(x, r_coefs, degree - 1);
-		/* imag */
-		c_outdata[j + base].cx_imag = ft_peval(x, i_coefs, degree - 1);
+	    for (j = k; j < length; j++)
+       {
+		  x = scale[j + base];
+		  /* real */
+		  c_outdata[j + base].cx_real = ft_peval(x, r_coefs, degree - 1);
+		  /* imag */
+		  c_outdata[j + base].cx_imag = ft_peval(x, i_coefs, degree - 1);
 	    }
-	}
+    }
 
 	tfree(r_coefs);
 	tfree(i_coefs);
 	tfree(scale);
 	return (void *) c_outdata;
 
-    } else {
+  }
+  else
+  {
 	/* all-real case */
 	double *coefs;
 
@@ -333,36 +344,155 @@ cx_deriv(void *data, short int type, int length, int *newlength, short int *newt
 	indata = (double *) data;
 	outdata = alloc_d(length);
 	scale = alloc_d(length);	/* XXX */
-	for (i = 0; i < length; i++)
-	    scale[i] = pl->pl_scale->v_realdata[i];
-	for (base = 0; base < length; base += grouping) {
+
+   /* Here I encountered a problem because when we issue an instruction like this:
+    * plot -deriv(vp(3)) to calculate something similar to the group delay, the code
+    * detects that vector vp(3) is real and it is believed that the frequency is also 
+    * real. The frequency is COMPLEX and the program aborts so I'm going to put the 
+    * check that the frequency is complex vector not to abort.
+    */
+
+
+   /* Original problematic code
+	* for (i = 0; i < length; i++)
+	*    scale[i] = pl->pl_scale->v_realdata[i];
+    */
+
+   /* Modified to deal with complex frequency vector */
+   if (pl->pl_scale->v_type == VF_COMPLEX)
+	  for (i = 0; i < length; i++)
+		   scale[i] = realpart(&pl->pl_scale->v_compdata[i]);
+	else
+	    for (i = 0; i < length; i++)
+		   scale[i] = pl->pl_scale->v_realdata[i];
+
+
+
+	for (base = 0; base < length; base += grouping)
+   {
 	    k = 0;
-	    for (i = degree; i < grouping; i += 1) {
-		if (!ft_polyfit(scale + i - degree + base,
+	    for (i = degree; i < grouping; i += 1)
+       {
+		  if (!ft_polyfit(scale + i - degree + base,
 		    indata + i - degree + base, coefs, degree, scratch))
-		{
+		   {
 		    fprintf(stderr, "ft_polyfit @ %d failed\n", i + base);
-		}
-		ft_polyderiv(coefs, degree);
+	    	}
+        ft_polyderiv(coefs, degree);
+
+		  /* for loop gets the beginning part */
+		  for (j = k; j <= i - degree / 2; j++)
+        {
+          /* Seems the same problem because the frequency vector is complex
+           * and the real part of the complex should be accessed because if we 
+           * run x = pl-> pl_scale-> v_realdata [base + j]; the execution will 
+           * abort.
+           */
+
+          if (pl->pl_scale->v_type == VF_COMPLEX)
+            x = realpart(&pl->pl_scale->v_compdata[j+base]);  /* For complex scale vector */
+          else
+            x = pl->pl_scale->v_realdata[j + base];           /* For real scale vector */
 
-		/* for loop gets the beginning part */
-		for (j = k; j <= i - degree / 2; j++) {
-		    x = pl->pl_scale->v_realdata[j + base];
 		    outdata[j + base] = ft_peval(x, coefs, degree - 1);
-		}
-		k = j;
+		  }
+	    	k = j;
 	    }
 
-	    for (j = k; j < length; j++) {
-		    x = pl->pl_scale->v_realdata[j + base];
+	    for (j = k; j < length; j++)
+       {
+          /* Again the same error */
+		  /* x = pl->pl_scale->v_realdata[j + base]; */
+          if (pl->pl_scale->v_type == VF_COMPLEX)
+            x = realpart(&pl->pl_scale->v_compdata[j+base]);  /* For complex scale vector */
+          else
+            x = pl->pl_scale->v_realdata[j + base];           /* For real scale vector */
+
 		    outdata[j + base] = ft_peval(x, coefs, degree - 1);
 	    }
-        }
+   }
 
 
 	tfree(coefs);
 	tfree(scale);	/* XXX */
-	return (void *) outdata;
+	return (char *) outdata;
+ }
+
+}
+
+
+void *
+cx_group_delay(void *data, short int type, int length, int *newlength, short int *newtype, struct plot *pl, struct plot *newpl, int grouping)
+{
+    complex *cc = (complex *) data;
+    double *v_phase = alloc_d(length);
+    double *datos,adjust_final;
+    double *group_delay = alloc_d(length);
+    int i;
+    /* char *datos_aux; */
+
+    /* Check to see if we have the frequency vector for the derivative */
+    if (!eq(pl->pl_scale->v_name, "frequency"))
+    {
+      fprintf(cp_err, "Internal error: cx_group_delay: need frequency based complex vector.\n");
+      return (NULL);
     }
 
+
+    if (type == VF_COMPLEX)
+     for (i = 0; i < length; i++)
+     {
+      v_phase[i] = radtodeg(cph(&cc[i]));
+     }
+    else
+    {
+      fprintf(cp_err, "Signal must be complex to calculate group delay\n");
+      return (NULL);
+    }
+
+
+    type = VF_REAL;
+
+    /* datos_aux = (char *)cx_deriv((char *)v_phase, type, length, newlength, newtype, pl, newpl, grouping);
+     * datos = (double *) datos_aux;
+     */
+    datos = (double *)cx_deriv((char *)v_phase, type, length, newlength, newtype, pl, newpl, grouping);
+    
+    /* With this global variable I will change how to obtain the group delay because
+     * it is defined as:
+     *
+     *  gd()=-dphase[rad]/dw[rad/s]
+     *
+     * if you have degrees in phase and frequency in Hz, must be taken into account
+     *
+     *  gd()=-dphase[deg]/df[Hz]/360
+     *  gd()=-dphase[rad]/df[Hz]/(2*pi)
+     */
+
+    if(cx_degrees)
+     {
+       adjust_final=1.0/360;
+     }
+     else
+     {
+       adjust_final=1.0/(2*M_PI);
+     }
+
+
+    for (i = 0; i < length; i++)
+    {
+     group_delay[i] = -datos[i]*adjust_final;
+    }
+
+    /* Adjust to Real because the result is Real */
+    *newtype = VF_REAL;
+
+    	
+    /* Set the type of Vector to "Time" because the speed of group units' s'
+     * The different types of vectors are INCLUDE \ Fte_cons.h
+     */
+    pl->pl_dvecs->v_type= SV_TIME;
+
+   return ((char *) group_delay);
+
 }
diff --git a/src/maths/cmaths/cmath4.h b/src/maths/cmaths/cmath4.h
index dcf55a18d..b08ef4258 100644
--- a/src/maths/cmaths/cmath4.h
+++ b/src/maths/cmaths/cmath4.h
@@ -13,6 +13,8 @@ void * cx_interpolate(void *data, short int type, int length, int *newlength,
 		      short int *newtype, struct plot *pl, struct plot *newpl, int grouping);
 void * cx_deriv(void *data, short int type, int length, int *newlength, short int *newtype, 
 		struct plot *pl, struct plot *newpl, int grouping);
+void * cx_group_delay(void *data, short int type, int length, int *newlength, short int *newtype, 
+		struct plot *pl, struct plot *newpl, int grouping);
 
 
 #endif