From ae80067fef4ac1794c5d51240ca60cc5baa161dc Mon Sep 17 00:00:00 2001
From: h_vogt <h_vogt>
Date: Sat, 17 Sep 2016 20:36:41 +0200
Subject: [PATCH] cmath2.c, new function stddev to calculate the standard
 deviation of all elements of a vector

---
 src/frontend/parse.c         |  1 +
 src/include/ngspice/fteext.h |  1 +
 src/maths/cmaths/cmath2.c    | 41 ++++++++++++++++++++++++++++++++++++
 src/maths/cmaths/cmath2.h    |  1 +
 4 files changed, 44 insertions(+)

diff --git a/src/frontend/parse.c b/src/frontend/parse.c
index 08c9c7e4d..691a660db 100644
--- a/src/frontend/parse.c
+++ b/src/frontend/parse.c
@@ -173,6 +173,7 @@ struct func ft_funcs[] = {
     { "floor",       cx_floor },
     { "ceil",        cx_ceil },
     { "mean",        cx_mean },
+    { "stddev",      cx_stddev },
     { "avg",         cx_avg }, /* A.Roldan 03/06/05 incremental average new function */
     { "group_delay", (cx_function_t*) cx_group_delay }, /* A.Roldan 10/06/05 group delay new function */
     { "vector",      cx_vector },
diff --git a/src/include/ngspice/fteext.h b/src/include/ngspice/fteext.h
index 7ef959c21..b59149c13 100644
--- a/src/include/ngspice/fteext.h
+++ b/src/include/ngspice/fteext.h
@@ -88,6 +88,7 @@ extern void *cx_sgauss(void *, short int , int , int *, short int *);
 extern void *cx_poisson(void *, short int , int , int *, short int *);
 extern void *cx_exponential(void *, short int , int , int *, short int *);
 extern void *cx_mean(void *, short int , int , int *, short int *);
+extern void *cx_stddev(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 *);
diff --git a/src/maths/cmaths/cmath2.c b/src/maths/cmaths/cmath2.c
index c27b3a5b6..e4ed30892 100644
--- a/src/maths/cmaths/cmath2.c
+++ b/src/maths/cmaths/cmath2.c
@@ -395,6 +395,47 @@ cx_mean(void *data, short int type, int length, int *newlength, short int *newty
 }
 
 
+/* Compute the standard deviation of all elements of a vector. */
+
+void *
+cx_stddev(void *data, short int type, int length, int *newlength, short int *newtype)
+{
+    *newlength = 1;
+    rcheck(length > 1, "stddev");
+    if (type == VF_REAL) {
+        double *mean = (double *)cx_mean(data, type, length, newlength, newtype);
+        double *d, sum = 0.;
+        double *dd = (double *)data;
+        int i;
+
+        d = alloc_d(1);
+        *newtype = VF_REAL;
+        for (i = 0; i < length; i++)
+            sum += (dd[i] - *mean) * (dd[i] - *mean);
+        *d = sqrt(sum / (length - 1));
+        tfree(mean);
+        return ((void *)d);
+    }
+    else {
+        ngcomplex_t *cmean = (ngcomplex_t *)cx_mean(data, type, length, newlength, newtype);
+        double *d, sum = 0., a, b;
+        ngcomplex_t *cc = (ngcomplex_t *)data;
+        int i;
+
+        d = alloc_d(1);
+        *newtype = VF_REAL;
+        for (i = 0; i < length; i++) {
+            a = realpart(cc[i]) - realpart(*cmean);
+            b = imagpart(cc[i]) - imagpart(*cmean);
+            sum += a * a + b * b;
+        }
+        *d = sqrt(sum / (length - 1));
+        tfree(cmean);
+        return ((void *)d);
+    }
+}
+
+
 void *
 cx_length(void *data, short int type, int length, int *newlength, short int *newtype)
 {
diff --git a/src/maths/cmaths/cmath2.h b/src/maths/cmaths/cmath2.h
index 4b71e2b12..54df8bbee 100644
--- a/src/maths/cmaths/cmath2.h
+++ b/src/maths/cmaths/cmath2.h
@@ -16,6 +16,7 @@ void * cx_sgauss(void *data, short int type, int length, int *newlength, short i
 void * cx_poisson(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_exponential(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_mean(void *data, short int type, int length, int *newlength, short int *newtype);
+void * cx_stddev(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_length(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_vector(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_unitvec(void *data, short int type, int length, int *newlength, short int *newtype);