diff --git a/include/eigen_polyfit.hpp b/include/eigen_polyfit.hpp
index 3b2d592..b7760a3 100644
--- a/include/eigen_polyfit.hpp
+++ b/include/eigen_polyfit.hpp
@@ -18,22 +18,20 @@ std::vector<T> polyfit(const std::vector<T> &yValues, const int degree) {
     int coeffs = degree + 1;
     size_t samples = yValues.size();
     
-    MatrixXf X(samples, coeffs);
-    MatrixXf Y(samples, 1);
+    MatrixXf A(samples, coeffs);
+    MatrixXf b(samples, 1);
     
-    // fill Y matrix
-    for (size_t i = 0; i < samples; i++) {
-        Y(i, 0) = yValues[i];
-    }
-    
-    // fill X matrix (Vandermonde matrix)
+    // fill b matrix
+    // fill A matrix (Vandermonde matrix)
     for (size_t row = 0; row < samples; row++) {
         for (int col = 0; col < coeffs; col++) {
-            X(row, col) = std::pow(row, col);
+            A(row, col) = std::pow(row, col);
         }
+        b(row) = yValues[row];
     }
     
+    // Solve Ax = b
     VectorXf coefficients;
-    coefficients = X.jacobiSvd(ComputeThinU | ComputeThinV).solve(Y);
+    coefficients = A.householderQr().solve(b);
     return std::vector<T>(coefficients.data(), coefficients.data() + coeffs);
 }
diff --git a/src/guess_tempo.cpp b/src/guess_tempo.cpp
index 3dfe0a3..308411d 100644
--- a/src/guess_tempo.cpp
+++ b/src/guess_tempo.cpp
@@ -76,7 +76,7 @@ std::vector<BPMFitness> estimate_bpm(const std::set<std::size_t>& onsets, const
     {
         std::ofstream broad_fitness_dump("broad_fitness.csv");
         broad_fitness_dump << "interval,fitness\n";
-        for (const auto& [interval, fitness, max_onset] : broad_fitness) {
+        for (const auto& [interval, fitness, _] : broad_fitness) {
             broad_fitness_dump << interval << "," << fitness << "\n";
         }
     }
@@ -212,6 +212,13 @@ Eigen::ArrayXf correct_bias(const std::vector<IntervalFitness>& fitness_results)
         y_values.begin(),
         [](const IntervalFitness& f){return f.fitness;}
     );
+    {
+        std::ofstream yvalues_dump("yvalues.csv");
+        yvalues_dump << "value\n";
+        for (const auto& value : y_values) {
+            yvalues_dump << value << "\n";
+        }     
+    }
     const auto coeffs = polyfit(y_values, 3);
     {
         std::ofstream coeffs_dump("coeffs.csv");
diff --git a/tests/meson.build b/tests/meson.build
index 4a3cc3c..5f00651 100644
--- a/tests/meson.build
+++ b/tests/meson.build
@@ -68,4 +68,11 @@ executable(
     'more_precise_music.cpp',
     dependencies: deps,
     include_directories: inc,
+)
+
+executable(
+    'test_polyfit',
+    'test_polyfit.cpp',
+    dependencies: deps,
+    include_directories: inc,
 )
\ No newline at end of file
diff --git a/tests/test_polyfit.cpp b/tests/test_polyfit.cpp
new file mode 100644
index 0000000..6402a4e
--- /dev/null
+++ b/tests/test_polyfit.cpp
@@ -0,0 +1,38 @@
+#include "Eigen/Dense"
+#include <iostream>
+
+template<typename T>
+std::vector<T> polyfit(const std::vector<T> &yValues, const int degree) {
+    using namespace Eigen;
+    
+    int coeffs = degree + 1;
+    size_t samples = yValues.size();
+    
+    MatrixXf X(samples, coeffs);
+    MatrixXf Y(samples, 1);
+    
+    // fill Y matrix
+    for (size_t i = 0; i < samples; i++) {
+        Y(i, 0) = yValues[i];
+    }
+    
+    // fill X matrix (Vandermonde matrix)
+    for (size_t row = 0; row < samples; row++) {
+        for (int col = 0; col < coeffs; col++) {
+            X(row, col) = std::pow(row, col);
+        }
+    }
+    
+    VectorXf coefficients;
+    coefficients = X.jacobiSvd(ComputeThinU | ComputeThinV).solve(Y);
+    return std::vector<T>(coefficients.data(), coefficients.data() + coeffs);
+}
+
+int main(void) {
+    // std::vector<float> x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
+    std::vector<float> y = {0.0, 0.8, 0.9, 0.1, -0.8, -1.0};
+    auto coeffs = polyfit(y, 3);
+    for (const auto& coeff : coeffs) {
+        std::cout << coeff << "\n";
+    }
+}
\ No newline at end of file