Oaconvolve

[eleroy]

Hi everyone,

I need to perform fast convolution to do a matched filter on an ESP32.

The signal I have is about 30 000 points with a 100 point filter signal. I can reach about 600 ms per convolution.

Here is my code:

from ulab import numpy as np
import time
import math
def oaconvolve(signal, filter_signal, window_size=2048):
    n=len(filter_signal) # Length of filter
    block_size=window_size-n
    n_blocks = math.ceil(len(signal)/block_size)
    padded_signal = np.zeros(n_blocks*window_size)
    padded_signal[:len(signal)] = signal
    padded_filter = np.zeros(window_size)
    padded_filter[:len(filter_signal)] = filter_signal[::-1]
    fft_filter_signal= np.fft.fft(padded_filter)
    final_signal = np.zeros(len(padded_signal)+n)
    subsignal = np.zeros(window_size)
    for i in range(n_blocks):
        subsignal[:block_size] = padded_signal[i*block_size:(i+1)*block_size]
        a = np.fft.fft(subsignal)
        convolved_block = np.fft.ifft(a*fft_filter_signal)
        final_signal[i*block_size:(i)*block_size+window_size]+=convolved_block
    return final_signal

Is there any chance I can make it faster by coding it in C ? Can I use directly the ulab routine in C ?

Thank you

[v923z]

One problem with the fft in numpy is that it doesn't allow you to re-use RAM. This is addressed in https://micropython-ulab.readthedocs.io/en/latest/ulab-utils.html#spectrogram. That alone is not going to help you, because that calculates the absolute value, and not the FFT. However, you can easily change that in

micropython-ulab/code/utils/utils.c

Lines 207 to 378 in 2b74236

	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);
	}

So, my suggestion would be to modify the spectrogram method of the utils module, and use the scratchpad and out keyword arguments in your implementation of the FFT. That'll save you the overhead of reserving space.

Is there any chance I can make it faster by coding it in C ?

ulab is written in C, and the FFT routine used is one of the fastest ones. The overhead of python is simply the argument inspection at the beginning of the function call, and the memory allocation. But the FFT code can be re-used, because that's already decoupled from python here:

micropython-ulab/code/numpy/fft/fft_tools.c

Lines 48 to 92 in 2b74236

	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;
	j = 0;
	for(size_t i = 0; i < n; i++) {
	if (j > i) {
	SWAP(mp_float_t, data[2*i], data[2*j]);
	SWAP(mp_float_t, data[2*i+1], data[2*j+1]);
	}
	m = n >> 1;
	while (j >= m && m > 0) {
	j -= m;
	m >>= 1;
	}
	j += m;
	}
	mmax = 1;
	while (n > mmax) {
	istep = mmax << 1;
	theta = MICROPY_FLOAT_CONST(-2.0)*isign*MP_PI/istep;
	wtemp = MICROPY_FLOAT_C_FUN(sin)(MICROPY_FLOAT_CONST(0.5) * theta);
	wpr = MICROPY_FLOAT_CONST(-2.0) * wtemp * wtemp;
	wpi = MICROPY_FLOAT_C_FUN(sin)(theta);
	wr = MICROPY_FLOAT_CONST(1.0);
	wi = MICROPY_FLOAT_CONST(0.0);
	for(m = 0; m < mmax; m++) {
	for(size_t i = m; i < n; i += istep) {
	j = i + mmax;
	tempr = wr * data[2*j] - wi * data[2*j+1];
	tempi = wr * data[2*j+1] + wi * data[2*j];
	data[2*j] = data[2*i] - tempr;
	data[2*j+1] = data[2*i+1] - tempi;
	data[2*i] += tempr;
	data[2*i+1] += tempi;
	}
	wtemp = wr;
	wr = wr*wpr - wi*wpi + wr;
	wi = wi*wpr + wtemp*wpi + wi;
	}
	mmax = istep;
	}
	}

[v923z]

Since fft_kernel is an in-place routine, if you don't need your original data, then you can simply call that, and be done with it. If you end up adding your method, I'd like to ask you to create a PR, because this could be useful in general.

[eleroy]

Thank you ! I did succeed to write it in C.

With a signal about 30 000 points and a filter about 200 points

I get this performance:

ulab C OAconvolve 502.704 ms
micropython+ulab OAconvolve 716.302 ms
ulab linear Convolve 3482.501 ms

I'll try to open a PR. I guess I will include the routine oaconvolve in signal.c.

[v923z]

Yes, oaconvolve is part of scipy.signal, so that would be the place for that.

[v923z]

I get this performance:

ulab C OAconvolve 502.704 ms micropython+ulab OAconvolve 716.302 ms ulab linear Convolve 3482.501 ms

Do you know where the time is used? Because it might actually make more sense to add the out keyword to the FFT, and do the rest in python, as in your original implementation.

[eleroy]

Hi,

I think for small kernel size and large signal size, the fft segments become small and the python iteration becomes slow.

For larger kernel size, the iteration is not prevalent and I guess it's the memory allocation that sets the gap.

[v923z]

What I was driving at is the question, where you burn time: you already have an implementation in python, and one of the issues with that is that the FFT can't re-use the RAM. The rest won't be significantly faster in C, and that's why I was trying to suggest to fix the FFT code, and leave the rest in python. But I would commit to that only, if I knew for sure that not too much time is lost with the python loop. Hence my question.

[eleroy]

Time of computation as a function of kernel size. As kernel gets small the advantage of oaconvolve reduces.

There has been discussion about a C++ OAConvolve implementation scipy/scipy#16802 that might be helpful