Use fastlog, and fix a discrepancy against Python implementation

classifier/ei_run_classifier.h

@@ -400,7 +400,7 @@ extern "C" EI_IMPULSE_ERROR run_classifier(
  * @param data_fn Function to retrieve data from sensors
  * @param debug Whether to log debug messages (default false)
  */
-EI_IMPULSE_ERROR run_impulse(
+__attribute__((unused)) EI_IMPULSE_ERROR run_impulse(
 #if defined(EI_CLASSIFIER_HAS_SAMPLER) && EI_CLASSIFIER_HAS_SAMPLER == 1
         EdgeSampler *sampler,
 #endif
@@ -457,7 +457,7 @@ EI_IMPULSE_ERROR run_impulse(
  * @param data_fn Function to retrieve data from sensors
  * @param debug Whether to log debug messages (default false)
  */
-EI_IMPULSE_ERROR run_impulse(
+__attribute__((unused)) EI_IMPULSE_ERROR run_impulse(
         ei_impulse_result_t *result,
 #ifdef __MBED__
         Callback<void(float*, size_t)> data_fn,

classifier/ei_run_dsp.h

@@ -35,7 +35,7 @@ namespace {
 
 using namespace ei;
 
-int extract_spectral_analysis_features(signal_t *signal, matrix_t *output_matrix, void *config_ptr) {
+__attribute__((unused)) int extract_spectral_analysis_features(signal_t *signal, matrix_t *output_matrix, void *config_ptr) {
     ei_dsp_config_spectral_analysis_t config = *((ei_dsp_config_spectral_analysis_t*)config_ptr);
 
     int ret;
@@ -102,7 +102,7 @@ int extract_spectral_analysis_features(signal_t *signal, matrix_t *output_matrix
     return EIDSP_OK;
 }
 
-int extract_raw_features(signal_t *signal, matrix_t *output_matrix, void *config_ptr) {
+__attribute__((unused)) int extract_raw_features(signal_t *signal, matrix_t *output_matrix, void *config_ptr) {
     ei_dsp_config_raw_t config = *((ei_dsp_config_raw_t*)config_ptr);
 
     // input matrix from the raw signal
@@ -120,7 +120,7 @@ int extract_raw_features(signal_t *signal, matrix_t *output_matrix, void *config
     return EIDSP_OK;
 }
 
-int extract_flatten_features(signal_t *signal, matrix_t *output_matrix, void *config_ptr) {
+__attribute__((unused)) int extract_flatten_features(signal_t *signal, matrix_t *output_matrix, void *config_ptr) {
     ei_dsp_config_flatten_t config = *((ei_dsp_config_flatten_t*)config_ptr);
 
     size_t expected_matrix_size = 0;
@@ -195,7 +195,7 @@ static int preemphasized_audio_signal_get_data(size_t offset, size_t length, flo
     return preemphasis->get_data(offset, length, out_ptr);
 }
 
-int extract_mfcc_features(signal_t *signal, matrix_t *output_matrix, void *config_ptr) {
+__attribute__((unused)) int extract_mfcc_features(signal_t *signal, matrix_t *output_matrix, void *config_ptr) {
     ei_dsp_config_mfcc_t config = *((ei_dsp_config_mfcc_t*)config_ptr);
 
     if (config.axes != 1) {

dsp/config.hpp

@@ -38,6 +38,9 @@
 #define EIDSP_i16                q15_t
 #define EIDSP_i8                 q7_t
 #define ARM_MATH_ROUNDING        1
+#elif EIDSP_USE_CMSIS_DSP == 0
+#define EIDSP_i16                int16_t
+#define EIDSP_i8                 int8_t
 #endif // EIDSP_USE_CMSIS_DSP
 
 #ifndef EIDSP_USE_ASSERTS

dsp/numpy.hpp

@@ -1149,6 +1149,51 @@ class numpy {
     }
 #endif
 
+    /**
+     * > 50% faster then the math.h log() function
+     * in return for a small loss in accuracy (0.00001 average diff with log())
+     * From: https://stackoverflow.com/questions/39821367/very-fast-approximate-logarithm-natural-log-function-in-c/39822314#39822314
+     * Licensed under the CC BY-SA 3.0
+     * @param a Input number
+     * @returns Natural log value of a
+     */
+    __attribute__((always_inline)) static inline float log(float a)
+    {
+        float m, r, s, t, i, f;
+        int32_t e, g;
+
+        g = (int32_t) * ((int32_t *)&a);
+        e = (g - 0x3f2aaaab) & 0xff800000;
+        g = g - e;
+        m = (float) * ((float *)&g);
+        i = (float)e * 1.19209290e-7f; // 0x1.0p-23
+        /* m in [2/3, 4/3] */
+        f = m - 1.0f;
+        s = f * f;
+        /* Compute log1p(f) for f in [-1/3, 1/3] */
+        r = fmaf(0.230836749f, f, -0.279208571f); // 0x1.d8c0f0p-3, -0x1.1de8dap-2
+        t = fmaf(0.331826031f, f, -0.498910338f); // 0x1.53ca34p-2, -0x1.fee25ap-2
+        r = fmaf(r, s, t);
+        r = fmaf(r, s, f);
+        r = fmaf(i, 0.693147182f, r); // 0x1.62e430p-1 // log(2)
+
+        return r;
+    }
+
+    /**
+     * Calculate the natural log value of a matrix. Does an in-place replacement.
+     * @param matrix Matrix (MxN)
+     * @returns 0 if OK
+     */
+    static int log(matrix_t *matrix)
+    {
+        for (size_t ix = 0; ix < matrix->rows * matrix->cols; ix++) {
+            matrix->buffer[ix] = numpy::log(matrix->buffer[ix]);
+        }
+
+        return EIDSP_OK;
+    }
+
 private:
     static int software_rfft(float *fft_input, float *output, size_t n_fft, size_t n_fft_out_features) {
         kiss_fft_cpx *fft_output = (kiss_fft_cpx*)ei_dsp_malloc(n_fft_out_features * sizeof(kiss_fft_cpx));

dsp/spectral/processing.hpp

@@ -85,7 +85,7 @@ namespace processing {
      * @param scale (float): The scaling factor (multiplies by this number).
      * @returns 0 when successful
      */
-    static int scale(float *signal, size_t signal_size, float scale = 1)
+    __attribute__((unused)) static int scale(float *signal, size_t signal_size, float scale = 1)
     {
         EI_DSP_MATRIX_B(temp, 1, signal_size, signal);
         return numpy::scale(&temp, scale);

dsp/speechpy/feature.hpp

@@ -111,7 +111,11 @@ class feature {
             EIDSP_ERR(EIDSP_OUT_OF_MEM);
         }
         for (uint16_t ix = 0; ix < num_filter + 2; ix++) {
-            freq_index[ix] = static_cast<int>(floor((coefficients + 1) * hertz[ix] / sampling_freq));
+            // here is a really annoying bug in Speechpy which calculates the frequency index wrong for the last bucket
+            // the last 'hertz' value is not 8,000 (with sampling rate 16,000) but 7,999.999999
+            // thus calculating the bucket to 64, not 65.
+            // we're adjusting this here a tiny bit to ensure we have the same result
+            freq_index[ix] = static_cast<int>(floor(((coefficients + 1) * hertz[ix] / sampling_freq) - 0.00001));
         }
         ei_dsp_free(hertz, hertz_mem_size);
 
@@ -214,7 +218,9 @@ class feature {
             EIDSP_ERR(EIDSP_MATRIX_SIZE_MISMATCH);
         }
 
-        memset(out_features->buffer, 0, out_features->rows * out_features->cols * sizeof(float));
+        for(int i = 0; i < out_features->rows * out_features->cols; i++) {
+            *(out_features->buffer + i) = 0;
+        }
 
         uint16_t coefficients = fft_length / 2 + 1;
 
@@ -393,8 +399,9 @@ class feature {
 
         // ok... now we need to calculate the MFCC from this...
         // first do log() over all features...
-        for (size_t ix = 0; ix < features_matrix.rows * features_matrix.cols; ix++) {
-            features_matrix.buffer[ix] = log(features_matrix.buffer[ix]);
+        ret = numpy::log(&features_matrix);
+        if (ret != EIDSP_OK) {
+            EIDSP_ERR(ret);
         }
 
         // now do DST type 2
@@ -406,15 +413,15 @@ class feature {
         // replace first cepstral coefficient with log of frame energy for DC elimination
         if (dc_elimination) {
             for (size_t row = 0; row < features_matrix.rows; row++) {
-                features_matrix.buffer[row * features_matrix.cols] = log(energy_matrix.buffer[row]);
+                features_matrix.buffer[row * features_matrix.cols] = numpy::log(energy_matrix.buffer[row]);
             }
         }
 
         // copy to the output...
         for (size_t row = 0; row < features_matrix.rows; row++) {
-            memcpy(out_features->buffer + (num_cepstral * row),
-                features_matrix.buffer + (features_matrix.cols * row),
-                num_cepstral * sizeof(float));
+            for(int i = 0; i < num_cepstral; i++) {
+                *(out_features->buffer + (num_cepstral * row) + i) = *(features_matrix.buffer + (features_matrix.cols * row) + i);
+            }
         }
 
         return EIDSP_OK;

dsp/speechpy/functions.hpp

@@ -39,7 +39,7 @@ class functions {
      * @returns The mel scale values(or a single mel).
      */
     static float frequency_to_mel(float f) {
-        return 1127.0 * log(1 + f / 700.0f);
+        return 1127.0 * numpy::log(1 + f / 700.0f);
     }
 
     /**

dsp/speechpy/processing.hpp

@@ -148,7 +148,7 @@ namespace processing {
      * @param cof (float): The preemphasising coefficient. 0 equals to no filtering.
      * @returns 0 when successful
      */
-    static int preemphasis(float *signal, size_t signal_size, int shift = 1, float cof = 0.98f)
+    __attribute__((unused)) static int preemphasis(float *signal, size_t signal_size, int shift = 1, float cof = 0.98f)
     {
         if (shift < 0) {
             shift = signal_size + shift;