Le Gall 5/3 Discrete Wavelet Transform

dwt.py

# Lifting scheme version of the discrete wavelet transform in JPEG 2000

# The approximation and detail coefficients are interleaved in the output of dwt (even 
# indices for approx, odd for detail). In practice you may need to separate them afterwards

# idwt is esentially reversing the order of the lifting scheme steps in dwt and
# flipping the sign of operations(i.e. addition becomes subtraction and vice versa)

# approx_len = approximation coefficients length = ceil(total_series_len / 2)
# detail_len = detail coefficients length = floor(total_series_len / 2)

# If you do out the math, the equivalent 1D convolution filters are:
#
# approx analysis: (-1/8, 2/8, 6/8, 2/8, -1/8)
# detail analysis: (-1/2, 1, -1/2)
#
# approx synthesis: (1/2, 1, 1/2)
# detail synthesis: (-1/8, -2/8, 6/8, -2/8, -1/8)
#
# These don't actually work due to the weird use of rounding though.
# Le Gall seems designed specifically for use in a lifting scheme.

import numpy as np

def slow_dwt(seq):
    seq = np.copy(seq)
    n = seq.shape[0]

    detail_len = n >> 1
    approx_len = n - detail_len

    for i in range(detail_len):
        seq[(i << 1) + 1] -= (seq[i << 1] + seq[min((i << 1) + 2, (approx_len - 1) << 1)]) >> 1
    
    for i in range(approx_len):
        seq[i << 1] += (seq[max((i << 1) - 1, 1)] + seq[min((i << 1) + 1, ((detail_len - 1) << 1) + 1)] + 2) >> 2

    return seq

def slow_idwt(seq):
    seq = np.copy(seq)
    n = seq.shape[0]

    detail_len = n >> 1
    approx_len = n - detail_len
    
    for i in range(approx_len):
        seq[i << 1] -= (seq[max((i << 1) - 1, 1)] + seq[min((i << 1) + 1, ((detail_len - 1) << 1) + 1)] + 2) >> 2

    for i in range(detail_len):
        seq[(i << 1) + 1] += (seq[i << 1] + seq[min((i << 1) + 2, (approx_len - 1) << 1)]) >> 1
    
    return seq

def dwt(seq):
    seq = np.copy(seq)
    n = seq.shape[0]

    detail_len = n >> 1
    approx_len = n - detail_len

    for i in range(1, detail_len):
        seq[(i << 1) - 1] -= (seq[(i - 1) << 1] + seq[i << 1]) >> 1
    seq[(detail_len << 1) - 1] -= (seq[(detail_len - 1) << 1] + seq[(approx_len - 1) << 1]) >> 1

    seq[0] += (seq[1] + seq[1] + 2) >> 2
    for i in range(1, approx_len - 1):
        seq[i << 1] += (seq[(i << 1) - 1] + seq[(i << 1) + 1] + 2) >> 2
    seq[(approx_len - 1) << 1] += (seq[(approx_len << 1) - 3] + seq[(detail_len << 1) - 1] + 2) >> 2

    return seq

def idwt(seq):
    seq = np.copy(seq)
    n = seq.shape[0]

    detail_len = n >> 1
    approx_len = n - detail_len
    
    seq[0] -= (seq[1] + seq[1] + 2) >> 2
    for i in range(1, approx_len - 1):
        seq[i << 1] -= (seq[(i << 1) - 1] + seq[(i << 1) + 1] + 2) >> 2
    seq[(approx_len - 1) << 1] -= (seq[(approx_len << 1) - 3] + seq[(detail_len << 1) - 1] + 2) >> 2

    for i in range(1, detail_len):
        seq[(i << 1) - 1] += (seq[(i - 1) << 1] + seq[i << 1]) >> 1
    seq[(detail_len << 1) - 1] += (seq[(detail_len - 1) << 1] + seq[(approx_len - 1) << 1]) >> 1

    return seq

seq = np.arange(11)

print(idwt(dwt(seq)))