Ondelettes

Fonction utilitaire pour afficher des courbes

In [1]: 
%matplotlib inline
from pylab import *

def dessine(*args):
    figure()
    for t in args:
        plot(range(len(t)), t, 'r')
    show()

Calcul de l'ondelette

Avec une base de 3 valeurs au lieu de 2 pour les ondelettes de Haar :

Au bord, on prolonge à gauche avec un miroir.

In [2]: 
def ondelette(x):
    if len(x) == 1:
        return x
    x = [x[1]] + x
    moyennes = []
    differences = []
    for i in range(1, len(x), 2):
        moyennes.append   ( 0.25 * x[i-1] + 0.5 * x[i] + 0.25 * x[i+1])
        differences.append(-0.25 * x[i-1] + 0.5 * x[i] - 0.25 * x[i+1])
    moyennes = ondelette(moyennes)
    # print x[1:], '→', moyennes + differences
    return moyennes + differences

On vérifie les cas simples :

In [3]: 
assert ondelette([1]) == [1]
assert ondelette([1, 1]) == [1, 0]
assert ondelette([1, 2]) == [1.5, -0.5]
assert ondelette([10, 11, 12, 13]) == [11.25, -0.75, -0.5, 0.0]
assert ondelette([10, 20, 30, 20]) == [20.0, -5.0, -5.0, 5.0]
assert ondelette([10, 11, 12, 13, 14, 10, 6, 2]) == [11.0625, 0.1875, -0.75, 1.875, -0.5, 0.0, 1.25, 0.0]

Affiche la courbe initiale et la transformée en ondelette. Elle est presque nulle partout.

In [4]: 
f = [256*(cos(i/50.)+cos(i/5.)/10.) for i in range(256)]
dessine(f, ondelette(f))

Ondelette inverse

On se rend compte de l'on obtient pour décompresser un système d'équation que l'on peut résoudre directement comme une jeu de dominos !

In [5]: 
def inverse(x):
    if len(x) <= 1:
        return x
    moyennes = inverse(x[:len(x)/2])
    differences = x[len(x)/2:]
    resultat = []
    for i in range(len(moyennes)):
        resultat.append(moyennes[i] + differences[i])
        if i != 0:
            resultat.append(-4*differences[i] + 2*resultat[-1] - resultat[-2])
        else:
            resultat.append(moyennes[i] - differences[i])
    # print x, '→', resultat
    return resultat

On vérifie que la décompression donne exactement la même valeur :

In [6]: 
for i in ([1, 1], [1, 2], [10, 11, 12, 13], [10, 20, 30, 20], [10, 11, 12, 13, 14, 10, 6, 2]):
    print inverse(ondelette(i)) == i, i

True [1, 1]
True [1, 2]
True [10, 11, 12, 13]
True [10, 20, 30, 20]
True [10, 11, 12, 13, 14, 10, 6, 2]

Malheureusement le résultat n'est pas terrible quand on modifie directement un des coefficients de l'ondelette.

In [80]: 
image = [cos(i/30.) for i in range(256)]

t = ondelette(image)
t[100] += 0.04
dessine(image, ondelette(image), inverse(t))

u = ondelette(image)
for i in range(64, len(t)):
    u[i] = int(u[i]/64+0.5) * 64   # Enleve 5 bits
dessine(ondelette(image), inverse(u))

Avec une base plus large

In [8]: 
def ondelette2(x):
    if len(x) == 1:
        return x
    if len(x) > 2:
        x = [x[2], x[1]] + x
    else:
        return ondelette(x)
    x += [x[-2]]
    moyennes = []
    differences = []
    for i in range(2, len(x)-2, 2):
        moyennes.append   ((x[i-2] + 4*x[i-1] + 6*x[i] + 4*x[i+1] + x[i+2]) / 16.)
        differences.append((x[i-2] - 4*x[i-1] + 6*x[i] - 4*x[i+1] + x[i+2]) / 16.)
    moyennes = ondelette2(moyennes)
    # print x[2:-1], '→', moyennes + differences
    return moyennes + differences
In [9]: 
assert ondelette2([1]) == [1]
assert ondelette2([1, 1]) == [1, 0]
assert ondelette2([1, 2]) == [1.5, -0.5]
assert ondelette2([10, 11, 12, 13]) == [11.3125, -0.5625, -0.25, -0.125]
assert ondelette2([10, 20, 30, 20]) == [20.625, -3.125, -2.5, 3.75]
assert ondelette2([10, 11, 12, 13, 14, 10, 6, 2]) == [11.07421875, 0.47265625, -0.453125, 1.3515625, -0.25, 0.0, 0.625, 0.5]
In [10]: 
f = [256*(cos(i/50.)+cos(i/5.)/10.) for i in range(256)]
dessine(f, ondelette2(f))
print "Moyenne des valeurs du signal", sum(abs(v) for v in f)/len(f)
print "Moyenne des valeurs de l'ondelette", sum(abs(v) for v in ondelette(f))/len(f)
print "Moyenne des valeurs de l'ondelette large", sum(abs(v) for v in ondelette2(f))/len(f)

Moyenne des valeurs du signal 155.622743271
Moyenne des valeurs de l'ondelette 2.78589055736
Moyenne des valeurs de l'ondelette large 1.0781068413

En JPEG2000 ce n'est pas la formule à utiliser car l'inversion est compliquée à faire numériquement. Voici ce qui est fait dans la norme. La première partie vérifie que l'inversion donne le bon résultat. Le codeur/décodeur fait les calculs sur place sans réserver de mémoire supplémentaire.

In [21]: 
# %pylab inline
import sympy

n = 16
i = [sympy.Symbol('i%s' % x) for x in range(n)]
o = [sympy.Symbol('o%s' % x) for x in range(n)]

s = i[::2]
d = i[1::2]

# Codeur jpeg2000

for x in range(len(d)):
    if x > 0:
        d[x] -= (s[x] + s[x-1]) / 2
    else:
        d[x] -= (s[x] + s[x]) / 2

for x in range(len(d)):
    if x != len(d)-1:
        s[x] += (d[x] + d[x+1]) / 4
    else:
        s[x] += (d[x] + d[x]) / 4

for x in range(len(d)-1):
    print "s=%-40s d=%s" % (s[x], d[x])

# Décodeur      

for x in range(len(d)):
    if x != len(d)-1:
        s[x] -= (d[x] + d[x+1])/4
    else:
        s[x] -= (d[x] + d[x])/4

for x in range(len(d)):
    if x != 0:
        d[x] += (s[x] + s[x-1])/2
    else:
        d[x] += (s[x] + s[x])/2

print "Ondelette inverse :"
for x in range(len(d)):
    print s[x], d[x],
print

    

s=5*i0/8 + i1/4 - i2/8 + i3/4              d=-i0 + i1
s=-i0/8 + 3*i2/4 + i3/4 - i4/8 + i5/4      d=-i0/2 - i2/2 + i3
s=-i2/8 + 3*i4/4 + i5/4 - i6/8 + i7/4      d=-i2/2 - i4/2 + i5
s=-i4/8 + 3*i6/4 + i7/4 - i8/8 + i9/4      d=-i4/2 - i6/2 + i7
s=-i10/8 + i11/4 - i6/8 + 3*i8/4 + i9/4    d=-i6/2 - i8/2 + i9
s=3*i10/4 + i11/4 - i12/8 + i13/4 - i8/8   d=-i10/2 + i11 - i8/2
s=-i10/8 + 3*i12/4 + i13/4 - i14/8 + i15/4 d=-i10/2 - i12/2 + i13
Ondelette inverse :
i0 i1 i2 i3 i4 i5 i6 i7 i8 i9 i10 i11 i12 i13 i14 i15
In [55]: 
def ondelette3(i):
    if len(i) == 1:
        return i
    s = i[::2]
    d = i[1::2]

    for x in range(len(d)):
        if x > 0:
            d[x] -= (s[x] + s[x-1]) / 2
        else:
            d[x] -= (s[x] + s[x]) / 2

    for x in range(len(d)):
        if x != len(d)-1:
            s[x] += (d[x] + d[x+1]) / 4
        else:
            s[x] += (d[x] + d[x]) / 4
    s = ondelette3(s)
    # print i, '→', s + d
    return s + d

f = [256*(cos(i/50.)+cos(i/5.)/10.) for i in range(256)]
# f = [10, 11, 12, 13]
dessine(f, ondelette3(f))
print "Moyenne des valeurs du signal", sum(abs(v) for v in f)/len(f)
print "Moyenne des valeurs de l'ondelette", sum(abs(v) for v in ondelette(f))/len(f)
print "Moyenne des valeurs de l'ondelette large", sum(abs(v) for v in ondelette2(f))/len(f)
print "Moyenne des valeurs de l'ondelette JPEG2000", sum(abs(v) for v in ondelette3(f))/len(f)

Moyenne des valeurs du signal 155.622743271
Moyenne des valeurs de l'ondelette 2.78589055736
Moyenne des valeurs de l'ondelette large 1.0781068413
Moyenne des valeurs de l'ondelette JPEG2000 24.2943704103
In [81]: 
def inverse3(i):
    if len(i) == 1:
        return i
    s = inverse3(i[:len(i)/2])
    d = i[len(i)/2:]
    for x in range(len(d)):
        if x != len(d)-1:
            s[x] -= (d[x]+d[x+1])/4
        else:
            s[x] -= (d[x]+d[x])/4

    for x in range(len(d)):
        if x != 0:
            d[x] += (s[x]+s[x-1])/2
        else:
            d[x] += (s[x]+s[x])/2
    t = []
    for x in range(len(d)):
        t.append(s[x])
        t.append(d[x])
    return t

f = [256*(cos(i/50.)+cos(i/5.)/10.) for i in range(256)]
t = ondelette3(f)
t[220] += 100

u = ondelette3(f)
for i in range(128, len(t)):
    u[i] = int(u[i]/16+0.5) * 16

dessine(f, inverse3(t), inverse3(u))