From 9cd41ba890c7da3e515d10e5a4cbb3a8b37f81ee Mon Sep 17 00:00:00 2001 From: eleroy Date: Tue, 26 Nov 2024 11:50:58 +0100 Subject: [PATCH] Move module to signal. oaconvolve now returns float when input are not complex. It should also work if there is no complex support. Revert filter+numpy changes --- code/numpy/filter.c | 208 +++++-------------------------------- code/numpy/filter.h | 1 - code/numpy/numpy.c | 5 - code/scipy/signal/signal.c | 171 ++++++++++++++++++++++++++++++ code/scipy/signal/signal.h | 1 + code/ulab.h | 4 + 6 files changed, 201 insertions(+), 189 deletions(-) diff --git a/code/numpy/filter.c b/code/numpy/filter.c index 7570ee59..79c1740a 100644 --- a/code/numpy/filter.c +++ b/code/numpy/filter.c @@ -10,7 +10,7 @@ * 2020 Scott Shawcroft for Adafruit Industries * 2020-2021 Zoltán Vörös * 2020 Taku Fukada - */ +*/ #include #include @@ -23,39 +23,33 @@ #include "../scipy/signal/signal.h" #include "carray/carray_tools.h" #include "filter.h" -#include "../numpy/fft/fft_tools.h" // For fft -#include "../ulab_tools.h" #if ULAB_NUMPY_HAS_CONVOLVE -mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) -{ +mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { static const mp_arg_t allowed_args[] = { - {MP_QSTR_a, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE}}, - {MP_QSTR_v, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE}}, + { MP_QSTR_a, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_v, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE } }, }; mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); - if (!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type) || !mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) - { + if(!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type) || !mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) { mp_raise_TypeError(MP_ERROR_TEXT("convolve arguments must be ndarrays")); } ndarray_obj_t *a = MP_OBJ_TO_PTR(args[0].u_obj); ndarray_obj_t *c = MP_OBJ_TO_PTR(args[1].u_obj); -// deal with linear arrays only -#if ULAB_MAX_DIMS > 1 - if ((a->ndim != 1) || (c->ndim != 1)) - { + // deal with linear arrays only + #if ULAB_MAX_DIMS > 1 + if((a->ndim != 1) || (c->ndim != 1)) { mp_raise_TypeError(MP_ERROR_TEXT("convolve arguments must be linear arrays")); } -#endif + #endif size_t len_a = a->len; size_t len_c = c->len; - if (len_a == 0 || len_c == 0) - { + if(len_a == 0 || len_c == 0) { mp_raise_TypeError(MP_ERROR_TEXT("convolve arguments must not be empty")); } @@ -63,12 +57,11 @@ mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a int32_t off = len_c - 1; uint8_t dtype = NDARRAY_FLOAT; -#if ULAB_SUPPORTS_COMPLEX - if ((a->dtype == NDARRAY_COMPLEX) || (c->dtype == NDARRAY_COMPLEX)) - { + #if ULAB_SUPPORTS_COMPLEX + if((a->dtype == NDARRAY_COMPLEX) || (c->dtype == NDARRAY_COMPLEX)) { dtype = NDARRAY_COMPLEX; } -#endif + #endif ndarray_obj_t *ndarray = ndarray_new_linear_array(len, dtype); mp_float_t *array = (mp_float_t *)ndarray->array; @@ -78,41 +71,33 @@ mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a int32_t as = a->strides[ULAB_MAX_DIMS - 1] / a->itemsize; int32_t cs = c->strides[ULAB_MAX_DIMS - 1] / c->itemsize; -#if ULAB_SUPPORTS_COMPLEX - if (dtype == NDARRAY_COMPLEX) - { + + #if ULAB_SUPPORTS_COMPLEX + if(dtype == NDARRAY_COMPLEX) { mp_float_t a_real, a_imag; mp_float_t c_real, c_imag = MICROPY_FLOAT_CONST(0.0); - for (int32_t k = -off; k < len - off; k++) - { + for(int32_t k = -off; k < len-off; k++) { mp_float_t accum_real = MICROPY_FLOAT_CONST(0.0); mp_float_t accum_imag = MICROPY_FLOAT_CONST(0.0); int32_t top_n = MIN(len_c, len_a - k); int32_t bot_n = MAX(-k, 0); - for (int32_t n = bot_n; n < top_n; n++) - { + for(int32_t n = bot_n; n < top_n; n++) { int32_t idx_c = (len_c - n - 1) * cs; int32_t idx_a = (n + k) * as; - if (a->dtype != NDARRAY_COMPLEX) - { + if(a->dtype != NDARRAY_COMPLEX) { a_real = ndarray_get_float_index(aarray, a->dtype, idx_a); a_imag = MICROPY_FLOAT_CONST(0.0); - } - else - { + } else { a_real = ndarray_get_float_index(aarray, NDARRAY_FLOAT, 2 * idx_a); a_imag = ndarray_get_float_index(aarray, NDARRAY_FLOAT, 2 * idx_a + 1); } - if (c->dtype != NDARRAY_COMPLEX) - { + if(c->dtype != NDARRAY_COMPLEX) { c_real = ndarray_get_float_index(carray, c->dtype, idx_c); c_imag = MICROPY_FLOAT_CONST(0.0); - } - else - { + } else { c_real = ndarray_get_float_index(carray, NDARRAY_FLOAT, 2 * idx_c); c_imag = ndarray_get_float_index(carray, NDARRAY_FLOAT, 2 * idx_c + 1); } @@ -124,15 +109,13 @@ mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a } return MP_OBJ_FROM_PTR(ndarray); } -#endif + #endif - for (int32_t k = -off; k < len - off; k++) - { + for(int32_t k = -off; k < len-off; k++) { mp_float_t accum = MICROPY_FLOAT_CONST(0.0); int32_t top_n = MIN(len_c, len_a - k); int32_t bot_n = MAX(-k, 0); - for (int32_t n = bot_n; n < top_n; n++) - { + for(int32_t n = bot_n; n < top_n; n++) { int32_t idx_c = (len_c - n - 1) * cs; int32_t idx_a = (n + k) * as; mp_float_t ai = ndarray_get_float_index(aarray, a->dtype, idx_a); @@ -145,146 +128,5 @@ mp_obj_t filter_convolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_a } MP_DEFINE_CONST_FUN_OBJ_KW(filter_convolve_obj, 2, filter_convolve); -#if ULAB_NUMPY_HAS_FFT_MODULE - -mp_obj_t filter_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) -{ - static const mp_arg_t allowed_args[] = { - {MP_QSTR_a, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE}}, - {MP_QSTR_v, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE}}, - }; - - mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; - mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); - - if (!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type) || !mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) - { - mp_raise_TypeError(MP_ERROR_TEXT("oa convolve arguments must be ndarrays")); - } - - ndarray_obj_t *a = MP_OBJ_TO_PTR(args[0].u_obj); - ndarray_obj_t *c = MP_OBJ_TO_PTR(args[1].u_obj); -// deal with linear arrays only -#if ULAB_MAX_DIMS > 1 - if ((a->ndim != 1) || (c->ndim != 1)) - { - mp_raise_TypeError(MP_ERROR_TEXT("convolve arguments must be linear arrays")); - } -#endif - size_t signal_len = a->len; - size_t kernel_len = c->len; - if (signal_len == 0 || kernel_len == 0) - { - mp_raise_TypeError(MP_ERROR_TEXT("convolve arguments must not be empty")); - } -#if !ULAB_SUPPORTS_COMPLEX - mp_raise_TypeError(MP_ERROR_TEXT("oaconvolve only implemented on port with complex support")); -#endif - size_t segment_len = kernel_len; // Min segment size, at least size of kernel - size_t fft_size = 1; - while (fft_size < segment_len + kernel_len - 1) - { - fft_size *= 2; // Find smallest power of 2 for fft size - } - segment_len = fft_size - kernel_len + 1; // Final segment size - size_t output_len = signal_len + kernel_len - 1; - uint8_t sz = 2 * sizeof(mp_float_t); - - // Buffer allocation - mp_float_t *kernel_fft_array = m_new0(mp_float_t, fft_size * 2); - mp_float_t *segment_fft_array = m_new0(mp_float_t, fft_size * 2); - ndarray_obj_t *output = ndarray_new_linear_array(output_len, NDARRAY_COMPLEX); - mp_float_t *output_array = (mp_float_t *)output->array; - uint8_t *kernel_array = (uint8_t *)c->array; - uint8_t *signal_array = (uint8_t *)a->array; - - // Copy kernel data - if (c->dtype == NDARRAY_COMPLEX) - { - for (size_t i = 0; i < kernel_len; i++) - { - memcpy(kernel_fft_array, kernel_array, sz); - kernel_fft_array += 2; - kernel_array += c->strides[ULAB_MAX_DIMS - 1]; - } - } - else - { - mp_float_t (*func)(void *) = ndarray_get_float_function(c->dtype); - for (size_t i = 0; i < kernel_len; i++) - { - // real part; the imaginary part is 0, no need to assign - *kernel_fft_array = func(kernel_array); - kernel_fft_array += 2; - kernel_array += c->strides[ULAB_MAX_DIMS - 1]; - } - } - kernel_fft_array -= 2 * kernel_len; - - // Compute kernel fft in place - fft_kernel(kernel_fft_array, fft_size, 1); - - mp_float_t real, imag; // For complex multiplication - size_t current_segment_size = segment_len; - mp_float_t (*funca)(void *) = ndarray_get_float_function(a->dtype); - for (size_t i = 0; i < signal_len; i += segment_len) - { - // Check if remaining data is less than segment length - if (i + segment_len > signal_len) - { - current_segment_size = signal_len - i; - } - // Load segment in buffer - memset(segment_fft_array, 0, fft_size * sz); // Reset buffer to zero - if (a->dtype == NDARRAY_COMPLEX) - { - for (size_t k = 0; k < current_segment_size; k++) - { - memcpy(segment_fft_array, signal_array, sz); - segment_fft_array += 2; - signal_array += a->strides[ULAB_MAX_DIMS - 1]; - } - } - else - { - for (size_t k = 0; k < current_segment_size; k++) - { - // real part; the imaginary part is 0, no need to assign - *segment_fft_array = funca(signal_array); - segment_fft_array += 2; - signal_array += a->strides[ULAB_MAX_DIMS - 1]; - } - } - segment_fft_array -= 2 * current_segment_size; - - // Compute segment fft in place - fft_kernel(segment_fft_array, fft_size, 1); - - // Product of kernel and segment fft (complex multiplication) - for (size_t j = 0; j < fft_size; j++) - { - real = segment_fft_array[j * 2] * kernel_fft_array[j * 2] - segment_fft_array[j * 2 + 1] * kernel_fft_array[j * 2 + 1]; - imag = segment_fft_array[j * 2] * kernel_fft_array[j * 2 + 1] + segment_fft_array[j * 2 + 1] * kernel_fft_array[j * 2]; - segment_fft_array[j * 2] = real; - segment_fft_array[j * 2 + 1] = imag; - } - // Inverse FFT in place - fft_kernel(segment_fft_array, fft_size, -1); - - // Overlap, Add and normalized inverse fft - for (size_t j = 0; j < fft_size * 2; j++) - { - if ((i * 2 + j) < (output_len * 2)) - { - output_array[i * 2 + j] += (segment_fft_array[j] / fft_size); - } - } - } - m_free(segment_fft_array); // Useless? - m_free(kernel_fft_array); // Useless? - return MP_OBJ_FROM_PTR(output); -} -MP_DEFINE_CONST_FUN_OBJ_KW(filter_oaconvolve_obj, 2, filter_oaconvolve); -#endif #endif diff --git a/code/numpy/filter.h b/code/numpy/filter.h index 479f926d..d6d0f172 100644 --- a/code/numpy/filter.h +++ b/code/numpy/filter.h @@ -17,5 +17,4 @@ #include "../ndarray.h" MP_DECLARE_CONST_FUN_OBJ_KW(filter_convolve_obj); -MP_DECLARE_CONST_FUN_OBJ_KW(filter_oaconvolve_obj); #endif diff --git a/code/numpy/numpy.c b/code/numpy/numpy.c index 639af7f9..eafd7728 100644 --- a/code/numpy/numpy.c +++ b/code/numpy/numpy.c @@ -216,11 +216,6 @@ static const mp_rom_map_elem_t ulab_numpy_globals_table[] = { #if ULAB_NUMPY_HAS_CONVOLVE { MP_ROM_QSTR(MP_QSTR_convolve), MP_ROM_PTR(&filter_convolve_obj) }, #endif - #if ULAB_NUMPY_HAS_FFT_MODULE - #if ULAB_NUMPY_HAS_CONVOLVE - { MP_ROM_QSTR(MP_QSTR_oaconvolve), MP_ROM_PTR(&filter_oaconvolve_obj) }, - #endif - #endif // functions of the numerical sub-module #if ULAB_NUMPY_HAS_ALL { MP_ROM_QSTR(MP_QSTR_all), MP_ROM_PTR(&numerical_all_obj) }, diff --git a/code/scipy/signal/signal.c b/code/scipy/signal/signal.c index f930a943..932a9e3b 100644 --- a/code/scipy/signal/signal.c +++ b/code/scipy/signal/signal.c @@ -19,6 +19,8 @@ #include "../../ulab.h" #include "../../ndarray.h" #include "../../numpy/carray/carray_tools.h" +#include "../../numpy/fft/fft_tools.h" // For fft +#include "../../ulab_tools.h" #if ULAB_SCIPY_SIGNAL_HAS_SOSFILT & ULAB_MAX_DIMS > 1 static void signal_sosfilt_array(mp_float_t *x, const mp_float_t *coeffs, mp_float_t *zf, const size_t len) { @@ -120,11 +122,180 @@ mp_obj_t signal_sosfilt(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar MP_DEFINE_CONST_FUN_OBJ_KW(signal_sosfilt_obj, 2, signal_sosfilt); #endif /* ULAB_SCIPY_SIGNAL_HAS_SOSFILT */ +#if ULAB_SCIPY_SIGNAL_HAS_OACONVOLVE + +mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) +{ + static const mp_arg_t allowed_args[] = { + {MP_QSTR_a, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE}}, + {MP_QSTR_v, MP_ARG_REQUIRED | MP_ARG_OBJ, {.u_rom_obj = MP_ROM_NONE}}, + }; + + mp_arg_val_t args[MP_ARRAY_SIZE(allowed_args)]; + mp_arg_parse_all(n_args, pos_args, kw_args, MP_ARRAY_SIZE(allowed_args), allowed_args, args); + + if (!mp_obj_is_type(args[0].u_obj, &ulab_ndarray_type) || !mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) + { + mp_raise_TypeError(MP_ERROR_TEXT("oa convolve arguments must be ndarrays")); + } + + ndarray_obj_t *a = MP_OBJ_TO_PTR(args[0].u_obj); + ndarray_obj_t *c = MP_OBJ_TO_PTR(args[1].u_obj); +// deal with linear arrays only +#if ULAB_MAX_DIMS > 1 + if ((a->ndim != 1) || (c->ndim != 1)) + { + mp_raise_TypeError(MP_ERROR_TEXT("convolve arguments must be linear arrays")); + } +#endif + size_t signal_len = a->len; + size_t kernel_len = c->len; + if (signal_len == 0 || kernel_len == 0) + { + mp_raise_TypeError(MP_ERROR_TEXT("convolve arguments must not be empty")); + } + + + + size_t segment_len = kernel_len; // Min segment size, at least size of kernel + size_t fft_size = 1; + while (fft_size < segment_len + kernel_len - 1) + { + fft_size *= 2; // Find smallest power of 2 for fft size + } + segment_len = fft_size - kernel_len + 1; // Final segment size + size_t output_len = signal_len + kernel_len - 1; + uint8_t dtype = NDARRAY_FLOAT; +#if ULAB_SUPPORTS_COMPLEX + if ((a->dtype == NDARRAY_COMPLEX) || (c->dtype == NDARRAY_COMPLEX)) + { + dtype = NDARRAY_COMPLEX; + } +#endif + + uint8_t sz = 2 * sizeof(mp_float_t); + + // Buffer allocation + mp_float_t *kernel_fft_array = m_new0(mp_float_t, fft_size * 2); + mp_float_t *segment_fft_array = m_new0(mp_float_t, fft_size * 2); + ndarray_obj_t *output = ndarray_new_linear_array(output_len, dtype); + mp_float_t *output_array = (mp_float_t *)output->array; + uint8_t *kernel_array = (uint8_t *)c->array; + uint8_t *signal_array = (uint8_t *)a->array; + + // Copy kernel data + if (c->dtype == NDARRAY_COMPLEX) + { + for (size_t i = 0; i < kernel_len; i++) + { + memcpy(kernel_fft_array, kernel_array, sz); + kernel_fft_array += 2; + kernel_array += c->strides[ULAB_MAX_DIMS - 1]; + } + } + else + { + mp_float_t (*func)(void *) = ndarray_get_float_function(c->dtype); + for (size_t i = 0; i < kernel_len; i++) + { + // real part; the imaginary part is 0, no need to assign + *kernel_fft_array = func(kernel_array); + kernel_fft_array += 2; + kernel_array += c->strides[ULAB_MAX_DIMS - 1]; + } + } + kernel_fft_array -= 2 * kernel_len; + + // Compute kernel fft in place + fft_kernel(kernel_fft_array, fft_size, 1); + + mp_float_t real, imag; // For complex multiplication + size_t current_segment_size = segment_len; + mp_float_t (*funca)(void *) = ndarray_get_float_function(a->dtype); + + for (size_t i = 0; i < signal_len; i += segment_len) + { + // Check if remaining data is less than segment length + if (i + segment_len > signal_len) + { + current_segment_size = signal_len - i; + } + // Load segment in buffer + memset(segment_fft_array, 0, fft_size * sz); // Reset buffer to zero + if (a->dtype == NDARRAY_COMPLEX) + { + for (size_t k = 0; k < current_segment_size; k++) + { + memcpy(segment_fft_array, signal_array, sz); + segment_fft_array += 2; + signal_array += a->strides[ULAB_MAX_DIMS - 1]; + } + } + else + { + for (size_t k = 0; k < current_segment_size; k++) + { + // real part; the imaginary part is 0, no need to assign + *segment_fft_array = funca(signal_array); + segment_fft_array += 2; + signal_array += a->strides[ULAB_MAX_DIMS - 1]; + } + } + segment_fft_array -= 2 * current_segment_size; + + // Compute segment fft in place + fft_kernel(segment_fft_array, fft_size, 1); + + // Product of kernel and segment fft (complex multiplication) + for (size_t j = 0; j < fft_size; j++) + { + real = segment_fft_array[j * 2] * kernel_fft_array[j * 2] - segment_fft_array[j * 2 + 1] * kernel_fft_array[j * 2 + 1]; + imag = segment_fft_array[j * 2] * kernel_fft_array[j * 2 + 1] + segment_fft_array[j * 2 + 1] * kernel_fft_array[j * 2]; + segment_fft_array[j * 2] = real; + segment_fft_array[j * 2 + 1] = imag; + } + // Inverse FFT in place + fft_kernel(segment_fft_array, fft_size, -1); + + #if ULAB_SUPPORTS_COMPLEX + if (dtype == NDARRAY_COMPLEX) + { + // Overlap, Add and normalized inverse fft + for (size_t j = 0; j < fft_size * 2; j++) + { + if ((i * 2 + j) < (output_len * 2)) + { + output_array[i * 2 + j] += (segment_fft_array[j] / fft_size); + } + } + continue; + } + #endif + + + for (size_t j = 0; j < fft_size; j++) + { + if ((i + j) < (output_len)) // adds only real part + { + output_array[i + j] += (segment_fft_array[j * 2] / fft_size); + } + } + + } + return MP_OBJ_FROM_PTR(output); +} +MP_DEFINE_CONST_FUN_OBJ_KW(signal_oaconvolve_obj, 2, signal_oaconvolve); +#endif /* ULAB_SCIPY_SIGNAL_HAS_OACONVOLVE */ + + static const mp_rom_map_elem_t ulab_scipy_signal_globals_table[] = { { MP_ROM_QSTR(MP_QSTR___name__), MP_ROM_QSTR(MP_QSTR_signal) }, #if ULAB_SCIPY_SIGNAL_HAS_SOSFILT & ULAB_MAX_DIMS > 1 { MP_ROM_QSTR(MP_QSTR_sosfilt), MP_ROM_PTR(&signal_sosfilt_obj) }, #endif + #if ULAB_SCIPY_SIGNAL_HAS_OACONVOLVE + { MP_ROM_QSTR(MP_QSTR_oaconvolve), MP_ROM_PTR(&signal_oaconvolve_obj) }, + #endif }; static MP_DEFINE_CONST_DICT(mp_module_ulab_scipy_signal_globals, ulab_scipy_signal_globals_table); diff --git a/code/scipy/signal/signal.h b/code/scipy/signal/signal.h index 033f6e4c..3b713187 100644 --- a/code/scipy/signal/signal.h +++ b/code/scipy/signal/signal.h @@ -19,5 +19,6 @@ extern const mp_obj_module_t ulab_scipy_signal_module; MP_DECLARE_CONST_FUN_OBJ_KW(signal_sosfilt_obj); +MP_DECLARE_CONST_FUN_OBJ_KW(signal_oaconvolve_obj); #endif /* _SCIPY_SIGNAL_ */ diff --git a/code/ulab.h b/code/ulab.h index 78b8a1ba..30c237b5 100644 --- a/code/ulab.h +++ b/code/ulab.h @@ -740,6 +740,10 @@ #define ULAB_SCIPY_SIGNAL_HAS_SOSFILT (1) #endif +#ifndef ULAB_SCIPY_SIGNAL_HAS_OACONVOLVE +#define ULAB_SCIPY_SIGNAL_HAS_OACONVOLVE (1) +#endif + #ifndef ULAB_SCIPY_HAS_OPTIMIZE_MODULE #define ULAB_SCIPY_HAS_OPTIMIZE_MODULE (1) #endif