Merge branch 'bugfix/mat_det' into 'master'

Bugfix for mat class in determinant calculation.

See merge request idf/esp-dsp!125

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased] 
 
 ### Changed
+- Bugfix for determinant calculation in mat.cpp 
 
 ### Added
 

modules/conv/test/test_dsps_conv_f32_ae32.c

@@ -14,6 +14,7 @@
 
 #include <string.h>
 #include <math.h>
+#include <malloc.h>
 #include "unity.h"
 #include "dsp_platform.h"
 #include "esp_log.h"
@@ -27,14 +28,15 @@ static const char *TAG = "dsps_conv";
 #define lenA  30
 #define lenB  30
 
-static float inputA[lenA];
-static float inputB[lenB];
-static float output_ref[lenA + lenB - 1 + 2];
-static float output_fwd[lenA + lenB - 1 + 2];
-static float output_back[lenA + lenB - 1 + 2];
-
 TEST_CASE("dsps_conv_f32 test output", "[dsps]")
 {
+    float *inputA = (float *)memalign(16, lenA * sizeof(float));
+    float *inputB = (float *)memalign(16, lenB * sizeof(float));
+
+    float *output_ref = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+    float *output_fwd = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+    float *output_back = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+
     int la = 3;
     int lb = 2;
 
@@ -58,14 +60,26 @@ TEST_CASE("dsps_conv_f32 test output", "[dsps]")
     float max_eps = 0.000001;
     for (size_t i = 0; i < (la + lb + 1); i++) {
         if (fabs(output_ref[i] - output_fwd[i]) > max_eps) {
-            ESP_LOGE(TAG, "la=%i, lb=%i, i=%i, ref=%2.3f, fwd=%2.3f", la, lb, i, output_ref[i], output_fwd[i]);
+            ESP_LOGI(TAG, "la=%i, lb=%i, i=%i, ref=%2.3f, fwd=%2.3f", la, lb, i, output_ref[i], output_fwd[i]);
         }
         TEST_ASSERT_EQUAL(output_ref[i], output_fwd[i]);
     }
+    free(inputA);
+    free(inputB);
+    free(output_ref);
+    free(output_fwd);
+    free(output_back);
 }
 
 TEST_CASE("dsps_conv_f32 functionality", "[dsps]")
 {
+    float *inputA = (float *)memalign(16, lenA * sizeof(float));
+    float *inputB = (float *)memalign(16, lenB * sizeof(float));
+
+    float *output_ref = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+    float *output_fwd = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+    float *output_back = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+
     for (size_t la = 2; la < lenA; la++) {
         for (size_t lb = 2; lb < lenB; lb++) {
             for (int i = 0 ; i < lenA ; i++) {
@@ -85,14 +99,19 @@ TEST_CASE("dsps_conv_f32 functionality", "[dsps]")
             float max_eps = 0.000001;
             for (size_t i = 0; i < (la + lb + 1); i++) {
                 if ((fabs(output_ref[i] - output_fwd[i]) > max_eps) || (fabs(output_ref[i] - output_back[i]) > max_eps) || (fabs(output_back[i] - output_fwd[i]) > max_eps)) {
-                    ESP_LOGE(TAG, "la=%i, lb=%i, i=%i, ref=%2.3f, fwd=%2.3f, back=%2.3f", la, lb, i, output_ref[i], output_fwd[i], output_back[i]);
+                    ESP_LOGI(TAG, "la=%i, lb=%i, i=%i, ref=%2.3f, fwd=%2.3f, back=%2.3f", la, lb, i, output_ref[i], output_fwd[i], output_back[i]);
                 }
                 TEST_ASSERT_EQUAL(output_ref[i], output_fwd[i]);
                 TEST_ASSERT_EQUAL(output_ref[i], output_back[i]);
                 TEST_ASSERT_EQUAL(output_back[i], output_fwd[i]);
             }
         }
     }
+    free(inputA);
+    free(inputB);
+    free(output_ref);
+    free(output_fwd);
+    free(output_back);
 }
 
 

modules/conv/test/test_dsps_conv_f32_ansi.c

@@ -14,10 +14,12 @@
 
 #include <string.h>
 #include <math.h>
+#include <malloc.h>
 #include "unity.h"
 #include "dsp_platform.h"
 #include "esp_log.h"
 
+#include "dsp_tests.h"
 #include "dsps_conv.h"
 #include "esp_attr.h"
 #include "esp_dsp.h"
@@ -27,13 +29,6 @@ static const char *TAG = "dsps_conv";
 #define lenA  20
 #define lenB  20
 
-static float inputA[lenA];
-static float inputB[lenB];
-static float output_ref[lenA + lenB - 1 + 2];
-static float output_fwd[lenA + lenB - 1 + 2];
-static float output_back[lenA + lenB - 1 + 2];
-
-
 esp_err_t dsps_conv_f32_ref(const float *Signal, const int siglen, const float *Kernel, const int kernlen, float *convout)
 {
     if (NULL == Signal) {
@@ -46,7 +41,6 @@ esp_err_t dsps_conv_f32_ref(const float *Signal, const int siglen, const float *
         return ESP_ERR_DSP_PARAM_OUTOFRANGE;
     }
 
-
     for (int n = 0; n < siglen + kernlen - 1; n++) {
         size_t kmin, kmax, k;
 
@@ -64,6 +58,13 @@ esp_err_t dsps_conv_f32_ref(const float *Signal, const int siglen, const float *
 
 TEST_CASE("dsps_conv_f32_ansi functionality", "[dsps]")
 {
+    float *inputA = (float *)memalign(16, lenA * sizeof(float));
+    float *inputB = (float *)memalign(16, lenB * sizeof(float));
+
+    float *output_ref = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+    float *output_fwd = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+    float *output_back = (float *)memalign(16, (lenA + lenB - 1 + 2) * sizeof(float));
+
     for (size_t la = 1; la < lenA; la++) {
         for (size_t lb = 1; lb < lenB; lb++) {
             for (int i = 0 ; i < lenA ; i++) {
@@ -83,14 +84,19 @@ TEST_CASE("dsps_conv_f32_ansi functionality", "[dsps]")
             float max_eps = 0.000001;
             for (size_t i = 0; i < (la + lb + 1); i++) {
                 if ((fabs(output_ref[i] - output_fwd[i]) > max_eps) || (fabs(output_ref[i] - output_back[i]) > max_eps) || (fabs(output_back[i] - output_fwd[i]) > max_eps)) {
-                    ESP_LOGE(TAG, "la=%i, lb=%i, i=%i, ref=%2.3f, fwd=%2.3f, back=%2.3f", la, lb, i, output_ref[i], output_fwd[i], output_back[i]);
+                    ESP_LOGI(TAG, "la=%i, lb=%i, i=%i, ref=%2.3f, fwd=%2.3f, back=%2.3f", la, lb, i, output_ref[i], output_fwd[i], output_back[i]);
                 }
                 TEST_ASSERT_EQUAL(output_ref[i], output_fwd[i]);
                 TEST_ASSERT_EQUAL(output_ref[i], output_back[i]);
                 TEST_ASSERT_EQUAL(output_back[i], output_fwd[i]);
             }
         }
     }
+    free(inputA);
+    free(inputB);
+    free(output_ref);
+    free(output_fwd);
+    free(output_back);
 }
 
 TEST_CASE("dsps_conv_f32_ansi draw", "[dsps]")

modules/matrix/mat/mat.cpp

@@ -743,28 +743,39 @@ Mat Mat::cofactor(int row, int col, int n)
 
 float Mat::det(int n)
 {
-    float D = 0; // Initialize result
-
     //  Base case : if matrix contains single element
     if (n == 1) {
         return (*this)(0, 0);
     }
 
-    Mat temp(this->rows, this->rows); // To store cofactors
-
-    int sign = 1;  // To store sign multiplier
-
-    // Iterate for each element of first row
-    for (int f = 0; f < n; f++) {
-        // Getting Cofactor of A[0][f]
-        Mat temp = this->cofactor(0, f, n);
-        D += (*this)(0, f) * temp.det(n - 1) * sign;
+    float det = 1.0;
+    Mat *temp = new Mat(n, n);
+    *temp = *this;
 
-        // terms are to be added with alternate sign
-        sign = -sign;
+    for (int i = 0; i < n; i++) {
+        int pivot = i;
+        for (int j = i + 1; j < n; j++) {
+            if (std::abs((*temp)(j, i)) > std::abs((*temp)(pivot, i))) {
+                pivot = j;
+            }
+        }
+        if (pivot != i) {
+            temp->swapRows(i, pivot);
+            det *= -1;
+        }
+        if ((*temp)(i, i) == 0) {
+            return 0;
+        }
+        det *= (*temp)(i, i);
+        for (int j = i + 1; j < n; j++) {
+            float factor = (*temp)(j, i) / (*temp)(i, i);
+            for (int k = i + 1; k < n; k++) {
+                (*temp)(j, k) -= factor * (*temp)(i, k);
+            }
+        }
     }
-
-    return D;
+    delete temp;
+    return det;
 }
 
 Mat Mat::adjoint()

modules/matrix/mul/test/test_mat_f32.cpp

@@ -238,7 +238,7 @@ TEST_CASE("Mat class operators", "[dspm]")
     std::cout << "inverse: " << std::endl;
     std::cout << result << std::endl;
     for (int i = 0 ; i < 3 * 3 ; i++) {
-        if (std::abs(result.data[i] - m_result[i]) > 1e-8) {
+        if (std::abs(result.data[i] - m_result[i]) > 1e-4) {
             printf("Error at[%i] = %f, expected= %f, calculated = %f \n", i, std::abs(result.data[i] - m_result[i]), m_result[i], result.data[i]);
             TEST_ASSERT_MESSAGE (false, "Error in inverse() operation!\n");
         }
@@ -257,3 +257,14 @@ TEST_CASE("Mat class operators", "[dspm]")
 
     delete[] check_array;
 }
+
+TEST_CASE("mat.cpp functionality", "[dsps]")
+{
+    int max_size = 10;
+    for (int i = 3 ; i < max_size ; i++) {
+        dspm::Mat A = dspm::Mat::eye(i);
+        float det = A.det(i);
+        printf("Det[%i] = %f\n", i, det);
+        TEST_ASSERT_EQUAL(det, 1);
+    }
+}

modules/matrix/mul/test/test_mat_sub_f32.cpp

@@ -821,8 +821,8 @@ static void test_mat_subset_inverse(void)
     std::cout << "inverse: " << std::endl;
     std::cout << result << std::endl;
     for (int i = 0; i < 3 * 3; i++) {
-        if (std::abs(result.data[i] - m_result[i]) > 1e-8) {
-            printf("Error at[%i] = %f, expected= %f, calculated = %f \n", i, std::abs(result.data[i] - m_result[i]), m_result[i], result.data[i]);
+        if (std::abs(result.data[i] - m_result[i]) > 1e-4) {
+            printf("Error at[%i] = %f, expected= %f, calculated = %f\n", i, std::abs(result.data[i] - m_result[i]), m_result[i], result.data[i]);
             TEST_ASSERT_MESSAGE (false, "Error in inverse() operation!\n");
         }
     }

modules/support/sfdr/test/test_dsps_sfdr_f32.c

@@ -24,11 +24,10 @@
 
 static const char *TAG = "dsps_sfdr_f32";
 
-static float data[1024];
-
 TEST_CASE("dsps_sfdr_f32 functionality", "[dsps]")
 {
-    int N = sizeof(data) / sizeof(float) / 2;
+    int N = 512;
+    float *data = (float *)malloc(N * 2 * sizeof(float));
     int check_bin = 32;
     float sfdr_exp = 4;
     for (int i = 0 ; i < N ; i++) {
@@ -40,4 +39,5 @@ TEST_CASE("dsps_sfdr_f32 functionality", "[dsps]")
     TEST_ASSERT_EQUAL( (int)20 * log10(sfdr_exp), (int)sfdr);
     ESP_LOGI(TAG, "dsps_sfdr_f32 = %f dB", sfdr);
     dsps_fft2r_deinit_fc32();
+    free(data);
 }

modules/support/snr/test/test_dsps_snr_f32.c

@@ -24,11 +24,11 @@
 
 static const char *TAG = "dsps_snr_f32";
 
-static float data[1024];
 
 TEST_CASE("dsps_snr_f32 functionality", "[dsps]")
 {
-    int N = sizeof(data) / sizeof(float) / 2;
+    int N = 512;
+    float *data = (float *)malloc(N * 2 * sizeof(float));
     int check_bin = 32;
     float snr_exp = 0.001;
     for (int i = 0 ; i < N ; i++) {
@@ -40,4 +40,5 @@ TEST_CASE("dsps_snr_f32 functionality", "[dsps]")
     TEST_ASSERT_EQUAL(-round(20 * log10(snr_exp) + 3), (int)round(snr));
     ESP_LOGI(TAG, "dsps_snr_f32 = %f dB", snr);
     dsps_fft2r_deinit_fc32();
+    free(data);
 }

modules/support/view/test/test_dsps_view.c

@@ -22,15 +22,15 @@
 
 static const char *TAG = "dsps_view";
 
-static float data[1024];
-
 TEST_CASE("dsps_view functionality", "[dsps]")
 {
+    float *data = (float *)malloc(1024 * sizeof(float));
     for (int i = 0 ; i < 1024 ; i++) {
         data[i] = -100;
     }
     data[256] = 0;
     dsps_view_spectrum(data, 1024, -100, 0);
 
     ESP_LOGI(TAG, "Just a check\n");
+    free(data);
 }

modules/windows/test/test_wind_f32.c

@@ -21,11 +21,11 @@
 #include "dsps_wind.h"
 
 static const int length = 1024;
-static float data[1024];
 
 // This test check if the window is symmetric
 TEST_CASE("dsps_wind_hann_f32: test Hann window for symmetry", "[dsps]")
 {
+    float *data = (float *)malloc(length * sizeof(float));
     dsps_wind_hann_f32(data, length);
     float hann_diff = 0;
     for (int i = 0 ; i < length / 2 ; i++) {
@@ -36,10 +36,12 @@ TEST_CASE("dsps_wind_hann_f32: test Hann window for symmetry", "[dsps]")
         TEST_ASSERT_EQUAL(0, hann_diff);
     }
     dsps_view(data, length, 64, 10, 0, 1, '.');
+    free(data);
 }
 
 TEST_CASE("dsps_wind_blackman_f32: test Blackman window for symmetry", "[dsps]")
 {
+    float *data = (float *)malloc(length * sizeof(float));
     dsps_wind_blackman_f32(data, length);
     float hann_diff = 0;
     for (int i = 0 ; i < length / 2 ; i++) {
@@ -49,10 +51,12 @@ TEST_CASE("dsps_wind_blackman_f32: test Blackman window for symmetry", "[dsps]")
         TEST_ASSERT_EQUAL(0, hann_diff);
     }
     dsps_view(data, length, 64, 10, 0, 1, '.');
+    free(data);
 }
 
 TEST_CASE("dsps_wind_blackman_harris_f32: test Blackman-Hariss window for symmetry", "[dsps]")
 {
+    float *data = (float *)malloc(length * sizeof(float));
     dsps_wind_blackman_harris_f32(data, length);
     float hann_diff = 0;
     for (int i = 0 ; i < length / 2 ; i++) {
@@ -63,10 +67,12 @@ TEST_CASE("dsps_wind_blackman_harris_f32: test Blackman-Hariss window for symmet
         TEST_ASSERT_EQUAL(0, hann_diff);
     }
     dsps_view(data, length, 64, 10, 0, 1, '.');
+    free(data);
 }
 
 TEST_CASE("dsps_wind_blackman_nuttall_f32: test Blackman-Nuttall window for symmetry", "[dsps]")
 {
+    float *data = (float *)malloc(length * sizeof(float));
     dsps_wind_blackman_nuttall_f32(data, length);
     float hann_diff = 0;
     for (int i = 0 ; i < length / 2 ; i++) {
@@ -77,10 +83,12 @@ TEST_CASE("dsps_wind_blackman_nuttall_f32: test Blackman-Nuttall window for symm
         TEST_ASSERT_EQUAL(0, hann_diff);
     }
     dsps_view(data, length, 64, 10, 0, 1, '.');
+    free(data);
 }
 
 TEST_CASE("dsps_wind_nuttall_f32: test Nuttall window for symmetry", "[dsps]")
 {
+    float *data = (float *)malloc(length * sizeof(float));
     dsps_wind_nuttall_f32(data, length);
     float hann_diff = 0;
     for (int i = 0 ; i < length / 2 ; i++) {
@@ -91,10 +99,12 @@ TEST_CASE("dsps_wind_nuttall_f32: test Nuttall window for symmetry", "[dsps]")
         TEST_ASSERT_EQUAL(0, hann_diff);
     }
     dsps_view(data, length, 64, 10, 0, 1, '.');
+    free(data);
 }
 
 TEST_CASE("dsps_wind_flat_top_f32: test Flat-Top window for symmetry", "[dsps]")
 {
+    float *data = (float *)malloc(length * sizeof(float));
     dsps_wind_flat_top_f32(data, length);
     float hann_diff = 0;
     for (int i = 0 ; i < length / 2 ; i++) {
@@ -105,4 +115,5 @@ TEST_CASE("dsps_wind_flat_top_f32: test Flat-Top window for symmetry", "[dsps]")
         TEST_ASSERT_EQUAL(0, hann_diff);
     }
     dsps_view(data, length, 64, 10, 0, 1, '.');
+    free(data);
 }