Daubchies d4 wavelet : transform et inverse transform

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;
}

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.