Lifting en ondelettes

Lifting en ondelettes

Un lifting en ondelettes est, en mathématiques, un schéma d’implantation d’une transformation en ondelettes un peu différent de celui plus habituel réalisé par les bancs de filtres.

Le lifting en ondelettes est l’expression retenue pour désigner le procédé d’amélioration des propriétés des ondelettes par utilisation réciproque des bandes passe-bas et passe-haut. On obtient ainsi une transformation qui, à partir d’une ondelette à (admettons) \ n\ moments nuls, produit une transformée dont les propriétés sont similaires à une ondelette à \ \tilde{n}\ moments nuls supérieur à \ n\ .

Avantages et applications

Le procédé d’élévation (lifting) décrit ici a pour avantage d’avoir une complexité très faible, de l’ordre de la moitié de celle de la transformée rapide de Fourier FFT. La transformation peut être effectuée in situ, c’est-à-dire sans utilisation de mémoire supplémentaire. On perd cependant la propriété d’invariance par translation.

Cette technique de transformation est notamment utilisée dans l’algorithme de compression d’images JPEG2000, et trouvera d’autres applications dans le traitement numérique du signal en général et la transmission ou le stockage de données échantillonnées (notamment la compression du son, ou de mesures physiques de précision).

Introduction

Le but de ce document est de synthétiser de façon concise les techniques et connaissances permettant la réalisation de la transformation version lifting en ondelettes, et en particulier, de réaliser la transformation d’entiers en coefficients d’ondelettes entières. Voir en fin d’article les références essentielles aux articles théoriques et tutoriels.

Sommaire


Transformée en ondelettes version lifting

Description générale de la transformée en ondelettes version lifting

La transformation en ondelettes version lifting est un processus permettant entre autres d’optimiser le nombre d’opérations à exécuter et l’occupation mémoire (le coût de calcul est la moitié du coût de calcul de la FFT !).

Le processus le plus courant pour obtenir une transformation en ondelettes est d’utiliser un banc de filtres (cf. Figure 1), cependant lorsqu’on regarde cette façon de procéder on constate que l’on effectue un sous-échantillonnage par deux après une opération de filtrage : on a donc dépensé en pure perte la moitié du coût de calcul effectué par le filtrage.

L’opération de lifting en ondelettes peut-être vu comme la transformation réalisée par le banc de filtres, mais en intervertissant les phases de filtrage et de sous-échantillonnage. On limite ainsi le nombre d’opérations à effectuer mais, nous perdons en revanche la propriété d’invariance par translation.

Une autre propriété intéressante est que le schéma de lifting est facilement inversible.

Le schéma de lifting est aussi lié aux processus d’échantillonnage et d’interpolation (non explicitement étudié ici) pour la transformation de signaux continus.

On désignera ici par \ \gamma\ les coefficients d’ondelettes et par \ \lambda\ les coefficients d’échelle obtenus par la transformation.

Pour une ondelette particulière (c.à.d. un couple de filtres (h,\ g) ou (\tilde{h},\ \tilde{g}) dans l’implémentation par banc de filtres) caractérisée en outre par (disons) \ n\ moments nuls (jouant un rôle dans le processus de décroissance des coefficients d’ondelettes à travers les résolutions) pour le filtre primaire et \ \tilde{n}\ moments nuls pour le filtre dual, le schéma d’implantation par lifting permet d’obtenir facilement des ondelettes de moments \ n\ et \ \tilde{n}\ plus élevé. On a donc élevé (lifté) l’ordre de cette ondelette par ce schéma (d’où la justification du nom lifting).

Dans la Figure 1, le signal original \left(\ X\ \right) passe par les deux filtres complémentaires \left[\ \tilde{h}\ \right] (passe-bas) et \left[\ \tilde{g}\ \right] (passe-haut), en sortie on peut sous-échantillonner par 2 \left[\ \downarrow 2\ \right] et on obtient alors respectivement des coefficients d’ondelettes \left(\ \gamma\ \right) et des coefficients d’échelle \left(\ \lambda\ \right).

Figure 1

Figure 1. Schéma d’une implémentation en banc de filtres d’une transformation en ondelettes.

La reconstruction du signal s’effectue par sur-échantillonnage par insertion de zéros \left[\ \uparrow 2\ \right] et passage par les filtres de synthèse \left[\ h\ \right] et \left[\ g\ \right], puis par simple addition pour obtenir le signal \left(\ Y\ \right).

On utilisera :

  • les ondelettes paresseuses (lazy wavelets) qui servent à séparer un vecteur en composantes paires et impaires,
  • ainsi qu’une matrice polyphase qui permet de travailler sélectivement sur les composantes paires ou impaires du signal.

On va factoriser la matrice polyphase et introduire alors deux opérations :

  • une opération de prédiction (predict) qui prédit les échantillons de rang pair à partir des échantillons de rang impair ;
  • une opération de mise à jour (update) qui permet de conserver sur une partie du signal la valeur moyenne de l’ensemble du signal.

Le formalisme utilisé est celui des articles de références [Sweldens, Valens]. (Attention cependant aux écarts de notations entre les différents articles).

Cette transformation va en outre permettre de réaliser une transformation sur des entiers qui donne des entiers. Cependant il faudra utiliser une étape supplémentaire utilisant aussi le lifting pour la mise à l’échelle.

Vers la réalisation de l’ondelette de Haar version lifting

On part des filtres d’analyse et de reconstruction de l’analyse multi-résolution classique par banc de filtres suivant le schéma classique de la Figure 1 :

  • \tilde{h}[n] = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ \underline{1} \end{bmatrix}
  • \tilde{g}[n] = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -\underline{1} \end{bmatrix}
  • \tilde{h}(z) = \frac{1}{\sqrt{2}} (1 + z^{-1})
  • \tilde{g}(z) = \frac{1}{\sqrt{2}} (1 - z^{-1})
  • h[n] = \frac{1}{\sqrt{2}} \begin{bmatrix} \underline{1} \\ 1 \end{bmatrix}
  • g[n] = \frac{1}{\sqrt{2}} \begin{bmatrix} \underline{1} \\ -1 \end{bmatrix}
  • h(z) = \frac{1}{\sqrt{2}} (1 + z^{-1})
  • g(z) = \frac{1}{\sqrt{2}} (1 - z^{-1})

Attention aux notations utilisées ci-dessus :

  • dans certaines références, on trouve souvent \tilde{g}(z^{-1}) et \tilde{h}(z^{-1}).
  • dans les expressions des filtres, le coefficient souligné correspond à l’indice n = 0.

Les filtres doivent en outre satisfaire les formules suivantes :

  • h(z) \bar{h}(z^{-1}) + g(z) \bar{g}(z^{-1}) = 2
  • h(z) \bar{h}(-z^{-1}) + g(z) \bar{g}(-z^{-1}) = 0

On a besoin de construire la matrice polyphase \tilde{P} ainsi que sa matrice duale P (le terme polyphase vient de la théorie du filtrage numérique où il est utilisé pour décrire le partitionnement d’une séquence d’échantillons en plusieurs sous-séquences qui peuvent être traitées en parallèle, les sous-séquences peuvent-être vues comme des versions d’elles-mêmes décalées en phase, d’où le nom). On utilise la formule suivante de décomposition polyphase sur les filtres précédents :

  • x(z) = x_e(z^2) + z^{-1}\ x_o(z^2), où xe désigne les échantillons pairs (even) et xo les échantillons impairs (odd) (moyen mnémotechnique : even a un nombre de lettres pair, odd a un nombre de lettres impair).
  • {P}(z) = \begin{pmatrix} {h_e}(z) & {g_e}(z) \\ {h_o}(z) & {g_o}(z) \end{pmatrix}
  • \tilde{P}(z) = \begin{pmatrix} \tilde{h_e}(z) & \tilde{h_o}(z) \\ \tilde{g_e}(z) & \tilde{g_o}(z) \end{pmatrix}

Il vient assez facilement :

  • \tilde{h}(z) = \begin{matrix} \underbrace{\tilde{h_e}(z^2)} \\ 1 \end{matrix} + z^{-1} \begin{matrix} \underbrace{\tilde{h}_o(z^2)} \\ 1 \end{matrix}
  • \tilde{g}(z) = \begin{matrix} \underbrace{\tilde{g_e}(z^2)} \\ 1 \end{matrix} + z^{-1} \begin{matrix} \underbrace{\tilde{g}_o(z^2)} \\ -1 \end{matrix}
  • \tilde{P}(z)=\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}

sachant que \det\tilde{P} = -1, et que pour \det\tilde{P} = 1, on a les propriétés suivantes :

  • \tilde{h}_e(z) = g_o(z^{-1})
  • \tilde{h}_o(z) = -g_e(z^{-1})
  • \tilde{g}_e(z) = -h_o(z^{-1})
  • \tilde{g}_o(z) = h_e(z^{-1})

Ces propriétés montrent que pour passer de l’analyse à la synthèse (de P à \tilde{P}, il suffit si \det\tilde{P} = 1, de prendre la matrice des cofacteurs en changeant de signe.

De plus sachant que P(z)\tilde{P}(z^{-1}) = I (l’inversion temporelle z − 1 est nécessaire pour compenser le délai introduit par le filtrage) :

  • {P}(z) = \frac{1}{\sqrt{2}}\ \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}

Lifting primaire

L’élévation primaire (primary lifting) est aussi appelé update. Le lifting à proprement parler consiste à modifier le filtre g en gardant la complémentarité avec h à l’aide de la formule :

  • g_{new}(z) = g(z) + h(z)\ s(z^2)

où s(z) est un polynôme de Laurent.

Soit la transformée en Z d’un filtre FIR :

  • h(z)=\sum\limits_{k = p}^q h[k]z^{-k}.

Cette somme est aussi nommée polynôme de Laurent ou encore série de Laurent.

La matrice polyphase devient alors :

  • P_{new}(z) = P(z)\ \begin{pmatrix} 1 & s(z) \\ 0 & 1 \end{pmatrix}.

De la même façon pour \tilde{h},

  • \tilde{h}_{new}(z) = \tilde{h}(z) + \tilde{g}(z)\ \tilde{s}(z^2),

où \tilde{s}(z) est un polynôme de Laurent.

Dès lors,

  • \tilde{P}_{new}(z) = \begin{pmatrix} 1 & \tilde{s}(z) \\ 0 & 1 \end{pmatrix}\ \tilde{P}(z).

Lifting dual

L’élévation duale (lifting dual) est aussi appelé prédiction.

Les formules précédentes modifient la sous-bande passe-bas à l’aide de la sous-bande passe-haut. Le lifting dual consiste en l’opération inverse : modifier la sous-bande passe-bas à l’aide de la sous-bande passe-haut.

  • h_{new}(z) = h(z) + g(z)\ t(z^2),
  • P_{new}(z) = P(z)\ \begin{pmatrix} 1 & 0 \\ t(z) & 1 \end{pmatrix},

et

  • \tilde{g}_{new}(z) = \tilde{g}(z) + \tilde{h}(z)\ \tilde{t}(z^2),
  • \tilde{P}_{new}(z) = \begin{pmatrix} 1 & 0 \\ t(z) & 1 \end{pmatrix}\ \tilde{P}(z),

où t est un polynôme de Laurent.

On a ainsi élevé (lifté) le niveau de sophistication de la transformation en ondelettes. On a de plus, pour avoir une transformation inversible :

  • t(z) = -\tilde{t}(z) et
  • s(z) = -\tilde{s}(z).

Factorisation des filtres

La démarche ici est en quelque sorte la démarche inverse afin de pouvoir passer de forme connue de pairs de filtres d’ondelettes à leur implémentation en termes de lifting d’ondelettes.

On peut réécrire :

  • \tilde{g}_{new}(z) = \tilde{g}(z) + \tilde{h}(z)\ \tilde{t}(z^2)

en  :

  • \tilde{g}(z) = \tilde{h}(z)\ \tilde{t}(z^2) + \tilde{g}_{new}(z),

qui est une forme de division avec reste (pour réaliser la division avec reste, on utilise l’algorithme d’Euclide) où \tilde{g}_{new} est le reste.

Alors :

  • \tilde{P}(z) = \begin{pmatrix} 1 & 0 \\ \tilde{t}(z) & 1 \end{pmatrix}\ \tilde{P}_{new}(z).

En itérant cet exemple, on peut arriver à obtenir une matrice polyphase qui est de la forme :

  • P(z) = \begin{pmatrix} K_1 & 0 \\ 0 & K_2 \end{pmatrix}\ \prod\limits_{i=1}^m \begin{pmatrix} 1 & s_i(z) \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ t_i(z) & 1 \end{pmatrix},

où K1 et K2 sont deux constantes (différentes de zéro) que l’on peut éventuellement factoriser en quatre étapes de lifting.

Figure 2

Figure 2. Lifting, chaque voie (sous-bande) est modifiée par l’autre.

Exemple de la transformée lifting en ondelettes de Haar

On est déjà parti des filtres de la transformation en ondelettes de Haar et on a obtenu la matrice polyphase et sa forme duale. Il reste à factoriser pour avoir l’implémentation en lifting.

  • \tilde{P}(z) = \frac{1}{\sqrt{2}}\ \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}

Lifting dual

On veut :

  • \tilde{P}(z) = \tilde{P}_{new}(z)\ \begin{pmatrix} 1 & 0 \\ \tilde{t}(z) & 1 \end{pmatrix} = \frac{1}{\sqrt{2}}\ \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix},
  • \tilde{P}(z) = \begin{pmatrix} \tilde{h_e}^{new}(z) & \tilde{h}_o^{new}(z) \\ \tilde{g_e}^{new}(z) & \tilde{g}_o^{new}(z) \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ \tilde{t}(z) & 1 \end{pmatrix} = \begin{pmatrix} \tilde{h}_e^{new}+\tilde{h}_o^{new}(z)\ \tilde{t}(z) & \tilde{h}_o^{new}(z) \\ \tilde{g}_e^{new}+\tilde{g}_o^{new}(z)\ \tilde{t}(z) & \tilde{g}_o^{new}(z) \end{pmatrix}.

Par identification, il vient :

  • \left\lbrace\begin{matrix} \frac{1}{\sqrt{2}} & = & \tilde{h_o}^{new}(z) \\ -\frac{1}{\sqrt{2}} & = & \tilde{g_o}^{new}(z) \\ \frac{1}{\sqrt{2}} & = & \frac{1}{\sqrt{2}}\ \tilde{t}(z) + \tilde{h_e}^{new}(z) \\ \frac{1}{\sqrt{2}} & = & -\frac{1}{\sqrt{2}}\ \tilde{t}(z) + \tilde{g_e}^{new}(z) \end{matrix}\right.

Dans notre cas très simple, il n’est pas nécessaire d’utiliser l’algorithme d’Euclide (pour un exemple détaillé de l’utilisation de l’algorithme d’Euclide, voir le cas de l’ondelette de Daubechies D4 ci-après) et on peut prendre

Rappel : soit

  • -\frac{1}{8} z^{-1} + \frac{3}{4} -\frac{1}{8} z = \tilde{t}(z)\left(\frac{1}{4}+\frac{1}{4}z\right)+\tilde{h}_e^{new}(z).

On applique l’algorithme d’Euclide, soit a(z) = b(z)\ q(z) + r(z).

Soient ici :

  • \left\lbrace\begin{matrix} a_0(z) & = & a(z) & = & -\frac{1}{8} z^{-1} + \frac{3}{4} - \frac{1}{8} z \\ b_0(z) & = & b(z) & = & \frac{1}{4} + \frac{1}{4} z \end{matrix}\right., et :
  • \left\lbrace\begin{matrix} a_1(z) & = & b_0(z) \\ b_1(z) & = & a_0(z)\ \%\ b_0(z) \\ q_1(z) & = & a_0(z)\ /\ b_0(z) \end{matrix}\right.,

où \ \%\ est l’opérateur modulo et \ /\ retourne le quotient entier. La procédure est itérative, plusieurs choix sont en général possible pour le choix de l’ordre du diviseur.

  • -\frac{1}{8} z^{-1} + \frac{3}{4} - \frac{1}{8} z = \left\lbrace\begin{matrix} (-\frac{1}{2} z^{-1} + \frac{7}{2}) & (\frac{1}{4} + \frac{1}{4} z) & -z \\ (-\frac{1}{2} z^{-1} - \frac{1}{2}) & (\frac{1}{4} + \frac{1}{4} z) & -1 \\ (\frac{7}{2} z^{-1} - \frac{1}{2}) & (\frac{1}{4} + \frac{1}{4} z) & -z^{-1} \end{matrix}\right.
  • \tilde{g}_e^{new}(z) = 0,

d’où :

  • \left\lbrace\begin{matrix} \tilde{t}(z) & = & -1 \\ \tilde{h}_e^{new}(z) & = & \frac{2}{\sqrt{2}} \end{matrix}\right..

Alors,

  • \tilde{P}(z) = \frac{1}{\sqrt{2}} \begin{pmatrix} 2 & 1 \\ 0 & -1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix}.

Lifting primaire

On désire réécrire \tilde{P}(z) sous la forme :

  • \tilde{P}(z) = \begin{pmatrix} \tilde{h}_e^{new}(z) & \tilde{h}_o^{new}(z) \\ \tilde{g}_e^{new}(z) & \tilde{g}_o^{new}(z) \end{pmatrix}\ \begin{pmatrix} 1 & \tilde{s}(z) \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix} = \frac{1}{\sqrt{2}}\ \begin{pmatrix} 2 & 1 \\ 0 & -1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix}

Par identification, il vient :

  • \left\lbrace\begin{matrix} \frac{2}{\sqrt{2}} & = &\tilde{h_e}^{new}(z) \\ 0 & = & \tilde{g_e}^{new}(z) \\ \frac{1}{\sqrt{2}} & = & \frac{2}{\sqrt{2}}\ \tilde{s}(z)+\tilde{h_o}^{new}(z) \\ -\frac{1}{\sqrt{2}} & = & 0\ \tilde{s}(z)+\tilde{g_o}^{new}(z) \end{matrix}\right.

On prend :

  • \tilde{s}(z) = \frac{1}{2} et \tilde{h}_o^{new}(z) = 0, et
  • \tilde{g}_o^{new}(z) = -\frac{1}{\sqrt{2}}

Alors :

  • \tilde{P}(z) = \frac{1}{\sqrt{2}}\ \begin{pmatrix} 2 & 0 \\ 0 & -1 \end{pmatrix}\ \begin{pmatrix} 1 & \frac{1}{2} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix}, et
  • {P}(z) = \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & -\frac{1}{2} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} \frac{1}{2} & 0 \\ 0 & -1 \end{pmatrix}\ \sqrt{2}.

Ces opérations purement matricielles s’implémentent très facilement sous Scilab ou Matlab.

On peut de façon plus condensée (que l’écriture matricielle) écrire le pseudo code pour l’algorithme de calcul in-place :

  • xe et xo représentent les coefficients pairs et impairs du signal (respectivement),
  • d représente les détails (c.à.d. les coefficients issus du filtrage passe-haut, donc les coefficients d’ondelettes),
  • s représente le signal grossier (smooth) issu du filtrage passe-bas et donc les coefficients d’échelle.

Alors :

\begin{matrix} d & = & x_o - x_e; \\ s & = & x_o + d / 2; \\ s & = & 2 / \sqrt{2} * s; \\ d & = & -1 / \sqrt{2} * d;\end{matrix}

Pour le passage à une résolution supérieure on prend le code précédent que l’on réinitialise en commençant par :

\begin{matrix} x_o & = & s(1:2:length(s)); \\ x_e & = & s(2:2:length(s)); \end{matrix}.

La reconstruction se fait en inversant l’ordre des opérations (et leurs signes, sauf pour la normalisation où l’on prend l’inverse).

\begin{matrix} d & = & -\sqrt{2} * d; \\ s & = & 1/ \sqrt{2} * s; \\ x_o & = & s - d / 2; \\ x_e & = & x_o - d; \end{matrix}.

Implémentation sous Scilab de la transformée lifting en ondelettes de Haar

Programme lifting_Haar.sce

// "lifting_Haar.sce"
// Ce programme implemente une transformation en ondelettes de Haar par le
// lifting scheme de Sweldens (à partir de lifting_int3.sce).
//
// Globalement par rapport à la transformée "classique", on
// procède d’abord au sous-échantillonnage avant d’appliquer le filtre.
// On sépare coefficient pairs et impair (''split'').

clear
xdel(winsid());
stacksize(60000000);

/// Le signal à tester est un signal ECG de fréquence
/// d’échantillonnage de 500Hz.
mtlb_load('donnees_rlp_sci.mat'); // d1, d2, d3, VR, VL, VF, V3, V4, V5, V6
sig = d2; // On choisit la dérivation d2.
N = 4096; // Nb de points à considérer.

sig = round(sig(1:N)); // On arrondit pour travailler avec des entiers.
sigold = sig;

//--------------------------------------------------------------------------------
//- Notations et données.
//--------------------------------------------------------------------------------

// On va prendre l’ondelette de Haar...
// On utilise la décomposition polyphase (Valens, Sweldens).
// * P est la matrice polyphase.
// * gam designe les coefficients d’ondelettes.
// * lam designe les coefficients d’échelle.
// * Pred est la matrice permettant de faire la prédiction.
// * Upda est la matrice permettant de faire la mise à jour. 
// * Mnor est la matrice de normalisation.
// * fnor est le facteur de normalisation.
// * xe est le vecteur de composantes paires.
// * xo est le vecteur de composantes impaires.
// * X est le vecteur xe xo.

Wsig = []; // Coefficients d’ondelettes.

reso = 2; // Nombre de résolution.

fnor = 1 / sqrt(2);
fnorinv = sqrt(2);

Mnor = [2 0; 0 -1];
Mnorinv = inv(Mnor);

Pred = [1 0; -1 1];
Upda = [1 (1 / 2); 0 1];

Predi = [1 0; 1 1];
Updai = [1 (-1 / 2); 0 1];

//-----------------------------------------------------------
//--- transformation duale (\tilde{P}) ----------------------
//-----------------------------------------------------------

for iter = 1:1:reso
 xe = sig(2:2:$);
 xo = sig(1:2:$);
 X = [xe; xo];
 
 Resu = fnor * Mnor * Upda * Pred * X; 
 
 gam = Resu(2, :);
 lam = Resu(1, :);
 
 //--------------- Affichage ------------------------------------------
 
 xset('window', max(winsid() + 1)),
  subplot(211),
   plot2d([[1:length(gam)]', [1:length(lam)]'] ,[gam', lam'], [2, 3], '181',
    "gamma (ondelette) @lambda (échelle)"),
   xtitle("coefficients d’ondelettes et d’échelle", "temps", "amplitude"),
  subplot(212),
   plot2d([[1:length(gam)]'], [gam'], [2], '181',
    "coefficients d’ondelettes"),
   xtitle("coefficients d’ondelettes", "temps", "amplitude");
 
 Wsig = [gam Wsig]; // stockage des coefficients d’ondelettes
 sig = lam;
 
 /// -- aparté ---------------------------------------------------------
 /// autre facon de l’ecrire "in-place" (plus compacte, moins matriciellement)
 
 d = xo - xe;
 s = xe + d / 2;
 s = 2 / sqrt(2) * s;
 d = -1 / sqrt(2) * d;
 
 // xset("window", max(winsid() + 1)),
 //  subplot(211),
 //   plot2d([[1:length(s)]'], [s'], [3], '081')
 //  subplot(212),
 //   plot2d([[1:length(d)]'], [d'], [2], '081')
 
 /// -- fin de l’aparté ----------------------------------------------
 
end // fin de l’itération duale

WSsig = [lam Wsig]; // stockage des coefficients d’échelle et d’ondelettes

//--------------- Affichage ------------------------------------------

xset('window', max(winsid() + 1)),
 plot2d([[1:length(WSsig)]'], [WSsig'], [3], '181',
  "coefficient d’ondelettes et d’échelle"),
 //xtitle("coefficient d’ondelettes et d’échelle", "temps", "tension"),
 xtitle("coefficient d’échelle et d’ondelettes", string(reso) + " résolutions", "amplitude");

xset("window", max(winsid() + 1)),
 plot2d([[1:length(lam)]' [1:length(gam)]'], [lam', gam'], [3, 2], '081');

//-----------------------------------------------------------------
//------ Transformation primaire (P) ------------------------------
//-----------------------------------------------------------------

dgam = gam; // Initialisation pour la reconstruction in-situ.
slam = lam;

for iter = 1:1:reso
 
 rx = zeros(1, length(Resu)); 
 Orig = fnorinv * Predi * Updai * (Mnorinv) * Resu; 
 
 rxe = Orig(1, :);
 rxo = Orig(2, :);
 
 rx(1:2:length(Resu)) = (rxo);
 rx(2:2:length(Resu)) = (rxe);
 
 if iter < reso
  gam = WSsig((length(rx) + 1):(2 * length(rx)));
  lam = rx;
 end
 
 /// -- aparté ------------------------------------------------------------------
 /// autre facon de l’ecrire "in-place" (moins matriciellement)
 rrx = zeros(1, length(Resu)); 
 d = dgam;
 s = slam; 
 
 d = -sqrt(2) * d;
 s = 1/sqrt(2) * s;
 rrxe = s - d / 2;
 rrxo = rrxe + d;
 rrx(1:2:length(Resu)) = (rrxo);
 rrx(2:2:length(Resu)) = (rrxe);
 
 xset('window', max(winsid() + 1));
 plot2d([[1:length(rrx)]'], [rrx'], [2], '181',
  "signal reconstruit in place");
 xtitle("signal reconstruit in place", "temps", "tension");
 
 if iter < reso
  dgam = WSsig((length(rrx) + 1):(2 * length(rrx)));
  slam = rrx;
 end
 /// -- fin de l’aparté -----------------------------------------------------
 
 Resu = [lam; gam];
 
end

//--------------- Affichage ------------------------------------------

xset('window', max(winsid() + 1)),
 subplot(211),
  plot2d([[1:length(rx)]', [1:length(sigold)]'], [rx', sigold'], [2, 3], '181',
   "signal reconstruit@signal original"),
  xtitle("signal reconstruit et signal original", "temps", "tension"),
 subplot(212),
  plot2d([[1:length(sigold)]'], [(rx - sigold)'], [2], '181',
   "erreur entre les courbes"),
  xtitle("différence entre le signal reconstruit et le signal original", "temps", "écart");

// Fin du programme.

Résultats du programme lifting_Haar.sce

On présente quelques illustrations du programme : le signal analysé est un signal ECG (fréquence d’échantillonnage 500 Hertz). Le résultat de l’analyse (i.e. coefficients d’ondelettes et coefficients d’échelle) sont présentés sur la Figure 3.

Figure 3

Figure 3. Coefficients d’ondelettes et d’échelle de la résolution 1.

L’erreur de reconstruction de la Figure 4 est due à la présence de l'irrationnel \sqrt{2} dans les calculs.

Figure 4

Figure 4. Signal original, signal reconstruit et erreur (différence entre les deux signaux.

Un exemple de coefficients d’échelle et d’ondelettes obtenus pour 3 niveaux de résolution est présenté sur la Figure 5.

Figure 5

Figure 5. Coefficients d’échelle et d’ondelettes sur 3 niveaux de résolutions.

Exemple de la transformée lifting en ondelettes de Haar entières

Description de la transformée lifting en ondelettes de Haar entières

Obtenir des coefficients d’ondelettes entiers peut-être un atout lorsque :

  • des valeurs entières sont recommandées (notamment pour les images),
  • pour certains processeurs embarqués ne traitant que les entiers,
  • l’occupation mémoire est prépondérante (un entier occupe moins de place qu’un flottant).

Lors d’une étape de lifting que l’on a étudié on a utilisé une relation de ce type :

  • x_{new}(z) \longleftarrow x(z) + s(z)\ y(z)

Puisque le signal y(z) n’est pas modifié par cette étape de lifting alors on peut arrondir et réécrire l’étape du lifting.

  • x_{new}(z) \longleftarrow x(z) + \lfloor s(z)\ y(z) \rfloor

où \lfloor . \rfloor note une opération d’arrondi.

Ce qui est très fort, c’est que, quelle que soit l’opération d’arrondi utilisée, la transformation est inversible :

  • x(z) \longleftarrow x_{new}(z) - \lfloor s(z)\ y(z) \rfloor

Il y a cependant une précaution à prendre qui concerne la mise à l’échelle (scaling), dans le cas où ces coefficents ne sont pas entiers. La solution dans ce cas consiste à transformer cette mise à l’échelle en 3 étapes supplémentaires de lifting suivant le modèle :

  • P(z) = \begin{pmatrix} K & 0 \\ 0 & \frac{1}{K} \end{pmatrix} = \begin{pmatrix} 1 & K-K^2 \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ \frac{1}{K} & 1 \end{pmatrix}\ \begin{pmatrix} 1 & K-1 \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}, ou
  • P(z) = \begin{pmatrix} K & 0 \\ 0 & \frac{1}{K} \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 1-\frac{1}{K} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ K & 1 \end{pmatrix}\ \begin{pmatrix} 1 & \frac{1}{K^2}-\frac{1}{K} \\ 0 & 1 \end{pmatrix}

On rappelle :

  • \tilde{P}(z) = \frac{1}{\sqrt{2}} \begin{pmatrix} 2 & 0 \\ 0 & -1 \end{pmatrix}\ \begin{pmatrix} 1 & \frac{1}{2} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix} = \begin{pmatrix} \sqrt{2} & 0 \\ 0 & -\frac{1}{\sqrt{2}} \end{pmatrix}\ \begin{pmatrix} 1 & \frac{1}{2} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix}

Pour notre ondelette de Haar, on a un signe différent dans la matrice de normalisation, il suffit juste d’adapter les expressions précédentes avec par exemple pour la première (attention cependant lors de l’inversion) :

  • P(z) = \begin{pmatrix} K & 0 \\ 0 & -\frac{1}{K} \end{pmatrix} = \begin{pmatrix} 1 & K-K^2 \\ 0 & -1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ \frac{1}{K} & 1 \end{pmatrix}\ \begin{pmatrix} 1 & K-1 \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix},

avec K = \sqrt{2}.

Implémentation sous Scilab de la transformée lifting en ondelettes de Haar entières

Programme lifting_Haar_int.sce

// "lifting_Haar_int.sce"
// Ce programme implemente une transformation en ondelettes de Haar entières par le
// ''lifting scheme'' de Sweldens.

// Globalement par rapport à la transformée "classique", on
// procède d’abord au sous-échantillonnage avant d’appliquer le filtre.
// On sépare coefficient pairs et impair (''split'').

clear
xdel(winsid());
stacksize(60000000);

/// Le signal à tester est un signal ECG de fréquence
/// d’échantillonnage de 500Hz.
mtlb_load('donnees_rlp_sci.mat'); // d1, d2, d3, VR, VL, VF, V3, V4, V5, V6
sig = d2; // On choisit la dérivation d2.
N = 4096; // Nb de points à considérer.

sig = round(sig(1:N)); // ATTENTION on arrondit au départ pour travailler sur des entiers.
sigold = sig;

//--------------------------------------------------------------------------------
//- Maintenant le lifting d’ondelettes.
//--------------------------------------------------------------------------------

// On va prendre l’ondelette de Haar...
// On utilise la décomposition polyphase (Valens).
// * P est la matrice polyphase.
// * gam désigne les coefficients d’ondelettes.
// * lam désigne les coefficients d’échelle.
// * Pred est la matrice permettant de faire la prédiction.
// * Upda est la matrice permettant de faire la mise à jour. 
// * Mnor est la matrice de normalisation.
// * xe est le vecteur de composantes paires.
// * xo est le vecteur de composantes impaires.
// * X est le vecteur xe xo.

Wsig = [];

reso = 3; // Nombre de résolutions.

K = 1 / sqrt(2);

Mnor = [K 0; 0 -1/K]; 

P1 = [1 (K - K^2); 0 -1]; 
P2 = [1 0; (-1 / K) 1];
P3 = [1 (K - 1); 0 1];
P4 = [1 0; 1 1];

P1i = [1 (K - K^2); 0 -1]; 
P2i = [1 0; (1 / K) 1];
P3i = [1 (1 - K); 0 1];
P4i = [1 0; -1 1];

Pred = [1 0; -1 1];
Upda = [1 (1 / 2); 0 1]; 

Predi = [1 0; 1 1];
Updai = [1 (-1 / 2); 0 1]; 

//-----------------------------------------------------------
//--- Transformation duale (\tilde{P}) ----------------------
//-----------------------------------------------------------

for iter = 1:1:reso
 
 xe = sig(2:2:$);
 xo = sig(1:2:$);
 
 X = [xe; xo];
 
 PX = Pred * X;
 
 Updax = Upda(1, 1) * PX(1, :) + round(Upda(1, 2) *PX(2, :));
 Upday = Upda(2, 1) * PX(1, :) + round(Upda(2, 2) *PX(2, :));
 
 R4x = P4(1, 1) * Updax + round(P4(1, 2) * Upday);
 R4y = P4(2, 1) * Updax + round(P4(2, 2) * Upday);
 
 R3x = P3(1, 1) * R4x + round(P3(1, 2) * R4y);
 R3y = P3(2, 1) * R4x + round(P3(2, 2) * R4y);
 
 R2x = P2(1, 1) * R3x + round(P2(1, 2) * R3y);
 R2y = round(P2(2, 1) * R3x) + round(P2(2, 2) * R3y);
 
 R1x = P1(1, 1) * R2x + round(P1(1, 2) * R2y);
 R1y = P1(2, 1) * R2x + round(P1(2, 2) * R2y);
 
 Resu = [R1x; R1y];
 
 gam = Resu(2, :);
 lam = Resu(1, :);
 
 //--------------- Affichage ------------------------------------------
 
 xset('window', max(winsid() + 1)),
  subplot(211),
   plot2d([[1:length(gam)]', [1:length(lam)]'], [gam', lam'], [2, 3], '181',
    "gamma (ondelette) @lambda (échelle)"),
   xtitle("coefficients d’ondelettes et d’échelle", "temps", "amplitude"),
  subplot(212),
   plot2d([[1:length(gam)]'] ,[gam'], [2], '181',
    "coefficients d’ondelettes"),
   xtitle("coefficients d’ondelettes", "temps","amplitude");
 
 Wsig = [gam Wsig]; // Stockage des coefficients d’ondelettes.
 sig = lam;
 
end // Fin de l’itération duale.

WSsig = [lam Wsig]; // Stockage des coefficients d’échelle et d’ondelettes.

//--------------- Affichage ------------------------------------------

xset('window', max(winsid() + 1));
plot2d([[1:length(WSsig)]'], [WSsig'], [3], '181',
 "coefficient d’ondelettes et d’échelle");
//xtitle("coefficient d’ondelettes et d’échelle", "temps", "tension");
xtitle("coefficient d’échelle et d’ondelettes", string(reso)+ " résolutions", "amplitude");

xset("window", max(winsid() + 1)),
 plot2d([[1:length(lam)]' [1:length(gam)]'], [lam', gam'], [3, 2], '081')

//-----------------------------------------------------------------
//------ Transformation primaire (P) ------------------------------
//-----------------------------------------------------------------

for iter = 1:1:reso
 
 rx = zeros(1, length(Resu));
 
 R1xi = P1i(1, 1) * Resu(1, :) + round(P1i(1, 2) * Resu(2, :));
 R1yi = P1i(2, 1) * Resu(1, :) + round(P1i(2, 2) * Resu(2, :));
 
 R2xi = P2i(1, 1) * R1xi + round(P2i(1, 2) * R1yi);
 R2yi = round(P2i(2, 1) * R1xi) + round(P2i(2, 2) * R1yi);
 
 R3xi = P3i(1, 1) * R2xi + round(P3i(1, 2) * R2yi);
 R3yi = P3i(2, 1) * R2xi + round(P3i(2, 2) * R2yi);
 
 R4xi = P4i(1, 1) * R3xi + round(P4i(1, 2) * R3yi);
 R4yi = P4i(2, 1) * R3xi + round(P4i(2, 2) * R3yi);
 
 MR = [R4xi; R4yi];
 
 Updaxi = Updai(1, 1) * MR(1, :) + round(Updai(1, 2) * MR(2, :));
 Updayi = Updai(2, 1) * MR(1, :) + round(Updai(2, 2) * MR(2, :));
 
 Orig = Predi * [Updaxi; Updayi];
 
 rxe = Orig(1, :);
 rxo = Orig(2, :);
 
 rx(1:2:length(Resu)) = (rxo);
 rx(2:2:length(Resu)) = (rxe);
 
 if iter < reso
  gam = WSsig((length(rx) + 1):(2 * length(rx)));
  lam = rx;
 end
 
 Resu = [lam; gam];
 
end

//--------------- Affichage ------------------------------------------

xset('window', max(winsid() + 1)),
 subplot(211),
  plot2d([[1:length(rx)]', [1:length(sigold)]'] ,[rx', sigold'], [2, 3], '181',
   "signal reconstruit@signal original"),
  xtitle("signal reconstruit et signal original", "temps", "tension"),
 subplot(212),
  plot2d([[1:length(sigold)]'], [(rx - sigold)'], [2], '181',
   "erreur entre les courbes"),
  xtitle("différence entre le signal reconstruit et le signal original", "temps", "écart");

// Fin du programme.

Résultats du programme lifting_Haar_int.sce

On a alors obtenu les résultats des Figures 6 et 7. On n'a pas les infimes d’erreurs de reconstruction précédentes grâce aux arrondis.


Figure 6. Coefficients d’ondelettes et d’échelle de la résolution 1.


Figure 7. Coefficients d’échelle et d’ondelettes sur 3 niveaux de résolutions.

On peut remarquer dans le programme précédent que l’on a perdu en partie le côté pratique du formalisme matriciel (à cause des opérations d’arrondis). L'inconvénient majeur de conserver des entiers est d’avoir augmenté la dynamique du signal de sortie d’un facteur 2 (voir la figure 7).

Exemple de transformée lifting en ondelettes Daubechies(4)

Description de la transformée lifting en ondelettes Daubechies(4)

Daubechies(4) signifie que les filtres de cette ondelette possèdent 4 coefficients et n\ =\ 2 moments nuls. On a aussi :

  • \tilde{h}(z) = -z^{-1}\ g(-z^{-1}), et
  • \tilde{g}(z) = z^{-1}\ h(-z^{-1}).

On trouve son expression dans les tables ou dans [Daubechies] :

  • \left\lbrace\begin{matrix} h_0 & = & h[0] & = & \frac{1 + \sqrt{3}}{4 \sqrt{2}} \\ h_1 & = & h[1] & = & \frac{3 + \sqrt{3}}{4 \sqrt{2}} \\ h_2 & = & h[2] & = & \frac{3 - \sqrt{3}}{4 \sqrt{2}} \\ h_3 & = & h[3] & = & \frac{1 - \sqrt{3}}{4 \sqrt{2}}\end{matrix}\right.

On a les propriétés intéressantes suivantes :

  • \left\lbrace\begin{matrix} h_0^2 + h_1^2 + h_2^2 + h_3^2 = 1 \\ h_3\ h_1 + h_2\ h_0 = 0 \\ h_1^2 = -\frac{h_2\ h_1\ h_0}{h_3} \end{matrix}\right., soit :
  • \left\lbrace\begin{matrix} h[0] & = & 0,482962913145 \\ h[1] & = & 0,836516303738 \\ h[2] & = & 0,224143868042 \\ h[3] & = & -0,129409522551 \end{matrix}\right.
  • h(z)\ =\ h_0 + h_1 z^{-1} + h_2 z^{-2} + h_3z^{-3}, et
  • g(z)\ =\ -h_3 z^2 + h_2 z^{1} - h_1 + h_0z^{-1}
  • P(z)\ =\ \begin{pmatrix} h_e(z) & g_e(z) \\ h_o(z) & g_o(z) \end{pmatrix}

On a alors pour h :

  • h(z)\ =\ h_e(z^2) + z^{-1} h_o(z^2)
  • h(z)\ =\ h_0 + h_2 z^{-2} + z^{-1} (h_1 + h_3 z^{-2}), d’où
  • h_e(z)\ =\ h_0 + h_2 z^{-1}, et
  • h_o(z)\ =\ h_1 + h_3 z^{-1}

De la même façon pour g :

  • g(z)\ =\ -h_3 z^2 - h_1 + z^{-1} (h_2 z^2 + h_0)
  • g_e(z)\ =\ -h_3 z - h_1, et :
  • g_o(z)\ =\ h_2 z + h_0.

D’où :

  • P(z) = \begin{pmatrix} h_0 + h_2 z^{-1} & -h_3 z - h_1 \\ h_1 + h_3 z^{-1} & h_2 z + h_0 \end{pmatrix}

Lifting primaire

On va écrire une étape de lifting primaire :

  • P(z) = \begin{pmatrix} 1 & s(z) \\ 0 & 1 \end{pmatrix} \begin{pmatrix} h_e^{new}(z) & g_e^{new}(z) \\ h_o^{new}(z) & g_o^{new}(z) \end{pmatrix}

Alors,

  • \left\lbrace\begin{matrix} h_o^{new}(z) & = & h_1 + h_3 z^{-1} \\ g_o^{new}(z) & = & h_2 z + h_0 \\ h_0 + h_2z^{-1} & = & s(z) (h_1 + h_3 z^{-1}) + h_e^{new}(z) \\ -h_3z - h_1 & = & s(z) (h_2z + h_0) + g_e^{new}(z) \end{matrix}\right.

On utilise alors la division euclidienne (plusieurs factorisations sont possibles suivant l’ordre dans lequel on choisit de prendre en compte diviseur et dividende) :

Le nombre de possibilités différentes est lié aux degrés (noté \vert{.}\vert) des polynômes de Laurent des diviseurs et dividendes. Le degré d’un polynôme de Laurent quelconque h(z) défini par :

  • h(z) = \sum\limits_{k = p}^q h[k] z^{-k}

est : \vert{h}\vert = q - p. Soit a = b\ q + r l’expression de la division euclienne des polynômes a par b, alors

  • si les degrés du diviseur et du dividende sont égaux, alors le nombre de possibilités est 4 . 3^{\vert{b}\vert - 1},
  • si maintenant \vert{a}\vert = \vert{b}\vert + 1, alors le nombre de possibilités est 3^{\vert{b}\vert}.

Par exemple :

Image article lifting math143.png

ou alors :

Image article lifting math144.png

etc.

On se limite à l’exploration de la première solution (les autres aboutiront aussi à des formulations différentes mais équivalentes).

On a alors :

  • s(z) = \frac{h_2}{h_3}, et
  • h_e^{new}(z) = \frac{h_3 h_0 - h_2 h_1}{h_3}

On explicite alors g_e^{new}(z) :

  • \begin{matrix} g_e^{new}(z) & = & -h_3 z - h1 - \frac{h_2}{h_3} (h_2 z + h_0) \\ & = & -(\frac{h_3^2 + h_2^2}{h_3}) z - \frac{h_1 h_3 + h_2 h_0}{h_3} \\ & = & -(\frac{h_3^2 + h_2^2}{h_3}) z \end{matrix}, d’où
  • P(z) = \begin{pmatrix} 1 & \frac{h2}{h3} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} \frac{h_3 h_0 - h_2 h_1}{h_3} & -\frac{h_3^2 + h_2^2}{h_3} z \\ h_1 + h_3 z^{-1} & h_2 z + h_0 \end{pmatrix}

Lifting dual

Maintenant, le lifting dual :

  • P(z) = \begin{pmatrix} 1 & \frac{h2}{h3} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ t(z) & 1 \end{pmatrix}\ \begin{pmatrix} h_e^{new}(z) & g_e^{new}(z) \\ h_o^{new}(z) & g_o^{new}(z) \end{pmatrix} = \begin{pmatrix} 1 & \frac{h2}{h3} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} \frac{h_3 h_0 - h_2 h_1}{h_3} & -\frac{h_3^2 + h_2^2}{h_3} z \\ h_1 + h_3 z^{-1} & h_2 z + h_0 \end{pmatrix}

Alors,

  • \left\lbrace\begin{matrix} h_e^{new}(z) & = & \frac{h_3 h_0 - h_2 h_1}{h_3} \\ g_e^{new}(z) & = & -\frac{h_3^2 + h_2^2}{h_3}z \\ h_1 + h_3 z^{-1} & = & t(z) (\frac{h_3 h_0 - h_2 h_1}{h_3}) + h_o^{new}(z) \\ h_2 z + h_0 & = & t(z) (-\frac{h_3^2 + h_2^2}{h_3} z) + g_o^{new}(z) \end{matrix}\right.

Effectuons la division euclidienne suivante :

  • nothumb.

Donc :

  • t(z) = \frac{h_3 h_1}{h_3 h_0 - h_2 h_1} + \frac{h_3^2}{h_3 h_0 - h_2 h_1} z^{-1} et
  • h_o^{new}(z) = 0, alors,
  • \begin{matrix} g_o^{new}(z) & = & h_2 z + h_0 - \left( \frac{h_3 h_1 + h_3^2 z^{-1}}{h_3 h_0 - h_2 h_1} \right)\ \left( -\frac{h_3^2 + h_2^2}{h_3} z \right) \\ & = & \left( h_2 + \frac{h_1 (h_3^2 + h_2^2)}{h_3 h_0 - h_2 h_1} \right) z + \frac{h_3 (h_3^2 + h_2^2)}{h_3 h_0 - h_2 h_1} + h_0 \end{matrix}.

Or,

  • h_2 + \frac{h_1 (h_3^2 + h_2^2)}{h_3 h_0 - h_2 h_1} = 0,

d’où

  • g_o^{new}(z) = \frac{h_3 (h_3^2 + h_2^2)}{h_3 h_0 - h_2 h_1} + h_0.

Voilà :

  • P(z) = \begin{pmatrix} 1 & \frac{h_2}{h_3} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ \frac{h_3h_1}{h_3 h_0 - h_2 h_1} + \frac{h_3^2}{h_3 h_0 - h_2 h_1} z^{-1} & 1 \end{pmatrix}\ \begin{pmatrix} \frac{h_3 h_0 - h_2 h_1}{h_3} & -\frac{h_3^2 + h_2^2}{h_3} z \\ 0 & \frac{h_3 (h_3^2 + h_2^2)}{h_3 h_0 - h_2 h_1} + h_0 \end{pmatrix}

On peut maintenant mettre en évidence la mise à l’échelle d’abord en réécrivant un peu les derniers coefficients obtenus à l’aide des formules sur les coefficients h_i, 0 \le i \le 3 :

  • -\frac{h_3^2 + h_2^2}{h_3} = \frac{h_3}{h_3 h_0 - h_2 h_1} et
  • \frac{h_3 (h_3^2 + h_2^2)}{h_3 h_0 - h_2 h_1} + h_0 = h_3 \frac{h_3^2 + h_2^2 + h_0^2 - \frac{h_2 h_1 h_0}{h_3}}{h_3 h_0 - h_2 h_1} = \frac{h_3}{h_3 h_0 - h_2 h_1},

donc,

  • P(z) = \begin{pmatrix} 1 & \frac{h2}{h3} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ \frac{h_3 h_1}{h_3 h_0 - h_2 h_1} + \frac{h_3^2}{h_3 h_0 - h_2 h_1} z^{-1} & 1 \end{pmatrix}\ \begin{pmatrix} \frac{h_3 h_0 - h_2 h_1}{h_3} & \frac{h_3}{h_3 h_0 - h_2 h_1} z \\ 0 & \frac{h_3}{h_3 h_0 - h_2 h_1} \end{pmatrix}

Mise à l’échelle

On peut donc introduire l’étape de mise à l’échelle (normalisation), en remarquant que :

  • \begin{pmatrix} K & \frac{1}{K} z \\ 0 & \frac{1}{K} \end{pmatrix} = \begin{pmatrix} 1 & z \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} K & 0 \\ 0 & \frac{1}{K} \end{pmatrix}

d’où finalement :

  • P(z) = \begin{pmatrix} 1 & \frac{h2}{h3} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ \frac{h_3 h_1}{h_3 h_0 - h_2 h_1} + \frac{h_3^2}{h_3 h_0 - h_2 h_1} z^{-1} & 1 \end{pmatrix}\ \begin{pmatrix} 1 & z \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} \frac{h_3 h_0 - h_2 h_1}{h_3} & 0 \\ 0 & \frac{h_3}{h_3 h_0 - h_2 h_1} \end{pmatrix},

c’est-à-dire :

  • P(z) = \begin{pmatrix} 1 & -\sqrt{3} \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ \frac{\sqrt{3}}{4} + \frac{\sqrt{3} - 2}{4} z^{-1} & 1 \end{pmatrix}\ \begin{pmatrix} 1 & z \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} \frac{\sqrt{3} + 1}{\sqrt{2}} & 0 \\ 0 & \frac{\sqrt{3} - 1}{\sqrt{2}} \end{pmatrix}, et alors :
  • \tilde{P}(z) = \begin{pmatrix} \frac{\sqrt{3}-1}{\sqrt{2}} & 0 \\ 0 & \frac{\sqrt{3}+1}{\sqrt{2}} \end{pmatrix}\ \begin{pmatrix} 1 & -z \\ 0 & 1 \end{pmatrix}\ \begin{pmatrix} 1 & 0 \\ -\frac{\sqrt{3}}{4} - \frac{\sqrt{3} - 2}{4} z^{-1} & 1 \end{pmatrix}\ \begin{pmatrix} 1 & \sqrt{3} \\ 0 & 1 \end{pmatrix}

Simplification de la mise en œuvre

On peut de façon plus condensée (que l’écriture matricielle) écrire le pseudo-code pour l’algorithme de calcul in situ :

  • xe et xo représente les coefficients pairs et impairs du signal (respectivement),
  • d représente les détails (c.à.d. les coefficients issu du filtrage passe-haut, donc les coefficients d’ondelettes),
  • s représente le signal grossier (smooth) issu du fitlrage passe-bas et donc les coefficients d’échelle.
\begin{matrix} s & = & x_o + \sqrt{3} x_e; \\ d & = & x_e - \sqrt{3} / 4 * s - (\sqrt{3} - 2) / 4 * [0\ s(1 : 1 : length(s) - 1)]; \\ s & = & s - [d(2:1:length(d))\ 0]; \\ d & = & (\sqrt{3} + 1) / \sqrt{2} * d; \\ s & = & (\sqrt{3} - 1) / \sqrt{2} * s; \end{matrix}

Pour le passage à une résolution supérieure, on prend le code précédent que l’on réinitialise en commençant par :

\begin{matrix} x_o & = & s(1 : 2 : length(s)); \\ x_e & = & s(2 : 2 : length(s)); \end{matrix}

L’implémentation peut se faire facilement suivant l’exemple de la Figure 8.

Figure 8

Figure 8. Schéma d’une implémentation du lifting

Implémentation sous Scilab de la transformée lifting en ondelettes Daubechies(4)

Programme lifting_D4.sce

// "lifting_D4.sce"
// Ce programme implémente une transformation en ondelettes par le
// ''lifting scheme'' de Sweldens pour les ondelettes Daubechies(4) (D4).
/// 
// Globalement par rapport à la transformée "classique", on
// procède d’abord au sous-échantillonnage avant d’appliquer le filtre
// et on sépare coefficient pairs et impairs (''split'').

clear
xdel(winsid());
stacksize(60000000);

/// Le signal à tester est un signal ECG de fréquence
/// d’échantillonnage de 500Hz.
mtlb_load('donnees_rlp_sci.mat'); // d1, d2, d3, VR, VL, VF, V3, V4, V5, V6
sig = d2; // On choisit la dérivation d2.
N = 4096; // Nb de points à considérer.

//////-----------------------------------------------------------------------------------
h0 = (1 + sqrt(3)) / (4 * sqrt(2));
h1 = (3 + sqrt(3)) / (4 * sqrt(2));
h2 = (3 - sqrt(3)) / (4 * sqrt(2));
h3 = (1 - sqrt(3)) / (4 * sqrt(2));
//////-----------------------------------------------------------------------------------

sig = round(sig(1:N));// ATTENTION on arrondit au départ pour travailler sur des entiers.
sigold = sig;
 
//--------------------------------------------------------------------------------
//- Notations et données.
//--------------------------------------------------------------------------------

// On va prendre l’ondelette de Daubechies D4...
// On utilise la décomposition polyphase (Valens).
// * P est la matrice polyphase.
// * gam désigne les coefficients d’ondelettes.
// * lam désigne les coefficients d’échelle.
// * Pred est la matrice permettant de faire la prédiction.
// * Upda est la matrice permettant de faire la mise à jour. 
// * Mnor est la matrice de normalisation.
// * xe est le vecteur de composantes paires.
// * xo est le vecteur de composantes impaires.
// * X est le vecteur xe xo.

Wsig = [];

reso = 3; // Nombre de résolutions.

K = (sqrt(3) + 1) / sqrt(2);

Mnor = [(1 / K) 0; 0 K]; 

P1 = [1 (K - K^2); 0 1]; 
P2 = [1 0; (-1 / K) 1];
P3 = [1 (K - 1); 0 1];
P4 = [1 0; 1 1];

P1i = [1 (-K + K^2); 0 1]; 
P2i = [1 0; (1 / K) 1];
P3i = [1 (-K + 1); 0 1];
P4i = [1 0; -1 1];

Pred = [1 0; sqrt(3) 1];
Upda = [1 (-sqrt(3) / 4); 0 1];
// Attention Updazm(1, 2) fait intervenir échantillon précédent.
Updazm = [0 ((2 - sqrt(3)) / 4); 0 0]; 
Pred2 = [1 0; 0 1];
// Attention Pred2zp(2, 1) fait intervenir échantillon suivant.
Pred2zp = [0 0; -1 0]; 

Predi = [1 0; -sqrt(3) 1];
Updai = [1 (sqrt(3) / 4); 0 1];
// Attention Updaizm(1, 2) fait intervenir échantillon précédent.
Updaizm = [0 ((sqrt(3) - 2) / 4); 0 0]; 
Pred2i = [1 0; 0 1];
// Attention Predi2zp(2, 1) fait intervenir échantillon suivant.
Pred2izp = [0 0; 1 0]; 

//-----------------------------------------------------------
//--- transformation duale (\tilde{P}) ----------------------
//-----------------------------------------------------------

for iter = 1:1:reso
 
 xe = sig(2:2:$);
 xo = sig(1:2:$);
 
 X = [xe; xo];
 
 // d (details) = passe-haut = coeff d’ondelettes.
 Px = Pred(1, 1) * X(1, :) + Pred(1, 2) * X(2, :);
 // s (smooth) = passe-bas = coeff d’échelle.
 Py = Pred(2, 1) * X(1, :) + Pred(2, 2) * X(2, :);
 
 Updax = Upda(1, 1) * Px + Upda(1, 2) * Py + Updazm(1, 2) * [0 Py(1:1:$-1)]; // d
 Upday = Upda(2, 1) * Px + Upda(2, 2) * Py; // s
 
 Pxx = Pred2(1, 1) * Updax + Pred2(1, 2) * Upday; // d
 Pyy = Pred2(2, 1) * Updax + Pred2zp(2, 1) * [Updax(2:1:$) 0] + Pred2(2, 2) * Upday; // s
 
 R4x = P4(1, 1) * Pxx + P4(1, 2) * Pyy; // d
 R4y = P4(2, 1) * Pxx + P4(2, 2) * Pyy; // s
 
 R3x = P3(1, 1) * R4x + P3(1, 2) * R4y; // d
 R3y = P3(2, 1) * R4x + P3(2, 2) * R4y; // s
 
 R2x = P2(1, 1) * R3x + P2(1, 2) * R3y; // d
 R2y = P2(2, 1) * R3x + P2(2, 2) * R3y; // s
 
 R1x = P1(1, 1) * R2x + P1(1, 2) * R2y; // d
 R1y = P1(2, 1) * R2x + P1(2, 2) * R2y; // s
 
 Resu = [R1x; R1y];
 
 gam = Resu(1, :); // d
 lam = Resu(2, :); // s
 
 //--------------- Affichage ------------------------------------------
 
 xset('window', max(winsid() + 1));
 subplot(211),
  plot2d([[1:length(gam)]', [1:length(lam)]'], [gam', lam'], [2, 3], '181',
   "gamma (ondelette) @lambda (echelle)");
 xtitle("coefficients d’ondelettes et d’échelle", "temps", "amplitude");
 subplot(212),
  plot2d([[2:length(gam)]'], [gam(2:1:$)'], [2], '181',
   "coefficients d’ondelettes");
 xtitle("coefficients d’ondelettes", "temps", "amplitude");
 Wsig = [gam Wsig]; // Stockage des coefficients d’ondelettes.
 sig = lam;
 
 ///-- aparté --------------------------------------------
 /// autre facon de l’écrire "in-situ" (plus compacte, moins matriciellement)
 
 s = xo + sqrt(3) * xe;
 d = xe - sqrt(3) / 4 * s - (sqrt(3) - 2) / 4 * [0 s(1:1:$-1)];
 s = s - [d(2:1:$) 0];
 d = (sqrt(3) + 1) / sqrt(2) * d;
 s = (sqrt(3) - 1) / sqrt(2) * s;
 
 // xset("window", max(winsid() + 1)),
 //  subplot(211),
 //   plot2d([[1:length(s)]'], [s'], [3], '081')
 //  subplot(212),
 //   plot2d([[1:length(d)]'], [d'], [2], '081')
 
 ///-- fin de l’aparté -----------------------------------
 
 /// comparaison avec les banc de filtres
 /// wt paquetage --- uvi_wave ---
 /// wcoef = wt(round(sigold), dH, dG, iter, 0, 0);
 /// xset("window", max(winsid() + 1)),
 /// plot2d([1:length(wcoef)]', wcoef, [2], '081')
 
end // Fin de l’itération duale.

WSsig = [lam Wsig]; // Stockage des coefficients d’échelle et d’ondelettes.

//--------------- Affichage ----------------------------------------

xset('window', max(winsid() + 1)),
 plot2d([[1:length(WSsig)]'], [WSsig'], [3], '181',
  "coefficient d’ondelettes et d’échelle"),
 xtitle("coefficient d’échelle et d’ondelettes",
  string(reso) + " résolutions", "amplitude");

xset("window", max(winsid() + 1)),
 plot2d([[1:length(lam)]' [1:length(gam)]'], [lam', gam'], [3, 2], '081')

//-----------------------------------------------------------------
//------ transformation primaire (P) ------------------------------
//-----------------------------------------------------------------

dgam = gam; // Initialisation pour le in-situ.
slam = lam;

for iter = 1:1:reso
 
 rx = zeros(1, length(Resu)); 
 
 R1xi = P1i(1, 1) * Resu(1, :) + P1i(1, 2) * Resu(2, :);
 R1yi = P1i(2, 1) * Resu(1, :) + P1i(2, 2) * Resu(2, :);
 
 R2xi = P2i(1, 1) * R1xi + P2i(1, 2) * R1yi;
 R2yi = P2i(2, 1) * R1xi + P2i(2, 2) * R1yi;
 
 R3xi = P3i(1, 1) * R2xi + P3i(1, 2) * R2yi;
 R3yi = P3i(2, 1) * R2xi + P3i(2, 2) * R2yi;
 
 R4xi = P4i(1, 1) * R3xi + P4i(1, 2) * R3yi;
 R4yi = P4i(2, 1) * R3xi + P4i(2, 2) * R3yi;
 
 Pxxi = Pred2i(1, 1) * R4xi + Pred2i(1, 2) * R4yi;
 Pyyi = Pred2i(2, 1) * R4xi + Pred2izp(2, 1) * [R4xi(2:1:$) 0] + Pred2i(2, 2) * R4yi; 
 
 Updaxi = Updai(1, 1) * Pxxi + Updai(1, 2) * Pyyi + Updaizm(1, 2) * [0 Pyyi(1:1:$-1)];
 Updayi = Updai(2, 1) * Pxxi + Updai(2, 2) * Pyyi;
 
 Pxi = Predi(1, 1) * Updaxi + Predi(1, 2) * Updayi;
 Pyi = Predi(2, 1) * Updaxi + Predi(2, 2) * Updayi;
 
 Orig = [Pxi; Pyi];
 
 rxe = Orig(1, :);
 rxo = Orig(2, :);
 
 rx(1:2:length(Resu)) = (rxo);
 rx(2:2:length(Resu)) = (rxe);
 
 if iter < reso
  gam = WSsig((length(rx) + 1):(2 * length(rx)));
  lam = rx;
 end
 
 /// -- aparté --------------------------------------------------------------
 /// autre facon de l’écrire "in-situ" (moins matriciellement)
 rrx = zeros(1, length(Resu));
 d = dgam;
 s = slam; 
 
 s = (sqrt(3) + 1) / sqrt(2) * s;
 d = (sqrt(3) - 1) / sqrt(2) * d;
 s = s + [d(2:1:$) 0];
 rrxe = d + sqrt(3) / 4 * s + (sqrt(3) - 2) / 4 * [0 s(1:1:$-1)];
 rrxo = s - sqrt(3) * rrxe;
 rrx(1:2:length(Resu)) = (rrxo);
 rrx(2:2:length(Resu)) = (rrxe);
 
 xset('window', max(winsid() + 1)),
  plot2d([[1:length(rrx)]'], [rrx'], [2], '181',
   "signal reconstruit in-place"),
  xtitle("signal reconstruit in-place", "temps", "tension");
 
 if iter < reso
  dgam = WSsig((length(rrx) + 1):(2 * length(rrx)));
  slam = rrx;
 end
 
 /// -- fin de l'aparté ----------------------------------------------------
 
 Resu = [gam; lam];
 
end

//--------------- Affichage ------------------------------------------

xset('window', max(winsid() + 1)),
 subplot(211),
  plot2d([[1:length(rx)]', [1:length(sigold)]'] ,[rx', sigold'], [2, 3], '181',
   "signal reconstruit@signal original"),
  xtitle("signal reconstruit et signal original", "temps", "tension"),
 subplot(212),
  plot2d([[1:length(sigold)]'] , [(rx - sigold)'], [2], '181',
   "erreur entre les courbes"),
  xtitle("difference entre le signal reconstruit et le signal original", "temps", "écart");

// Fin du programme.

Résultats du programme lifting_D4.sce

On présente encore quelques illustrations du programme : le signal analysé est un signal ECG (fréquence d’échantillonnage 500 hertz). Le résultat de l’analyse (c.à.d. coefficients d’ondelettes et coefficients d’échelle) sont présentés sur la Figure 9.

Figure 9

Figure 9. Coefficients d’ondelettes et d’échelle de la résolution 1.

Les infimes erreurs de reconstruction de la Figure 10 sont dues à la présence des rationnels \sqrt{3} et \sqrt{2} dans les calculs.

Figure 10

Figure 10. Signal original, signal reconstruit et erreur (différence entre les deux signaux).

Un exemple de coefficients d’échelle et d’ondelettes obtenus pour 3 niveaux de résolution est présenté sur la Figure 11.

Figure 11

Figure 11. Coefficients d’échelle et d’ondelettes sur 3 niveaux de résolutions.

(On peut également, comme pour l’ondelette de Haar, effectuer avec l’ondelette de Daubechies la transformation d’entiers en coefficients d’ondelettes entiers.)

Voir aussi

Références externes


Wikimedia Foundation. 2010.

Contenu soumis à la licence CC-BY-SA. Source : Article Lifting en ondelettes de Wikipédia en français (auteurs)

Игры ⚽ Поможем сделать НИР

Regardez d'autres dictionnaires:

  • Lifting En Ondelettes — Un lifting en ondelettes est, en mathématiques, un schéma d’implantation d’une transformation en ondelettes un peu différent de celui plus habituel réalisé par les bancs de filtres. Le lifting en ondelettes est l’expression retenue pour désigner… …   Wikipédia en Français

  • Ondelette — En ondelette de Daubechies 2 Une ondelette est une fonction à la base de la décomposition en ondelettes, décomposition similaire à la transformée de Fourier à court terme, utilisée dans le traitement du signal. Elle correspond à l idée intuitive… …   Wikipédia en Français

  • Projet:Mathématiques/Liste des articles de mathématiques — Cette page n est plus mise à jour depuis l arrêt de DumZiBoT. Pour demander sa remise en service, faire une requête sur WP:RBOT Cette page recense les articles relatifs aux mathématiques, qui sont liés aux portails de mathématiques, géométrie ou… …   Wikipédia en Français

  • .JP2 — JPEG 2000 Comparaison du JPEG 2000 avec d autres formats JPEG 2000 ou ISO/CEI 15444 1 est une norme commune à l’ISO, la CEI et l’UIT T. C’est une norme de compression d’images produite par le groupe de travail Joint Photographic Experts Group.… …   Wikipédia en Français

  • .jp2 — JPEG 2000 Comparaison du JPEG 2000 avec d autres formats JPEG 2000 ou ISO/CEI 15444 1 est une norme commune à l’ISO, la CEI et l’UIT T. C’est une norme de compression d’images produite par le groupe de travail Joint Photographic Experts Group.… …   Wikipédia en Français

  • Codestream — JPEG 2000 Comparaison du JPEG 2000 avec d autres formats JPEG 2000 ou ISO/CEI 15444 1 est une norme commune à l’ISO, la CEI et l’UIT T. C’est une norme de compression d’images produite par le groupe de travail Joint Photographic Experts Group.… …   Wikipédia en Français

  • ISO 15444 — JPEG 2000 Comparaison du JPEG 2000 avec d autres formats JPEG 2000 ou ISO/CEI 15444 1 est une norme commune à l’ISO, la CEI et l’UIT T. C’est une norme de compression d’images produite par le groupe de travail Joint Photographic Experts Group.… …   Wikipédia en Français

  • JPEG2000 — JPEG 2000 Comparaison du JPEG 2000 avec d autres formats JPEG 2000 ou ISO/CEI 15444 1 est une norme commune à l’ISO, la CEI et l’UIT T. C’est une norme de compression d’images produite par le groupe de travail Joint Photographic Experts Group.… …   Wikipédia en Français

  • JPEG 2000 — Comparaison du JPEG 2000 avec d autres formats JPEG 2000 ou ISO/CEI 15444 1 est une norme commune à l’ISO, la CEI et l’UIT T. C’est une norme de compression d’images produite par le groupe de travail Joint Photographic Experts Group. JPEG 2000… …   Wikipédia en Français

  • JPG 2000 — JPEG 2000 Comparaison du JPEG 2000 avec d autres formats JPEG 2000 ou ISO/CEI 15444 1 est une norme commune à l’ISO, la CEI et l’UIT T. C’est une norme de compression d’images produite par le groupe de travail Joint Photographic Experts Group.… …   Wikipédia en Français

Share the article and excerpts

Direct link
… Do a right-click on the link above
and select “Copy Link”