Implementing Haar wavelet in python without packages

Question

I am trying to write a code to implement discrete wavelet transform (haar wavelet dwt) without using packages in python.

So far I've found a link where they implemented something similar, the link Is this wavelet transform implementation correct?. It doesn't give any errors while running, but the end result isn't correct. The code I ran is :

def discreteHaarWaveletTransform(x):
    N = len(x)
    output = [0.0]*N

    length = N >> 1
    while True:
        for i in range(0,length):
            summ = x[i * 2] + x[i * 2 + 1]
            difference = x[i * 2] - x[i * 2 + 1]
            output[i] = summ
            output[length + i] = difference

        if length == 1:
            return output

        #Swap arrays to do next iteration
        x = output[:length << 1]
        length >>= 1

Input :

list=[56, 40, 8, 24, 48, 48, 40, 16]

Current output :

[280, -24, 64, 40, 16, -16, 0, 24]

Expected output :

[35, -3, 16, 10, 8, -8, 0, 12]

Is there something obvious I can't see?

Answer 1

Something like this should do the trick -- it's almost a literal translation of this answer for C. It probably means that this is not very Python-y code, but if you're not using numpy for this that's not very Python-y anyway.

The main reason you did not get the output you expected is that you forgot to scale the output after filtering. This makes the coefficients at each next level approximately twice as high.

Mind that a scaling of ½ gives you the expected output, but a scaling of ½√2 is more commonly used, to preserve the L2-norm of signal under the wavelet transform.

def haarFWT ( signal, level ):

    s = .5;                  # scaling -- try 1 or ( .5 ** .5 )

    h = [ 1,  1 ];           # lowpass filter
    g = [ 1, -1 ];           # highpass filter        
    f = len ( h );           # length of the filter

    t = signal;              # 'workspace' array
    l = len ( t );           # length of the current signal
    y = [0] * l;             # initialise output

    t = t + [ 0, 0 ];        # padding for the workspace

    for i in range ( level ):

        y [ 0:l ] = [0] * l; # initialise the next level 
        l2 = l // 2;         # half approximation, half detail

        for j in range ( l2 ):            
            for k in range ( f ):                
                y [j]    += t [ 2*j + k ] * h [ k ] * s;
                y [j+l2] += t [ 2*j + k ] * g [ k ] * s;

        l = l2;              # continue with the approximation
        t [ 0:l ] = y [ 0:l ] ;

    return y

def main():
    s0 = [ 56, 40, 8, 24, 48, 48, 40, 16 ];

    print( "level 0" );
    print( s0 );

    print( "level 1" );
    print( haarFWT (s0, 1 ) );

    print( "level 2" );
    print( haarFWT (s0, 2 ) );

    print( "level 3" );
    print( haarFWT (s0, 3 ) );

if __name__ == "__main__":
    main()
# run with: >>> execfile ( "haarwavelet.py" )

Answer 2

A numpy-based approach to perform a multi-level 1D DWT signal decomposition using the Haar wavelet can be implemented this way. The results are consistent with the pywavelets implementation of pywt.wavedec(signal, "haar", mode="zero").

def scratch_haar_wavedec(signal, level=None, scale=np.sqrt(0.5)):
    """
    A numpy-based implementation of a multi-level 1D wavelet decomposition
    using the Haar wavelet. The input signal must be at least 2 samples long,
    and is always zero-padded to a length of (2^(level+1)) samples.

    Parameters
    ----------
    signal : np.ndarray
        A numpy array containing the input signal.
    level : int or None, optional
        A number of levels of the decomposition or 
        None for automatic estimation.
        The default is None.
    scale : float, optional
        A scaling factor for the decomposition. 
        The default value preserves L2 norm of the transformed signal.
        The default is np.sqrt(0.5).

    Returns
    -------
    res : List[np.ndarray]
        A list of numpy arrays containing the result in pywavelets format:
        [cA(n), cD(n), cD(n-1), ..., cD(1)], i.e., the approximate coefficients
        from the n-th level and all the detail coefficients ordered 
        from the n-th level to the 1st.

    """
    if len(signal) < 2: return None
    
    # get the level n (by maximizing n while 2^n <= signal length) if not specified
    n = int(np.log2(len(signal))) if level is None else level
    
    h = np.array([1.0,  1.0]) # low-pass filter (sum)
    g = np.array([1.0, -1.0]) # high-pass filter (difference)
    # zero-pad the signal to 2^(n+1) length
    t = np.pad(signal, (0,2**(n+1)-len(signal)), constant_values=0)
    
    res = list()
    for i in range(n):
        # calculate the low-pass and high-pass coefficients using correlation
        # the next iteration will process the new low-pass coefficients
        d = np.correlate(t, g)[::2] * scale
        t = np.correlate(t, h)[::2] * scale
        # store the high-pass (detail) coefficients (stripped zeros from right)
        res.append(d[:np.argwhere(d).max()+1])
    
    # store the final low-pass (approximate) coefficients
    res.append(t[:np.argwhere(t).max()+1])
    
    return res[::-1]

One may replace the correlation with convolution, i.e.:

d = np.convolve(t[::-1], g)[-2::-2] * scale
t = np.convolve(t[::-1], h)[-2::-2] * scale