From e38546f54e4fdb3cbc90ef14b43b2b2e526e46df Mon Sep 17 00:00:00 2001 From: eleroy Date: Mon, 25 Nov 2024 23:54:01 +0100 Subject: [PATCH 1/4] - Add oaconvolve method to signal module --- code/numpy/filter.c | 208 ++++++++++++++++++++++++++++++++++++++------ code/numpy/filter.h | 1 + code/numpy/numpy.c | 5 ++ 3 files changed, 189 insertions(+), 25 deletions(-) diff --git a/code/numpy/filter.c b/code/numpy/filter.c index 79c1740a..7570ee59 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,33 +23,39 @@ #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")); } @@ -57,11 +63,12 @@ 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; @@ -71,33 +78,41 @@ 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); } @@ -109,13 +124,15 @@ 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); @@ -128,5 +145,146 @@ 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 d6d0f172..479f926d 100644 --- a/code/numpy/filter.h +++ b/code/numpy/filter.h @@ -17,4 +17,5 @@ #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 eafd7728..639af7f9 100644 --- a/code/numpy/numpy.c +++ b/code/numpy/numpy.c @@ -216,6 +216,11 @@ 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) }, From 9cd41ba890c7da3e515d10e5a4cbb3a8b37f81ee Mon Sep 17 00:00:00 2001 From: eleroy Date: Tue, 26 Nov 2024 11:50:58 +0100 Subject: [PATCH 2/4] 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 From 8571d9eb69b6f602544a31d4bfd1c7eceea6aa50 Mon Sep 17 00:00:00 2001 From: eleroy Date: Tue, 26 Nov 2024 20:12:36 +0100 Subject: [PATCH 3/4] Bump ulab version. Add documentation and test for oaconvolve --- code/ulab.c | 2 +- docs/manual/source/scipy-signal.rst | 36 ++++++++++++++++++++++++++++- tests/1d/scipy/oaconvolve.py | 16 +++++++++++++ tests/1d/scipy/oaconvolve.py.exp | 1 + 4 files changed, 53 insertions(+), 2 deletions(-) create mode 100644 tests/1d/scipy/oaconvolve.py create mode 100644 tests/1d/scipy/oaconvolve.py.exp diff --git a/code/ulab.c b/code/ulab.c index 12f37027..17058454 100644 --- a/code/ulab.c +++ b/code/ulab.c @@ -33,7 +33,7 @@ #include "user/user.h" #include "utils/utils.h" -#define ULAB_VERSION 6.6.1 +#define ULAB_VERSION 6.6.2 #define xstr(s) str(s) #define str(s) #s diff --git a/docs/manual/source/scipy-signal.rst b/docs/manual/source/scipy-signal.rst index d1f34818..0f078eba 100644 --- a/docs/manual/source/scipy-signal.rst +++ b/docs/manual/source/scipy-signal.rst @@ -2,9 +2,10 @@ scipy.signal ============ -This module defines the single function: +This module defines the functions: 1. `scipy.signal.sosfilt <#sosfilt>`__ +2. `scipy.signal.oaconvolve <#oaconvolve>`__ sosfilt ------- @@ -64,6 +65,39 @@ initial values are assumed to be 0. ======================================== zf: array([[37242.0, 74835.0], [1026187.0, 1936542.0]], dtype=float) + +oaconvolve +------- + +``scipy``: +https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.oaconvolve.html + +Convolve two N-dimensional arrays using the overlap-add method. + +Convolve in1 and in2 using the overlap-add method. Similarly to numpy.convolve, +this method works in full mode and other modes can be obtained by slicing the result? + +This is generally much faster than linear convolve for large arrays (n > ~500), +and generally much faster than fftconvolve (not implemented yet in ulab) when one array is much larger +than the other (e.g. for a matched filter algorithm), but can be slower when only a few output values +are needed or when the arrays are very similar in shape, and can only output float arrays (int or object array inputs will be cast to float). + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + from ulab import scipy as spy + + x = np.array((1,2,3)) + y = np.array((1,10,100,1000)) + result = spy.signal.oaconvolve(x, y) + print('result: ', result) + +.. parsed-literal:: + + result: array([1.0, 12.00024, 123.0001, 1230.0, 2300.0, 3000.0], dtype=float32) + diff --git a/tests/1d/scipy/oaconvolve.py b/tests/1d/scipy/oaconvolve.py new file mode 100644 index 00000000..48e8492a --- /dev/null +++ b/tests/1d/scipy/oaconvolve.py @@ -0,0 +1,16 @@ +import math + +try: + from ulab import scipy, numpy as np +except ImportError: + import scipy + import numpy as np + +x = np.array((1,2,3)) +y = np.array((1,10,100,1000)) +result = (scipy.signal.oaconvolve(x, y)) +ref_result = np.array([1, 12, 123, 1230, 2300, 3000],dtype=np.float) +cmp_result = [] +for p,q in zip(list(result), list(ref_result)): + cmp_result.append(math.isclose(p, q, rel_tol=1e-06, abs_tol=5e-04)) +print(cmp_result) diff --git a/tests/1d/scipy/oaconvolve.py.exp b/tests/1d/scipy/oaconvolve.py.exp new file mode 100644 index 00000000..63a3ac63 --- /dev/null +++ b/tests/1d/scipy/oaconvolve.py.exp @@ -0,0 +1 @@ +[True, True, True, True, True, True] From 1ec776c8350d1504758463485eeb77e7eecba1af Mon Sep 17 00:00:00 2001 From: eleroy Date: Wed, 27 Nov 2024 08:39:51 +0100 Subject: [PATCH 4/4] Fix formatting. Add docs in jupyter notebook. --- code/scipy/signal/signal.c | 72 ++++++++++------------------ docs/scipy-signal.ipynb | 96 +++++++++++++++++++++++++++++++------- 2 files changed, 105 insertions(+), 63 deletions(-) diff --git a/code/scipy/signal/signal.c b/code/scipy/signal/signal.c index 932a9e3b..53554e2d 100644 --- a/code/scipy/signal/signal.c +++ b/code/scipy/signal/signal.c @@ -10,6 +10,7 @@ * 2020 Scott Shawcroft for Adafruit Industries * 2020-2021 Zoltán Vörös * 2020 Taku Fukada + * 2024 Edouard Leroy (oaconvolve implementation) */ #include @@ -19,7 +20,7 @@ #include "../../ulab.h" #include "../../ndarray.h" #include "../../numpy/carray/carray_tools.h" -#include "../../numpy/fft/fft_tools.h" // For fft +#include "../../numpy/fft/fft_tools.h" #include "../../ulab_tools.h" #if ULAB_SCIPY_SIGNAL_HAS_SOSFILT & ULAB_MAX_DIMS > 1 @@ -124,8 +125,7 @@ MP_DEFINE_CONST_FUN_OBJ_KW(signal_sosfilt_obj, 2, signal_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) -{ +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}}, @@ -134,8 +134,7 @@ mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw 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("oa convolve arguments must be ndarrays")); } @@ -143,15 +142,13 @@ mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw 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)) - { + 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) - { + if (signal_len == 0 || kernel_len == 0) { mp_raise_TypeError(MP_ERROR_TEXT("convolve arguments must not be empty")); } @@ -159,16 +156,14 @@ mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw 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) - { + 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)) - { + if ((a->dtype == NDARRAY_COMPLEX) || (c->dtype == NDARRAY_COMPLEX)) { dtype = NDARRAY_COMPLEX; } #endif @@ -184,20 +179,16 @@ mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw 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++) - { + 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 - { + else { mp_float_t (*func)(void *) = ndarray_get_float_function(c->dtype); - for (size_t i = 0; i < kernel_len; i++) - { + 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; @@ -213,28 +204,22 @@ mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw 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) - { + 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) - { + 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++) - { + 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++) - { + 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; @@ -247,8 +232,7 @@ mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw 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++) - { + 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; @@ -258,13 +242,10 @@ mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw fft_kernel(segment_fft_array, fft_size, -1); #if ULAB_SUPPORTS_COMPLEX - if (dtype == NDARRAY_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)) - { + 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); } } @@ -273,14 +254,11 @@ mp_obj_t signal_oaconvolve(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw #endif - for (size_t j = 0; j < fft_size; j++) - { - if ((i + j) < (output_len)) // adds only real part - { + 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); } diff --git a/docs/scipy-signal.ipynb b/docs/scipy-signal.ipynb index ec10d2e6..8abfcde5 100644 --- a/docs/scipy-signal.ipynb +++ b/docs/scipy-signal.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2021-01-12T16:11:12.111639Z", @@ -11,10 +11,19 @@ }, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Populating the interactive namespace from numpy and matplotlib\n" + "ename": "ModuleNotFoundError", + "evalue": "No module named 'matplotlib'", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[1], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m \u001b[43mget_ipython\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun_line_magic\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mpylab\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43minline\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\Users\\EL221079\\Documents\\MICROPYTHON\\oaconvolve\\.venv\\Lib\\site-packages\\IPython\\core\\interactiveshell.py:2480\u001b[0m, in \u001b[0;36mInteractiveShell.run_line_magic\u001b[1;34m(self, magic_name, line, _stack_depth)\u001b[0m\n\u001b[0;32m 2478\u001b[0m kwargs[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlocal_ns\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mget_local_scope(stack_depth)\n\u001b[0;32m 2479\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbuiltin_trap:\n\u001b[1;32m-> 2480\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mfn\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 2482\u001b[0m \u001b[38;5;66;03m# The code below prevents the output from being displayed\u001b[39;00m\n\u001b[0;32m 2483\u001b[0m \u001b[38;5;66;03m# when using magics with decorator @output_can_be_silenced\u001b[39;00m\n\u001b[0;32m 2484\u001b[0m \u001b[38;5;66;03m# when the last Python token in the expression is a ';'.\u001b[39;00m\n\u001b[0;32m 2485\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mgetattr\u001b[39m(fn, magic\u001b[38;5;241m.\u001b[39mMAGIC_OUTPUT_CAN_BE_SILENCED, \u001b[38;5;28;01mFalse\u001b[39;00m):\n", + "File \u001b[1;32mc:\\Users\\EL221079\\Documents\\MICROPYTHON\\oaconvolve\\.venv\\Lib\\site-packages\\IPython\\core\\magics\\pylab.py:159\u001b[0m, in \u001b[0;36mPylabMagics.pylab\u001b[1;34m(self, line)\u001b[0m\n\u001b[0;32m 155\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m 156\u001b[0m \u001b[38;5;66;03m# invert no-import flag\u001b[39;00m\n\u001b[0;32m 157\u001b[0m import_all \u001b[38;5;241m=\u001b[39m \u001b[38;5;129;01mnot\u001b[39;00m args\u001b[38;5;241m.\u001b[39mno_import_all\n\u001b[1;32m--> 159\u001b[0m gui, backend, clobbered \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshell\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43menable_pylab\u001b[49m\u001b[43m(\u001b[49m\u001b[43margs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mgui\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mimport_all\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mimport_all\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 160\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_show_matplotlib_backend(args\u001b[38;5;241m.\u001b[39mgui, backend)\n\u001b[0;32m 161\u001b[0m \u001b[38;5;28mprint\u001b[39m(\n\u001b[0;32m 162\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m%\u001b[39m\u001b[38;5;124mpylab is deprecated, use \u001b[39m\u001b[38;5;124m%\u001b[39m\u001b[38;5;124mmatplotlib inline and import the required libraries.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 163\u001b[0m )\n", + "File \u001b[1;32mc:\\Users\\EL221079\\Documents\\MICROPYTHON\\oaconvolve\\.venv\\Lib\\site-packages\\IPython\\core\\interactiveshell.py:3719\u001b[0m, in \u001b[0;36mInteractiveShell.enable_pylab\u001b[1;34m(self, gui, import_all, welcome_message)\u001b[0m\n\u001b[0;32m 3692\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Activate pylab support at runtime.\u001b[39;00m\n\u001b[0;32m 3693\u001b[0m \n\u001b[0;32m 3694\u001b[0m \u001b[38;5;124;03mThis turns on support for matplotlib, preloads into the interactive\u001b[39;00m\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 3715\u001b[0m \u001b[38;5;124;03m This argument is ignored, no welcome message will be displayed.\u001b[39;00m\n\u001b[0;32m 3716\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 3717\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mIPython\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcore\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpylabtools\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m import_pylab\n\u001b[1;32m-> 3719\u001b[0m gui, backend \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43menable_matplotlib\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgui\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 3721\u001b[0m \u001b[38;5;66;03m# We want to prevent the loading of pylab to pollute the user's\u001b[39;00m\n\u001b[0;32m 3722\u001b[0m \u001b[38;5;66;03m# namespace as shown by the %who* magics, so we execute the activation\u001b[39;00m\n\u001b[0;32m 3723\u001b[0m \u001b[38;5;66;03m# code in an empty namespace, and we update *both* user_ns and\u001b[39;00m\n\u001b[0;32m 3724\u001b[0m \u001b[38;5;66;03m# user_ns_hidden with this information.\u001b[39;00m\n\u001b[0;32m 3725\u001b[0m ns \u001b[38;5;241m=\u001b[39m {}\n", + "File \u001b[1;32mc:\\Users\\EL221079\\Documents\\MICROPYTHON\\oaconvolve\\.venv\\Lib\\site-packages\\IPython\\core\\interactiveshell.py:3665\u001b[0m, in \u001b[0;36mInteractiveShell.enable_matplotlib\u001b[1;34m(self, gui)\u001b[0m\n\u001b[0;32m 3662\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib_inline\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mbackend_inline\u001b[39;00m\n\u001b[0;32m 3664\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mIPython\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcore\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m pylabtools \u001b[38;5;28;01mas\u001b[39;00m pt\n\u001b[1;32m-> 3665\u001b[0m gui, backend \u001b[38;5;241m=\u001b[39m \u001b[43mpt\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfind_gui_and_backend\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgui\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpylab_gui_select\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 3667\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m gui \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m 3668\u001b[0m \u001b[38;5;66;03m# If we have our first gui selection, store it\u001b[39;00m\n\u001b[0;32m 3669\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpylab_gui_select \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n", + "File \u001b[1;32mc:\\Users\\EL221079\\Documents\\MICROPYTHON\\oaconvolve\\.venv\\Lib\\site-packages\\IPython\\core\\pylabtools.py:338\u001b[0m, in \u001b[0;36mfind_gui_and_backend\u001b[1;34m(gui, gui_select)\u001b[0m\n\u001b[0;32m 321\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfind_gui_and_backend\u001b[39m(gui\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, gui_select\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m 322\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Given a gui string return the gui and mpl backend.\u001b[39;00m\n\u001b[0;32m 323\u001b[0m \n\u001b[0;32m 324\u001b[0m \u001b[38;5;124;03m Parameters\u001b[39;00m\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 335\u001b[0m \u001b[38;5;124;03m 'WXAgg','Qt4Agg','module://matplotlib_inline.backend_inline','agg').\u001b[39;00m\n\u001b[0;32m 336\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[1;32m--> 338\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\n\u001b[0;32m 340\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m _matplotlib_manages_backends():\n\u001b[0;32m 341\u001b[0m backend_registry \u001b[38;5;241m=\u001b[39m matplotlib\u001b[38;5;241m.\u001b[39mbackends\u001b[38;5;241m.\u001b[39mregistry\u001b[38;5;241m.\u001b[39mbackend_registry\n", + "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'matplotlib'" ] } ], @@ -22,6 +31,11 @@ "%pylab inline" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, { "cell_type": "markdown", "metadata": {}, @@ -31,7 +45,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-01-29T20:50:20.813162Z", @@ -49,7 +63,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2022-01-29T20:50:21.613220Z", @@ -131,7 +145,7 @@ }, { "cell_type": "code", - "execution_count": 57, + "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-05-07T07:35:35.126401Z", @@ -147,7 +161,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-05-19T19:11:18.145548Z", @@ -162,7 +176,7 @@ }, { "cell_type": "code", - "execution_count": 58, + "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-05-07T07:35:38.725924Z", @@ -225,9 +239,10 @@ "source": [ "# scipy.signal\n", "\n", - "This module defines the single function:\n", + "This module defines the functions:\n", "\n", - "1. [scipy.signal.sosfilt](#sosfilt)" + "1. [scipy.signal.sosfilt](#sosfilt)\n", + "2. [scipy.signal.oaconvolve](#oaconvolve)" ] }, { @@ -245,7 +260,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-06-19T20:24:10.529668Z", @@ -277,7 +292,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-06-19T20:27:39.508508Z", @@ -314,11 +329,60 @@ "print('y: ', y)\n", "print('\\n' + '='*40 + '\\nzf: ', zf)" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## oaconvolve\n", + "\n", + "``scipy``:\n", + "https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.oaconvolve.html\n", + "\n", + "Convolve two N-dimensional arrays using the overlap-add method.\n", + "\n", + "Convolve in1 and in2 using the overlap-add method. Similarly to numpy.convolve, \n", + "this method works in full mode and other modes can be obtained by slicing the result?\n", + "\n", + "This is generally much faster than linear convolve for large arrays (n > ~500), \n", + "and generally much faster than fftconvolve (not implemented yet in ulab) when one array is much larger \n", + "than the other (e.g. for a matched filter algorithm), but can be slower when only a few output values \n", + "are needed or when the arrays are very similar in shape, and can only output float arrays (int or object array inputs will be cast to float).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "ename": "ModuleNotFoundError", + "evalue": "No module named 'ulab'", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[2], line 3\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[38;5;66;03m# code to be run in micropython\u001b[39;00m\n\u001b[1;32m----> 3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mulab\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m numpy \u001b[38;5;28;01mas\u001b[39;00m np\n\u001b[0;32m 4\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mulab\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m scipy \u001b[38;5;28;01mas\u001b[39;00m spy\n\u001b[0;32m 6\u001b[0m x \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray((\u001b[38;5;241m1\u001b[39m,\u001b[38;5;241m2\u001b[39m,\u001b[38;5;241m3\u001b[39m))\n", + "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'ulab'" + ] + } + ], + "source": [ + " # code to be run in micropython\n", + " \n", + "from ulab import numpy as np\n", + "from ulab import scipy as spy\n", + "\n", + "x = np.array((1,2,3))\n", + "y = np.array((1,10,100,1000))\n", + "result = spy.signal.oaconvolve(x, y)\n", + "print('result: ', result)" + ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": ".venv", "language": "python", "name": "python3" }, @@ -332,7 +396,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.11.9" }, "toc": { "base_numbering": 1,