diff --git a/code/numpy/fft/fft.c b/code/numpy/fft/fft.c index e0893c2f..d4cab9e5 100644 --- a/code/numpy/fft/fft.c +++ b/code/numpy/fft/fft.c @@ -5,7 +5,7 @@ * * The MIT License (MIT) * - * Copyright (c) 2019-2021 Zoltán Vörös + * Copyright (c) 2019-2024 Zoltán Vörös * 2020 Scott Shawcroft for Adafruit Industries * 2020 Taku Fukada */ @@ -43,16 +43,16 @@ //| #if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE static mp_obj_t fft_fft(mp_obj_t arg) { - return fft_fft_ifft_spectrogram(arg, FFT_FFT); + return fft_fft_ifft(arg, FFT_FFT); } MP_DEFINE_CONST_FUN_OBJ_1(fft_fft_obj, fft_fft); #else static mp_obj_t fft_fft(size_t n_args, const mp_obj_t *args) { if(n_args == 2) { - return fft_fft_ifft_spectrogram(n_args, args[0], args[1], FFT_FFT); + return fft_fft_ifft(n_args, args[0], args[1], FFT_FFT); } else { - return fft_fft_ifft_spectrogram(n_args, args[0], mp_const_none, FFT_FFT); + return fft_fft_ifft(n_args, args[0], mp_const_none, FFT_FFT); } } @@ -71,7 +71,7 @@ MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(fft_fft_obj, 1, 2, fft_fft); #if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE static mp_obj_t fft_ifft(mp_obj_t arg) { - return fft_fft_ifft_spectrogram(arg, FFT_IFFT); + return fft_fft_ifft(arg, FFT_IFFT); } MP_DEFINE_CONST_FUN_OBJ_1(fft_ifft_obj, fft_ifft); @@ -79,9 +79,9 @@ MP_DEFINE_CONST_FUN_OBJ_1(fft_ifft_obj, fft_ifft); static mp_obj_t fft_ifft(size_t n_args, const mp_obj_t *args) { NOT_IMPLEMENTED_FOR_COMPLEX() if(n_args == 2) { - return fft_fft_ifft_spectrogram(n_args, args[0], args[1], FFT_IFFT); + return fft_fft_ifft(n_args, args[0], args[1], FFT_IFFT); } else { - return fft_fft_ifft_spectrogram(n_args, args[0], mp_const_none, FFT_IFFT); + return fft_fft_ifft(n_args, args[0], mp_const_none, FFT_IFFT); } } diff --git a/code/numpy/fft/fft_tools.c b/code/numpy/fft/fft_tools.c index 04c0b2f7..bc98b3d3 100644 --- a/code/numpy/fft/fft_tools.c +++ b/code/numpy/fft/fft_tools.c @@ -5,7 +5,7 @@ * * The MIT License (MIT) * - * Copyright (c) 2019-2021 Zoltán Vörös + * Copyright (c) 2019-2024 Zoltán Vörös */ #include @@ -45,7 +45,7 @@ imag[i] = data[2i+1] */ -void fft_kernel_complex(mp_float_t *data, size_t n, int isign) { +void fft_kernel(mp_float_t *data, size_t n, int isign) { size_t j, m, mmax, istep; mp_float_t tempr, tempi; mp_float_t wtemp, wr, wpr, wpi, wi, theta; @@ -94,9 +94,9 @@ void fft_kernel_complex(mp_float_t *data, size_t n, int isign) { /* * The following function is a helper interface to the python side. * It has been factored out from fft.c, so that the same argument parsing - * routine can be called from scipy.signal.spectrogram. + * routine can be called from utils.spectrogram. */ -mp_obj_t fft_fft_ifft_spectrogram(mp_obj_t data_in, uint8_t type) { +mp_obj_t fft_fft_ifft(mp_obj_t data_in, uint8_t type) { if(!mp_obj_is_type(data_in, &ulab_ndarray_type)) { mp_raise_NotImplementedError(MP_ERROR_TEXT("FFT is defined for ndarrays only")); } @@ -134,20 +134,10 @@ mp_obj_t fft_fft_ifft_spectrogram(mp_obj_t data_in, uint8_t type) { } data -= 2 * len; - if((type == FFT_FFT) || (type == FFT_SPECTROGRAM)) { - fft_kernel_complex(data, len, 1); - if(type == FFT_SPECTROGRAM) { - ndarray_obj_t *spectrum = ndarray_new_linear_array(len, NDARRAY_FLOAT); - mp_float_t *sarray = (mp_float_t *)spectrum->array; - for(size_t i = 0; i < len; i++) { - *sarray++ = MICROPY_FLOAT_C_FUN(sqrt)(data[0] * data[0] + data[1] * data[1]); - data += 2; - } - m_del(mp_float_t, data, 2 * len); - return MP_OBJ_FROM_PTR(spectrum); - } + if(type == FFT_FFT) { + fft_kernel(data, len, 1); } else { // inverse transform - fft_kernel_complex(data, len, -1); + fft_kernel(data, len, -1); // TODO: numpy accepts the norm keyword argument for(size_t i = 0; i < 2 * len; i++) { *data++ /= len; @@ -202,7 +192,7 @@ void fft_kernel(mp_float_t *real, mp_float_t *imag, size_t n, int isign) { } } -mp_obj_t fft_fft_ifft_spectrogram(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im, uint8_t type) { +mp_obj_t fft_fft_ifft(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im, uint8_t type) { if(!mp_obj_is_type(arg_re, &ulab_ndarray_type)) { mp_raise_NotImplementedError(MP_ERROR_TEXT("FFT is defined for ndarrays only")); } @@ -258,15 +248,8 @@ mp_obj_t fft_fft_ifft_spectrogram(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_i data_im -= len; } - if((type == FFT_FFT) || (type == FFT_SPECTROGRAM)) { + if(type == FFT_FFT) { fft_kernel(data_re, data_im, len, 1); - if(type == FFT_SPECTROGRAM) { - for(size_t i=0; i < len; i++) { - *data_re = MICROPY_FLOAT_C_FUN(sqrt)(*data_re * *data_re + *data_im * *data_im); - data_re++; - data_im++; - } - } } else { // inverse transform fft_kernel(data_re, data_im, len, -1); // TODO: numpy accepts the norm keyword argument @@ -275,13 +258,9 @@ mp_obj_t fft_fft_ifft_spectrogram(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_i *data_im++ /= len; } } - if(type == FFT_SPECTROGRAM) { - return MP_OBJ_FROM_PTR(out_re); - } else { - mp_obj_t tuple[2]; - tuple[0] = MP_OBJ_FROM_PTR(out_re); - tuple[1] = MP_OBJ_FROM_PTR(out_im); - return mp_obj_new_tuple(2, tuple); - } + mp_obj_t tuple[2]; + tuple[0] = MP_OBJ_FROM_PTR(out_re); + tuple[1] = MP_OBJ_FROM_PTR(out_im); + return mp_obj_new_tuple(2, tuple); } #endif /* ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE */ diff --git a/code/numpy/fft/fft_tools.h b/code/numpy/fft/fft_tools.h index 9444232f..aa598201 100644 --- a/code/numpy/fft/fft_tools.h +++ b/code/numpy/fft/fft_tools.h @@ -14,15 +14,14 @@ enum FFT_TYPE { FFT_FFT, FFT_IFFT, - FFT_SPECTROGRAM, }; #if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE void fft_kernel(mp_float_t *, size_t , int ); -mp_obj_t fft_fft_ifft_spectrogram(mp_obj_t , uint8_t ); +mp_obj_t fft_fft_ifft(mp_obj_t , uint8_t ); #else void fft_kernel(mp_float_t *, mp_float_t *, size_t , int ); -mp_obj_t fft_fft_ifft_spectrogram(size_t , mp_obj_t , mp_obj_t , uint8_t ); +mp_obj_t fft_fft_ifft(size_t , mp_obj_t , mp_obj_t , uint8_t ); #endif /* ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE */ #endif /* _FFT_TOOLS_ */ diff --git a/code/numpy/poly.c b/code/numpy/poly.c index ca5aa638..8b8e1435 100644 --- a/code/numpy/poly.c +++ b/code/numpy/poly.c @@ -121,7 +121,7 @@ mp_obj_t poly_polyfit(size_t n_args, const mp_obj_t *args) { ndarray_obj_t *beta = ndarray_new_linear_array(deg+1, NDARRAY_FLOAT); mp_float_t *betav = (mp_float_t *)beta->array; // x[0..(deg+1)] contains now the product X^T * y; we can get rid of y - m_del(float, y, leny); + m_del(mp_float_t, y, leny); // now, we calculate beta, i.e., we apply prod = (X^T * X)^(-1) on x = X^T * y; x is a column vector now for(uint8_t i=0; i < deg+1; i++) { diff --git a/code/ulab.c b/code/ulab.c index 948d621f..66582693 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.5.4 +#define ULAB_VERSION 6.5.5 #define xstr(s) str(s) #define str(s) #s diff --git a/code/ulab.h b/code/ulab.h index a384bfe2..36462e5b 100644 --- a/code/ulab.h +++ b/code/ulab.h @@ -440,7 +440,7 @@ // Note that in this case, the input also must be numpythonic, // i.e., the real an imaginary parts cannot be passed as two arguments #ifndef ULAB_FFT_IS_NUMPY_COMPATIBLE -#define ULAB_FFT_IS_NUMPY_COMPATIBLE (0) +#define ULAB_FFT_IS_NUMPY_COMPATIBLE (1) #endif #ifndef ULAB_FFT_HAS_FFT diff --git a/code/utils/utils.c b/code/utils/utils.c index 1d282d86..8477f5de 100644 --- a/code/utils/utils.c +++ b/code/utils/utils.c @@ -5,7 +5,7 @@ * * The MIT License (MIT) * - * Copyright (c) 2020-2021 Zoltán Vörös + * Copyright (c) 2020-2024 Zoltán Vörös */ #include @@ -16,6 +16,7 @@ #include "py/misc.h" #include "utils.h" +#include "../ulab_tools.h" #include "../numpy/fft/fft_tools.h" #if ULAB_HAS_UTILS_MODULE @@ -203,23 +204,180 @@ MP_DEFINE_CONST_FUN_OBJ_KW(utils_from_uint32_buffer_obj, 1, utils_from_uint32_bu //| ... //| -mp_obj_t utils_spectrogram(size_t n_args, const mp_obj_t *args) { - #if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE - return fft_fft_ifft_spectrogram(args[0], FFT_SPECTROGRAM); +mp_obj_t utils_spectrogram(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args) { + static const mp_arg_t allowed_args[] = { + { MP_QSTR_, MP_ARG_REQUIRED | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_NONE }} , + #if !ULAB_FFT_IS_NUMPY_COMPATIBLE + { MP_QSTR_, MP_ARG_OBJ, { .u_rom_obj = MP_ROM_NONE } }, + #endif + { MP_QSTR_scratchpad, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_out, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_NONE } }, + { MP_QSTR_log, MP_ARG_KW_ONLY | MP_ARG_OBJ, { .u_rom_obj = MP_ROM_FALSE } }, + }; + + 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_raise_NotImplementedError(MP_ERROR_TEXT("spectrogram is defined for ndarrays only")); + } + ndarray_obj_t *in = MP_OBJ_TO_PTR(args[0].u_obj); + + #if ULAB_MAX_DIMS > 1 + if(in->ndim != 1) { + mp_raise_TypeError(MP_ERROR_TEXT("spectrogram is implemented for 1D arrays only")); + } + #endif + + size_t len = in->len; + // Check if input is of length of power of 2 + if((len & (len-1)) != 0) { + mp_raise_ValueError(MP_ERROR_TEXT("input array length must be power of 2")); + } + + ndarray_obj_t *out = NULL; + + #if ULAB_FFT_IS_NUMPY_COMPATIBLE + mp_obj_t scratchpad_object = args[1].u_obj; + mp_obj_t out_object = args[2].u_obj; + mp_obj_t log_object = args[3].u_obj; #else + mp_obj_t scratchpad_object = args[2].u_obj; + mp_obj_t out_object = args[3].u_obj; + mp_obj_t log_object = args[4].u_obj; + #endif + + if(out_object != mp_const_none) { + if(!mp_obj_is_type(out_object, &ulab_ndarray_type)) { + mp_raise_TypeError(MP_ERROR_TEXT("out must be an ndarray")); + } + + out = MP_OBJ_TO_PTR(out_object); + if((out->dtype != NDARRAY_FLOAT) || (out->ndim != 1)){ + mp_raise_TypeError(MP_ERROR_TEXT("out array must be a 1D array of float type")); + } + if(len != out->len) { + mp_raise_ValueError(MP_ERROR_TEXT("input and out arrays must have same length")); + } + } else { + out = ndarray_new_linear_array(len, NDARRAY_FLOAT); + } + + ndarray_obj_t *scratchpad = NULL; + mp_float_t *tmp = NULL; + + if(scratchpad_object != mp_const_none) { + if(!mp_obj_is_type(scratchpad_object, &ulab_ndarray_type)) { + mp_raise_TypeError(MP_ERROR_TEXT("scratchpad must be an ndarray")); + } + + scratchpad = MP_OBJ_TO_PTR(scratchpad_object); + if(!ndarray_is_dense(scratchpad) || (scratchpad->ndim != 1) || (scratchpad->dtype != NDARRAY_FLOAT)) { + mp_raise_ValueError(MP_ERROR_TEXT("scratchpad must be a 1D dense float array")); + } + if(scratchpad->len != 2 * len) { + mp_raise_ValueError(MP_ERROR_TEXT("scratchpad must be twice as long as input")); + } + + tmp = (mp_float_t *)scratchpad->array; + } else { + tmp = m_new0(mp_float_t, 2 * len); + } + + uint8_t *array = (uint8_t *)in->array; + + #if ULAB_FFT_IS_NUMPY_COMPATIBLE // ULAB_SUPPORTS_COMPLEX is automatically true + if(in->dtype == NDARRAY_COMPLEX) { + uint8_t sz = 2 * sizeof(mp_float_t); + for(size_t i = 0; i < len; i++) { + memcpy(tmp, array, sz); + tmp += 2; + array += in->strides[ULAB_MAX_DIMS - 1]; + } + } else { + mp_float_t (*func)(void *) = ndarray_get_float_function(in->dtype); + for(size_t i = 0; i < len; i++) { + *tmp++ = func(array); // real part + *tmp++ = 0; // imaginary part, clear + array += in->strides[ULAB_MAX_DIMS - 1]; + } + } + + tmp -= 2 * len; + fft_kernel(tmp, len, 1); + #else // we might have two real input vectors + + ndarray_obj_t *in2 = NULL; + + if(n_args == 2) { + if(!mp_obj_is_type(args[1].u_obj, &ulab_ndarray_type)) { + mp_raise_NotImplementedError(MP_ERROR_TEXT("spectrogram is defined for ndarrays only")); + } + in2 = MP_OBJ_TO_PTR(args[1].u_obj); + + #if ULAB_MAX_DIMS > 1 + if(in2->ndim != 1) { + mp_raise_TypeError(MP_ERROR_TEXT("spectrogram is implemented for 1D arrays only")); + } + #endif + if(len != in2->len) { + mp_raise_TypeError(MP_ERROR_TEXT("input arrays are not compatible")); + } + } + + mp_float_t (*func)(void *) = ndarray_get_float_function(in->dtype); + + for(size_t i = 0; i < len; i++) { + *tmp++ = func(array); // real part; imageinary will be cleared later + array += in->strides[ULAB_MAX_DIMS - 1]; + } + if(n_args == 2) { - return fft_fft_ifft_spectrogram(n_args, args[0], args[1], FFT_SPECTROGRAM); + mp_float_t (*func2)(void *) = ndarray_get_float_function(in2->dtype); + array = (uint8_t *)in2->array; + for(size_t i = 0; i < len; i++) { + *tmp++ = func2(array); + array += in2->strides[ULAB_MAX_DIMS - 1]; + } + tmp -= len; } else { - return fft_fft_ifft_spectrogram(n_args, args[0], mp_const_none, FFT_SPECTROGRAM); + // if there is only one input argument, clear the imaginary part + memset(tmp, 0, len * sizeof(mp_float_t)); } - #endif + + tmp -= len; + + fft_kernel(tmp, tmp + len, len, 1); + #endif /* ULAB_FFT_IS_NUMPY_COMPATIBLE */ + + mp_float_t *spectrum = (mp_float_t *)out->array; + uint8_t spectrum_sz = out->strides[ULAB_MAX_DIMS - 1] / sizeof(mp_float_t); + + for(size_t i = 0; i < len; i++) { + #if ULAB_FFT_IS_NUMPY_COMPATIBLE + *spectrum = MICROPY_FLOAT_C_FUN(sqrt)(*tmp * *tmp + *(tmp + 1) * *(tmp + 1)); + tmp += 2; + #else + *spectrum = MICROPY_FLOAT_C_FUN(sqrt)(*tmp * *tmp + *(tmp + len) * *(tmp + len)); + tmp++; + #endif + if(log_object == mp_const_true) { + *spectrum = MICROPY_FLOAT_C_FUN(log)(*spectrum); + } + spectrum += spectrum_sz; + } + + if(scratchpad_object == mp_const_none) { + tmp -= len; + #if ULAB_FFT_IS_NUMPY_COMPATIBLE + tmp -= len; + #endif + m_del(mp_float_t, tmp, 2 * len); + } + return MP_OBJ_FROM_PTR(out); } -#if ULAB_SUPPORTS_COMPLEX & ULAB_FFT_IS_NUMPY_COMPATIBLE -MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(utils_spectrogram_obj, 1, 1, utils_spectrogram); -#else -MP_DEFINE_CONST_FUN_OBJ_VAR_BETWEEN(utils_spectrogram_obj, 1, 2, utils_spectrogram); -#endif +MP_DEFINE_CONST_FUN_OBJ_KW(utils_spectrogram_obj, 1, utils_spectrogram); #endif /* ULAB_UTILS_HAS_SPECTROGRAM */ diff --git a/docs/manual/source/conf.py b/docs/manual/source/conf.py index 4a2cea6e..6bbff154 100644 --- a/docs/manual/source/conf.py +++ b/docs/manual/source/conf.py @@ -27,7 +27,7 @@ author = 'Zoltán Vörös' # The full version, including alpha/beta/rc tags -release = '6.5.0' +release = '6.5.5' # -- General configuration --------------------------------------------------- diff --git a/docs/manual/source/ulab-intro.rst b/docs/manual/source/ulab-intro.rst index 81019e33..42e5b260 100644 --- a/docs/manual/source/ulab-intro.rst +++ b/docs/manual/source/ulab-intro.rst @@ -8,9 +8,10 @@ Enter ulab ``ulab`` is a ``numpy``-like module for ``micropython`` and its derivatives, meant to simplify and speed up common mathematical operations on arrays. ``ulab`` implements a small subset of ``numpy`` -and ``scipy``. The functions were chosen such that they might be useful -in the context of a microcontroller. However, the project is a living -one, and suggestions for new features are always welcome. +and ``scipy``, as well as a number of functions manipulating byte +arrays. The functions were chosen such that they might be useful in the +context of a microcontroller. However, the project is a living one, and +suggestions for new features are always welcome. This document discusses how you can use the library, starting from building your own firmware, through questions like what affects the @@ -265,9 +266,9 @@ functions that are part of ``numpy``, you have to import ``numpy`` as p = np.array([1, 2, 3]) np.polyval(p, x) -There are a couple of exceptions to this rule, namely ``fft``, and -``linalg``, which are sub-modules even in ``numpy``, thus you have to -write them out as +There are a couple of exceptions to this rule, namely ``fft``, +``linalg``, and ``random``, which are sub-modules even in ``numpy``, +thus you have to write them out as .. code:: python diff --git a/docs/manual/source/ulab-ndarray.rst b/docs/manual/source/ulab-ndarray.rst index ef54a242..04036338 100644 --- a/docs/manual/source/ulab-ndarray.rst +++ b/docs/manual/source/ulab-ndarray.rst @@ -2330,12 +2330,12 @@ future version of ``ulab``. a = np.array(range(9), dtype=np.float) print("a:\t", a) - print("a < 5:\t", a[a < 5]) + print("a[a < 5]:\t", a[a < 5]) .. parsed-literal:: a: array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0], dtype=float) - a < 5: array([0.0, 1.0, 2.0, 3.0, 4.0], dtype=float) + a[a < 5]: array([0.0, 1.0, 2.0, 3.0, 4.0], dtype=float) diff --git a/docs/manual/source/ulab-utils.rst b/docs/manual/source/ulab-utils.rst index 82476bad..4305dfb9 100644 --- a/docs/manual/source/ulab-utils.rst +++ b/docs/manual/source/ulab-utils.rst @@ -52,7 +52,9 @@ Here is an example without keyword arguments a = bytearray([1, 1, 0, 0, 0, 0, 0, 255]) print('a: ', a) print() - print('unsigned integers: ', utils.from_uint32_buffer(a)) + print('unsigned integers: ', utils.from_uint32_buffe + print('original vector:\n', y) + print('\nspectrum:\n', a)r(a)) b = bytearray([1, 1, 0, 0, 0, 0, 0, 255]) print('\nb: ', b) @@ -144,9 +146,53 @@ In addition to the Fourier transform and its inverse, ``ulab`` also sports a function called ``spectrogram``, which returns the absolute value of the Fourier transform, also known as the power spectrum. This could be used to find the dominant spectral component in a time series. -The arguments are treated in the same way as in ``fft``, and ``ifft``. -This means that, if the firmware was compiled with complex support, the -input can also be a complex array. +The positional arguments are treated in the same way as in ``fft``, and +``ifft``. This means that, if the firmware was compiled with complex +support and ``ULAB_FFT_IS_NUMPY_COMPATIBLE`` is defined to be 1 in +``ulab.h``, the input can also be a complex array. + +And easy way to find out if the FFT is ``numpy``-compatible is to check +the number of values ``fft.fft`` returns, when called with a single real +argument of length other than 2: + +.. code:: + + # code to be run in micropython + + from ulab import numpy as np + + if len(np.fft.fft(np.zeros(4))) == 2: + print('FFT is NOT numpy compatible (real and imaginary parts are treated separately)') + else: + print('FFT is numpy compatible (complex inputs/outputs)') + +.. parsed-literal:: + + FFT is numpy compatible (complex inputs/outputs) + + + + +Depending on the ``numpy``-compatibility of the FFT, the ``spectrogram`` +function takes one or two positional arguments, and three keyword +arguments. If the FFT is ``numpy`` compatible, one positional argument +is allowed, and it is a 1D real or complex ``ndarray``. If the FFT is +not ``numpy``-compatible, if a single argument is supplied, it will be +treated as the real part of the input, and if two positional arguments +are supplied, they are treated as the real and imaginary parts of the +signal. + +The keyword arguments are as follows: + +1. ``scratchpad = None``: must be a 1D, dense, floating point array, + twice as long as the input array; the ``scratchpad`` will be used as + a temporary internal buffer to perform the Fourier transform; the + ``scratchpad`` can repeatedly be re-used. +2. ``out = None``: must be a 1D, not necessarily dense, floating point + array that will store the results +3. ``log = False``: must be either ``True``, or ``False``; if ``True``, + the ``spectrogram`` returns the logarithm of the absolute values of + the Fourier transform. .. code:: @@ -169,17 +215,24 @@ input can also be a complex array. array([0.0, 0.009775015390171337, 0.01954909674625918, ..., -0.5275140569487312, -0.5357931822978732, -0.5440211108893697], dtype=float64) spectrum: - array([187.8635087634579, 315.3112063607119, 347.8814873399374, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float64) + array([187.8635087634578, 315.3112063607119, 347.8814873399375, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float64) As such, ``spectrogram`` is really just a shorthand for -``np.sqrt(a*a + b*b)``, however, it saves significant amounts of RAM: -the expression ``a*a + b*b`` has to allocate memory for ``a*a``, -``b*b``, and finally, their sum. In contrast, ``spectrogram`` calculates -the spectrum internally, and stores it in the memory segment that was -reserved for the real part of the Fourier transform. +``np.abs(np.fft.fft(signal))``, if the FFT is ``numpy``-compatible, or +``np.sqrt(a*a + b*b)`` if the FFT returns the real (``a``) and imaginary +(``b``) parts separately. However, ``spectrogram`` saves significant +amounts of RAM: the expression ``a*a + b*b`` has to allocate memory for +``a*a``, ``b*b``, and finally, their sum. Similarly, ``np.abs`` returns +a new array. This issue is compounded even more, if ``np.log()`` is used +on the absolute value. + +In contrast, ``spectrogram`` handles all calculations in the same +internal arrays, and allows one to re-use previously reserved RAM. This +can be especially useful in cases, when ``spectogram`` is called +repeatedly, as in the snippet below. .. code:: @@ -188,25 +241,34 @@ reserved for the real part of the Fourier transform. from ulab import numpy as np from ulab import utils as utils - x = np.linspace(0, 10, num=1024) - y = np.sin(x) + n = 1024 + t = np.linspace(0, 2 * np.pi, num=1024) + scratchpad = np.zeros(2 * n) - a, b = np.fft.fft(y) + for _ in range(10): + signal = np.sin(t) + utils.spectrogram(signal, out=signal, scratchpad=scratchpad, log=True) - print('\nspectrum calculated the hard way:\n', np.sqrt(a*a + b*b)) + print('signal: ', signal) - a = utils.spectrogram(y) + for _ in range(10): + signal = np.sin(t) + out = np.log(utils.spectrogram(signal)) - print('\nspectrum calculated the lazy way:\n', a) + print('out: ', out) .. parsed-literal:: - - spectrum calculated the hard way: - array([187.8635087634579, 315.3112063607119, 347.8814873399374, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float64) - - spectrum calculated the lazy way: - array([187.8635087634579, 315.3112063607119, 347.8814873399374, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float64) + signal: array([-27.38260169844543, 6.237834411021073, -0.4038327279002965, ..., -0.9795967096969854, -0.4038327279002969, 6.237834411021073], dtype=float64) + out: array([-27.38260169844543, 6.237834411021073, -0.4038327279002965, ..., -0.9795967096969854, -0.4038327279002969, 6.237834411021073], dtype=float64) + +Note that ``scratchpad`` is reserved only once, and then is re-used in +the first loop. By assigning ``signal`` to the output, we save +additional RAM. This approach avoids the usual problem of memory +fragmentation, which would happen in the second loop, where both +``spectrogram``, and ``np.log`` must reserve RAM in each iteration. + + diff --git a/docs/ulab-change-log.md b/docs/ulab-change-log.md index 90b86307..4600c7c4 100644 --- a/docs/ulab-change-log.md +++ b/docs/ulab-change-log.md @@ -1,5 +1,11 @@ Sat, 14 Sep 2024 +version 6.5.5 + + add scratchpad, out, log keyword arguments to spectrum + +Sat, 14 Sep 2024 + version 6.5.4 fix roll, when shift is 0 diff --git a/docs/ulab-convert.ipynb b/docs/ulab-convert.ipynb index bd587912..ef5e24f5 100644 --- a/docs/ulab-convert.ipynb +++ b/docs/ulab-convert.ipynb @@ -57,11 +57,11 @@ "# -- Project information -----------------------------------------------------\n", "\n", "project = 'The ulab book'\n", - "copyright = '2019-2022, Zoltán Vörös and contributors'\n", + "copyright = '2019-2024, Zoltán Vörös and contributors'\n", "author = 'Zoltán Vörös'\n", "\n", "# The full version, including alpha/beta/rc tags\n", - "release = '5.1.0'\n", + "release = '6.5.5'\n", "\n", "\n", "# -- General configuration ---------------------------------------------------\n", @@ -190,6 +190,7 @@ " numpy-universal\n", " numpy-fft\n", " numpy-linalg\n", + " numpy-random\n", " scipy-linalg\n", " scipy-optimize\n", " scipy-signal\n", @@ -215,7 +216,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": { "ExecuteTime": { "end_time": "2022-02-09T06:27:21.647179Z", @@ -256,14 +257,49 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2022-02-09T06:27:42.024028Z", "start_time": "2022-02-09T06:27:36.109093Z" } }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n", + "/home/v923z/anaconda3/lib/python3.11/site-packages/nbconvert/exporters/exporter.py:349: MissingIDFieldWarning: Code cell is missing an id field, this will become a hard error in future nbformat versions. You may want to use `normalize()` on your notebooks before validations (available since nbformat 5.1.4). Previous versions of nbformat are fixing this issue transparently, and will stop doing so in the future.\n", + " _, nbc = validator.normalize(nbc)\n" + ] + } + ], "source": [ "files = ['ulab-intro',\n", " 'ulab-ndarray',\n", @@ -271,6 +307,7 @@ " 'numpy-universal',\n", " 'numpy-fft',\n", " 'numpy-linalg',\n", + " 'numpy-random',\n", " 'scipy-linalg',\n", " 'scipy-optimize',\n", " 'scipy-signal',\n", @@ -449,7 +486,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.11.7" }, "toc": { "base_numbering": 1, diff --git a/docs/ulab-intro.ipynb b/docs/ulab-intro.ipynb index 67d6b608..30266418 100644 --- a/docs/ulab-intro.ipynb +++ b/docs/ulab-intro.ipynb @@ -10,13 +10,6 @@ } }, "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Matplotlib is building the font cache; this may take a moment.\n" - ] - }, { "name": "stdout", "output_type": "stream", @@ -38,7 +31,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2022-01-07T18:13:14.590799Z", @@ -239,7 +232,7 @@ "source": [ "## Enter ulab\n", "\n", - "`ulab` is a `numpy`-like module for `micropython` and its derivatives, meant to simplify and speed up common mathematical operations on arrays. `ulab` implements a small subset of `numpy` and `scipy`. The functions were chosen such that they might be useful in the context of a microcontroller. However, the project is a living one, and suggestions for new features are always welcome. \n", + "`ulab` is a `numpy`-like module for `micropython` and its derivatives, meant to simplify and speed up common mathematical operations on arrays. `ulab` implements a small subset of `numpy` and `scipy`, as well as a number of functions manipulating byte arrays. The functions were chosen such that they might be useful in the context of a microcontroller. However, the project is a living one, and suggestions for new features are always welcome. \n", "\n", "This document discusses how you can use the library, starting from building your own firmware, through questions like what affects the firmware size, what are the trade-offs, and what are the most important differences to `numpy` and `scipy`, respectively. The document is organised as follows:\n", "\n", @@ -404,7 +397,7 @@ "np.polyval(p, x)\n", "```\n", "\n", - "There are a couple of exceptions to this rule, namely `fft`, and `linalg`, which are sub-modules even in `numpy`, thus you have to write them out as \n", + "There are a couple of exceptions to this rule, namely `fft`, `linalg`, and `random`, which are sub-modules even in `numpy`, thus you have to write them out as \n", "\n", "```python\n", "from ulab import numpy as np\n", @@ -842,7 +835,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.9.13" }, "toc": { "base_numbering": 1, diff --git a/docs/ulab-utils.ipynb b/docs/ulab-utils.ipynb index 3c5e5c66..8dce72c0 100644 --- a/docs/ulab-utils.ipynb +++ b/docs/ulab-utils.ipynb @@ -31,7 +31,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2022-01-29T16:53:11.972661Z", @@ -49,7 +49,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2022-01-29T16:59:24.652277Z", @@ -77,7 +77,7 @@ " if args.unix: # tests the code on the unix port. Note that this works on unix only\n", " with open('/dev/shm/micropython.py', 'w') as fout:\n", " fout.write(cell)\n", - " proc = subprocess.Popen([\"../micropython/ports/unix/micropython-2\", \"/dev/shm/micropython.py\"], \n", + " proc = subprocess.Popen([\"../micropython/ports/unix/build-2/micropython-2\", \"/dev/shm/micropython.py\"], \n", " stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n", " print(proc.stdout.read().decode(\"utf-8\"))\n", " print(proc.stderr.read().decode(\"utf-8\"))\n", @@ -291,7 +291,9 @@ "a = bytearray([1, 1, 0, 0, 0, 0, 0, 255])\n", "print('a: ', a)\n", "print()\n", - "print('unsigned integers: ', utils.from_uint32_buffer(a))\n", + "print('unsigned integers: ', utils.from_uint32_buffe\n", + "print('original vector:\\n', y)\n", + "print('\\nspectrum:\\n', a)r(a))\n", "\n", "b = bytearray([1, 1, 0, 0, 0, 0, 0, 255])\n", "print('\\nb: ', b)\n", @@ -398,12 +400,53 @@ "source": [ "## spectrogram\n", "\n", - "In addition to the Fourier transform and its inverse, `ulab` also sports a function called `spectrogram`, which returns the absolute value of the Fourier transform, also known as the power spectrum. This could be used to find the dominant spectral component in a time series. The arguments are treated in the same way as in `fft`, and `ifft`. This means that, if the firmware was compiled with complex support, the input can also be a complex array." + "In addition to the Fourier transform and its inverse, `ulab` also sports a function called `spectrogram`, which returns the absolute value of the Fourier transform, also known as the power spectrum. This could be used to find the dominant spectral component in a time series. The positional arguments are treated in the same way as in `fft`, and `ifft`. This means that, if the firmware was compiled with complex support and `ULAB_FFT_IS_NUMPY_COMPATIBLE` is defined to be 1 in `ulab.h`, the input can also be a complex array. \n", + "\n", + "And easy way to find out if the FFT is `numpy`-compatible is to check the number of values `fft.fft` returns, when called with a single real argument of length other than 2: " ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "FFT is numpy compatible (complex inputs/outputs)\n", + "\n", + "\n" + ] + } + ], + "source": [ + "%%micropython -unix 1\n", + "\n", + "from ulab import numpy as np\n", + "\n", + "if len(np.fft.fft(np.zeros(4))) == 2:\n", + " print('FFT is NOT numpy compatible (real and imaginary parts are treated separately)')\n", + "else:\n", + " print('FFT is numpy compatible (complex inputs/outputs)')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Depending on the `numpy`-compatibility of the FFT, the `spectrogram` function takes one or two positional arguments, and three keyword arguments. If the FFT is `numpy` compatible, one positional argument is allowed, and it is a 1D real or complex `ndarray`. If the FFT is not `numpy`-compatible, if a single argument is supplied, it will be treated as the real part of the input, and if two positional arguments are supplied, they are treated as the real and imaginary parts of the signal.\n", + "\n", + "The keyword arguments are as follows:\n", + "\n", + "1. `scratchpad = None`: must be a 1D, dense, floating point array, twice as long as the input array; the `scratchpad` will be used as a temporary internal buffer to perform the Fourier transform; the `scratchpad` can repeatedly be re-used.\n", + "1. `out = None`: must be a 1D, not necessarily dense, floating point array that will store the results\n", + "1. `log = False`: must be either `True`, or `False`; if `True`, the `spectrogram` returns the logarithm of the absolute values of the Fourier transform." + ] + }, + { + "cell_type": "code", + "execution_count": 15, "metadata": { "ExecuteTime": { "end_time": "2022-01-29T16:59:56.400603Z", @@ -419,7 +462,7 @@ " array([0.0, 0.009775015390171337, 0.01954909674625918, ..., -0.5275140569487312, -0.5357931822978732, -0.5440211108893697], dtype=float64)\n", "\n", "spectrum:\n", - " array([187.8635087634579, 315.3112063607119, 347.8814873399374, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float64)\n", + " array([187.8635087634578, 315.3112063607119, 347.8814873399375, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float64)\n", "\n", "\n" ] @@ -444,12 +487,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "As such, `spectrogram` is really just a shorthand for `np.sqrt(a*a + b*b)`, however, it saves significant amounts of RAM: the expression `a*a + b*b` has to allocate memory for `a*a`, `b*b`, and finally, their sum. In contrast, `spectrogram` calculates the spectrum internally, and stores it in the memory segment that was reserved for the real part of the Fourier transform." + "As such, `spectrogram` is really just a shorthand for `np.abs(np.fft.fft(signal))`, if the FFT is `numpy`-compatible, or `np.sqrt(a*a + b*b)` if the FFT returns the real (`a`) and imaginary (`b`) parts separately. However, `spectrogram` saves significant amounts of RAM: the expression `a*a + b*b` has to allocate memory for `a*a`, `b*b`, and finally, their sum. Similarly, `np.abs` returns a new array. This issue is compounded even more, if `np.log()` is used on the absolute value. \n", + "\n", + "In contrast, `spectrogram` handles all calculations in the same internal arrays, and allows one to re-use previously reserved RAM. This can be especially useful in cases, when `spectogram` is called repeatedly, as in the snippet below." ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 34, "metadata": { "ExecuteTime": { "end_time": "2022-01-29T16:59:48.485610Z", @@ -461,12 +506,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "\n", - "spectrum calculated the hard way:\n", - " array([187.8635087634579, 315.3112063607119, 347.8814873399374, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float64)\n", - "\n", - "spectrum calculated the lazy way:\n", - " array([187.8635087634579, 315.3112063607119, 347.8814873399374, ..., 84.45888934298905, 347.8814873399374, 315.3112063607118], dtype=float64)\n", + "signal: array([-27.38260169844543, 6.237834411021073, -0.4038327279002965, ..., -0.9795967096969854, -0.4038327279002969, 6.237834411021073], dtype=float64)\n", + "out: array([-27.38260169844543, 6.237834411021073, -0.4038327279002965, ..., -0.9795967096969854, -0.4038327279002969, 6.237834411021073], dtype=float64)\n", "\n", "\n" ] @@ -478,17 +519,34 @@ "from ulab import numpy as np\n", "from ulab import utils as utils\n", "\n", - "x = np.linspace(0, 10, num=1024)\n", - "y = np.sin(x)\n", + "n = 1024\n", + "t = np.linspace(0, 2 * np.pi, num=1024)\n", + "scratchpad = np.zeros(2 * n)\n", "\n", - "a, b = np.fft.fft(y)\n", + "for _ in range(10):\n", + " signal = np.sin(t)\n", + " utils.spectrogram(signal, out=signal, scratchpad=scratchpad, log=True)\n", "\n", - "print('\\nspectrum calculated the hard way:\\n', np.sqrt(a*a + b*b))\n", + "print('signal: ', signal)\n", "\n", - "a = utils.spectrogram(y)\n", + "for _ in range(10):\n", + " signal = np.sin(t)\n", + " out = np.log(utils.spectrogram(signal))\n", "\n", - "print('\\nspectrum calculated the lazy way:\\n', a)" + "print('out: ', out)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that `scratchpad` is reserved only once, and then is re-used in the first loop. By assigning `signal` to the output, we save additional RAM. This approach avoids the usual problem of memory fragmentation, which would happen in the second loop, where both `spectrogram`, and `np.log` must reserve RAM in each iteration." ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] } ], "metadata": { @@ -507,7 +565,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.5" + "version": "3.11.7" }, "toc": { "base_numbering": 1, diff --git a/tests/1d/numpy/fft.py b/tests/1d/numpy/fft.py index 1a1dee75..6b79f74c 100644 --- a/tests/1d/numpy/fft.py +++ b/tests/1d/numpy/fft.py @@ -10,8 +10,12 @@ y = np.sin(x) if use_ulab: - a, b = np.fft.fft(y) - c, d = np.fft.ifft(a, b) + if 'real' in dir(np): + a = np.fft.fft(y) + c = np.real(np.fft.ifft(a)) + else: + a, b = np.fft.fft(y) + c, d = np.fft.ifft(a, b) # c should be equal to y cmp_result = [] for p,q in zip(list(y), list(c)): @@ -19,8 +23,12 @@ print(cmp_result) z = np.zeros(len(x)) - a, b = np.fft.fft(y, z) - c, d = np.fft.ifft(a, b) + if 'real' in dir(np): + a = np.fft.fft(y) + c = np.real(np.fft.ifft(a)) + else: + a, b = np.fft.fft(y, z) + c, d = np.fft.ifft(a, b) # c should be equal to y cmp_result = [] for p,q in zip(list(y), list(c)):