intdwt97.cpp

小波变换去噪和97提升小波和离散小波变换的matlab程序

// LiftDwt.cpp: implementation of the CLiftDwt class.
//  作者：国防科技大学理学院数学系 张增辉    2001.11
//  QQ: 15105473   email: [e-mail]
//
//  使用带参数的9-7小波滤波器进行的整数小波变换，采用了参数9-7小波滤波器的提升格式，但是当参数为0.5时候，提升格式不同，因此，此程序中参数选择范围不能包括0.5
//  
//
//
//
//
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "IntDwt97.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CIntDwt97::CIntDwt97()
{
	p_t = 0.75;
	shift = 4;
	filterset(p_t);  // 设置滤波器的系数
}

CIntDwt97::~CIntDwt97()
{
	delete []h;
	delete []h1;
}

void CIntDwt97::Dwt2D(int buffer[], int height, int width, int lv)
{
	int i,j,k;
	int nheight,nwidth;
 
	
	for ( k=0; k<lv; k++){

       nheight=height>>k;
	   nwidth=width>>k;

	   int *pdata;
	   pdata = new int [nwidth];
		for (i=0; i<nheight; i++){
			for(j=0; j<nwidth; j++)
				*(pdata+j) = *(buffer+i*width+j);
			Dwt1D(pdata,nwidth);
			for(j=0; j<nwidth; j++)
				*(buffer+i*width+j) = *(pdata+j);
		}
		delete []pdata;

        int *p1data;
		p1data = new int [nheight];
		for(j=0; j<nwidth; j++){
			for (i=0; i<nheight; i++)
				*(p1data+i)=*(buffer+i*width+j);
			Dwt1D(p1data,nheight);
			for(i=0; i<nwidth; i++)
				*(buffer+i*width+j) = *(p1data+i);
		}
		delete []p1data;
	}
}

void CIntDwt97::IDwt2D(int buffer[], int height, int width, int lv)
{

	int i,j,k;
 	int nheight,nwidth;
 
   	for (k=0; k<lv; k++){
 		nheight=height>>(lv-k-1);
 		nwidth=width>>(lv-k-1);
 
		int *pdata;
		pdata = new int [nheight];
 		for(j=0; j<nwidth; j++){
 			for (i=0; i<nheight; i++)
 				*(pdata+i) = *(buffer+i*width+j);
 			IDwt1D(pdata,nheight);
 			for(i=0; i<nwidth; i++)
 				*(buffer+i*width+j) = *(pdata+i);
 		}
		delete []pdata;

		int *p1data;
		p1data  = new int [nwidth]; 		
 		for (i=0; i<nheight; i++){
 			for(j=0; j<nwidth; j++)
 				*(p1data+j)=*(buffer+i*width+j);
 			IDwt1D(p1data,nwidth);
 			for(j=0; j<nwidth; j++)
 				*(buffer+i*width+j) = *(p1data+j);
 		}
		delete []p1data;
 	}
}

void CIntDwt97::Dwt1D(int *buffer, int buflen)
{
	int i;
	int itemp;

	int *d, *s, *p;
	p = new int [buflen+2*shift];
	d = new int [(buflen>>1)+shift];
	s = new int [(buflen>>1)+shift];

	/*    使用周期延拓方式	/*
	for (i=0; i<shift; i++){
		*(p+i) = *(buffer+buflen-shift+i);
		*(p+buflen+shift+i) = *(buffer+i);
	}
	for (i=0; i<buflen; i++) *(p+i+shift) = *(buffer+i);
	/* 数据延拓结束  计算p1 = d1          p2 = s1 */

   /*  对称延拓 */
	for (i=0; i<shift; i++){
		*(p+i) = *(buffer+shift-i);
		*(p+buflen+shift+i) = *(buffer+buflen-i-2);
	}
	for (i=0; i<buflen; i++) *(p+i+shift) = *(buffer+i);	
	/* 数据延拓结束  计算p1 = d1          p2 = s1 */

	for (i=0; i<(buflen>>1)+shift; i++){
		itemp = i-(shift>>1);
		*(d+i) = *(p+shift+2*itemp+1) + int(floor(p_alfa*( *(p+shift+2*itemp)+*(p+shift+2*itemp+2))+0.5));
	}

	for (i=0; i<(buflen>>1)+shift-1; i++){
		itemp = i+1-(shift>>1);
		*(s+i+1) = *(p+shift+2*itemp) + int(floor( p_beta*( *(d+i+1)+*(d+i+1-1) )+0.5));
	}
   /*  d1 s1 计算结束 下计算 d2  */
	for (i=0; i<(buflen>>1)+shift-1; i++){
		//itemp = i-(shift>>1);
		*(p+shift+(buflen>>1)+i) =  *(d+i) + int(floor(p_gama*( *(s+i)+*(s+i+1))+0.5)); 
	}
	/*  s2 */
    for (i=0; i<(buflen>>1)+shift-1; i++){
	//	itemp = i-(shift>>1);
		*(p+i+1) = *(s+i+1) + int(floor(p_delta*(*(p+(buflen>>1)+shift+i+1) + *(p+(buflen>>1)+shift+i+1-1) )+0.5));
	}

	/* d3 */
	for (i=0; i<(buflen>>1); i++)
		*(d+i) = *(p+i+(buflen>>1)+shift+(shift>>1)) + int(floor(p_kesa*(1-p_kesa)*(*(p+i+(shift>>1)))+0.5));
	/* s3 */
	for (i=0; i<(buflen>>1); i++)
		*(s+i) = *(p+i+(shift>>1)) -int(floor(1/p_kesa*(*(d+i))+0.5));
    /* d4 */
	for (i=0; i<(buflen>>1); i++)
		*(buffer+i+(buflen>>1)) = *(d+i) +int(floor((p_kesa-1)*(*(s+i))+0.5));

	/* s4 */
	for (i=0; i<(buflen>>1); i++)
		*(buffer+i) = *(s+i) + *(buffer+i+(buflen>>1));

/*	for (i=0; i<(buflen>>1); i++){
   *(buffer+i) = p_kesa*(*(p+i+(shift>>1)));  
 	 *(buffer+i+(buflen>>1)) = 1/p_kesa*(*(p+i+(buflen>>1)+shift+(shift>>1)));	
	}
*/
    delete []p;
	delete []d;
	delete []s;
    
}

void CIntDwt97::IDwt1D(int *buffer, int buflen)
{
	int i;
    int *p1, *p2, *s, *d;
	p1 = new int [(buflen>>1)+shift];
	p2 = new int [(buflen>>1)+shift];
	s  = new int [(buflen>>1)+shift];
	d  = new int [(buflen>>1)+shift];

  
   /* s3 */
   for (i=0; i<(buflen>>1); i++)
	   *(s+i) = *(buffer+i) - *(buffer+i+(buflen>>1));
   /* d3 */
   for (i=0; i<(buflen>>1); i++)
	   *(d+i) =*(buffer+i+(buflen>>1))- int(floor((p_kesa-1)*(*(s+i))+0.5));
   /*  s2 */
   for (i=0; i<(buflen>>1); i++)
	   *(buffer+i) = *(s+i) + int(floor((1/p_kesa)*(*(d+i))+0.5));
   /*  d2 */
   for (i=0; i<(buflen>>1); i++)
	   *(buffer+i+(buflen>>1)) = *(d+i) - int(floor(p_kesa*(1-p_kesa)*(*(buffer+i))+0.5));


   /*    延拓方式为周期延拓     /*
	for (i=0; i<(shift>>1); i++){
		p1[i] = buffer[i+(buflen>>1)-(shift>>1)];
		p1[i+(shift>>1)+(buflen>>1)] = buffer[i];
		p2[i] = buffer[i+buflen-(shift>>1)];
		p2[i+(shift>>1)+(buflen>>1)] = buffer[i+(buflen>>1)];
	}
	for (i=0; i<(buflen>>1); i++){
		p1[i+(shift>>1)] = buffer[i];
		p2[i+(shift>>1)] = buffer[i+(buflen>>1)];
	}

	/*  延拓结束 */	

     	/*  对称延拓	*/
   	for (i=0; i<(shift>>1); i++){
		p1[i] = buffer[(shift>>1)-i];
		p1[i+(shift>>1)+(buflen>>1)] = buffer[(buflen>>1)-i-1];
		p2[i] = buffer[(buflen>>1)+(shift>>1)-i-1];
		p2[i+(shift>>1)+(buflen>>1)] = buffer[buflen-i-2];
	}
	for (i=0; i<(buflen>>1); i++){
		p1[i+(shift>>1)] = buffer[i];
		p2[i+(shift>>1)] = buffer[i+(buflen>>1)];
	}
/*   对称延拓结束  */


   /*  s1  */
   for (i=0; i<(buflen>>1)+shift-1; i++)
	   *(s+i+1) = *(p1+i+1) - int(floor(p_delta*( *(p2+i+1) + *(p2+i+1-1))+0.5));
   /*  d1 */
   for (i=0; i<(buflen>>1)+shift-1; i++)
	   *(d+i) = *(p2+i) - int(floor(p_gama*(*(s+i) + *(s+i+1))+0.5));
    /* p1 = s0 */
   for (i=0; i<(buflen>>1)+shift-1; i++)
	   *(p1+i+1) = *(s+i+1) - int(floor(p_beta*( *(d+i+1)+*(d+i+1-1))+0.5));
   /*  p2 = d0 */
   for (i=0; i<(buflen>>1)+shift-1; i++)
	   *(p2+i) = *(d+i) - int(floor(p_alfa*(*(p1+i) + *(p1+i+1))+0.5));

   for (i=0; i<(buflen>>1); i++){
	   *(buffer+2*i) = *(p1+i+(shift>>1)) ;
	   *(buffer+2*i+1) = *(p2+i+(shift>>1));
   }


	delete []p1;
	delete []p2;
	delete []s;
	delete []d;
}

void CIntDwt97::filterset(double t)
{
	h = new double [9];
	h1 = new double [7];
	*h = (8*t*t*t-6*t*t+3*t)/(1+2*t)*(1/32.)*sqrt(2);
	*(h+1) = (-16*t*t*t+20*t*t-12*t+3)/(1+2*t)*(1/32.)*sqrt(2);
	*(h+2) = (2*t-3)/(1+2*t)*(1/8.)*sqrt(2);
	*(h+3) = (16*t*t*t-20*t*t+28*t+5)/(1+2*t)*(1/32.)*sqrt(2);
	*(h+4) = (-8*t*t*t+6*t*t+5*t+20)/(1+2*t)*(1/16.)*sqrt(2);
	*(h+5) = *(h+3);
	*(h+6) = *(h+2);
	*(h+7) = *(h+1);
	*(h+8) = *(h+0);

	double r0, r1, s0, t0;
	r0 = (*(h+4))-2*(*h)*(*(h+3))/(*(h+1));
	r1 = (*(h+2))-(*h)-(*h)*(*(h+3))/(*(h+1));
	s0 = (*(h+3))-(*(h+1))-(*(h+1))*r0/r1;
	t0 = r0-2*r1;

	p_alfa = (*h)/(*(h+1));
	p_beta = (*(h+1))/r1;
	p_gama = r1/s0;
	p_delta = s0/t0;
	p_kesa = t0;


}