Daubchies d4 wavelet : transform et inverse transform

Soyez le premier à donner votre avis sur cette source.

Snippet vu 6 133 fois - Téléchargée 37 fois

Contenu du snippet

Cette class permet d'effectuer la transformation (et l'inverse) en utilisant l'ondelette de daubechies a 4 moments.

Remarque: pour traiter l'effet bord du a ce type d'ondelette, j'ai utiliser la technique de "wrap around" tel que T+1=T-1 et donc T+2=T-2, au lieu du T+1=0 et T+2=1

Source / Exemple :


#ifndef DAUB4_H
#define DAUB4_H

#include <cmath>

	static const double c0=0.6830127;
	static const double c1=1.1830127;
	static const double c2=0.3169873;
	static const double c3=-0.1830127;
	static const double b0=c3;
	static const double b1=-c2;
	static const double b2=c1;
	static const double b3=-c0;
	static const double sqr= pow(2,-0.5);
 
class Daub
{
	private:
	int T;

	public:
	
	Daub(int n, double data[]):T(n){}
	static void InvTrans(double* data, int n);
	static void Transform(double* data, int n);
	~Daub(){}
};

void Daub::InvTrans(double* data, int n)
{

	int i,j;
	double* dt=new double[2*n];

	for(i=0; i<2*n; i++){dt[i]=0;}

	for(i=0, j=0; i<2*n; i+=2, j+=1){dt[i]=data[j];}

	for(i=0, j=0; i<n; i+=2, j+=1)
	{
		if(i<2)
		{
		data[0]=sqr*(c0*dt[0]+c2*dt[1]+b0*dt[n]+b2*dt[n+1]);
		data[1]=sqr*(c1*dt[0]+c3*dt[1]+b1*dt[n]+b3*dt[n+1]);
		}
		else
		{
		data[i]  =sqr*(c0*dt[i]+c2*dt[i-2]+b0*dt[n+i]+b2*dt[n+i-2]);
		data[i+1]=sqr*(c1*dt[i]+c3*dt[i-2]+b1*dt[n+i]+b3*dt[n+i-2]);
		}
	}
delete [] dt;
}

void Daub::Transform(double* data, int n)
{
	int i,g, wrp1, wrp2;
	double* dt= new double[n];

	if(n!=2){wrp1=n-2;wrp2=n-3;}
	else{wrp1=0;wrp2=1;}

	for(i=0, g=0 ; i<n ;i+=2, g++)
	{
		if(g==n/2-1)
		{
	dt[g]=sqr*(c0*data[i]+c1*data[i+1]+c2*data[wrp1]+c3*data[wrp2]);
	dt[g+n/2]=sqr*(b0*data[i]+b1*data[i+1]+b2*data[wrp1]+b3*data[wrp2]);
		}
		else
		{
	dt[g]=sqr*(c0*data[i]+c1*data[i+1]+c2*data[i+2]+c3*data[i+3]);
	dt[g+n/2]=sqr*(b0*data[i]+b1*data[i+1]+b2*data[i+2]+b3*data[i+3]);
		}
	}
	for(i=0; i<n ;i++){data[i]=dt[i];}
	delete [] dt;
}

#endif

Conclusion :


#include <iostream.h>
#include "Daub4.h"
//example de débruitage en ondelette

int main()
{
const int N=32;
int m,n;
//série bruité
double tabint[N] = {0.61739 ,2.96222,
3.26528 ,2.27141,3.82775,4.88492,7.23847,7.38413,
11.1395 ,8.95568,10.6950,11.3998,12.0592,13.6894,
14.6041 ,15.9120,17.1826,17.1258,18.5656,18.7819,
20.3416 ,20.4525,21.5375,23.8455,26.9481,26.8730,
26.1749 ,27.1218,28.9285,31.3041,31.1834,31.5983,
};

//calcule du nombre d'itération nécessaire pour effectuer la transformation à toutes les échelles

m=static_cast<int>(log(N)/log(2)-1);

//on effectue la transformation compléte

for(int i=0;i<=m;i++)
{
n=static_cast<int>(N/pow(2,i));
Daub::Transform(tabint, n);
}

//on annule le premier et certain des derniers de transformation
tabint[1]=0;
for(i=15;i<=31;i++){tabint[i]=0;}

//on reconstruit le signal avec inverse transform
for(i=m;i>=0;i--)
{
n=static_cast<int>(N/pow(2,i));
Daub::InvTrans(tabint, n);
}

//affichage des résultats
for(i=0 ; i<N;i++){cout<< tabint[i]<<endl;}
return 0;
}

A voir également

Ajouter un commentaire

Commentaires

coucou747
Messages postés
12336
Date d'inscription
mardi 10 février 2004
Statut
Modérateur
Dernière intervention
30 juillet 2012
26 -
je ne comrpends pas grand chose, a ton explication, tu pourrais détailler stp ?
plus_plus_fab
Messages postés
232
Date d'inscription
vendredi 9 janvier 2004
Statut
Membre
Dernière intervention
8 janvier 2005
-
salut,

quelques commentaires en vrac :

tu fais du C++ :
-#include <cmath>, met le en dessous des #ifndef , #define. Et utilises l'espace de nomage std.
- tu devrais mettre tes variables static const dans la classe (private).
- Le constructeur est faux
Daub(int n, double data[]) : T(n) {}
- indata est inutile, et T doit etre private.
- Invtrans et transform n'ont pas a etre static.
- pour l'efficacité, il faut à tout prix éviter les temporaires
dt[g]=sqr*(c0*data[i]+c1*data[i+1]+c2*data[wrp1]+c3*data[wrp2]) en crée un bon paquet par exemple !
- et pour finir, utlilise double*& ou lieu de double* dans la déclaration de transform, ça t'évitera de copier à la fin.

et pour aller + loin, utilise la généricité au lieu d'imposer un type double.

bonne continuation ...
cs_hitcher
Messages postés
21
Date d'inscription
jeudi 8 avril 2004
Statut
Membre
Dernière intervention
15 septembre 2009
-
Pour explication cet opérateurs mathématiques transforme un signal de longueur 2^i , avec i appartenant à Nen 2 série de longeurs 2^(i-1)
L'une contenant les coefficients d'echelle et l'autre les coefficients de transformation.
Cet opérateur est très utilisé un débruitage, reconnaissance de forme et filtrage,....
pour plus d'info sur les ondelettes: http://www.bearcave.com/misl/misl_tech/wavelets/


Pour ce qui est du code en lui même, effecitvement je suis totalement autodidacte et je n'arrive pas à mettre en oeuvre certaine des suggestions. Par contre la non utilisation des templates est voulu car l'output ne peut en aucun cas être int.
plus_plus_fab
Messages postés
232
Date d'inscription
vendredi 9 janvier 2004
Statut
Membre
Dernière intervention
8 janvier 2005
-
avec un exemple, c'est plus clair ! et ça laisse apparaitre la fragilité de la conception ...
En gros, tu donnes cette classe a quelqu'un qui connait quand meme le principe, il est incapable de faire quoi que ce soit ! (mon cas)

mon idee serait de faire une conception totalement différente qui ressemblerait à ceci :

template <class Alg>
class Signal
{
private:
typedef typename Alg::value_type value_type;
typedef int integer_type;
typedef typename Alg::seuillage_type seuillage_type;
typedef void (typename Alg::*ptr_to_fseuil)(seuillage_type);

public:
Signal(integer_type size, const value_type* signal_input);
~Signal();
void transform(const Alg&);
void restore(ptr_to_fseuil);
const T* coef_begin()const;
const T* coef_end()const;
const T* sig_res_begin()const;
const T* sig_res_end()const;
const T* sig_init_begin()const;
const T* sig_init_end()const;

private:
T* coefs;
T* sig_res;
const T* sig_init;

};

tu n'aurais plus qu'a ecrire la class Alg, et ça fonctionnerai que ce soit de la transformée de fourrier, de l'ondelettes de Haar, ou à 4 moments, ...
cette class "Alg" (Daub), devra fournir les types et méthodes suivantes :
types : value_type, seuillage_type.
methodes :
void transform(const value_type* input, value_type* output)const;
void restore(const value_type* input, value_type* output, ptr_to_fseuil)const;
Le constructeur sera vide, il n'y aura pas de membres pour cette classe.

pour moi, restore applique une fonction de seuillage (présente dans la classe Daub (Alg)) sur les coefficients calculées (input), puis restaure le signal à partir des coefficients seuillés (output).

Tu te sens de le faire ? Je peux t'aider en tout cas, ça m'intéresse, contacte moi en privé.

Quand ce sera fonctionnel, il faudra optimiser les calculs à coup d'élimination de temporaires et de valarray<>, on en est pas la ...

@+
plus_plus_fab
Messages postés
232
Date d'inscription
vendredi 9 janvier 2004
Statut
Membre
Dernière intervention
8 janvier 2005
-
Si on veut separer le seuillage de la "restoration", on peut faire comme ça : (je pense que c'est mieux)

void restore(const value_type* input, value_type* output)const;

void seuillage(value_type* input_output, seuillage_type seuil)const;
j'ignore pour l'instant quel parametre invoquer pour effecteur un seuillage, donc c'est un peu au pif ça ...

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.