Add a function m3avg(vector) for filtering of trap ringing.

Moving average with D(n) = (C(n-1)/2 + C(n) + C(n+1)/2)/2

src/frontend/parse.c

@@ -363,8 +363,9 @@ struct func ft_funcs[] = {
     { "ceil",        cx_ceil },
     { "mean",        cx_mean },
     { "stddev",      cx_stddev },
-    { "avg",         cx_avg }, /* A.Roldan 03/06/05 incremental average new function */
+    { "avg",         cx_avg },
-    { "group_delay", (cx_function_t*)(void *) cx_group_delay }, /* A.Roldan 10/06/05 group delay new function */
+    { "m3avg",       cx_m3avg },
+    { "group_delay", (cx_function_t*)(void *) cx_group_delay },
     { "vector",      cx_vector },
     { "cvector",     cx_cvector },
     { "unitvec",     cx_unitvec },

src/include/ngspice/fteext.h

@@ -103,7 +103,7 @@ extern void *cx_max(void *, short int , int , int *, short int *);
 extern void *cx_min(void *, short int , int , int *, short int *);
 extern void *cx_d(void *, short int , int , int *, short int *);
 extern void *cx_avg(void *, short int , int , int *, short int *);
-
+extern void *cx_m3avg(void *, short int , int , int *, short int *);
 
 /* cmath3.c */
 

src/maths/cmaths/cmath2.c

@@ -361,6 +361,62 @@ cx_avg(void *data, short int type, int length, int *newlength, short int *newtyp
     }
 }
 
+/* Compute the moving average of a vector.
+ * Using three elements ((n-1)/2 + n + (n+1)/2)/2
+ * Useful for removing trap ringing.
+ */
+
+void*
+cx_m3avg(void* data, short int type, int length, int* newlength, short int* newtype)
+{
+    double sum_real = 0.0, sum_imag = 0.0;
+    int i;
+
+    if (type == VF_REAL) {
+
+        double* d = alloc_d(length);
+        double* dd = (double*)data;
+
+        int nlen = length - 1;
+
+        *newtype = VF_REAL;
+        *newlength = length;
+
+        d[0] = dd[0];
+
+        for (i = 1; i < nlen; i++) {
+            d[i] = (dd[i-1] + dd[i+1]) / 4. + dd[i] / 2.;
+        }
+
+        d[nlen] = dd[nlen];
+
+        return ((void*)d);
+
+    }
+    else {
+
+        ngcomplex_t* c = alloc_c(length);
+        ngcomplex_t* cc = (ngcomplex_t*)data;
+
+        int nlen = length - 1;
+
+        *newtype = VF_COMPLEX;
+        *newlength = length;
+
+        realpart(c[0]) = realpart(cc[0]);
+
+        for (i = 0; i < length; i++) {
+            realpart(c[i]) = (realpart(cc[i - 1]) + realpart(cc[i + 1])) / 4. + realpart(cc[i]) / 2.;
+            imagpart(c[i]) = (imagpart(cc[i - 1]) + imagpart(cc[i + 1])) / 4. + imagpart(cc[i]) / 2.;
+        }
+
+        realpart(c[nlen]) = realpart(cc[nlen]);
+        imagpart(c[nlen]) = imagpart(cc[nlen]);
+
+        return ((void*)c);
+
+    }
+}
 
 /* Compute the mean of a vector. */
 

src/maths/cmaths/cmath2.h

@@ -33,6 +33,7 @@ void * cx_max(void *data, short int type, int length, int *newlength, short int
 void * cx_min(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_d(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_avg(void *data, short int type, int length, int *newlength, short int *newtype);
+void * cx_m3avg(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_floor(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_ceil(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_nint(void *data, short int type, int length, int *newlength, short int *newtype);