diff --git a/src/CoeffComponents.cpp b/src/CoeffComponents.cpp
index a17b99f..b17ed63 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.
@@ -114,13 +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
+            updateFilterStateOnCoefEdit(row, col, value);
         };
     }
 
@@ -138,6 +134,96 @@ 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
+
+    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::QR(this->ffcoeffs);
+
+        while (!processor->filterState->zeros.isEmpty())
+        {
+            auto *r = processor->filterState->zeros.getFirst();
+            processor->filterState->remove(r);
+        }
+
+        for (size_t i=0 ; i< zeros.size(); i++)
+        {
+            auto &r = zeros[i].first;
+            auto &order = zeros[i].second;
+            processor->filterState->add(order, r);
+        }
+    }
+    else if (col == 3)
+    {
+        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++)
+        {
+            
+            auto &r = poles[i].first;
+            auto &order = poles[i].second;
+            processor->filterState->add(-order, r);
+        }
+
+        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};
+        };
+        
+        for (size_t i=0 ; i< prev_sz; i++)
+        {
+            auto *r = processor->filterState->poles[i];
+            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);
+                i--;
+                prev_sz--;
+            } 
+        }
+    }    
+}
+
 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 e896f17..8e1ae06 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>
 
@@ -28,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
new file mode 100644
index 0000000..3894a4c
--- /dev/null
+++ b/src/CoefficientsToRoots.cpp
@@ -0,0 +1,243 @@
+// #pragma once
+// #include <juce_audio_processors/juce_audio_processors.h>
+// #include <juce_dsp/juce_dsp.h>
+
+#include "CoefficientsToRoots.h"
+
+#include <algorithm>
+#include <cmath>
+
+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.
+    size_t degree = 0;
+    while( degree < coefs.size() && coefs[degree] != 1.0)
+        degree++;
+    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( int i = static_cast<int>(coefs.size()-1); i >= 0 && std::abs(coefs[(size_t)i]) == 0.0 ; --i)
+    {
+        ++orderAtZero;
+        --degree;
+    }
+
+    // append root at zero of "orderAtZero" order
+    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)
+    {
+        // 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;
+    }
+
+    // build companion Matrix
+    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 + (degree-1)] = -coefs[coef_idx];  // fill last column with negative vals of coefs
+        if (i!=degree-1)
+        {
+            A[j * degree + i] = 1.0;  // fill subdiagonal entries         
+        }
+    }
+
+    // QR algorithm
+    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)
+    {
+
+        if (++iter > degree * MaxIterations) break;
+
+
+        // Reset R,Q
+        for (size_t r = 0; r < shift_idx; ++r)
+        {
+            for (size_t c = 0; c < shift_idx; ++c)
+            {
+                R[r * degree + c] = 0.0;
+                Q[r * degree + c] = 0.0;
+            }
+        }
+
+        // subtract rayleigh quotient shift
+        const double shift = A[(shift_idx - 1) * degree + (shift_idx - 1)];
+        for (size_t i = 0; i < shift_idx; ++i) 
+            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]
+            for (size_t row = 0 ; row< shift_idx; row++)
+                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++)
+            {
+                // compute dot product ...
+                double dot_product {0.0};
+                for (size_t curr_row = 0; curr_row < shift_idx; ++curr_row)
+                {
+                    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 * degree + curr_col];
+                }
+            }
+
+            // normalize column of q
+            double norm = 0.0;
+            for (size_t k = 0; k < shift_idx; ++k)
+            {
+                norm += v[k] * v[k];
+            }
+            
+            double normReciprocal = 1.0 / std::sqrt(norm);
+            
+            for (size_t k = 0; k < shift_idx; ++k)
+            {
+                Q[k * degree + col] = v[k] * normReciprocal;
+            }
+        }
+
+        // step 2 - calculate R :  R = Q A
+        for (size_t row = 0; row < shift_idx; row++)
+        {
+            for (size_t col = 0; col < shift_idx; col++)
+            {
+                double sum {0.0};
+                for (size_t inner = 0; inner < shift_idx; inner++)
+                {
+                    sum += (Q[inner * degree + row] * A[inner * degree + col]);     // Q transposed
+                }
+                R[row* degree + col] = sum;
+            }
+        }
+
+        // step 3 - A = R Q
+        for (size_t row = 0; row < shift_idx; row++)
+        {
+            for (size_t col = 0; col < shift_idx; col++)
+            {
+                double sum {0.0};
+                for (size_t inner = 0; inner< shift_idx; ++inner)
+                {
+                    sum += R[row * degree + inner] * Q[inner * degree + 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 * degree + i] += shift;
+
+        // check only the last subdiagonal entry (real eigenvalue)
+        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) * degree + (shift_idx-3)]) < Epsilon)
+        {
+            shift_idx -= 2;
+        }
+
+    }
+
+    // Extract roots — read diagonal
+    extractRoots(roots, A, degree);
+    return roots;
+}
+
+void CoefficientsToRoots::extractRoots(std::vector<std::pair<c128, int>> & roots, const std::vector<double>& M, size_t degree)
+{
+        
+    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);
+    };
+
+    size_t i = 0;
+    while (i < degree)
+    {
+        
+        if ( i == degree - 1 || std::abs(M[(i+1) * degree + i]) < Epsilon )
+        {
+            // Real eigenvalue on diagonal
+            c128 newRoot (M[i * degree + i], 0.0);
+            addRoot(newRoot);
+            ++i;
+        }
+        else
+        {
+            // 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 * 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;
+
+            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(halfSum + halfSqrt, 0.0));
+                addRoot(c128(halfSum - halfSqrt, 0.0));
+            }
+            else
+            {
+                // Complex conjugate pair
+                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.
+            }
+            i += 2;
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/CoefficientsToRoots.h b/src/CoefficientsToRoots.h
new file mode 100644
index 0000000..05b2891
--- /dev/null
+++ b/src/CoefficientsToRoots.h
@@ -0,0 +1,53 @@
+#pragma once
+
+#include "FilterState.h"
+#include <vector>
+#include <utility>
+
+class CoefficientsToRoots
+{
+	/*	
+		Class to convert polynomial coefficients into roots using a QR method.
+	*/
+
+	public:
+		/*	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);
+	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;
+
+		/* 	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; 
+		
+		/*	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<double>&, size_t);
+};
\ No newline at end of file
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[])