single dwt algorithm location (source code)

[pramodbutte]

I have an existing wavelet denoising function written in python using pywt. Due to some performance requirements we have moved most of the code over in Rust. One of the last remaining code is the discreet wavelent transform (bior3.9).

My issue is. I am noticing a small difference in my implementation of the wavelets compared to pywavelets results. I have been trying to understand how the wavelet decomposition is implemented in pywavelet. but I cannot seem to find the source code. can any please guide me to the source code where the actual convolution of the wavelet and signal is computed.

thanks in advance.

P.S. I was able to follow till

with nogil:
            retval_a = c_wt.double_dec_a(&data[0], data_size, wavelet.w,
                                <double *>cA.data, output_len, mode)
            retval_d = c_wt.double_dec_d(&data[0], data_size, wavelet.w,
                            <double *>cD.data, output_len, mode)

but I am unable to locate c_wt.double_dec_a implementation.

in case anyone is curious here is the rust implementation

fn wavelet_denoising(data: &mut PyArray1<f64>)-> PyResult<PyObject>{
    let gil = Python::acquire_gil();
    let py = gil.python();

    let dec_lo = vec![-0.0006797443727836989, 0.002039233118351097, 0.005060319219611981, -0.020618912641105536, -0.014112787930175844, 0.09913478249423216, 0.012300136269419315, -0.32019196836077857, 0.0020500227115698858, 0.9421257006782068, 0.9421257006782068, 0.0020500227115698858, -0.32019196836077857, 0.012300136269419315, 0.09913478249423216, -0.014112787930175844, -0.020618912641105536, 0.005060319219611981, 0.002039233118351097, -0.0006797443727836989];
    let dec_hi = vec![-0.0, 0.0, -0.0, 0.0, -0.0, 0.0, -0.0, 0.0, -0.1767766952966369, 0.5303300858899106, -0.5303300858899106, 0.1767766952966369, -0.0, 0.0, -0.0, 0.0, -0.0, 0.0, -0.0, 0.0];
    let rec_lo = vec![0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1767766952966369, 0.5303300858899106, 0.5303300858899106, 0.1767766952966369, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
    let rec_hi = vec![-0.0006797443727836989, -0.002039233118351097, 0.005060319219611981, 0.020618912641105536, -0.014112787930175844, -0.09913478249423216, 0.012300136269419315, 0.32019196836077857, 0.0020500227115698858, -0.9421257006782068, 0.9421257006782068, -0.0020500227115698858, -0.32019196836077857, -0.012300136269419315, 0.09913478249423216, 0.014112787930175844, -0.020618912641105536, -0.005060319219611981, 0.002039233118351097, 0.0006797443727836989];

    let dec_len = 20;
    let dec_offset = 0;
    let mut data = data.as_array_mut().to_vec();
    let data_len = data.len() as f64;
    let dwt_max_level : usize = (data_len/(dec_len as f64 -1.0)).log2().floor() as usize;
    let wvlt = wavelet::Wavelet{
        length: dec_len, offset: dec_offset , dec_lo: dec_lo, dec_hi: dec_hi, rec_lo: rec_lo,rec_hi: rec_hi
    };
    data.transform(Operation::Forward,&wvlt,5);
    let (data_bool, sigma) = mad_based_outlier(&data[480..].to_vec());
    let uthresh = (sigma/0.6744897501960817) * ((2.0* (data.len() as f64).ln()).sqrt());
    let mut d1: Vec<f64> = data[0..60].to_vec();
    let mut d2: Vec<f64> = data[60..].iter().map(|x| (x/x.abs()) * (x.abs()- uthresh).max(0.0)).collect::<Vec<f64>>();
    let mut data: Vec<f64> = d1.into_iter().chain(d2.into_iter()).collect(); 

    data.transform(Operation::Inverse,&wvlt,5);
    Ok(data.to_object(py))
}

and equivalent python code

def wavelet_smooth(x, wavelet="bior3.9", level=1, title=None):
    coeff = wavedec(x, wavelet)
    sigma = mad(coeff[-level])
    uthresh = sigma * np.sqrt(2 * np.log(len(x)))
    coeff[1:] = (threshold(k, value=uthresh, mode="soft") for k in coeff[1:])
    y = waverec(coeff, wavelet)
    return y

[grlee77]

Is this for a proprietary implementation? (If so that is fine, just curious to know what all PyWavelets is being used for).

double_dec_a gets created by some templating. The code underlying float_dec_a and double_dec_a is here: https://github.com/PyWavelets/pywt/blob/b65662ff4f7c33064d584e4b5ffec92bcffd19b2/pywt/_extensions/c/wt.template.c#L318 whereas the actual templating is done via #define / #undefine statements as in: https://github.com/PyWavelets/pywt/blob/master/pywt/_extensions/c/wt.c

double_dec_a ends up calling another templated function to do the actual convolution: https://github.com/PyWavelets/pywt/blob/b65662ff4f7c33064d584e4b5ffec92bcffd19b2/pywt/_extensions/c/convolution.template.c#L127

The convolution code ends up being a little hard to read mainly do to all the various boundary handling modes and the fact that it handles edge cases where the filter may be longer than the data, etc.

[pramodbutte]

This is being used in cancer research in a prototype being built in my lab. nothing proprietary.

If you are just curious regarding use case: My lab is building a prototype which uses fluorescence generated by ultra-fast lasers to detect brain tumor during surgery. We use pywavelets for denoising the signal from the tissue (mostly getting rid of thermal and shot noise prior to averaging multiple pulses)

[pramodbutte]

Thanks for the reply!

[grlee77]

My lab is building a prototype which uses fluorescence generated by ultra-fast lasers to detect brain tumor during surgery.

Okay, so you need real-time processing during surgery to guide the procedure? What is the typical image size and frame rate you need to process (I assume this is 2D + time data)? Contact me via email if you want to discuss further, as I don't want to clutter the issues page with an extended discussion here. (I have some prototype GPU-based implementations of various transforms and denoising algorithms that could potentially be of use. I have started looking at means to obtain funding to turn this in to a better maintained public package like PyWavelets, but that is still at the planning stage and part of why I was interested in potential downstream uses for such a package.)