Compare commits
3 commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
30c7f4f49e | ||
|
|
feb9ad9ac7 | ||
|
|
decc6013c2 |
10 changed files with 324 additions and 100 deletions
|
|
@ -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);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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 <math.h>
|
||||
|
|
@ -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 */
|
||||
|
|
|
|||
|
|
@ -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_ */
|
||||
|
|
|
|||
|
|
@ -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++) {
|
||||
|
|
|
|||
|
|
@ -33,7 +33,7 @@
|
|||
#include "user/user.h"
|
||||
#include "utils/utils.h"
|
||||
|
||||
#define ULAB_VERSION 6.5.0
|
||||
#define ULAB_VERSION 6.5.1
|
||||
#define xstr(s) str(s)
|
||||
#define str(s) #s
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 <math.h>
|
||||
|
|
@ -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);
|
||||
#else
|
||||
if(n_args == 2) {
|
||||
return fft_fft_ifft_spectrogram(n_args, args[0], args[1], FFT_SPECTROGRAM);
|
||||
} else {
|
||||
return fft_fft_ifft_spectrogram(n_args, args[0], mp_const_none, 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) {
|
||||
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 {
|
||||
// if there is only one input argument, clear the imaginary part
|
||||
memset(tmp, 0, len * sizeof(mp_float_t));
|
||||
}
|
||||
|
||||
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 */
|
||||
|
||||
|
|
|
|||
|
|
@ -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.0'\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": 5,
|
||||
"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.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\n",
|
||||
"/home/v923z/anaconda3/lib/python3.9/site-packages/nbconvert/exporters/exporter.py:305: 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",
|
||||
" nbformat.validate(nbc, relax_add_props=True)\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.9.13"
|
||||
},
|
||||
"toc": {
|
||||
"base_numbering": 1,
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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\n",
|
||||
"1. `out = None`: must be a 1D, not necessarily dense, floating point array that will store the results\n",
|
||||
"1. `log = False`: must 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.9.13"
|
||||
},
|
||||
"toc": {
|
||||
"base_numbering": 1,
|
||||
|
|
|
|||
Loading…
Reference in a new issue