From ebc5334c2ef4586850546d1bafd0c9fd16a93ae5 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Fri, 22 May 2026 10:16:50 +0300
Subject: [PATCH 01/16] add prime implementation for calculating roots from
 coefficients based on the GramSchmidt method.

Computes both real and complex roots

Cons: Inefficient as number of coefficients scale up. Optimizations and alternative methods that provide a more efficient way to handle input should be considered for providing a practical application of this feature
---
 src/CoeffComponents.cpp     |  45 +++++++-
 src/CoeffComponents.h       |   1 -
 src/CoefficientsToRoots.cpp | 205 ++++++++++++++++++++++++++++++++++++
 src/CoefficientsToRoots.h   |  18 ++++
 4 files changed, 267 insertions(+), 2 deletions(-)
 create mode 100644 src/CoefficientsToRoots.cpp
 create mode 100644 src/CoefficientsToRoots.h

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index a17b99f..5323169 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -1,4 +1,6 @@
 #include "CoeffComponents.h"
+#include "RootsToCoefficients.h"
+#include "CoefficientsToRoots.cpp"
 
 #include <iostream>
 #define COL_WIDTH 100 // TODO make this value relate to the global UI.
@@ -118,9 +120,50 @@ juce::Component* CoefficientsComponent::refreshComponentForCell(int row, int col
                 ffcoeffs[static_cast<size_t>(row)] = value;
             else if (col == 3)
                 fbcoeffs[static_cast<size_t>(row)] = value;
-
+            
             // TODO: calculate roots through the coefficient2roots function.
             // TODO: notify other listeners about this change
+            if (col == 2)
+            {
+
+                auto zeros = CoefficientsToRoots::GramSchmidt(this->ffcoeffs);
+
+                std::cout<<"Existing Zeros"<<std::endl;
+                for (int i=0 ; i< processor->filterState->zeros.size(); i++)
+                {
+                    auto *r = processor->filterState->zeros[i];
+                    std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+                    // processor->filterState->remove(r);
+                }
+                std::cout<<std::endl;
+
+                // for (int i=0 ; i< zeros.size(); i++)
+                // {
+                //     auto &r = zeros[i];
+                //     processor->filterState->add(i, r);
+                // }
+
+            }
+            else if (col == 3)
+            {
+                auto poles =  CoefficientsToRoots::GramSchmidt(this->fbcoeffs);
+
+                std::cout<<"Existing Poles"<<std::endl;
+                for (int i=0 ; i< processor->filterState->poles.size(); i++)
+                {
+                    auto *r = processor->filterState->poles[i];
+                    std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+                    // processor->filterState->remove(r);
+                }
+                std::cout<<std::endl;
+
+                // for (int i=0 ; i< poles.size(); i++)
+                // {
+                //     auto &r = poles[i];
+                //     processor->filterState->add(i, r);
+                // }
+            }
+            
         };
     }
 
diff --git a/src/CoeffComponents.h b/src/CoeffComponents.h
index e896f17..b1aa6c3 100644
--- a/src/CoeffComponents.h
+++ b/src/CoeffComponents.h
@@ -1,7 +1,6 @@
 #pragma once
 
 #include "FilterState.h"
-#include "RootsToCoefficients.h"
 #include <vector>
 #include <string>
 
diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
new file mode 100644
index 0000000..c07d208
--- /dev/null
+++ b/src/CoefficientsToRoots.cpp
@@ -0,0 +1,205 @@
+// #pragma once
+// #include <juce_audio_processors/juce_audio_processors.h>
+// #include <juce_dsp/juce_dsp.h>
+
+#include "CoefficientsToRoots.h"
+
+#include <algorithm>
+#include <iostream>
+#include <cmath>
+
+std::vector<c128> CoefficientsToRoots::GramSchmidt(std::vector<double> coefs)
+{
+    // filter out leading zeros and leading 1.0
+    int degree = 0;
+    while(coefs[degree] ==0 || coefs[degree] == 1.0)
+        degree++;
+    degree = coefs.size() - degree;
+    jassert( degree!=0 );
+
+    if (degree==1)
+    {
+        std::cout<<"Returing 1 root --> "<<coefs[degree]<<std::endl;
+        return {static_cast<c128>(coefs[degree])};
+    }
+
+    // build companion Matrix
+    std::cout<<"Building companion matrix with degree "<<degree<<std::endl;
+    std::vector<std::vector<double>> A(degree, std::vector<double>(degree, 0.0));
+    // Build companion matrix (Hessenberg form)
+    for (int i=0; i< degree; i++)
+    {
+        int j = i+1;
+        int coef_idx = degree-i;
+        A[i][degree-1] = -coefs[coef_idx];  // fill last column with negative vals of coefs
+        if (i!=degree-1)
+        {
+            A[j][i] = 1.0;  // fill subdiagonal entries         
+        }
+    }
+
+        
+    std::cout<<"// Check companion matrix"<<std::endl;
+    printCheck(A);
+
+
+    // QR algorithm
+    std::vector<std::vector<double>> Q(degree, std::vector<double>(degree));
+    std::vector<std::vector<double>> R(degree, std::vector<double>(degree));
+    
+    int iter = 0;
+    while(true) // could alternatively check the convergence of the values of the sub-diagonal elements
+    {
+
+        if (++iter > degree * MaxIterations) break;
+
+
+        // Reset R,Q
+        for (int r = 0; r < degree; ++r)
+        {
+            for (int c = 0; c < degree; ++c)
+            {
+                R[r][c] = 0.0;
+                Q[r][c] = 0.0;
+            }
+        }
+
+        // step 1 - calculate Q : Q = A[col] - projections onto the previous Q[col]
+        for (int col = 0 ; col < degree; col++)
+        {
+            // set v with column A[col]
+            std::vector<double> v(degree);
+            for (int row = 0 ; row< degree; row++)
+                v[row] = A[row][col];
+
+            // subtract projection : v[col] = A[col] - proj onto the previous (<col) orthonormal colums of Q.
+            for (int curr_col = 0; curr_col < col; curr_col++)
+            {
+                // compute dot product ...
+                double dot_product {0.0};
+                for (int curr_row = 0; curr_row < degree; ++curr_row)
+                {
+                    dot_product += Q[curr_row][curr_col] * A[curr_row][col];
+                }
+
+                // .. and subtract it from the current column vector
+                for (int curr_row=0; curr_row < degree; ++curr_row)
+                {
+                    v[curr_row] -= dot_product * Q[curr_row][curr_col];
+                }
+            }
+
+            // normalize column of q
+            double norm = 0.0;
+            for (int k = 0; k < degree; ++k)
+            {
+                norm += v[k] * v[k];
+            }
+            norm = std::sqrt(norm);
+            
+            for (int k = 0; k < degree; ++k)
+            {
+                Q[k][col] = v[k] / norm;
+            }
+        }
+
+        // step 2 - calculate R :  R = Q A
+        for (int row = 0; row < degree; row++)
+        {
+            for (int col = 0; col < degree; col++)
+            {
+                for (int inner = 0; inner < degree; inner++)
+                {
+                    R[row][col] += (Q[inner][row] * A[inner][col]);     // Q transposed
+                }
+            }
+        }
+
+        // step 3 - A = Q R
+        for (int row = 0; row < degree; row++)
+        {
+            for (int col = 0; col < degree; col++)
+            {
+                double sum {0.0};
+                for (int inner = 0; inner< degree; ++inner)
+                {
+                    sum += Q[inner][col] * R[row][inner];
+                }
+                A[row][col] = sum; // resetting A[row][col]
+            }
+        }
+        
+        // step 4 sum all sub-diagonal entries
+        double subdiag_sum {0.0};
+        for (int i = 1; i < degree; ++i)
+        {
+            subdiag_sum += std::abs(A[i][i-1]);
+        }
+
+        if (subdiag_sum < Epsilon)
+            break;
+    }
+    
+    // Extract roots — read diagonal
+    return extractRoots(A, degree);
+}
+ 
+void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &matrix)
+{
+    int n = matrix.size();
+    for (int i=0; i< n; i++)
+    {
+        for (int j=0; j<n; j++)
+        {
+            std::cout<<matrix[i][j]<<" ";
+        }
+        std::cout<<std::endl;
+    }
+}
+
+
+std::vector<c128> CoefficientsToRoots::extractRoots(const std::vector<std::vector<double>>& M, int degree)
+{
+    // https://www.mosismath.com/Eigenvalues/EigenvalsBasics.html --> public void CalcEigenvals()
+    std::vector<c128> roots;
+ 
+    int i = 0;
+    while (i < degree)
+    {
+        double a, b, c;
+        a = 1.0;
+        b = -M[i][i] - M[i+1][i+1];
+        c = M[i][i] * M[i+1][i+1] - M[i][i+1] * M[i+1][i];
+
+        if ((b * b) >= (4.0 * a * c))
+        {
+            roots.emplace_back(
+                (-b + std::sqrt((b * b) - (4.0 * a * c))) / 2.0 / a, 
+                0.0);
+            roots.emplace_back(
+                (-b - std::sqrt((b * b) - (4.0 * a * c))) / 2.0 / a,
+                0.0);
+            ++i;
+        }
+        else
+        {
+            double realPart = -b / 2.0 / a;
+            double imagPart =
+                std::sqrt((4.0 * a * c) - (b * b))
+                / 2.0 / a;
+
+            roots.emplace_back(realPart,  imagPart);
+            roots.emplace_back(realPart, -imagPart);
+
+            i += 2;
+        }
+    }
+
+    std::cout<<"Calculated Roots: ";
+    for (auto root : roots)
+        std::cout<<"("<<root.real()<<","<<root.imag()<<")  ";
+    std::cout<<std::endl;
+
+    return roots;
+
+}
\ No newline at end of file
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
new file mode 100644
index 0000000..3b6519f
--- /dev/null
+++ b/src/CoefficientsToRoots.h
@@ -0,0 +1,18 @@
+#pragma once
+
+#include "FilterState.h"
+#include <vector>
+
+class CoefficientsToRoots
+{
+	public:
+		static std::vector<c128> GramSchmidt(std::vector<double> coefs);
+	private:
+
+		// TODO Finetune these parameters
+	    static constexpr double Epsilon = 1e-3;
+	    static constexpr int MaxIterations = 100000;
+
+		static void printCheck(const std::vector<std::vector<double>> &matrix);
+		static std::vector<c128> extractRoots(const std::vector<std::vector<double>>& M, int n);
+};
\ No newline at end of file

From 1241d5292cd0ee01adf8aa52261fecffe707de7d Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Fri, 22 May 2026 10:22:35 +0300
Subject: [PATCH 02/16] improve extract roots method, symplifying it and making
 it more readable and maintainable. + added a todo as a comment for improving
 the data structure that stores the extracted roots in the upcoming updates

---
 src/CoefficientsToRoots.cpp | 66 +++++++++++++++++++++++--------------
 1 file changed, 42 insertions(+), 24 deletions(-)

diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index c07d208..cc7ae62 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -160,37 +160,56 @@ void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &mat
 
 std::vector<c128> CoefficientsToRoots::extractRoots(const std::vector<std::vector<double>>& M, int degree)
 {
-    // https://www.mosismath.com/Eigenvalues/EigenvalsBasics.html --> public void CalcEigenvals()
+
+    std::cout<<"// Check Roots"<<std::endl;
+    CoefficientsToRoots::printCheck(M);
+    
+    // TODO : Consider making this a map of root key and count index. 
+    // requires defining a resolution thresshold that unites keys of close values
     std::vector<c128> roots;
- 
     int i = 0;
     while (i < degree)
     {
-        double a, b, c;
-        a = 1.0;
-        b = -M[i][i] - M[i+1][i+1];
-        c = M[i][i] * M[i+1][i+1] - M[i][i+1] * M[i+1][i];
-
-        if ((b * b) >= (4.0 * a * c))
+        
+        if ( i == degree - 1 || std::abs(M[i+1][i] < Epsilon) )
         {
-            roots.emplace_back(
-                (-b + std::sqrt((b * b) - (4.0 * a * c))) / 2.0 / a, 
-                0.0);
-            roots.emplace_back(
-                (-b - std::sqrt((b * b) - (4.0 * a * c))) / 2.0 / a,
-                0.0);
+            // Real eigenvalue on diagonal
+            roots.emplace_back(M[i][i], 0.0);
             ++i;
         }
         else
         {
-            double realPart = -b / 2.0 / a;
-            double imagPart =
-                std::sqrt((4.0 * a * c) - (b * b))
-                / 2.0 / a;
-
-            roots.emplace_back(realPart,  imagPart);
-            roots.emplace_back(realPart, -imagPart);
-
+            // check if conjugate by solving the equation which instantiates by the 2x2 cell:
+            // https://www.physicsforums.com/threads/how-do-i-estimate-complex-eigenvalues.170108/post-1330547
+            // https://www.mosismath.com/Eigenvalues/EigenvalsQR.html
+            // det(A - λI) = 0
+            // =>  det| a-λ    b |  
+            //        |   c  d-λ | = 0 
+            // => (a-λ)(d-λ) - bc = 0
+            // => λ^2 - (a+d)λ + (ad-bc) = 0
+            // => λ^2 - (a+d)λ + detA = 0 --> solve with discriminant
+            const double a = M[i][i];
+            const double b = M[i][i+1];
+            const double c = M[i+1][i];
+            const double d = M[i+1][i+1];
+            const double det  = a*d - b*c;
+            // solve with discriminant
+            const double discriminant = (a+d)*(a+d) - 4.0*det;
+
+            if ( discriminant >= 0.0)
+            {
+                // tow real roots
+                roots.emplace_back(((a+d) + std::sqrt(discriminant)) / 2.0, 0.0);
+                roots.emplace_back(((a+d) - std::sqrt(discriminant)) / 2.0, 0.0);
+            }
+            else
+            {
+                // Complex conjugate pair
+                const double re = (a+d) / 2.0;
+                const double im = std::sqrt(-discriminant) / 2.0;
+                roots.emplace_back(re,  im);
+                roots.emplace_back(re, -im);
+            }
             i += 2;
         }
     }
@@ -201,5 +220,4 @@ std::vector<c128> CoefficientsToRoots::extractRoots(const std::vector<std::vecto
     std::cout<<std::endl;
 
     return roots;
-
-}
\ No newline at end of file
+}

From 6f7856160b0e60243bd9cbd657c4a676c9f26bcd Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Tue, 26 May 2026 09:37:28 +0300
Subject: [PATCH 03/16] fix typo in CoefficientsToRoots::extractRoots method

---
 src/CoefficientsToRoots.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index cc7ae62..0191f52 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -171,7 +171,7 @@ std::vector<c128> CoefficientsToRoots::extractRoots(const std::vector<std::vecto
     while (i < degree)
     {
         
-        if ( i == degree - 1 || std::abs(M[i+1][i] < Epsilon) )
+        if ( i == degree - 1 || std::abs(M[i+1][i]) < Epsilon )
         {
             // Real eigenvalue on diagonal
             roots.emplace_back(M[i][i], 0.0);

From e1fcfe86bfdf77634bf653b89318af385caf033b Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Tue, 26 May 2026 19:04:14 +0300
Subject: [PATCH 04/16] implement updating of the filter state on coefficient
 edits

---
 src/CoeffComponents.cpp     | 151 ++++++++++++++++++++++++------------
 src/CoeffComponents.h       |   1 +
 src/CoefficientsToRoots.cpp |  46 +++++++----
 src/CoefficientsToRoots.h   |   8 +-
 4 files changed, 141 insertions(+), 65 deletions(-)

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index 5323169..f95b03a 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -116,54 +116,7 @@ juce::Component* CoefficientsComponent::refreshComponentForCell(int row, int col
         // update data when editing finishes
         label->onTextChange = [this, row, col, label]{
             double value = label->getText().getDoubleValue();
-            if (col == 2)
-                ffcoeffs[static_cast<size_t>(row)] = value;
-            else if (col == 3)
-                fbcoeffs[static_cast<size_t>(row)] = value;
-            
-            // TODO: calculate roots through the coefficient2roots function.
-            // TODO: notify other listeners about this change
-            if (col == 2)
-            {
-
-                auto zeros = CoefficientsToRoots::GramSchmidt(this->ffcoeffs);
-
-                std::cout<<"Existing Zeros"<<std::endl;
-                for (int i=0 ; i< processor->filterState->zeros.size(); i++)
-                {
-                    auto *r = processor->filterState->zeros[i];
-                    std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
-                    // processor->filterState->remove(r);
-                }
-                std::cout<<std::endl;
-
-                // for (int i=0 ; i< zeros.size(); i++)
-                // {
-                //     auto &r = zeros[i];
-                //     processor->filterState->add(i, r);
-                // }
-
-            }
-            else if (col == 3)
-            {
-                auto poles =  CoefficientsToRoots::GramSchmidt(this->fbcoeffs);
-
-                std::cout<<"Existing Poles"<<std::endl;
-                for (int i=0 ; i< processor->filterState->poles.size(); i++)
-                {
-                    auto *r = processor->filterState->poles[i];
-                    std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
-                    // processor->filterState->remove(r);
-                }
-                std::cout<<std::endl;
-
-                // for (int i=0 ; i< poles.size(); i++)
-                // {
-                //     auto &r = poles[i];
-                //     processor->filterState->add(i, r);
-                // }
-            }
-            
+            updateFilterStateOnCoefEdit(row, col, value);
         };
     }
 
@@ -181,6 +134,108 @@ juce::Component* CoefficientsComponent::refreshComponentForCell(int row, int col
     return label;
 }
 
+
+void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double value)
+{
+    // calculate roots through the coefficient2roots function && notify other listeners about this change
+
+    // @TODO Clean print statements (after ensuring robust implementation)...
+
+    if (col == 2)
+        ffcoeffs[static_cast<size_t>(row)] = value;
+    else if (col == 3)
+        fbcoeffs[static_cast<size_t>(row)] = value;
+    
+    if (col == 2)
+    {
+
+        auto zeros = CoefficientsToRoots::GramSchmidt(this->ffcoeffs);
+
+#ifdef DEBUG_C2R
+        std::cout<<"Previous Zeros"<<std::endl;
+        for (int i=0 ; i< processor->filterState->zeros.size(); i++)
+        {
+            auto *r = processor->filterState->zeros[i];
+            std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+        }std::cout<<std::endl;
+#endif
+
+        std::cout<<"Removing Zeroes"<<std::endl;
+        while (!processor->filterState->zeros.isEmpty())
+        {
+            auto *r = processor->filterState->zeros.getFirst();
+            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+            processor->filterState->remove(r);
+        }
+        std::cout<<std::endl;
+
+        std::cout<<"Adding Zeros"<<std::endl;
+        for (int i=0 ; i< zeros.size(); i++)
+        {
+            auto &r = zeros[i].first;
+            auto &order = zeros[i].second;
+            std::cout<<"\t("<<r.real()<<","<< r.imag()<<") ";
+            processor->filterState->add(order, r);
+        }
+        std::cout<<std::endl;
+
+#ifdef DEBUG_C2R
+        std::cout<<"Final Zeros"<<std::endl;
+        for (int i=0 ; i< processor->filterState->zeros.size(); i++)
+        {
+            auto *r = processor->filterState->zeros[i];
+
+            std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+        }std::cout<<std::endl;
+#endif
+    }
+    else if (col == 3)
+    {
+        auto poles =  CoefficientsToRoots::GramSchmidt(this->fbcoeffs);
+
+#ifdef DEBUG_C2R
+        std::cout<<"Previous Poles"<<std::endl;
+        for (int i=0 ; i< processor->filterState->poles.size(); i++)
+        {
+            auto *r = processor->filterState->poles[i];
+            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+        }std::cout<<std::endl;
+#endif
+        std::cout<<"Adding Poles"<<std::endl;
+        for (int i=0 ; i< poles.size(); i++)
+        {
+            
+            auto &r = poles[i].first;
+            auto &order = poles[i].second;
+            processor->filterState->add(-order, r);
+            std::cout<<"\t("<<r.real()<<","<< r.imag()<<") ";
+        }
+        std::cout<<std::endl;
+
+        std::cout<<"Removing Poles"<<std::endl;
+        size_t sz = processor->filterState->poles.size() - poles.size();
+        for (int i=0 ; i< sz; i++)
+        {
+            auto *r = processor->filterState->poles.getFirst();
+            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+            processor->filterState->remove(r);
+        }
+        std::cout<<std::endl;
+#ifdef DEBUG_C2R
+        std::cout<<"Final Poles"<<std::endl;
+        for (int i=0 ; i< processor->filterState->poles.size(); i++)
+        {
+            auto *r = processor->filterState->poles[i];
+            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+        }std::cout<<std::endl;
+#endif
+    }
+#ifdef DEBUG_C2R
+    std::cout<<"Root order "<<processor->filterState->totalOrder <<" finiteZerosOrder "<< processor->filterState->finiteZerosOrder<<std::endl;
+#endif
+    
+}
+
 void CoefficientsComponent::valueTreePropertyChanged (juce::ValueTree& node, const juce::Identifier& property)
 {
     if(property == IDs::ValueRe || property == IDs::ValueIm)    // when dragging roots
diff --git a/src/CoeffComponents.h b/src/CoeffComponents.h
index b1aa6c3..8e1ae06 100644
--- a/src/CoeffComponents.h
+++ b/src/CoeffComponents.h
@@ -27,6 +27,7 @@ class CoefficientsComponent final
 
         void toggleCollapseExpand();
         void updateCoeffTable();
+        void updateFilterStateOnCoefEdit(int row, int col, double value);
 
         // override juce::Component
         // void paint(juce::Graphics &g) override; // TODO is this trully needed?
diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index 0191f52..42b9856 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -8,7 +8,7 @@
 #include <iostream>
 #include <cmath>
 
-std::vector<c128> CoefficientsToRoots::GramSchmidt(std::vector<double> coefs)
+std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<double> coefs)
 {
     // filter out leading zeros and leading 1.0
     int degree = 0;
@@ -20,7 +20,7 @@ std::vector<c128> CoefficientsToRoots::GramSchmidt(std::vector<double> coefs)
     if (degree==1)
     {
         std::cout<<"Returing 1 root --> "<<coefs[degree]<<std::endl;
-        return {static_cast<c128>(coefs[degree])};
+        return {std::make_pair(static_cast<c128>(coefs[degree]), 1)};
     }
 
     // build companion Matrix
@@ -157,16 +157,30 @@ void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &mat
     }
 }
 
-
-std::vector<c128> CoefficientsToRoots::extractRoots(const std::vector<std::vector<double>>& M, int degree)
+std::vector<std::pair<c128, int>> CoefficientsToRoots::extractRoots(const std::vector<std::vector<double>>& M, int degree)
 {
 
     std::cout<<"// Check Roots"<<std::endl;
     CoefficientsToRoots::printCheck(M);
     
-    // TODO : Consider making this a map of root key and count index. 
-    // requires defining a resolution thresshold that unites keys of close values
-    std::vector<c128> roots;
+    std::vector<std::pair<c128, int>> roots;
+
+    auto addRoot = [&](c128 newVal) {
+        for (auto& [val, order] : roots)
+        {
+            double diff_re = std::abs(val.real() - newVal.real());
+            double diff_im = std::abs(val.imag() - newVal.imag());
+            double scale = std::abs(val) + std::abs(newVal) + 1e-7; // + 1e-7 to avoid division with zero
+
+            if (diff_re / scale < tolerance && diff_im / scale < tolerance)
+            {
+                order++;
+                return;
+            }
+        }
+        roots.emplace_back(newVal, 1);
+    };
+
     int i = 0;
     while (i < degree)
     {
@@ -174,7 +188,8 @@ std::vector<c128> CoefficientsToRoots::extractRoots(const std::vector<std::vecto
         if ( i == degree - 1 || std::abs(M[i+1][i]) < Epsilon )
         {
             // Real eigenvalue on diagonal
-            roots.emplace_back(M[i][i], 0.0);
+            c128 newRoot (M[i][i], 0.0);
+            addRoot(newRoot);
             ++i;
         }
         else
@@ -199,25 +214,26 @@ std::vector<c128> CoefficientsToRoots::extractRoots(const std::vector<std::vecto
             if ( discriminant >= 0.0)
             {
                 // tow real roots
-                roots.emplace_back(((a+d) + std::sqrt(discriminant)) / 2.0, 0.0);
-                roots.emplace_back(((a+d) - std::sqrt(discriminant)) / 2.0, 0.0);
+                addRoot(c128(((a+d) + std::sqrt(discriminant)) / 2.0, 0.0));
+                addRoot(c128(((a+d) - std::sqrt(discriminant)) / 2.0, 0.0));
             }
             else
             {
                 // Complex conjugate pair
                 const double re = (a+d) / 2.0;
                 const double im = std::sqrt(-discriminant) / 2.0;
-                roots.emplace_back(re,  im);
-                roots.emplace_back(re, -im);
+                addRoot(c128(re,im));
+                // addRoot(c128(re,-im)); // Note: this is added automatically using addRoot. Commenting this, removes the bug of overlapping roots.
             }
             i += 2;
         }
     }
 
+#ifdef DEBUG_C2R
     std::cout<<"Calculated Roots: ";
-    for (auto root : roots)
-        std::cout<<"("<<root.real()<<","<<root.imag()<<")  ";
+    for (auto [root,order] : roots)
+        std::cout<<"("<<root.real()<<","<<root.imag()<<") - "<<order<<", ";
     std::cout<<std::endl;
-
+#endif
     return roots;
 }
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index 3b6519f..2a45225 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -2,17 +2,21 @@
 
 #include "FilterState.h"
 #include <vector>
+#include <utility>
+
+#define DEBUG_C2R
 
 class CoefficientsToRoots
 {
 	public:
-		static std::vector<c128> GramSchmidt(std::vector<double> coefs);
+		static std::vector<std::pair<c128, int>> GramSchmidt(std::vector<double> coefs); // TODO change return type to bool --> for when there is a pole outside the unit. Check will be done outside of this.
 	private:
 
 		// TODO Finetune these parameters
 	    static constexpr double Epsilon = 1e-3;
 	    static constexpr int MaxIterations = 100000;
+		static constexpr double tolerance = 1e-4; // threshold for considering two roots with negligible diff the same. Expressed in % after scaling differences, since zeros can much more than 1.0.
 
 		static void printCheck(const std::vector<std::vector<double>> &matrix);
-		static std::vector<c128> extractRoots(const std::vector<std::vector<double>>& M, int n);
+		static std::vector<std::pair<c128, int>> extractRoots(const std::vector<std::vector<double>>& M, int n);
 };
\ No newline at end of file

From fba0a7146b7e0693b06bf0bb1d1aa26ec157d68b Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Wed, 27 May 2026 01:16:46 +0300
Subject: [PATCH 05/16] fix 4 bugs all at once.

1. fix bug in updateFilterStateOnCoefEdit when adding a pole that already exists (increasing its order), and then removing it (removing the previously added new pole)

2. fix index of returned root when degree equals 1.0

3. fix sign of returned root when degree equals 1.0

4. fix coef_index when constructing the companion matrix. Now it always initializes the matrix with the last degree-lengthed coefficients (correctly ignoring the leading zeros and 1.0)
---
 src/CoeffComponents.cpp     | 22 +++++++++++++++++++++-
 src/CoefficientsToRoots.cpp |  4 ++--
 2 files changed, 23 insertions(+), 3 deletions(-)

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index f95b03a..c1b22eb 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -213,12 +213,32 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         std::cout<<std::endl;
 
         std::cout<<"Removing Poles"<<std::endl;
+
+        auto isNewPole = [&](c128 val) -> std::pair<bool, int> {
+            for (auto& [r, order] : poles) 
+            {
+                if (val.real() == r.real() && val.imag() == r.imag())
+                    return {true, order};
+            }
+            return {false, 0};
+        };
+        
         size_t sz = processor->filterState->poles.size() - poles.size();
         for (int i=0 ; i< sz; i++)
         {
             auto *r = processor->filterState->poles.getFirst();
             std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
-            processor->filterState->remove(r);
+            c128 key(c128(r->value.re.get(),r->value.im.get()));
+            auto [found, order] = isNewPole(key);
+
+            if ( found )
+            {
+                r->order= order;
+            }
+            else
+            {
+                processor->filterState->remove(r);
+            } 
         }
         std::cout<<std::endl;
 #ifdef DEBUG_C2R
diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index 42b9856..f608670 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -20,7 +20,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
     if (degree==1)
     {
         std::cout<<"Returing 1 root --> "<<coefs[degree]<<std::endl;
-        return {std::make_pair(static_cast<c128>(coefs[degree]), 1)};
+        return {std::make_pair(static_cast<c128>(-coefs[coefs.size() - 1]), 1)};
     }
 
     // build companion Matrix
@@ -30,7 +30,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
     for (int i=0; i< degree; i++)
     {
         int j = i+1;
-        int coef_idx = degree-i;
+        int coef_idx = coefs.size()-1-i;
         A[i][degree-1] = -coefs[coef_idx];  // fill last column with negative vals of coefs
         if (i!=degree-1)
         {

From f0dc2f2571cbc7361852bfe742b008d7e09a4234 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Thu, 28 May 2026 16:17:29 +0300
Subject: [PATCH 06/16] minor improvements and fixes:

change termination condition, now all subdiagonal values should be below epsilon (not their sum)

fix typo : step 3 - A = R Q

change tolerance threshold to 5e-2

add order in prints for debugging
---
 src/CoeffComponents.cpp     | 16 ++++++++--------
 src/CoefficientsToRoots.cpp | 26 +++++++++++++++-----------
 src/CoefficientsToRoots.h   |  2 +-
 3 files changed, 24 insertions(+), 20 deletions(-)

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index c1b22eb..fa58b50 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -156,7 +156,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         for (int i=0 ; i< processor->filterState->zeros.size(); i++)
         {
             auto *r = processor->filterState->zeros[i];
-            std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+            std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
         }std::cout<<std::endl;
 #endif
 
@@ -164,7 +164,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         while (!processor->filterState->zeros.isEmpty())
         {
             auto *r = processor->filterState->zeros.getFirst();
-            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
             processor->filterState->remove(r);
         }
         std::cout<<std::endl;
@@ -174,7 +174,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         {
             auto &r = zeros[i].first;
             auto &order = zeros[i].second;
-            std::cout<<"\t("<<r.real()<<","<< r.imag()<<") ";
+            std::cout<<"\t("<<r.real()<<","<< r.imag()<<") - "<<order<<", ";
             processor->filterState->add(order, r);
         }
         std::cout<<std::endl;
@@ -185,7 +185,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         {
             auto *r = processor->filterState->zeros[i];
 
-            std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+            std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
         }std::cout<<std::endl;
 #endif
     }
@@ -198,7 +198,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         for (int i=0 ; i< processor->filterState->poles.size(); i++)
         {
             auto *r = processor->filterState->poles[i];
-            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
         }std::cout<<std::endl;
 #endif
         std::cout<<"Adding Poles"<<std::endl;
@@ -208,7 +208,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
             auto &r = poles[i].first;
             auto &order = poles[i].second;
             processor->filterState->add(-order, r);
-            std::cout<<"\t("<<r.real()<<","<< r.imag()<<") ";
+            std::cout<<"\t("<<r.real()<<","<< r.imag()<<") - "<<order<<", ";
         }
         std::cout<<std::endl;
 
@@ -227,7 +227,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         for (int i=0 ; i< sz; i++)
         {
             auto *r = processor->filterState->poles.getFirst();
-            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
             c128 key(c128(r->value.re.get(),r->value.im.get()));
             auto [found, order] = isNewPole(key);
 
@@ -246,7 +246,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         for (int i=0 ; i< processor->filterState->poles.size(); i++)
         {
             auto *r = processor->filterState->poles[i];
-            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") ";
+            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
         }std::cout<<std::endl;
 #endif
     }
diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index f608670..3ce0c0c 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -115,7 +115,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
             }
         }
 
-        // step 3 - A = Q R
+        // step 3 - A = R Q
         for (int row = 0; row < degree; row++)
         {
             for (int col = 0; col < degree; col++)
@@ -123,27 +123,31 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
                 double sum {0.0};
                 for (int inner = 0; inner< degree; ++inner)
                 {
-                    sum += Q[inner][col] * R[row][inner];
+                    sum += R[row][inner] * Q[inner][col];
                 }
                 A[row][col] = sum; // resetting A[row][col]
             }
         }
         
-        // step 4 sum all sub-diagonal entries
-        double subdiag_sum {0.0};
+        // step 4 check all sub-diagonal entries
+        bool subdiag_low {true};
         for (int i = 1; i < degree; ++i)
         {
-            subdiag_sum += std::abs(A[i][i-1]);
+            if (std::abs(A[i][i-1]) > Epsilon)
+            {
+                subdiag_low = false;
+                break;
+            }
         }
 
-        if (subdiag_sum < Epsilon)
+        if (subdiag_low)
             break;
     }
-    
+
     // Extract roots — read diagonal
     return extractRoots(A, degree);
 }
- 
+
 void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &matrix)
 {
     int n = matrix.size();
@@ -231,9 +235,9 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::extractRoots(const std::v
 
 #ifdef DEBUG_C2R
     std::cout<<"Calculated Roots: ";
-    for (auto [root,order] : roots)
-        std::cout<<"("<<root.real()<<","<<root.imag()<<") - "<<order<<", ";
+    for (auto root : roots)
+        std::cout<<"("<<root.first.real()<<","<<root.first.imag()<<") - "<<root.second<<", ";
     std::cout<<std::endl;
 #endif
     return roots;
-}
+}
\ No newline at end of file
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index 2a45225..f0a244d 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -15,7 +15,7 @@ class CoefficientsToRoots
 		// TODO Finetune these parameters
 	    static constexpr double Epsilon = 1e-3;
 	    static constexpr int MaxIterations = 100000;
-		static constexpr double tolerance = 1e-4; // threshold for considering two roots with negligible diff the same. Expressed in % after scaling differences, since zeros can much more than 1.0.
+		static constexpr double tolerance = 5e-2; // threshold for considering two roots with negligible diff the same. Expressed in % after scaling differences, since zeros can much more than 1.0.
 
 		static void printCheck(const std::vector<std::vector<double>> &matrix);
 		static std::vector<std::pair<c128, int>> extractRoots(const std::vector<std::vector<double>>& M, int n);

From 2e9d98b80dd8dddb52517fdaabbd9fbf84c6a626 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Fri, 29 May 2026 18:26:36 +0300
Subject: [PATCH 07/16] supress gcc warnings

---
 src/CoefficientsToRoots.cpp | 54 ++++++++++++++++++-------------------
 src/CoefficientsToRoots.h   |  6 ++---
 2 files changed, 30 insertions(+), 30 deletions(-)

diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index 3ce0c0c..c780a4d 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -11,8 +11,8 @@
 std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<double> coefs)
 {
     // filter out leading zeros and leading 1.0
-    int degree = 0;
-    while(coefs[degree] ==0 || coefs[degree] == 1.0)
+    size_t degree = 0;
+    while(coefs[degree] ==0.0 || coefs[degree] == 1.0)
         degree++;
     degree = coefs.size() - degree;
     jassert( degree!=0 );
@@ -27,10 +27,10 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
     std::cout<<"Building companion matrix with degree "<<degree<<std::endl;
     std::vector<std::vector<double>> A(degree, std::vector<double>(degree, 0.0));
     // Build companion matrix (Hessenberg form)
-    for (int i=0; i< degree; i++)
+    for (size_t i=0; i< degree; i++)
     {
-        int j = i+1;
-        int coef_idx = coefs.size()-1-i;
+        size_t j = i+1;
+        size_t coef_idx = coefs.size()-1-i;
         A[i][degree-1] = -coefs[coef_idx];  // fill last column with negative vals of coefs
         if (i!=degree-1)
         {
@@ -47,7 +47,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
     std::vector<std::vector<double>> Q(degree, std::vector<double>(degree));
     std::vector<std::vector<double>> R(degree, std::vector<double>(degree));
     
-    int iter = 0;
+    size_t iter = 0;
     while(true) // could alternatively check the convergence of the values of the sub-diagonal elements
     {
 
@@ -55,9 +55,9 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
 
 
         // Reset R,Q
-        for (int r = 0; r < degree; ++r)
+        for (size_t r = 0; r < degree; ++r)
         {
-            for (int c = 0; c < degree; ++c)
+            for (size_t c = 0; c < degree; ++c)
             {
                 R[r][c] = 0.0;
                 Q[r][c] = 0.0;
@@ -65,25 +65,25 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
         }
 
         // step 1 - calculate Q : Q = A[col] - projections onto the previous Q[col]
-        for (int col = 0 ; col < degree; col++)
+        for (size_t col = 0 ; col < degree; col++)
         {
             // set v with column A[col]
             std::vector<double> v(degree);
-            for (int row = 0 ; row< degree; row++)
+            for (size_t row = 0 ; row< degree; row++)
                 v[row] = A[row][col];
 
             // subtract projection : v[col] = A[col] - proj onto the previous (<col) orthonormal colums of Q.
-            for (int curr_col = 0; curr_col < col; curr_col++)
+            for (size_t curr_col = 0; curr_col < col; curr_col++)
             {
                 // compute dot product ...
                 double dot_product {0.0};
-                for (int curr_row = 0; curr_row < degree; ++curr_row)
+                for (size_t curr_row = 0; curr_row < degree; ++curr_row)
                 {
                     dot_product += Q[curr_row][curr_col] * A[curr_row][col];
                 }
 
                 // .. and subtract it from the current column vector
-                for (int curr_row=0; curr_row < degree; ++curr_row)
+                for (size_t curr_row=0; curr_row < degree; ++curr_row)
                 {
                     v[curr_row] -= dot_product * Q[curr_row][curr_col];
                 }
@@ -91,24 +91,24 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
 
             // normalize column of q
             double norm = 0.0;
-            for (int k = 0; k < degree; ++k)
+            for (size_t k = 0; k < degree; ++k)
             {
                 norm += v[k] * v[k];
             }
             norm = std::sqrt(norm);
             
-            for (int k = 0; k < degree; ++k)
+            for (size_t k = 0; k < degree; ++k)
             {
                 Q[k][col] = v[k] / norm;
             }
         }
 
         // step 2 - calculate R :  R = Q A
-        for (int row = 0; row < degree; row++)
+        for (size_t row = 0; row < degree; row++)
         {
-            for (int col = 0; col < degree; col++)
+            for (size_t col = 0; col < degree; col++)
             {
-                for (int inner = 0; inner < degree; inner++)
+                for (size_t inner = 0; inner < degree; inner++)
                 {
                     R[row][col] += (Q[inner][row] * A[inner][col]);     // Q transposed
                 }
@@ -116,12 +116,12 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
         }
 
         // step 3 - A = R Q
-        for (int row = 0; row < degree; row++)
+        for (size_t row = 0; row < degree; row++)
         {
-            for (int col = 0; col < degree; col++)
+            for (size_t col = 0; col < degree; col++)
             {
                 double sum {0.0};
-                for (int inner = 0; inner< degree; ++inner)
+                for (size_t inner = 0; inner< degree; ++inner)
                 {
                     sum += R[row][inner] * Q[inner][col];
                 }
@@ -131,7 +131,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
         
         // step 4 check all sub-diagonal entries
         bool subdiag_low {true};
-        for (int i = 1; i < degree; ++i)
+        for (size_t i = 1; i < degree; ++i)
         {
             if (std::abs(A[i][i-1]) > Epsilon)
             {
@@ -150,10 +150,10 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
 
 void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &matrix)
 {
-    int n = matrix.size();
-    for (int i=0; i< n; i++)
+    size_t n = matrix.size();
+    for (size_t i=0; i< n; i++)
     {
-        for (int j=0; j<n; j++)
+        for (size_t j=0; j<n; j++)
         {
             std::cout<<matrix[i][j]<<" ";
         }
@@ -161,7 +161,7 @@ void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &mat
     }
 }
 
-std::vector<std::pair<c128, int>> CoefficientsToRoots::extractRoots(const std::vector<std::vector<double>>& M, int degree)
+std::vector<std::pair<c128, int>> CoefficientsToRoots::extractRoots(const std::vector<std::vector<double>>& M, size_t degree)
 {
 
     std::cout<<"// Check Roots"<<std::endl;
@@ -185,7 +185,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::extractRoots(const std::v
         roots.emplace_back(newVal, 1);
     };
 
-    int i = 0;
+    size_t i = 0;
     while (i < degree)
     {
         
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index f0a244d..d4441cf 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -14,9 +14,9 @@ class CoefficientsToRoots
 
 		// TODO Finetune these parameters
 	    static constexpr double Epsilon = 1e-3;
-	    static constexpr int MaxIterations = 100000;
+	    static constexpr size_t MaxIterations = 100000;
 		static constexpr double tolerance = 5e-2; // threshold for considering two roots with negligible diff the same. Expressed in % after scaling differences, since zeros can much more than 1.0.
 
-		static void printCheck(const std::vector<std::vector<double>> &matrix);
-		static std::vector<std::pair<c128, int>> extractRoots(const std::vector<std::vector<double>>& M, int n);
+		static void printCheck(const std::vector<std::vector<double>> &);
+		static std::vector<std::pair<c128, int>> extractRoots(const std::vector<std::vector<double>>&, size_t);
 };
\ No newline at end of file

From 6cd4ac6d51999d522addcdc9c6e8d6ae4c505f81 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Wed, 3 Jun 2026 17:02:32 +0300
Subject: [PATCH 08/16] implement rayleigh quotient shift to speed up QR
 convergence

---
 src/CoeffComponents.cpp     | 14 ++++----
 src/CoefficientsToRoots.cpp | 64 +++++++++++++++++++++----------------
 src/CoefficientsToRoots.h   |  4 +--
 3 files changed, 45 insertions(+), 37 deletions(-)

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index fa58b50..7367b56 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -170,7 +170,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         std::cout<<std::endl;
 
         std::cout<<"Adding Zeros"<<std::endl;
-        for (int i=0 ; i< zeros.size(); i++)
+        for (size_t i=0 ; i< zeros.size(); i++)
         {
             auto &r = zeros[i].first;
             auto &order = zeros[i].second;
@@ -181,7 +181,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
 
 #ifdef DEBUG_C2R
         std::cout<<"Final Zeros"<<std::endl;
-        for (int i=0 ; i< processor->filterState->zeros.size(); i++)
+        for (size_t i=0 ; i< processor->filterState->zeros.size(); i++)
         {
             auto *r = processor->filterState->zeros[i];
 
@@ -195,14 +195,14 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
 
 #ifdef DEBUG_C2R
         std::cout<<"Previous Poles"<<std::endl;
-        for (int i=0 ; i< processor->filterState->poles.size(); i++)
+        for (size_t i=0 ; i< processor->filterState->poles.size(); i++)
         {
             auto *r = processor->filterState->poles[i];
             std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
         }std::cout<<std::endl;
 #endif
         std::cout<<"Adding Poles"<<std::endl;
-        for (int i=0 ; i< poles.size(); i++)
+        for (size_t i=0 ; i< poles.size(); i++)
         {
             
             auto &r = poles[i].first;
@@ -223,8 +223,8 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
             return {false, 0};
         };
         
-        size_t sz = processor->filterState->poles.size() - poles.size();
-        for (int i=0 ; i< sz; i++)
+        size_t sz = static_cast<size_t>(processor->filterState->poles.size()) - poles.size();
+        for (size_t i=0 ; i< sz; i++)
         {
             auto *r = processor->filterState->poles.getFirst();
             std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
@@ -243,7 +243,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         std::cout<<std::endl;
 #ifdef DEBUG_C2R
         std::cout<<"Final Poles"<<std::endl;
-        for (int i=0 ; i< processor->filterState->poles.size(); i++)
+        for (size_t i=0 ; i< processor->filterState->poles.size(); i++)
         {
             auto *r = processor->filterState->poles[i];
             std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index c780a4d..d366a18 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -47,29 +47,34 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
     std::vector<std::vector<double>> Q(degree, std::vector<double>(degree));
     std::vector<std::vector<double>> R(degree, std::vector<double>(degree));
     
-    size_t iter = 0;
-    while(true) // could alternatively check the convergence of the values of the sub-diagonal elements
+    size_t shift_idx {degree}, iter {0};
+    while(shift_idx > 1)
     {
 
         if (++iter > degree * MaxIterations) break;
 
 
         // Reset R,Q
-        for (size_t r = 0; r < degree; ++r)
+        for (size_t r = 0; r < shift_idx; ++r)
         {
-            for (size_t c = 0; c < degree; ++c)
+            for (size_t c = 0; c < shift_idx; ++c)
             {
                 R[r][c] = 0.0;
                 Q[r][c] = 0.0;
             }
         }
 
+        // subtract rayleigh quotient shift
+        const double shift = A[shift_idx - 1][shift_idx - 1];
+        for (size_t i = 0; i < shift_idx; ++i) 
+            A[i][i] -= shift; // Decompose (Ak - sI)
+
         // step 1 - calculate Q : Q = A[col] - projections onto the previous Q[col]
-        for (size_t col = 0 ; col < degree; col++)
+        for (size_t col = 0 ; col < shift_idx; col++)
         {
             // set v with column A[col]
-            std::vector<double> v(degree);
-            for (size_t row = 0 ; row< degree; row++)
+            std::vector<double> v(shift_idx);
+            for (size_t row = 0 ; row< shift_idx; row++)
                 v[row] = A[row][col];
 
             // subtract projection : v[col] = A[col] - proj onto the previous (<col) orthonormal colums of Q.
@@ -77,13 +82,13 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
             {
                 // compute dot product ...
                 double dot_product {0.0};
-                for (size_t curr_row = 0; curr_row < degree; ++curr_row)
+                for (size_t curr_row = 0; curr_row < shift_idx; ++curr_row)
                 {
                     dot_product += Q[curr_row][curr_col] * A[curr_row][col];
                 }
 
                 // .. and subtract it from the current column vector
-                for (size_t curr_row=0; curr_row < degree; ++curr_row)
+                for (size_t curr_row=0; curr_row < shift_idx; ++curr_row)
                 {
                     v[curr_row] -= dot_product * Q[curr_row][curr_col];
                 }
@@ -91,24 +96,24 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
 
             // normalize column of q
             double norm = 0.0;
-            for (size_t k = 0; k < degree; ++k)
+            for (size_t k = 0; k < shift_idx; ++k)
             {
                 norm += v[k] * v[k];
             }
             norm = std::sqrt(norm);
             
-            for (size_t k = 0; k < degree; ++k)
+            for (size_t k = 0; k < shift_idx; ++k)
             {
                 Q[k][col] = v[k] / norm;
             }
         }
 
         // step 2 - calculate R :  R = Q A
-        for (size_t row = 0; row < degree; row++)
+        for (size_t row = 0; row < shift_idx; row++)
         {
-            for (size_t col = 0; col < degree; col++)
+            for (size_t col = 0; col < shift_idx; col++)
             {
-                for (size_t inner = 0; inner < degree; inner++)
+                for (size_t inner = 0; inner < shift_idx; inner++)
                 {
                     R[row][col] += (Q[inner][row] * A[inner][col]);     // Q transposed
                 }
@@ -116,32 +121,35 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
         }
 
         // step 3 - A = R Q
-        for (size_t row = 0; row < degree; row++)
+        for (size_t row = 0; row < shift_idx; row++)
         {
-            for (size_t col = 0; col < degree; col++)
+            for (size_t col = 0; col < shift_idx; col++)
             {
                 double sum {0.0};
-                for (size_t inner = 0; inner< degree; ++inner)
+                for (size_t inner = 0; inner< shift_idx; ++inner)
                 {
                     sum += R[row][inner] * Q[inner][col];
                 }
                 A[row][col] = sum; // resetting A[row][col]
             }
         }
-        
-        // step 4 check all sub-diagonal entries
-        bool subdiag_low {true};
-        for (size_t i = 1; i < degree; ++i)
+
+
+        // step 4 add shift back in A and check sub-diagonal entries
+        for (size_t i = 0; i < shift_idx; ++i) 
+            A[i][i] += shift;
+
+        // check only the last subdiagonal entry (real eigenvalue)
+        if (std::abs(A[shift_idx-1][shift_idx-2]) < Epsilon)
         {
-            if (std::abs(A[i][i-1]) > Epsilon)
-            {
-                subdiag_low = false;
-                break;
-            }
+            shift_idx--;
+        }
+        // check the entry above the 2x2 block in search of complex conjugate pairs
+        else if (shift_idx > 2 && std::abs(A[shift_idx-2][shift_idx-3]) < Epsilon)
+        {
+            shift_idx -= 2;
         }
 
-        if (subdiag_low)
-            break;
     }
 
     // Extract roots — read diagonal
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index d4441cf..2089276 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -4,7 +4,7 @@
 #include <vector>
 #include <utility>
 
-#define DEBUG_C2R
+// #define DEBUG_C2R
 
 class CoefficientsToRoots
 {
@@ -14,7 +14,7 @@ class CoefficientsToRoots
 
 		// TODO Finetune these parameters
 	    static constexpr double Epsilon = 1e-3;
-	    static constexpr size_t MaxIterations = 100000;
+	    static constexpr size_t MaxIterations = 100;
 		static constexpr double tolerance = 5e-2; // threshold for considering two roots with negligible diff the same. Expressed in % after scaling differences, since zeros can much more than 1.0.
 
 		static void printCheck(const std::vector<std::vector<double>> &);

From bfb04e3a1ba66e1bb59499a9e28cf5abf6554097 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Tue, 9 Jun 2026 12:25:55 +0300
Subject: [PATCH 09/16] fix 3 bugs:

trailing zeros in (feedback) coefficients (when adding default poles) are now parsed out of the companion matrix and roots at 0.,0. (of order equal to #numOfTrailingZeros) are appended directly into the root list before building the companion matrix and applying the QR

prev_sz is now calculated directly before appending new poles on the filterstate, taking into account that some updated roots may not result in removing previous poles, but changing their order instead

leading zeros if (feedforward) coefficients (when total order of poles exceeds total order of zeros) are now tracked using std::numeric_limits instead of comparing doubles with the == operator.
---
 src/CoeffComponents.cpp     | 12 +++++++-----
 src/CoefficientsToRoots.cpp | 33 +++++++++++++++++++++++++--------
 src/CoefficientsToRoots.h   |  2 +-
 3 files changed, 33 insertions(+), 14 deletions(-)

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index 7367b56..fbd29a1 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -193,6 +193,8 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
     {
         auto poles =  CoefficientsToRoots::GramSchmidt(this->fbcoeffs);
 
+        size_t prev_sz = static_cast<size_t>(processor->filterState->poles.size());
+
 #ifdef DEBUG_C2R
         std::cout<<"Previous Poles"<<std::endl;
         for (size_t i=0 ; i< processor->filterState->poles.size(); i++)
@@ -213,7 +215,6 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         std::cout<<std::endl;
 
         std::cout<<"Removing Poles"<<std::endl;
-
         auto isNewPole = [&](c128 val) -> std::pair<bool, int> {
             for (auto& [r, order] : poles) 
             {
@@ -223,21 +224,22 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
             return {false, 0};
         };
         
-        size_t sz = static_cast<size_t>(processor->filterState->poles.size()) - poles.size();
-        for (size_t i=0 ; i< sz; i++)
+        for (size_t i=0 ; i< prev_sz; i++)
         {
-            auto *r = processor->filterState->poles.getFirst();
+            auto *r = processor->filterState->poles[i];
             std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
             c128 key(c128(r->value.re.get(),r->value.im.get()));
             auto [found, order] = isNewPole(key);
 
             if ( found )
             {
-                r->order= order;
+                r->order= -order;
             }
             else
             {
                 processor->filterState->remove(r);
+                i--;
+                prev_sz--;
             } 
         }
         std::cout<<std::endl;
diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index d366a18..5812599 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -7,20 +7,37 @@
 #include <algorithm>
 #include <iostream>
 #include <cmath>
+#include <limits>
 
 std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<double> coefs)
 {
+
     // filter out leading zeros and leading 1.0
     size_t degree = 0;
-    while(coefs[degree] ==0.0 || coefs[degree] == 1.0)
+    while( coefs[degree] < std::numeric_limits<double>::min() || coefs[degree] - 1.0 < std::numeric_limits<double>::min())
         degree++;
     degree = coefs.size() - degree;
+ std::cout<<"Degree "<<degree<<std::endl;
+   // filter out trailing zeros increasing the order of root at Zero (usually for default poles)
+    size_t orderAtZero = 0;
+    for( size_t i = coefs.size()-1; i >= 0 && std::abs(coefs[i]) < std::numeric_limits<double>::min() ; --i)
+    {
+        ++orderAtZero;
+        --degree;
+    }
     jassert( degree!=0 );
+    
+    // append root at zero of "orderAtZero" order
+    // size_t numRoots = (orderAtZero>0) ? degree+1 : degree;
+    std::vector<std::pair<c128, int>> roots;
+    if (orderAtZero>0)
+        roots.emplace_back(std::make_pair(0.0, orderAtZero));
 
     if (degree==1)
     {
-        std::cout<<"Returing 1 root --> "<<coefs[degree]<<std::endl;
-        return {std::make_pair(static_cast<c128>(-coefs[coefs.size() - 1]), 1)};
+        std::cout<<"Returing 1 non-zero root --> "<<coefs[degree]<<std::endl;
+        roots.emplace_back( std::make_pair(static_cast<c128>(-coefs[coefs.size() - 1]), 1) );
+        return roots;
     }
 
     // build companion Matrix
@@ -30,7 +47,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
     for (size_t i=0; i< degree; i++)
     {
         size_t j = i+1;
-        size_t coef_idx = coefs.size()-1-i;
+        size_t coef_idx = coefs.size()-1 - orderAtZero -i;
         A[i][degree-1] = -coefs[coef_idx];  // fill last column with negative vals of coefs
         if (i!=degree-1)
         {
@@ -153,7 +170,8 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<d
     }
 
     // Extract roots — read diagonal
-    return extractRoots(A, degree);
+    extractRoots(roots, A, degree);
+    return roots;
 }
 
 void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &matrix)
@@ -169,13 +187,13 @@ void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &mat
     }
 }
 
-std::vector<std::pair<c128, int>> CoefficientsToRoots::extractRoots(const std::vector<std::vector<double>>& M, size_t degree)
+void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots, const std::vector<std::vector<double>>& M, size_t degree)
 {
 
     std::cout<<"// Check Roots"<<std::endl;
     CoefficientsToRoots::printCheck(M);
     
-    std::vector<std::pair<c128, int>> roots;
+    // std::vector<std::pair<c128, int>> roots;
 
     auto addRoot = [&](c128 newVal) {
         for (auto& [val, order] : roots)
@@ -247,5 +265,4 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::extractRoots(const std::v
         std::cout<<"("<<root.first.real()<<","<<root.first.imag()<<") - "<<root.second<<", ";
     std::cout<<std::endl;
 #endif
-    return roots;
 }
\ No newline at end of file
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index 2089276..4dde842 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -18,5 +18,5 @@ class CoefficientsToRoots
 		static constexpr double tolerance = 5e-2; // threshold for considering two roots with negligible diff the same. Expressed in % after scaling differences, since zeros can much more than 1.0.
 
 		static void printCheck(const std::vector<std::vector<double>> &);
-		static std::vector<std::pair<c128, int>> extractRoots(const std::vector<std::vector<double>>&, size_t);
+		static void extractRoots(std::vector<std::pair<c128, int>> &, const std::vector<std::vector<double>>&, size_t);
 };
\ No newline at end of file

From 307dd97b9bb9dba2cbc2530e14cb19241f757eaa Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Wed, 10 Jun 2026 12:58:04 +0300
Subject: [PATCH 10/16] 2 fixes:

Fix bug when degree equals to zero. Now the roots (including potential roots at (0,0)), are returned

reverse last fix on parsing leading zeros by using std::numeric_limits instead of direct comparison. Now we compare directly using the == operator.
---
 src/CoefficientsToRoots.cpp | 21 +++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index 5812599..ba6b1f1 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -7,35 +7,36 @@
 #include <algorithm>
 #include <iostream>
 #include <cmath>
-#include <limits>
 
 std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<double> coefs)
 {
 
-    // filter out leading zeros and leading 1.0
+    // filter out leading zeros, increasing counter until the leading 1.0 is found.
     size_t degree = 0;
-    while( coefs[degree] < std::numeric_limits<double>::min() || coefs[degree] - 1.0 < std::numeric_limits<double>::min())
+    while( degree < coefs.size() && coefs[degree] != 1.0)
         degree++;
-    degree = coefs.size() - degree;
- std::cout<<"Degree "<<degree<<std::endl;
+    degree = coefs.size() - degree -1; // subtract 1 for the leading 1.0
+ 
    // filter out trailing zeros increasing the order of root at Zero (usually for default poles)
     size_t orderAtZero = 0;
-    for( size_t i = coefs.size()-1; i >= 0 && std::abs(coefs[i]) < std::numeric_limits<double>::min() ; --i)
+    for( int i = static_cast<int>(coefs.size()-1); i >= 0 && std::abs(coefs[(size_t)i]) == 0.0 ; --i)
     {
         ++orderAtZero;
         --degree;
     }
-    jassert( degree!=0 );
-    
+
     // append root at zero of "orderAtZero" order
-    // size_t numRoots = (orderAtZero>0) ? degree+1 : degree;
     std::vector<std::pair<c128, int>> roots;
     if (orderAtZero>0)
         roots.emplace_back(std::make_pair(0.0, orderAtZero));
 
+    if (degree == 0 )
+        // Returing root at (0,0) if any
+        return roots;
+
     if (degree==1)
     {
-        std::cout<<"Returing 1 non-zero root --> "<<coefs[degree]<<std::endl;
+        // Returing 1 non-zero root + root at (0,0) if any
         roots.emplace_back( std::make_pair(static_cast<c128>(-coefs[coefs.size() - 1]), 1) );
         return roots;
     }

From 5332c6db6dbba1cc738dc6750091196cc2fe9789 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Wed, 10 Jun 2026 14:39:23 +0300
Subject: [PATCH 11/16] add documentation on CoefficientsToRoots class and
 change method name from GramSchmidt to QR

---
 src/CoeffComponents.cpp     |  4 ++--
 src/CoefficientsToRoots.cpp |  4 ++--
 src/CoefficientsToRoots.h   | 40 +++++++++++++++++++++++++++++++++++--
 3 files changed, 42 insertions(+), 6 deletions(-)

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index fbd29a1..2ba61b5 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -149,7 +149,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
     if (col == 2)
     {
 
-        auto zeros = CoefficientsToRoots::GramSchmidt(this->ffcoeffs);
+        auto zeros = CoefficientsToRoots::QR(this->ffcoeffs);
 
 #ifdef DEBUG_C2R
         std::cout<<"Previous Zeros"<<std::endl;
@@ -191,7 +191,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
     }
     else if (col == 3)
     {
-        auto poles =  CoefficientsToRoots::GramSchmidt(this->fbcoeffs);
+        auto poles =  CoefficientsToRoots::QR(this->fbcoeffs);
 
         size_t prev_sz = static_cast<size_t>(processor->filterState->poles.size());
 
diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index ba6b1f1..c1993ee 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -8,7 +8,7 @@
 #include <iostream>
 #include <cmath>
 
-std::vector<std::pair<c128, int>> CoefficientsToRoots::GramSchmidt(std::vector<double> coefs)
+std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> coefs)
 {
 
     // filter out leading zeros, increasing counter until the leading 1.0 is found.
@@ -254,7 +254,7 @@ void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots
                 const double re = (a+d) / 2.0;
                 const double im = std::sqrt(-discriminant) / 2.0;
                 addRoot(c128(re,im));
-                // addRoot(c128(re,-im)); // Note: this is added automatically using addRoot. Commenting this, removes the bug of overlapping roots.
+                // addRoot(c128(re,-im)); // Note: this is added automatically later using FilterState::add method. Commenting this, removes the bug of overlapping roots.
             }
             i += 2;
         }
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index 4dde842..65fcc1a 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -8,15 +8,51 @@
 
 class CoefficientsToRoots
 {
+	/*	
+		Class to convert polynomial coefficients into roots using a QR method.
+	*/
+
 	public:
-		static std::vector<std::pair<c128, int>> GramSchmidt(std::vector<double> coefs); // TODO change return type to bool --> for when there is a pole outside the unit. Check will be done outside of this.
+		/*	Method used to convert polynomial coefficients into roots using a QR method.
+			Consists of the following steps:
+				- Preprocessing:
+					- Filter out leading zeros if any
+					- Filter out trailing zeros if any
+					- Build companion matrix of Hessenberg form from given coefficients
+					- Calculate rayleigh quotient shift 
+				- apply shifted QR iterations (starting from bottom right) with deflation when subdiagonal entries become ~0
+						- Step 1 calculate Q0 via Gram Schmidt orthogonalization method
+						- Step 2 calculate R0 :  R0 = Q0 A0
+						- Step 3 calculate A1 : A1 = R0 Q0
+						- Step 4 add shift back in A1 and check sub-diagonal entries for deflation
+				- Extract eigenvalues:
+					- merge roots that have smaller scaled difference than tolerance
+					- scan subdiagonal:
+						- if values on the subdiagonal smaller than Epsilon, then real root
+						- else 2 solutions. Solve det(A - λI) = 0 for 2x2 square matrix:
+							- if discriminant >= 0, then two real solutions
+							- else then root has a complex conjugate (conjugate is not appended on returned vector, as this is done automatically when using FilterState::add to add that root in the Value tree)
+			Returns complex roots paired with their corresponding order.
+		 */
+		static std::vector<std::pair<c128, int>> QR(std::vector<double> coefs); // TODO change return type to bool --> for when there is a pole outside the unit. Check will be done outside of this.
 	private:
 
 		// TODO Finetune these parameters
+
+		/*	Threshold for detecting convergence (near-zero) of the subdiagonal elements in QR iteration.*/
 	    static constexpr double Epsilon = 1e-3;
+
+		/*	Maximum QR iterations per eigenvalue block to prevent infinite loops.*/
 	    static constexpr size_t MaxIterations = 100;
-		static constexpr double tolerance = 5e-2; // threshold for considering two roots with negligible diff the same. Expressed in % after scaling differences, since zeros can much more than 1.0.
 
+		/* 	Threshold for considering two roots with negligible diff the same. 
+			Expressed in % after scaling differences, since zeros may lie outside the unit circle. */
+		static constexpr double tolerance = 5e-2; 
+
+		/*	Method added to debug - will be removed when cleaning out the files */
 		static void printCheck(const std::vector<std::vector<double>> &);
+		
+		/*	Extracts roots from the eigenvalues of the converged quasi-triangular QR matrix and merges duplicates. 
+			For more details see description of QR method */
 		static void extractRoots(std::vector<std::pair<c128, int>> &, const std::vector<std::vector<double>>&, size_t);
 };
\ No newline at end of file

From b2d6f3a493c7614184cc2c3434770db6112148ca Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Wed, 10 Jun 2026 14:49:49 +0300
Subject: [PATCH 12/16] clean print statements

---
 src/CoeffComponents.cpp     | 56 +------------------------------------
 src/CoefficientsToRoots.cpp | 34 +---------------------
 src/CoefficientsToRoots.h   |  5 ----
 3 files changed, 2 insertions(+), 93 deletions(-)

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index 2ba61b5..c400deb 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -139,8 +139,6 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
 {
     // calculate roots through the coefficient2roots function && notify other listeners about this change
 
-    // @TODO Clean print statements (after ensuring robust implementation)...
-
     if (col == 2)
         ffcoeffs[static_cast<size_t>(row)] = value;
     else if (col == 3)
@@ -148,46 +146,20 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
     
     if (col == 2)
     {
-
         auto zeros = CoefficientsToRoots::QR(this->ffcoeffs);
 
-#ifdef DEBUG_C2R
-        std::cout<<"Previous Zeros"<<std::endl;
-        for (int i=0 ; i< processor->filterState->zeros.size(); i++)
-        {
-            auto *r = processor->filterState->zeros[i];
-            std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
-        }std::cout<<std::endl;
-#endif
-
-        std::cout<<"Removing Zeroes"<<std::endl;
         while (!processor->filterState->zeros.isEmpty())
         {
             auto *r = processor->filterState->zeros.getFirst();
-            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
             processor->filterState->remove(r);
         }
-        std::cout<<std::endl;
 
-        std::cout<<"Adding Zeros"<<std::endl;
         for (size_t i=0 ; i< zeros.size(); i++)
         {
             auto &r = zeros[i].first;
             auto &order = zeros[i].second;
-            std::cout<<"\t("<<r.real()<<","<< r.imag()<<") - "<<order<<", ";
             processor->filterState->add(order, r);
         }
-        std::cout<<std::endl;
-
-#ifdef DEBUG_C2R
-        std::cout<<"Final Zeros"<<std::endl;
-        for (size_t i=0 ; i< processor->filterState->zeros.size(); i++)
-        {
-            auto *r = processor->filterState->zeros[i];
-
-            std::cout<<"("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
-        }std::cout<<std::endl;
-#endif
     }
     else if (col == 3)
     {
@@ -195,26 +167,14 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
 
         size_t prev_sz = static_cast<size_t>(processor->filterState->poles.size());
 
-#ifdef DEBUG_C2R
-        std::cout<<"Previous Poles"<<std::endl;
-        for (size_t i=0 ; i< processor->filterState->poles.size(); i++)
-        {
-            auto *r = processor->filterState->poles[i];
-            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
-        }std::cout<<std::endl;
-#endif
-        std::cout<<"Adding Poles"<<std::endl;
         for (size_t i=0 ; i< poles.size(); i++)
         {
             
             auto &r = poles[i].first;
             auto &order = poles[i].second;
             processor->filterState->add(-order, r);
-            std::cout<<"\t("<<r.real()<<","<< r.imag()<<") - "<<order<<", ";
         }
-        std::cout<<std::endl;
 
-        std::cout<<"Removing Poles"<<std::endl;
         auto isNewPole = [&](c128 val) -> std::pair<bool, int> {
             for (auto& [r, order] : poles) 
             {
@@ -227,7 +187,6 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
         for (size_t i=0 ; i< prev_sz; i++)
         {
             auto *r = processor->filterState->poles[i];
-            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
             c128 key(c128(r->value.re.get(),r->value.im.get()));
             auto [found, order] = isNewPole(key);
 
@@ -242,20 +201,7 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
                 prev_sz--;
             } 
         }
-        std::cout<<std::endl;
-#ifdef DEBUG_C2R
-        std::cout<<"Final Poles"<<std::endl;
-        for (size_t i=0 ; i< processor->filterState->poles.size(); i++)
-        {
-            auto *r = processor->filterState->poles[i];
-            std::cout<<"\t("<<r->value.re.get()<<","<< r->value.im.get()<<") - "<<r->order<<", ";
-        }std::cout<<std::endl;
-#endif
-    }
-#ifdef DEBUG_C2R
-    std::cout<<"Root order "<<processor->filterState->totalOrder <<" finiteZerosOrder "<< processor->filterState->finiteZerosOrder<<std::endl;
-#endif
-    
+    }    
 }
 
 void CoefficientsComponent::valueTreePropertyChanged (juce::ValueTree& node, const juce::Identifier& property)
diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index c1993ee..c3eab15 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -5,7 +5,6 @@
 #include "CoefficientsToRoots.h"
 
 #include <algorithm>
-#include <iostream>
 #include <cmath>
 
 std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> coefs)
@@ -42,7 +41,6 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
     }
 
     // build companion Matrix
-    std::cout<<"Building companion matrix with degree "<<degree<<std::endl;
     std::vector<std::vector<double>> A(degree, std::vector<double>(degree, 0.0));
     // Build companion matrix (Hessenberg form)
     for (size_t i=0; i< degree; i++)
@@ -56,11 +54,6 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
         }
     }
 
-        
-    std::cout<<"// Check companion matrix"<<std::endl;
-    printCheck(A);
-
-
     // QR algorithm
     std::vector<std::vector<double>> Q(degree, std::vector<double>(degree));
     std::vector<std::vector<double>> R(degree, std::vector<double>(degree));
@@ -175,27 +168,9 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
     return roots;
 }
 
-void CoefficientsToRoots::printCheck(const std::vector<std::vector<double>> &matrix)
-{
-    size_t n = matrix.size();
-    for (size_t i=0; i< n; i++)
-    {
-        for (size_t j=0; j<n; j++)
-        {
-            std::cout<<matrix[i][j]<<" ";
-        }
-        std::cout<<std::endl;
-    }
-}
-
 void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots, const std::vector<std::vector<double>>& M, size_t degree)
 {
-
-    std::cout<<"// Check Roots"<<std::endl;
-    CoefficientsToRoots::printCheck(M);
-    
-    // std::vector<std::pair<c128, int>> roots;
-
+        
     auto addRoot = [&](c128 newVal) {
         for (auto& [val, order] : roots)
         {
@@ -259,11 +234,4 @@ void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots
             i += 2;
         }
     }
-
-#ifdef DEBUG_C2R
-    std::cout<<"Calculated Roots: ";
-    for (auto root : roots)
-        std::cout<<"("<<root.first.real()<<","<<root.first.imag()<<") - "<<root.second<<", ";
-    std::cout<<std::endl;
-#endif
 }
\ No newline at end of file
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index 65fcc1a..3caf706 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -4,8 +4,6 @@
 #include <vector>
 #include <utility>
 
-// #define DEBUG_C2R
-
 class CoefficientsToRoots
 {
 	/*	
@@ -48,9 +46,6 @@ class CoefficientsToRoots
 		/* 	Threshold for considering two roots with negligible diff the same. 
 			Expressed in % after scaling differences, since zeros may lie outside the unit circle. */
 		static constexpr double tolerance = 5e-2; 
-
-		/*	Method added to debug - will be removed when cleaning out the files */
-		static void printCheck(const std::vector<std::vector<double>> &);
 		
 		/*	Extracts roots from the eigenvalues of the converged quasi-triangular QR matrix and merges duplicates. 
 			For more details see description of QR method */

From 50f60fecd5325bb04640d5d9b7f2ab5873768341 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Sat, 13 Jun 2026 13:15:06 +0300
Subject: [PATCH 13/16] add tests for Coefficients2Roots::QR function (order
 computation only). An empirical tolerance variable is used allowing errors of
 1 order per 4 orders (scaling per root)

---
 tests/CoefficientsToRootsTest.h | 190 ++++++++++++++++++++++++++++++++
 tests/Main.cpp                  |   1 +
 2 files changed, 191 insertions(+)
 create mode 100644 tests/CoefficientsToRootsTest.h

diff --git a/tests/CoefficientsToRootsTest.h b/tests/CoefficientsToRootsTest.h
new file mode 100644
index 0000000..277b931
--- /dev/null
+++ b/tests/CoefficientsToRootsTest.h
@@ -0,0 +1,190 @@
+#pragma once
+#include <juce_core/juce_core.h>
+#include <juce_dsp/juce_dsp.h>
+#include <juce_data_structures/juce_data_structures.h>
+#include "TestHelper.h"
+#include "../src/PluginProcessor.h"
+#include "../src/CoefficientsToRoots.h"
+#include "../src/RootsToCoefficients.h"
+
+
+class CoefficientsToRootsTest : public juce::UnitTest
+{
+public:
+    CoefficientsToRootsTest() : UnitTest("CoefficientsToRootsTest", "Math") 
+	{
+	}
+
+    void runTest() override
+    {
+		AudioPluginAudioProcessor processor; // shouldn't be a field of this class as its static object is defined.
+
+		performTest("1:1 real pole order 1", 							processor, { {-1, -0.5, 0} } );
+		performTest("2:1 real pole order 2", 							processor, { {-2, -0.5, 0} } );
+		performTest("2:2 real poles order 1",  							processor, { {-1, -0.5, 0}, {-1, 0.3, 0} } );
+		performTest("2:1 complex pole order 1 a", 						processor, { {-1, 0.3, 0.4} } );
+		performTest("2:1 complex pole order 1 b",						processor, { {-1, 0.3, 0.4} } );
+		performTest("2:1 complex pole order 1 c",						processor, { {-1, 0.3, 0.4} } );
+		performTest("3:1 real pole order 2 + order 1",					processor, { {-2, -0.01, 0}, {-1, 0.3, 0} } );
+		performTest("3:1 real pole order 3 a",							processor, { {-3, -0.5, 0} } );
+		performTest("3:1 real pole order 3 b",							processor, { {-3, -0.9, 0} } );
+		performTest("3:3 real poles order 1", 							processor, { {-1, -0.5, 0}, {-1, 0.3, 0}, {-1, 0.2, 0} } );
+		performTest("3:1 complex order 1 + 1 real order 1", 			processor, { {-1, 0.3, 0.4}, {-1, -0.5, 0} } );
+		performTest("3:1 complex order 1 + 1 real order 1 b",			processor, { {-1, 0.3, 0.4}, {-1, -0.5, 0} } );
+		performTest("3:1 complex order 1 + 1 real order 2",				processor, { {-1, 0.3, 0.4}, {-2, -0.5, 0} } );
+		performTest("3:2 real order 1 + 1 real order 2",				processor, { {-1,-0.5,0}, {-1,0.3,0}, {-2,0.2,0} });
+		performTest("4:4 real poles order 1",							processor, { {-1,-0.5,0}, {-1,0.3,0}, {-1,0.2,0}, {-1,0.1,0} } );
+		performTest("4:2 real poles order 2",							processor, { {-2,-0.5,0}, {-2,0.3,0} } );
+		performTest("4:2 complex poles order 1 a",						processor, { {-1,0.3,0.4}, {-1,-0.1,0.3} } );
+		performTest("4:2 complex poles order 1 b",						processor, { {-1,0.3,0.4}, {-1,-0.1,0.3} } );
+		performTest("4:1 real order 3 + 1 real order 1",				processor, { {-3,-0.5,0}, {-1,0.3,0} } );
+		performTest("5:1 real pole order 5",							processor, { {-5,-0.5,0} } );
+		performTest("5:1 real order 4 + 1 real order 1 a",				processor, { {-4,-0.5,0}, {-1,0.1,0} } );
+		performTest("5:1 real order 4 + 1 real order 1 b",				processor, { {-4,-0.5,0}, {-1,0.2,0} } );
+		performTest("5:1 real order 4 + 1 real order 1 c",				processor, { {-4,-0.5,0}, {-1,0.3,0} } );
+		performTest("5:1 real order 4 + 1 real order 1 d",				processor, { {-4,-0.5,0}, {-1,0.3,0} } );
+		performTest("5:1 real order 3 + 1 real order 2 a",				processor, { {-3,-0.5,0}, {-2,0.3,0} } );
+		performTest("5:1 real order 3 + 1 real order 2 b",				processor, { {-3,-0.5,0}, {-2,0.3,0} } );
+		performTest("5:1 real order 3 + 1 complex order 1 a",			processor, { {-3,-0.5,0}, {-1,0.1,0.4} } );
+		performTest("5:1 real order 3 + 1 complex order 1 b",			processor, { {-3,-0.5,0}, {-1,0.2,0.3} } );
+		performTest("5:1 real order 1 + 1 complex order 2 a",			processor, { {-1,-0.5,0}, {-2,0.3,0.1} } );
+		performTest("5:1 real order 1 + 1 complex order 2 b",			processor, { {-1,-0.5,0}, {-2,0.4,0.44} } );
+		performTest("5:2 real order 2 + 1 real order 1",				processor, { {-2,-0.5,0}, {-2,0.3,0}, {-1,0.2,0} } );
+		performTest("5:1 real order 2 + 3 real order 1",				processor, { {-2,-0.5,0}, {-1,0.3,0}, {-1,0.2,0}, {-1,0.1,0} } );
+		performTest("5:5 real poles order 1 a",							processor, { {-1,-0.5,0}, {-1,0.3,0}, {-1,0.2,0}, {-1,0.1,0}, {-1,-0.3,0} } );
+		performTest("5:5 real poles order 1 b", 						processor, { {-1,-0.1,0}, {-1,0.2,0}, {-1,0.3,0}, {-1,0.4,0}, {-1,-0.5,0} } );
+		performTest("5:1 real order 1 + 1 complex + 2 real order 1 a", 	processor, { {-1,-0.5,0}, {-1,0.3,0.4}, {-1,0.2,0}, {-1,0.1,0} } );
+		performTest("5:1 real order 1 + 1 complex + 2 real order 1 b", 	processor, { {-1,-0.2,0}, {-1,0.1,0.4}, {-1,0.4,0}, {-1,0.8,0} } );
+		performTest("5:1 real order 1 + 2 complex order 1 a",			processor, { {-1,-0.1,0}, {-1,0.3,0.4}, {-1,-0.1,0.3} } );
+		performTest("5:1 real order 1 + 2 complex order 1 b",			processor, { {-1,-0.9,0}, {-1,0.13,0.4}, {-1,-0.13,0.53} } );
+		performTest("6:1 real pole order 6 a",							processor, { {-6,-0.01,0} } );
+		performTest("6:1 real pole order 6 b",							processor, { {-6,-0.99,0} } );
+		performTest("6:1 real order 5 + 1 real order 1",				processor, { {-5,-0.5,0}, {-1,0.3,0} } );
+		performTest("6:1 real order 4 + 1 complex order 1 a",			processor, { {-4,-0.05,0}, {-1,0.13,0.14} } );
+		performTest("6:1 real order 4 + 1 complex order 1 b",			processor, { {-4,-0.95,0}, {-1,0.31,0.41} } );
+		performTest("6:1 real order 4 + 2 real order 1 a",				processor, { {-4,-0.1,0}, {-1,0.35,0}, {-1,0.12,0} } );
+		performTest("6:1 real order 4 + 2 real order 1 b",				processor, { {-4,-0.6,0}, {-1,0.63,0}, {-1,0.21,0} } );
+		performTest("6:1 real order 3 + 3 real order 1",				processor, { {-3,-0.5,0}, {-1,0.3,0}, {-1,0.2,0}, {-1,0.1,0} } );
+		performTest("6:1 real order 3 + 1 complex + 1 real order 1",	processor, { {-3,-0.5,0}, {-1,0.3,0.4}, {-1,0.2,0} } );
+		performTest("6:1 complex pole order 3",							processor, { {-3,0.99,0.4} } );
+		performTest("6:1 complex pole order 3 b",						processor, { {-3,0.1,0.4} } );
+		performTest("6:3 real poles order 2",							processor, { {-2,-0.5,0}, {-2,0.3,0}, {-2,0.2,0} } );
+		performTest("6:2 real order 2 + 1 complex order 1",				processor, { {-2,-0.5,0}, {-2,0.3,0}, {-1,0.2,0.3} } );
+		performTest("6:1 real order 2 + 2 complex order 1",				processor, { {-2,-0.5,0}, {-1,0.3,0.4}, {-1,-0.1,0.3} } );
+		performTest("6:3 complex poles order 1",						processor, { {-1,0.3,0.4}, {-1,-0.1,0.3}, {-1,0.1,0.2} } );
+		performTest("7:1 real pole order 7 a",							processor, { {-7,-0.9,0} } );
+		performTest("7:1 real pole order 7 b",							processor, { {-7,-0.1,0} } );
+		performTest("7:1 real order 6 + 1 real order 1",				processor, { {-6,-0.5,0}, {-1,0.3,0} } );
+		performTest("7:1 real order 5 + 1 complex order 1 a",			processor, { {-5,-0.5,0}, {-1,0.9,0.1} } );
+		performTest("7:1 real order 5 + 1 complex order 1 b",			processor, { {-5,-0.5,0}, {-1,0.3,0.4} } );
+		performTest("7:1 real order 4 + 1 complex + 1 real order 1",	processor, { {-4,-0.5,0}, {-1,0.3,0.4}, {-1,0.2,0} } );
+		performTest("7:1 real order 4 + 3 real order 1",				processor, { {-4,-0.5,0}, {-1,0.3,0}, {-1,0.2,0}, {-1,0.1,0} } );
+		performTest("7:1 real order 3 + 1 complex order 2",				processor, { {-3,-0.5,0}, {-2,0.3,0.4} } );
+		performTest("7:1 complex order 3 + 1 real order 1",	 			processor, { {-3,0.3,0.4}, {-1,-0.5,0} } );
+		performTest("7:3 complex + 1 real order 1",						processor, { {-1,0.3,0.4}, {-1,-0.1,0.3}, {-1,0.1,0.2}, {-1,-0.5,0} } );
+		performTest("7:1 complex + 5 real order 1",						processor, { {-1,0.3,0.4}, {-1,-0.5,0}, {-1,0.2,0}, {-1,0.1,0}, {-1,-0.3,0}, {-1,-0.2,0} } );
+		performTest("7:7 real poles order 1",							processor, { {-1,-0.5,0},{-1,0.3,0},{-1,0.2,0},{-1,0.1,0},{-1,-0.3,0},{-1,-0.2,0},{-1,0.4,0} } );
+		performTest("8:1 real pole order 8",							processor, { {-8,-0.5,0} } );
+		performTest("8:1 real order 7 + 1 real order 1",				processor, { {-7,-0.5,0}, {-1,0.3,0} } );
+		performTest("8:1 real order 6 + 1 complex order 1",				processor, { {-6,-0.5,0}, {-1,0.3,0.4} } );
+		performTest("8:8 real poles order 1",							processor, { {-1,-0.5,0},{-1,0.3,0},{-1,0.2,0},{-1,0.1,0},{-1,-0.3,0},{-1,-0.2,0},{-1,0.4,0},{-1,-0.4,0} } );
+		performTest("8:1 complex + 6 real order 1",						processor, { {-1,0.3,0.4},{-1,-0.5,0},{-1,0.2,0},{-1,0.1,0},{-1,-0.3,0},{-1,-0.2,0},{-1,0.4,0} } );
+		performTest("8:2 complex + 4 real order 1",						processor, { {-1,0.3,0.4},{-1,-0.1,0.3},{-1,-0.5,0},{-1,0.2,0},{-1,0.1,0},{-1,-0.3,0} } );
+		performTest("8:4 complex poles order 1", 						processor, { {-1,0.3,0.4},{-1,-0.1,0.3},{-1,0.1,0.2},{-1,-0.2,0.1} } );
+		performTest("9:1 real pole order 9",							processor, { {-9,-0.5,0} } );
+		performTest("9 real poles order 1",								processor, { {-1,-0.5,0},{-1,0.3,0},{-1,0.2,0},{-1,0.1,0},{-1,-0.3,0},{-1,-0.2,0},{-1,0.4,0},{-1,-0.4,0},{-1,0.45,0} } );
+		performTest("9:1 complex + 7 real order 1", 					processor, { {-1,0.3,0.4},{-1,-0.5,0},{-1,0.2,0},{-1,0.1,0},{-1,-0.3,0},{-1,-0.2,0},{-1,0.4,0},{-1,-0.4,0} } );
+		performTest("9:1 real order 5 + 1 real order 4",  				processor, { {-5,-0.5,0}, {-4,0.3,0} } );
+		performTest("9:1 real order 5 + 1 complex order 2", 			processor, { {-5,-0.5,0}, {-2,0.3,0.4} } );
+		performTest("10:1 real pole order 10",							processor, { {-10,-0.5,0} } );
+		performTest("10:2 real poles order 5",							processor, { {-5,-0.5,0}, {-5,0.3,0} } );
+		performTest("10:1 real order 5 + order 4 + order 1", 			processor, { {-5,-0.5,0}, {-4,0.3,0}, {-1,0.2,0} } );
+		performTest("10:1 real order 4 + order 3 + order 2 + order 1", 	processor, { {-4,-0.5,0},{-3,0.3,0},{-2,0.2,0},{-1,0.1,0} } );
+		performTest("10:1 real order 4 + order 3 + complex + order 1", 	processor, { {-4,-0.5,0},{-3,0.3,0},{-1,0.2,0.3},{-1,0.1,0} } );
+		performTest("14:1 real pole order 14",							processor, { {-14,-0.5,0} } );
+		performTest("20:1 real pole order 20",							processor, { {-20,-0.5,0} } );
+		performTest("28:1 real pole order 28",							processor, { {-28,-0.5,0} } );
+		performTest("32:1 real pole order 32",							processor, { {-32,-0.5,0} } );
+	}
+
+private:
+	void performTest(
+		const juce::String testName,
+		AudioPluginAudioProcessor& processor,
+		std::vector<std::array<double, 3>> givenRoots)
+	{
+		beginTest(testName);
+		
+		int zeros_order{0}, poles_order{0};
+		size_t zeros_orderTolerance{0}, poles_orderTolerance{0};
+		for (auto r : givenRoots)
+		{
+			int order = static_cast<int>(r[0]) * (r[2] == 0 ? 1 : 2);
+			jassert(std::abs(order) > 0);
+			if (order>0)
+			{
+				zeros_order+=order;
+				zeros_orderTolerance+= (r[0] / RootOrderToleranceStep);	// TODO: consider log scale growth instead of linear
+			}
+			else
+			{
+				poles_order+=order;
+				poles_orderTolerance+= (std::abs(r[0]) / RootOrderToleranceStep);	// TODO: consider log scale growth instead of linear
+			}
+		}
+		// causality
+		if (zeros_order > -poles_order)
+			poles_order = -zeros_order;
+
+		std::cout<<"zeros_orderTolerance "<<zeros_orderTolerance<<std::endl;
+		std::cout<<"poles_orderTolerance "<<poles_orderTolerance<<std::endl;
+		
+		auto* state = processor.filterState.get();
+		TestHelper::makeFilterState(state, givenRoots, 1);
+
+        if (!state->zeros.isEmpty())
+        {
+            auto ffcoeffs = RootsToCoefficients::CalculatePolynomialCoefficientsFrom(state->zeros);
+            auto newZeros = CoefficientsToRoots::QR(ffcoeffs);
+
+            // Expect to have the same order
+			int curr_order{0};
+			for (auto r : newZeros)
+			{
+				curr_order += r.second * (r.first.imag() == 0 ? 1 : 2);
+			}
+			
+			// expectEquals( static_cast<int>(zeros_order), curr_order );
+			const size_t absDiff {static_cast<size_t>(std::abs(static_cast<int>(zeros_order) - curr_order))};
+		    expectLessOrEqual( absDiff, zeros_orderTolerance);
+			if (absDiff)
+				std::cout<<"Warning : computation error { total order = "<<zeros_order<<" != "<<curr_order<< " by "<<absDiff<<"}"<<std::endl;
+        }
+		
+		if (!state->poles.isEmpty())
+        {
+            auto fbcoeffs = RootsToCoefficients::CalculatePolynomialCoefficientsFrom(state->poles);
+			auto newPoles = CoefficientsToRoots::QR(fbcoeffs);
+
+            // Expect to have the same order
+			int curr_order{0};
+			for (auto r : newPoles)
+			{
+				curr_order -= r.second * (r.first.imag() == 0 ? 1 : 2); // have to reverse sign, since pole-logic is outside GramSchmidt function
+			}
+
+			// expectEquals( static_cast<int>(poles_order), curr_order );
+			const size_t absDiff {static_cast<size_t>(std::abs(static_cast<int>(poles_order) - curr_order))};
+		    expectLessOrEqual( absDiff , poles_orderTolerance);
+			if (absDiff)
+				std::cout<<"Warning : computation error { total order = "<<poles_order<<" != "<<curr_order<< " by "<<absDiff<<"}"<<std::endl;
+			
+
+        }
+	}
+
+	/* Empirical value : Root computatation accuracy degrades for roots of higher-order as their order increases. This step allows 1 order of error per 4 orders.*/
+	static constexpr int RootOrderToleranceStep {4};
+};
+
+static CoefficientsToRootsTest CoefficientsToRootsTest;
\ No newline at end of file
diff --git a/tests/Main.cpp b/tests/Main.cpp
index b279ab0..01274f0 100644
--- a/tests/Main.cpp
+++ b/tests/Main.cpp
@@ -2,6 +2,7 @@
 #include "RootsToCoefficientsTest.h"
 #include "ProcessorChainModifierTest.h"
 #include "PhaseFrequencyResponseTest.h"
+#include "CoefficientsToRootsTest.h"
 
 //==============================================================================
 int main (int argc, char* argv[])

From e75f2ae3f6cc5c427abfde5e99b54feac13143dd Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Sun, 14 Jun 2026 10:43:08 +0300
Subject: [PATCH 14/16] 2 significant performance improvements applied
 (optimization and micro-optimization) :

Optimization change: replace 2D vectors for A,Q,R with 1D vectors, applying manual indexing ([row][col] now becomes [row * degree + col]).

Micro-optimization change : Avoid repeated memory allocation by moving the v vector for storing current columns of A outside loops, and setting its size to {max shift_idx}, that equals to degree of the polynomial
---
 src/CoefficientsToRoots.cpp | 56 +++++++++++++++++++------------------
 src/CoefficientsToRoots.h   |  2 +-
 2 files changed, 30 insertions(+), 28 deletions(-)

diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index c3eab15..75090f0 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -41,24 +41,25 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
     }
 
     // build companion Matrix
-    std::vector<std::vector<double>> A(degree, std::vector<double>(degree, 0.0));
+    std::vector<double> A(degree * degree, 0.0);
     // Build companion matrix (Hessenberg form)
     for (size_t i=0; i< degree; i++)
     {
         size_t j = i+1;
         size_t coef_idx = coefs.size()-1 - orderAtZero -i;
-        A[i][degree-1] = -coefs[coef_idx];  // fill last column with negative vals of coefs
+        A[i * degree + (degree-1)] = -coefs[coef_idx];  // fill last column with negative vals of coefs
         if (i!=degree-1)
         {
-            A[j][i] = 1.0;  // fill subdiagonal entries         
+            A[j * degree + i] = 1.0;  // fill subdiagonal entries         
         }
     }
 
     // QR algorithm
-    std::vector<std::vector<double>> Q(degree, std::vector<double>(degree));
-    std::vector<std::vector<double>> R(degree, std::vector<double>(degree));
+    std::vector<double> Q(degree * degree);
+    std::vector<double> R(degree * degree);
     
     size_t shift_idx {degree}, iter {0};
+    std::vector<double> v(shift_idx); // vector for storing current column of A. Note: only the first {0 to (curr val of shift_idx - 1) } indexes are used per iteration.
     while(shift_idx > 1)
     {
 
@@ -70,23 +71,22 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
         {
             for (size_t c = 0; c < shift_idx; ++c)
             {
-                R[r][c] = 0.0;
-                Q[r][c] = 0.0;
+                R[r * degree + c] = 0.0;
+                Q[r * degree + c] = 0.0;
             }
         }
 
         // subtract rayleigh quotient shift
-        const double shift = A[shift_idx - 1][shift_idx - 1];
+        const double shift = A[(shift_idx - 1) * degree + (shift_idx - 1)];
         for (size_t i = 0; i < shift_idx; ++i) 
-            A[i][i] -= shift; // Decompose (Ak - sI)
+            A[i * degree + i] -= shift; // Decompose (Ak - sI)
 
         // step 1 - calculate Q : Q = A[col] - projections onto the previous Q[col]
         for (size_t col = 0 ; col < shift_idx; col++)
         {
             // set v with column A[col]
-            std::vector<double> v(shift_idx);
             for (size_t row = 0 ; row< shift_idx; row++)
-                v[row] = A[row][col];
+                v[row] = A[row * degree + col];
 
             // subtract projection : v[col] = A[col] - proj onto the previous (<col) orthonormal colums of Q.
             for (size_t curr_col = 0; curr_col < col; curr_col++)
@@ -95,13 +95,13 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
                 double dot_product {0.0};
                 for (size_t curr_row = 0; curr_row < shift_idx; ++curr_row)
                 {
-                    dot_product += Q[curr_row][curr_col] * A[curr_row][col];
+                    dot_product += Q[curr_row * degree + curr_col] * A[curr_row * degree + col];
                 }
 
                 // .. and subtract it from the current column vector
                 for (size_t curr_row=0; curr_row < shift_idx; ++curr_row)
                 {
-                    v[curr_row] -= dot_product * Q[curr_row][curr_col];
+                    v[curr_row] -= dot_product * Q[curr_row * degree + curr_col];
                 }
             }
 
@@ -115,7 +115,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
             
             for (size_t k = 0; k < shift_idx; ++k)
             {
-                Q[k][col] = v[k] / norm;
+                Q[k * degree + col] = v[k] / norm;
             }
         }
 
@@ -124,10 +124,12 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
         {
             for (size_t col = 0; col < shift_idx; col++)
             {
+                double sum {0.0};
                 for (size_t inner = 0; inner < shift_idx; inner++)
                 {
-                    R[row][col] += (Q[inner][row] * A[inner][col]);     // Q transposed
+                    sum += (Q[inner * degree + row] * A[inner * degree + col]);     // Q transposed
                 }
+                R[row* degree + col] = sum;
             }
         }
 
@@ -139,24 +141,24 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
                 double sum {0.0};
                 for (size_t inner = 0; inner< shift_idx; ++inner)
                 {
-                    sum += R[row][inner] * Q[inner][col];
+                    sum += R[row * degree + inner] * Q[inner * degree + col];
                 }
-                A[row][col] = sum; // resetting A[row][col]
+                A[row* degree + col] = sum; // resetting A[row][col]
             }
         }
 
 
         // step 4 add shift back in A and check sub-diagonal entries
         for (size_t i = 0; i < shift_idx; ++i) 
-            A[i][i] += shift;
+            A[i * degree + i] += shift;
 
         // check only the last subdiagonal entry (real eigenvalue)
-        if (std::abs(A[shift_idx-1][shift_idx-2]) < Epsilon)
+        if (std::abs(A[(shift_idx-1)* degree + (shift_idx-2)]) < Epsilon)
         {
             shift_idx--;
         }
         // check the entry above the 2x2 block in search of complex conjugate pairs
-        else if (shift_idx > 2 && std::abs(A[shift_idx-2][shift_idx-3]) < Epsilon)
+        else if (shift_idx > 2 && std::abs(A[(shift_idx-2) * degree + (shift_idx-3)]) < Epsilon)
         {
             shift_idx -= 2;
         }
@@ -168,7 +170,7 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
     return roots;
 }
 
-void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots, const std::vector<std::vector<double>>& M, size_t degree)
+void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots, const std::vector<double>& M, size_t degree)
 {
         
     auto addRoot = [&](c128 newVal) {
@@ -191,10 +193,10 @@ void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots
     while (i < degree)
     {
         
-        if ( i == degree - 1 || std::abs(M[i+1][i]) < Epsilon )
+        if ( i == degree - 1 || std::abs(M[(i+1) * degree + i]) < Epsilon )
         {
             // Real eigenvalue on diagonal
-            c128 newRoot (M[i][i], 0.0);
+            c128 newRoot (M[i * degree + i], 0.0);
             addRoot(newRoot);
             ++i;
         }
@@ -209,10 +211,10 @@ void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots
             // => (a-λ)(d-λ) - bc = 0
             // => λ^2 - (a+d)λ + (ad-bc) = 0
             // => λ^2 - (a+d)λ + detA = 0 --> solve with discriminant
-            const double a = M[i][i];
-            const double b = M[i][i+1];
-            const double c = M[i+1][i];
-            const double d = M[i+1][i+1];
+            const double a = M[i * degree + i];
+            const double b = M[i * degree + i+1];
+            const double c = M[(i+1) * degree + i];
+            const double d = M[(i+1) * degree + i+1];
             const double det  = a*d - b*c;
             // solve with discriminant
             const double discriminant = (a+d)*(a+d) - 4.0*det;
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index 3caf706..7ff6be4 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -49,5 +49,5 @@ class CoefficientsToRoots
 		
 		/*	Extracts roots from the eigenvalues of the converged quasi-triangular QR matrix and merges duplicates. 
 			For more details see description of QR method */
-		static void extractRoots(std::vector<std::pair<c128, int>> &, const std::vector<std::vector<double>>&, size_t);
+		static void extractRoots(std::vector<std::pair<c128, int>> &, const std::vector<double>&, size_t);
 };
\ No newline at end of file

From f4bb67303c032bf165f7031495fff769a5977145 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Sun, 14 Jun 2026 14:12:31 +0300
Subject: [PATCH 15/16] apply micro-opt by replacing divisions with
 multiplications a) on the loop over subdiabonal elements and b) by avoiding
 double call of sqrt + replacing division with 2.0 by multipying with 0.5
 instead on the extractRoots method

---
 src/CoefficientsToRoots.cpp | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/CoefficientsToRoots.cpp b/src/CoefficientsToRoots.cpp
index 75090f0..3894a4c 100644
--- a/src/CoefficientsToRoots.cpp
+++ b/src/CoefficientsToRoots.cpp
@@ -111,11 +111,12 @@ std::vector<std::pair<c128, int>> CoefficientsToRoots::QR(std::vector<double> co
             {
                 norm += v[k] * v[k];
             }
-            norm = std::sqrt(norm);
+            
+            double normReciprocal = 1.0 / std::sqrt(norm);
             
             for (size_t k = 0; k < shift_idx; ++k)
             {
-                Q[k * degree + col] = v[k] / norm;
+                Q[k * degree + col] = v[k] * normReciprocal;
             }
         }
 
@@ -219,17 +220,20 @@ void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots
             // solve with discriminant
             const double discriminant = (a+d)*(a+d) - 4.0*det;
 
+            const double halfSum = 0.5 * (a + d);
+            const double halfSqrt = 0.5 * std::sqrt(std::abs(discriminant));
+
             if ( discriminant >= 0.0)
             {
                 // tow real roots
-                addRoot(c128(((a+d) + std::sqrt(discriminant)) / 2.0, 0.0));
-                addRoot(c128(((a+d) - std::sqrt(discriminant)) / 2.0, 0.0));
+                addRoot(c128(halfSum + halfSqrt, 0.0));
+                addRoot(c128(halfSum - halfSqrt, 0.0));
             }
             else
             {
                 // Complex conjugate pair
-                const double re = (a+d) / 2.0;
-                const double im = std::sqrt(-discriminant) / 2.0;
+                const double re = halfSum;
+                const double im = halfSqrt;
                 addRoot(c128(re,im));
                 // addRoot(c128(re,-im)); // Note: this is added automatically later using FilterState::add method. Commenting this, removes the bug of overlapping roots.
             }

From 445fdd635e83ab1f0064d0f45ecb24a1fc94c044 Mon Sep 17 00:00:00 2001
From: Paschalis Melissas <melissaspaschalis@gmail.com>
Date: Sun, 14 Jun 2026 16:09:51 +0300
Subject: [PATCH 16/16] skip updating of the filterstate in case that computed
 poles are violating causality constraint (at least one pole outside the unit
 circle - magnitude >= 1.0)

---
 src/CoeffComponents.cpp   | 22 +++++++++++++++++++++-
 src/CoefficientsToRoots.h |  2 +-
 2 files changed, 22 insertions(+), 2 deletions(-)

diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index c400deb..b17ed63 100644
--- a/src/CoeffComponents.cpp
+++ b/src/CoeffComponents.cpp
@@ -134,7 +134,6 @@ juce::Component* CoefficientsComponent::refreshComponentForCell(int row, int col
     return label;
 }
 
-
 void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double value)
 {
     // calculate roots through the coefficient2roots function && notify other listeners about this change
@@ -165,6 +164,27 @@ void CoefficientsComponent::updateFilterStateOnCoefEdit(int row, int col, double
     {
         auto poles =  CoefficientsToRoots::QR(this->fbcoeffs);
 
+        // check stability 
+        auto filterStable = [&]() -> bool {
+            for (auto& [r, order] : poles)
+            {
+                const double re = r.real();
+                const double im = r.imag();
+                const double magnitude { std::sqrt ( re * re + im * im) };
+                if ( magnitude >= processor->filterState->maxPoleMagnitude)
+                    return false;
+            }
+            return true;
+        };
+
+        if (!filterStable())
+        {
+            // @TODO handle that case, reporting back to the user a message about the causality violation
+            // unstable filter, do nothing & return
+            return;
+        }
+
+
         size_t prev_sz = static_cast<size_t>(processor->filterState->poles.size());
 
         for (size_t i=0 ; i< poles.size(); i++)
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
index 7ff6be4..05b2891 100644
--- a/src/CoefficientsToRoots.h
+++ b/src/CoefficientsToRoots.h
@@ -32,7 +32,7 @@ class CoefficientsToRoots
 							- else then root has a complex conjugate (conjugate is not appended on returned vector, as this is done automatically when using FilterState::add to add that root in the Value tree)
 			Returns complex roots paired with their corresponding order.
 		 */
-		static std::vector<std::pair<c128, int>> QR(std::vector<double> coefs); // TODO change return type to bool --> for when there is a pole outside the unit. Check will be done outside of this.
+		static std::vector<std::pair<c128, int>> QR(std::vector<double> coefs);
 	private:
 
 		// TODO Finetune these parameters