I'm looking for a way to partially reconstruct branches of a wavelet decomposition, such that the sum would recreate the original signal. This could be achieved using Matlab using:
DATA = [0,1,2,3,4,5,6,7,8,9]
N_LEVELS = 2;
WAVELET_NAME = 'db4';
[C,L] = wavedec(DATA, N_LEVELS, WAVELET_NAME);
A2 = wrcoef('a', C, L, WAVELET_NAME, 2);
D2 = wrcoef('d', C, L, WAVELET_NAME, 2);
D1 = wrcoef('d', C, L, WAVELET_NAME, 1);
A2+D2+D1
ans =
0.0000 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000
I'd like to achieve the same using pywt, but I'm not sure how to go about this.
The pywt.waverec
function creates a full reconstruction, but doesn't have a level parameter for partial reconstruction.
The pywt.upcoef
function does what I need for a single level but I'm not sure how to expand this for multiple levels:
>>> import pywt
>>> data = [1,2,3,4,5,6]
>>> (cA, cD) = pywt.dwt(data, 'db2', 'smooth')
>>> n = len(data)
>>> pywt.upcoef('a', cA, 'db2', take=n) + pywt.upcoef('d', cD, 'db2', take=n)
array([ 1., 2., 3., 4., 5., 6.])
I managed to write my own version of the wrcoef
function which appears to work as expected:
import pywt
import numpy as np
def wrcoef(X, coef_type, coeffs, wavename, level):
N = np.array(X).size
a, ds = coeffs[0], list(reversed(coeffs[1:]))
if coef_type =='a':
return pywt.upcoef('a', a, wavename, level=level)[:N]
elif coef_type == 'd':
return pywt.upcoef('d', ds[level-1], wavename, level=level)[:N]
else:
raise ValueError("Invalid coefficient type: {}".format(coef_type))
level = 4
X = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]
coeffs = pywt.wavedec(X, 'db1', level=level)
A4 = wrcoef(X, 'a', coeffs, 'db1', level)
D4 = wrcoef(X, 'd', coeffs, 'db1', level)
D3 = wrcoef(X, 'd', coeffs, 'db1', 3)
D2 = wrcoef(X, 'd', coeffs, 'db1', 2)
D1 = wrcoef(X, 'd', coeffs, 'db1', 1)
print A4 + D4 + D3 + D2 + D1
# Results:
[ 9.99200722e-16 1.00000000e+00 2.00000000e+00 3.00000000e+00
4.00000000e+00 5.00000000e+00 6.00000000e+00 7.00000000e+00
8.00000000e+00 9.00000000e+00 1.00000000e+01 1.10000000e+01
1.20000000e+01 1.30000000e+01 1.40000000e+01 1.50000000e+01
1.60000000e+01 1.70000000e+01]