diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 8900bc71..c0477efb 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -16,6 +16,7 @@ on:
 
 jobs:
   micropython:
+    continue-on-error: true
     strategy:
         matrix:
             os:
@@ -28,10 +29,10 @@ jobs:
       env:
         GITHUB_CONTEXT: ${{ toJson(github) }}
       run: echo "$GITHUB_CONTEXT"
-    - name: Set up Python 3.10
-      uses: actions/setup-python@v1
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v5
       with:
-        python-version: "3.10"
+        python-version: "3.12"
 
     - name: Install requirements
       run: |
@@ -44,10 +45,10 @@ jobs:
         gcc --version
         python3 --version
     - name: Checkout ulab
-      uses: actions/checkout@v1
+      uses: actions/checkout@v4
 
     - name: Checkout micropython repo
-      uses: actions/checkout@v2
+      uses: actions/checkout@v4
       with:
         repository: micropython/micropython
         path: micropython
@@ -56,6 +57,7 @@ jobs:
       run: ./build.sh ${{ matrix.dims }}
 
   circuitpython:
+    continue-on-error: true
     strategy:
         matrix:
             os:
@@ -68,10 +70,10 @@ jobs:
       env:
         GITHUB_CONTEXT: ${{ toJson(github) }}
       run: echo "$GITHUB_CONTEXT"
-    - name: Set up Python 3.10
-      uses: actions/setup-python@v1
+    - name: Set up Python 3.12
+      uses: actions/setup-python@v5
       with:
-        python-version: "3.10"
+        python-version: "3.12"
 
     - name: Versions
       run: |
@@ -79,7 +81,7 @@ jobs:
         python3 --version
 
     - name: Checkout ulab
-      uses: actions/checkout@v1
+      uses: actions/checkout@v4
 
     - name: Install requirements
       run: |
diff --git a/README.md b/README.md
index e7feeca6..5c91d671 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # ulab
 
-[![Documentation Status](https://readthedocs.org/projects/micropython-ulab-robert/badge/?version=latest)](https://micropython-ulab-robert.readthedocs.io/en/latest/?badge=latest)
+[![Documentation Status](https://readthedocs.org/projects/micropython-ulab/badge/?version=latest)](https://micropython-ulab.readthedocs.io/en/latest/index.html)
 
 `ulab` is a `numpy`-like array manipulation library for [micropython](http://micropython.org/) and [CircuitPython](https://circuitpython.org/).
 The module is written in C, defines compact containers (`ndarray`s) for numerical data of one to four
@@ -112,7 +112,6 @@ of the user manual.
 1. `MaixPy` https://github.com/sipeed/MaixPy
 1. `OpenMV` https://github.com/openmv/openmv
 1. `pimoroni-pico` https://github.com/pimoroni/pimoroni-pico
-3. `pycom` https://pycom.io/
 
 ## Compiling
 
diff --git a/build-cp.sh b/build-cp.sh
index ee6a95ef..9081a618 100755
--- a/build-cp.sh
+++ b/build-cp.sh
@@ -41,8 +41,7 @@ HERE="$(dirname -- "$(readlinkf_posix -- "${0}")" )"
 rm -rf circuitpython/extmod/ulab; ln -s "$HERE" circuitpython/extmod/ulab
 dims=${1-2}
 make -C circuitpython/mpy-cross -j$NPROC
-sed -e '/MICROPY_PY_UHASHLIB/s/1/0/' < circuitpython/ports/unix/mpconfigport.h > circuitpython/ports/unix/mpconfigport_ulab.h
-make -k -C circuitpython/ports/unix -j$NPROC DEBUG=1 MICROPY_PY_FFI=0 MICROPY_PY_BTREE=0 MICROPY_SSL_AXTLS=0 MICROPY_PY_USSL=0 CFLAGS_EXTRA="-DMP_CONFIGFILE=\"<mpconfigport_ulab.h>\" -Wno-tautological-constant-out-of-range-compare -Wno-unknown-pragmas -DULAB_MAX_DIMS=$dims" BUILD=build-$dims PROG=micropython-$dims
+make -k -C circuitpython/ports/unix -j$NPROC DEBUG=1 MICROPY_PY_FFI=0 MICROPY_PY_BTREE=0 MICROPY_SSL_AXTLS=0 MICROPY_PY_USSL=0 CFLAGS_EXTRA="-Wno-tautological-constant-out-of-range-compare -Wno-unknown-pragmas -DULAB_MAX_DIMS=$dims" BUILD=build-$dims PROG=micropython-$dims
 
 # bash test-common.sh "${dims}" "circuitpython/ports/unix/micropython-$dims"
 
diff --git a/code/ndarray.c b/code/ndarray.c
index 5e662851..26ced6fc 100644
--- a/code/ndarray.c
+++ b/code/ndarray.c
@@ -6,7 +6,7 @@
  *
  * The MIT License (MIT)
  *
- * Copyright (c) 2019-2022 Zoltán Vörös
+ * Copyright (c) 2019-2024 Zoltán Vörös
  *               2020 Jeff Epler for Adafruit Industries
  *               2020 Taku Fukada
 */
@@ -509,8 +509,9 @@ static size_t multiply_size(size_t a, size_t b) {
     return result;
 }
 
-ndarray_obj_t *ndarray_new_ndarray(uint8_t ndim, size_t *shape, int32_t *strides, uint8_t dtype) {
+ndarray_obj_t *ndarray_new_ndarray(uint8_t ndim, size_t *shape, int32_t *strides, uint8_t dtype, uint8_t *buffer) {
     // Creates the base ndarray with shape, and initialises the values to straight 0s
+    // optionally, values can be supplied via the last argument
     ndarray_obj_t *ndarray = m_new_obj(ndarray_obj_t);
     ndarray->base.type = &ulab_ndarray_type;
     ndarray->dtype = dtype == NDARRAY_BOOL ? NDARRAY_UINT8 : dtype;
@@ -536,9 +537,13 @@ ndarray_obj_t *ndarray_new_ndarray(uint8_t ndim, size_t *shape, int32_t *strides
 
     // if the length is 0, still allocate a single item, so that contractions can be handled
     size_t len = multiply_size(ndarray->itemsize, MAX(1, ndarray->len));
-    uint8_t *array = m_new0(byte, len);
-    // this should set all elements to 0, irrespective of the of the dtype (all bits are zero)
-    // we could, perhaps, leave this step out, and initialise the array only, when needed
+    uint8_t *array;
+    array = buffer;
+    if(array == NULL) {
+        // this should set all elements to 0, irrespective of the of the dtype (all bits are zero)
+        // we could, perhaps, leave this step out, and initialise the array only, when needed
+        array = m_new0(byte, len);
+    }
     ndarray->array = array;
     ndarray->origin = array;
     return ndarray;
@@ -547,12 +552,12 @@ ndarray_obj_t *ndarray_new_ndarray(uint8_t ndim, size_t *shape, int32_t *strides
 ndarray_obj_t *ndarray_new_dense_ndarray(uint8_t ndim, size_t *shape, uint8_t dtype) {
     // creates a dense array, i.e., one, where the strides are derived directly from the shapes
     // the function should work in the general n-dimensional case
-    int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
-    strides[ULAB_MAX_DIMS-1] = (int32_t)ulab_binary_get_size(dtype);
-    for(size_t i=ULAB_MAX_DIMS; i > 1; i--) {
-        strides[i-2] = strides[i-1] * MAX(1, shape[i-1]);
-    }
-    return ndarray_new_ndarray(ndim, shape, strides, dtype);
+    // int32_t *strides = m_new(int32_t, ULAB_MAX_DIMS);
+    // strides[ULAB_MAX_DIMS - 1] = (int32_t)ulab_binary_get_size(dtype);
+    // for(size_t i = ULAB_MAX_DIMS; i > 1; i--) {
+    //     strides[i-2] = strides[i-1] * MAX(1, shape[i-1]);
+    // }
+    return ndarray_new_ndarray(ndim, shape, NULL, dtype, NULL);
 }
 
 ndarray_obj_t *ndarray_new_ndarray_from_tuple(mp_obj_tuple_t *_shape, uint8_t dtype) {
@@ -650,7 +655,7 @@ ndarray_obj_t *ndarray_copy_view(ndarray_obj_t *source) {
     if(source->boolean) {
         dtype = NDARRAY_BOOL;
     }
-    ndarray_obj_t *ndarray = ndarray_new_ndarray(source->ndim, source->shape, strides, dtype);
+    ndarray_obj_t *ndarray = ndarray_new_ndarray(source->ndim, source->shape, strides, dtype, NULL);
     ndarray_copy_array(source, ndarray, 0);
     return ndarray;
 }
@@ -1884,7 +1889,7 @@ mp_obj_t ndarray_unary_op(mp_unary_op_t op, mp_obj_t self_in) {
             #if ULAB_SUPPORTS_COMPLEX
             if(self->dtype == NDARRAY_COMPLEX) {
                 int32_t *strides = strides_from_shape(self->shape, NDARRAY_FLOAT);
-                ndarray_obj_t *target = ndarray_new_ndarray(self->ndim, self->shape, strides, NDARRAY_FLOAT);
+                ndarray_obj_t *target = ndarray_new_ndarray(self->ndim, self->shape, strides, NDARRAY_FLOAT, NULL);
                 ndarray = MP_OBJ_TO_PTR(carray_abs(self, target));
             } else {
             #endif
diff --git a/code/ndarray.h b/code/ndarray.h
index 2ce84bed..3e82b385 100644
--- a/code/ndarray.h
+++ b/code/ndarray.h
@@ -188,7 +188,7 @@ int32_t *ndarray_contract_strides(ndarray_obj_t *, uint8_t );
 ndarray_obj_t *ndarray_from_iterable(mp_obj_t , uint8_t );
 ndarray_obj_t *ndarray_new_dense_ndarray(uint8_t , size_t *, uint8_t );
 ndarray_obj_t *ndarray_new_ndarray_from_tuple(mp_obj_tuple_t *, uint8_t );
-ndarray_obj_t *ndarray_new_ndarray(uint8_t , size_t *, int32_t *, uint8_t );
+ndarray_obj_t *ndarray_new_ndarray(uint8_t , size_t *, int32_t *, uint8_t , uint8_t *);
 ndarray_obj_t *ndarray_new_linear_array(size_t , uint8_t );
 ndarray_obj_t *ndarray_new_view(ndarray_obj_t *, uint8_t , size_t *, int32_t *, int32_t );
 bool ndarray_is_dense(ndarray_obj_t *);
diff --git a/code/numpy/carray/carray.c b/code/numpy/carray/carray.c
index 159c191d..bf0d7166 100644
--- a/code/numpy/carray/carray.c
+++ b/code/numpy/carray/carray.c
@@ -25,9 +25,11 @@
 
 #if ULAB_SUPPORTS_COMPLEX
 
+//| import builtins
+//|
 //| import ulab.numpy
 
-//| def real(val):
+//| def real(val: ulab.numpy.ndarray) -> ulab.numpy.ndarray:
 //|     """
 //|     Return the real part of the complex argument, which can be
 //|     either an ndarray, or a scalar."""
@@ -54,7 +56,7 @@ mp_obj_t carray_real(mp_obj_t _source) {
 
 MP_DEFINE_CONST_FUN_OBJ_1(carray_real_obj, carray_real);
 
-//| def imag(val):
+//| def imag(val: ulab.numpy.ndarray) -> ulab.numpy.ndarray:
 //|     """
 //|     Return the imaginary part of the complex argument, which can be
 //|     either an ndarray, or a scalar."""
@@ -82,7 +84,9 @@ MP_DEFINE_CONST_FUN_OBJ_1(carray_imag_obj, carray_imag);
 
 #if ULAB_NUMPY_HAS_CONJUGATE
 
-//| def conjugate(val):
+//| def conjugate(
+//|     val: builtins.complex | ulab.numpy.ndarray
+//| ) -> builtins.complex | ulab.numpy.ndarray:
 //|     """
 //|     Return the conjugate of the complex argument, which can be
 //|     either an ndarray, or a scalar."""
diff --git a/code/numpy/create.c b/code/numpy/create.c
index a2e0fce3..e8be7c24 100644
--- a/code/numpy/create.c
+++ b/code/numpy/create.c
@@ -1073,19 +1073,10 @@ mp_obj_t create_frombuffer(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
                 len = count;
             }
         }
-        ndarray_obj_t *ndarray = m_new_obj(ndarray_obj_t);
-        ndarray->base.type = &ulab_ndarray_type;
-        ndarray->dtype = dtype == NDARRAY_BOOL ? NDARRAY_UINT8 : dtype;
-        ndarray->boolean = dtype == NDARRAY_BOOL ? NDARRAY_BOOLEAN : NDARRAY_NUMERIC;
-        ndarray->ndim = 1;
-        ndarray->len = len;
-        ndarray->itemsize = sz;
-        ndarray->shape[ULAB_MAX_DIMS - 1] = len;
-        ndarray->strides[ULAB_MAX_DIMS - 1] = sz;
 
+        size_t *shape = ndarray_shape_vector(0, 0, 0, len);
         uint8_t *buffer = bufinfo.buf;
-        ndarray->array = buffer + offset;
-        return MP_OBJ_FROM_PTR(ndarray);
+        return ndarray_new_ndarray(1, shape, NULL, dtype, buffer + offset);
     }
     return mp_const_none;
 }
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 <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 */
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/numerical.c b/code/numpy/numerical.c
index e3b42525..0961e3c0 100644
--- a/code/numpy/numerical.c
+++ b/code/numpy/numerical.c
@@ -746,7 +746,7 @@ mp_obj_t numerical_argsort(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw
     numerical_reduce_axes(ndarray, ax, shape, strides);
 
     // We could return an NDARRAY_UINT8 array, if all lengths are shorter than 256
-    ndarray_obj_t *indices = ndarray_new_ndarray(ndarray->ndim, ndarray->shape, NULL, NDARRAY_UINT16);
+    ndarray_obj_t *indices = ndarray_new_ndarray(ndarray->ndim, ndarray->shape, NULL, NDARRAY_UINT16, NULL);
     int32_t *istrides = m_new0(int32_t, ULAB_MAX_DIMS);
     numerical_reduce_axes(indices, ax, shape, istrides);
 
@@ -1186,13 +1186,19 @@ mp_obj_t numerical_roll(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_ar
         mp_raise_TypeError(MP_ERROR_TEXT("roll argument must be an ndarray"));
     }
     ndarray_obj_t *ndarray = MP_OBJ_TO_PTR(args[0].u_obj);
-    uint8_t *array = ndarray->array;
     ndarray_obj_t *results = ndarray_new_dense_ndarray(ndarray->ndim, ndarray->shape, ndarray->dtype);
 
     int32_t shift = mp_obj_get_int(args[1].u_obj);
+
+    if(shift == 0) {
+        ndarray_copy_array(ndarray, results, 0);
+        return MP_OBJ_FROM_PTR(results);
+    }
+
     int32_t _shift = shift < 0 ? -shift : shift;
 
     size_t counter;
+    uint8_t *array = ndarray->array;
     uint8_t *rarray = (uint8_t *)results->array;
 
     if(args[2].u_obj == mp_const_none) { // roll the flattened array
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/numpy/random/random.c b/code/numpy/random/random.c
index 9be600fa..165f11b5 100644
--- a/code/numpy/random/random.c
+++ b/code/numpy/random/random.c
@@ -41,7 +41,7 @@ MP_DEFINE_CONST_OBJ_TYPE(
     MP_QSTR_generator,
     MP_TYPE_FLAG_NONE,
     print, random_generator_print,
-    make_new, random_generator_make_new, 
+    make_new, random_generator_make_new,
     locals_dict, &random_generator_locals_dict
 );
 #else
@@ -76,11 +76,12 @@ mp_obj_t random_generator_make_new(const mp_obj_type_t *type, size_t n_args, siz
     if(args[0] == mp_const_none) {
         #ifndef MICROPY_PY_RANDOM_SEED_INIT_FUNC
         mp_raise_ValueError(MP_ERROR_TEXT("no default seed"));
-        #endif
+        #else
         random_generator_obj_t *generator = m_new_obj(random_generator_obj_t);
         generator->base.type = &random_generator_type;
         generator->state = MICROPY_PY_RANDOM_SEED_INIT_FUNC;
         return MP_OBJ_FROM_PTR(generator);
+        #endif
     } else if(mp_obj_is_int(args[0])) {
         random_generator_obj_t *generator = m_new_obj(random_generator_obj_t);
         generator->base.type = &random_generator_type;
@@ -89,7 +90,7 @@ mp_obj_t random_generator_make_new(const mp_obj_type_t *type, size_t n_args, siz
     } else if(mp_obj_is_type(args[0], &mp_type_tuple)){
         mp_obj_tuple_t *seeds = MP_OBJ_TO_PTR(args[0]);
         mp_obj_t *items = m_new(mp_obj_t, seeds->len);
-        
+
         for(uint8_t i = 0; i < seeds->len; i++) {
             random_generator_obj_t *generator = m_new_obj(random_generator_obj_t);
             generator->base.type = &random_generator_type;
@@ -175,6 +176,8 @@ static mp_obj_t random_normal(size_t n_args, const mp_obj_t *pos_args, mp_map_t
 
     mp_float_t *array = (mp_float_t *)ndarray->array;
 
+    // numpy's random supports only dense output arrays, so we can simply
+    // loop through the elements in a linear fashion
     for(size_t i = 0; i < ndarray->len; i = i + 2) {
         #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
         uint32_t x = pcg32_next(&self->state);
@@ -246,7 +249,7 @@ static mp_obj_t random_random(size_t n_args, const mp_obj_t *pos_args, mp_map_t
         if(!mp_obj_is_type(out, &ulab_ndarray_type)) {
             mp_raise_TypeError(MP_ERROR_TEXT("out has wrong type"));
         }
-        
+
         ndarray = MP_OBJ_TO_PTR(out);
 
         if(ndarray->dtype != NDARRAY_FLOAT) {
@@ -281,10 +284,10 @@ static mp_obj_t random_random(size_t n_args, const mp_obj_t *pos_args, mp_map_t
 
     mp_float_t *array = (mp_float_t *)ndarray->array;
 
-    // numpy's random supports only dense output arrays, so we can simply 
+    // numpy's random supports only dense output arrays, so we can simply
     // loop through the elements in a linear fashion
     for(size_t i = 0; i < ndarray->len; i++) {
-        
+
         #if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
         uint32_t x = pcg32_next(&self->state);
         *array = (float)(int32_t)(x >> 8) * 0x1.0p-24f;
@@ -373,4 +376,3 @@ const mp_obj_module_t ulab_numpy_random_module = {
     .base = { &mp_type_module },
     .globals = (mp_obj_dict_t*)&mp_module_ulab_numpy_random_globals,
 };
-
diff --git a/code/ulab.h b/code/ulab.h
index 8d06fdc2..78b8a1ba 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 <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);
+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 83f03d43..c4a27e26 100644
--- a/docs/manual/source/conf.py
+++ b/docs/manual/source/conf.py
@@ -29,7 +29,6 @@
 # The full version, including alpha/beta/rc tags
 release = '6.6.0'
 
-
 # -- General configuration ---------------------------------------------------
 
 # Add any Sphinx extension module names here, as strings. They can be
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 91f86c28..d3e7587d 100644
--- a/docs/ulab-change-log.md
+++ b/docs/ulab-change-log.md
@@ -4,6 +4,18 @@ version 6.6.0
 
     add numpy.take
 
+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
+
 Wed, 6 Mar 2024
 
 version 6.5.2
diff --git a/docs/ulab-convert.ipynb b/docs/ulab-convert.ipynb
index 53eba42c..4ce30e5a 100644
--- a/docs/ulab-convert.ipynb
+++ b/docs/ulab-convert.ipynb
@@ -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,7 +257,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {
     "ExecuteTime": {
      "end_time": "2022-02-09T06:27:42.024028Z",
@@ -304,6 +305,7 @@
     "         'numpy-universal',\n",
     "         'numpy-fft',\n",
     "         'numpy-linalg',\n",
+    "        'numpy-random',\n",
     "         'scipy-linalg',\n",
     "         'scipy-optimize',\n",
     "         'scipy-signal',\n",
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-ndarray.ipynb b/docs/ulab-ndarray.ipynb
index 0e2bb1b0..8d67ed9b 100644
--- a/docs/ulab-ndarray.ipynb
+++ b/docs/ulab-ndarray.ipynb
@@ -3270,7 +3270,7 @@
      "output_type": "stream",
      "text": [
       "a:\t array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0], dtype=float)\n",
-      "a < 5:\t array([0.0, 1.0, 2.0, 3.0, 4.0], dtype=float)\n",
+      "a[a < 5]:\t array([0.0, 1.0, 2.0, 3.0, 4.0], dtype=float)\n",
       "\n",
       "\n"
      ]
@@ -3283,7 +3283,7 @@
     "\n",
     "a = np.array(range(9), dtype=np.float)\n",
     "print(\"a:\\t\", a)\n",
-    "print(\"a < 5:\\t\", a[a < 5])"
+    "print(\"a[a < 5]:\\t\", a[a < 5])"
    ]
   },
   {
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)):