1

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?

Siddhant Samal
  • 23
  • 2
  • 5
  • see the answer below - does that solve your problem? – alle_meije Sep 10 '19 at 21:44
  • Here is another implementation you can refer to. It uses numpy, but you can replace the numpy operations with pure Python if desired: https://tim.cogan.dev/wavelet/ – Tim Mar 11 '22 at 15:19

1 Answers1

3

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" )
alle_meije
  • 2,424
  • 1
  • 19
  • 40