wavelet.h

东南大学版本

#include "math.h"
#define FW double
#define TB 0
double hb[9],hr[7];
void flwt(FW* pa,int Len)
{
	FW* p=new FW [Len];
	for(int i=0;i<Len;i++)
	{
		p[i]=pa[i];
	}
	for(i=0;i<Len;i++)
	{
		if(i<(Len/2))
		{
			pa[i]=p[i<<1];
		}
		else
		{
			pa[i]=p[2*(i-Len/2)+1];
		}
	}
	FW* c=pa;
	FW* d=pa+Len/2;
/////////////////////////////
	for(i=0;i<Len/2;i++)
	{
		int Lenc=Len/2;
        d[i]=d[i]+((1*c[(i-1+Lenc)%Lenc]+   
			      -9*c[i]+
				  -9*c[(i+1+Lenc)%Lenc]+
				   1*c[(i+2+Lenc)%Lenc]+2)/16.0);
	}
////////////////////////////
	for(i=0;i<Len/2;i++)
	{
		int Lenc=Len/2;
        c[i]=c[i]+(1.0/4)*(d[i]+2)+(1.0/4)*(d[(i-1+Lenc)%Lenc]+2);
	}
/////////////////////////////
	for(i=0;i<Len/2;i++)
	{
        d[i]=d[i]/sqrt(2);
		c[i]=c[i]*sqrt(2);
	}
}

void iflwt(FW* pa,int Len)
{
//////////////////////////sa
	FW* c=pa;
	FW* d=pa+Len/2;
	for(int i=0;i<Len/2;i++)
	{
        d[i]=d[i]*sqrt(2);
		c[i]=c[i]/sqrt(2);
	}
////////////////////////////
	for(i=0;i<Len/2;i++)
	{
		int Lenc=Len/2;
        c[i]=c[i]-(1.0/4)*(d[i]+2)-(1.0/4)*(d[(i-1+Lenc)%Lenc]+2);
	}
///////////////////////////
	for(i=0;i<Len/2;i++)
	{
		int Lenc=Len/2;
        d[i]=d[i]+((-1*c[(i-1+Lenc)%Lenc]+   
			      9*c[i]+
				  9*c[(i+1+Lenc)%Lenc]-
				   1*c[(i+2+Lenc)%Lenc]+2)/16.0);
	}
	FW* p=new FW [Len];
	for(i=0;i<Len;i++)
	{
		p[i]=pa[i];
	}
	for(i=0;i<Len/2;i++)
	{
		pa[2*i]=p[i];
		pa[2*i+1]=p[i+Len/2];
	}
	delete p;
}
void flwt2D(FW** pd,int height,int width)
{
	double t=1.24087;
	for(int i=0;i<height;i++)
	{
		flwt(pd[i],width);
	}

	FW** pd1=new FW* [width];
	for(i=0;i<width;i++)
	{
		pd1[i]=new FW [height];
	}
	for(i=0;i<width;i++)
	{
		for(int j=0;j<height;j++)
		{
			pd1[i][j]=pd[j][i];
		}
	}

	for(i=0;i<width;i++)
	{
	    flwt(pd1[i],height);
	}

	for(i=0;i<width;i++)
	{
		for(int j=0;j<height;j++)
		{
			pd[j][i]=pd1[i][j];
		}
	}

	for(i=0;i<width;i++)
	{
		delete pd1[i];
	}
	delete pd1;
}

void iflwt2D(FW** pd,int height,int width)
{
	double t=1.24087;
	FW** pd1=new FW* [width];
	for(int i=0;i<width;i++)
	{
		pd1[i]=new FW [height];
	}
	for(i=0;i<width;i++)
	{
		for(int j=0;j<height;j++)
		{
			pd1[i][j]=pd[j][i];
		}
	}

	for(i=0;i<width;i++)
	{
		iflwt(pd1[i],height);
	}

	for(i=0;i<width;i++)
	{
		for(int j=0;j<height;j++)
		{
			pd[j][i]=pd1[i][j];
		}
	}

	for(i=0;i<width;i++)
	{
		delete pd1[i];
	}
	delete pd1;

	for(i=0;i<height;i++)
	{
		iflwt(pd[i],width);
	}
}
void flwt2D(FW** pd,int height,int width,int N)
{
	for(int i=0;i<N;i++)
	{
		flwt2D(pd,height,width);
		height>>=1;
		width>>=1;
	}
}
void iflwt2D(FW** pd,int height,int width,int N)
{
	height>>=(N-1);
	width>>=(N-1);
	for(int i=0;i<N;i++)
	{
		iflwt2D(pd,height,width);
		height<<=1;
		width<<=1;
	}
}