From 474a1e47b85b4ee6fabdd71bb9faede18fb53d83 Mon Sep 17 00:00:00 2001
From: Holger Vogt <holger.vogt@uni-due.de>
Date: Sun, 3 Jul 2022 11:28:00 +0200
Subject: [PATCH] =?UTF-8?q?Avoid=20jump=20in=20phase=20when=20crossing=209?=
 =?UTF-8?q?0=C2=B0=20boundaries.=20Code=20similar=20to=20cx=5Fcph=20in=20c?=
 =?UTF-8?q?math1.c?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/maths/cmaths/cmath4.c | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/maths/cmaths/cmath4.c b/src/maths/cmaths/cmath4.c
index afdc25d55..475559456 100644
--- a/src/maths/cmaths/cmath4.c
+++ b/src/maths/cmaths/cmath4.c
@@ -496,19 +496,23 @@ cx_group_delay(void *data, short int type, int length, int *newlength, short int
         return (NULL);
     }
 
-
-    if (type == VF_COMPLEX)
-        for (i = 0; i < length; i++)
-        {
-            v_phase[i] = radtodeg(cph(cc[i]));
+    if (type == VF_COMPLEX) {
+        /*  accept continuous phase over 90° boundaries */
+        double last_ph = cph(cc[0]);
+        v_phase[0] = radtodeg(last_ph);
+        for (i = 1; i < length; i++) {
+            double ph = cph(cc[i]);
+            last_ph = ph - (2 * M_PI) * floor((ph - last_ph) / (2 * M_PI) + 0.5);
+            v_phase[i] = radtodeg(last_ph);
+//            fprintf(stderr, "v_phase %e, cc %e %e\n", v_phase[i], cc[i].cx_real, cc[i].cx_imag);
         }
+    }
     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);