wavelet.cpp

信号消燥

					}
					if (mod!=0)
					{
				    vector X_temp2;
				    X_temp2.GetLeft(X,mod);
			        Y.Link(X_temp2);
					}

				}
				else
				{
                     wextend(X,Y,"ppd",l,'b');
				}
			 }
			 break;
	case 'l':		     
		     if ( mode=="zpd")
		     {
				 Y.Set(n+l);
				 for (i=1;i<=n;i++)
				 {
					 Y[l+i]=X[i];
				 }
		     }
			 else if (mode=="sym")
			 {
				 	int a=l/n,mod=l%n;
				 if (a%2==0)
				 {   
					 if (mod!=0)
					 {
                     Y.GetLeft(X,mod);
					 Y.Invers();
					 }

					 Y.Link(X);
					 for (i=2;i<=a+1;i++)
					 {
						 X.Invers();
						 Y.Link(X);
					 }
				 }
				 else if (a%2==1)
				 {
					 if (mod!=0)
					 {
						 Y.GetRight(X,mod);
					 }
					 
					 for (i=1;i<=a+1;i++)
					 {
						 X.Invers();
						 Y.Link(X);
					 }
				 }
			 }
			 else if (mode=="ppd")
			 {
				 int a=l/n,mod=l%n;
				 if (mod!=0)
				 {
					  Y.GetRight(X,mod);
				 }
				
				 for (i=1;i<=a+1;i++)
				 {
					 Y.Link(X);
				 }
			 }
			 else if (mode=="per")
			 {
				if (n%2==1)
				{
					X.Add(X[n]);
					int nn=X.Size();
                    int a=l/nn,mod=l%nn;
					if (mod!=0)
					{
                     Y.GetRight(X,mod);
					}
					
				    for (i=1;i<=a+1;i++)
					{
					   Y.Link(X);
					}
				}
				else
				{
                     wextend(X,Y,"ppd",l,'l');
				}
			 }
             break;
	case 'r':
		    if ( mode=="zpd")
		     {
				 Y.Set(n+l);
				 for (i=1;i<=n;i++)
				 {
					 Y[i]=X[i];
				 }
		     }
			 else if (mode=="sym")
			 {
				 	int a=l/n,mod=l%n;
				 if (a%2==0)
				 {
                     Y=X;
					 for (i=2;i<=a+1;i++)
					 {
						 X.Invers();
						 Y.Link(X);
					 }
					 if (mod!=0)
					 {
						 vector X_temp;
					 X_temp.GetRight(X,mod);
					 X_temp.Invers();
					 Y.Link(X_temp);
					 }
					 
				 }
				 else if (a%2==1)
				 {
					 Y=X;
					 for (i=2;i<=a+1;i++)
					 {
						 X.Invers();
						 Y.Link(X);
					 }
					 if (mod!=0)
					 {
						 X.Invers();
					 vector X_temp1;
					 X_temp1.GetLeft(X,mod);
					 Y.Link(X_temp1);
					 }
					 
				 }
			 }
			 else if (mode=="ppd")
			 {
				 int a=l/n,mod=l%n;
				 Y=X;
				 for (i=2;i<=a+1;i++)
				 {
					 Y.Link(X);
				 }
				 if (mod!=0)
				 {
                 vector X_temp2;
				 X_temp2.GetLeft(X,mod);
			     Y.Link(X_temp2);
				 }
				 
			 }
			 else if (mode=="per")
			 {
				if (n%2==1)
				{
					X.Add(X[n]);
					int nn=X.Size();
                    int a=l/nn,mod=l%nn;
				    Y=X;
				    for (i=2;i<=a+1;i++)
					{
					   Y.Link(X);
					}
					if (mod!=0)
					{
                    vector X_temp2;
				    X_temp2.GetLeft(X,mod);
			        Y.Link(X_temp2);
					}
				    
				}
				else
				{
                     wextend(X,Y,"ppd",l,'r');
				}
			 }
			 break;
	}
}

void Cwavelet::dwt(vector X,CString wavetype,int n,CString mode)
{
	//单尺度一维离散小波变换;
	//wavetype为小波基函数的名称；n为基函数的长度，如db2，wavetype="db",n=2;
	//mode为延拓模式，其具体的取值见wextend函数。一般取"per";
	//编写时间:2006.09.20

	//计算滤波器系数
	if (wavetype=="db")
	{
			dbaux(n);
	        orthfilt(db_w);
	}

	//判断延拓模式，并根据此计算扩展长度等；
	int lf=Lo_D.Size(),lx=X.Size();
	int LenEXT,LenKEPT;//lenEXT为延拓的长度，LenKEPT为中间截取的长度；
	if (mode=="per")
	{
		LenEXT=lf/2;
		if (lx%2==0)
		{
		LenKEPT=lx;
		}
	    else
		{
		LenKEPT=2*(lx/2+1);
		}
	}
	else
	{
		LenEXT=lf-1;
		LenKEPT=lx+lf-1;
	}
	
   vector y;
   wextend(X,y,mode,LenEXT,'b');
   vector cc1,cc2;
   cc1.LConvolute(y,Lo_D);  //卷积计算；
   int dd=((cc1.Size()-LenKEPT)/2)+1;
   cc2.GetMid(cc1,dd,LenKEPT);
   dyaddown(cc2,cA,0);
   // cA.Print(16);
   cc1.LConvolute(y,Hi_D);
   dd=((cc1.Size()-LenKEPT)/2)+1;
   cc2.GetMid(cc1,dd,LenKEPT);
   dyaddown(cc2,cD,0);
   //cD.Print(16);

}


void Cwavelet::wnoise(int num, int len, double mean,double variance,matrix & data)
{
	//此函数产生仿真数据(包括原始信号,噪声信号,以及含噪信号);
	//num表示产生数据类型,num=1 产生Blocks信号,num=2产生Heavy sine信号;
	//注:如果想增加信号的类型,只需增加num即可;
	//len表示要产生的数据的长度;
	//mean,variance分别表示要产生的噪声的均值和方差;
	//data表示最后输出矩阵;
	//编写时间2006.09.27

    vector t(11),h(11);
    t.Init(0.1,0.13,0.15,0.23,0.25,0.40,0.44,0.65,0.76,0.78,0.81);
    h.Init(4.0,-5.0,3.0,-4.0,5.0,-4.2,2.1,4.3,-3.1,5.1,-4.2);
    vector signal; //用于存储信号;

	if (num==1)  //Blocks函数;
    {
		 vector tt,x(len);
	    tt.Cut(0.0,1.0,len-1);
	    for (int i=1;i<=11;i++)
		{
		   vector temp(len);
		   temp.Const(t[i]);
		   temp-=tt;
		   for (int j=1;j<=len;j++)
		   {
			    if (temp[j]>0.0)
				{
				temp[j]=0.0;
				}
			    else if (temp[j]<0.0)
				{
				temp[j]=1.0;
				}
			    else
				{
                temp[j]=0.5;
				}
		   }
		temp*=h[i];
		x+=temp;
		}
		signal=x;
    }
	else if (num==2)   //Heavy sine函数;
	{
		const pi=3.14159265358979;
 		vector x;
	    x.Cut(0.0,1.0,len-1);

		vector temp1(len),temp2(len);
		int j;
		temp1.Const(-0.3);
		temp2.Const(0.72);
		temp1+=x;
		temp2-=x;
		for ( j=1;j<=len;j++)
		   {
			    if (temp1[j]>0.0)
				{
				temp1[j]=1.0;
				}
			    else if (temp1[j]<0.0)
				{
				temp1[j]=-1.0;
				}
			    else
				{
                temp1[j]=0.0;
				}

				if (temp2[j]>0.0)
				{
				temp2[j]=1.0;
				}
			    else if (temp2[j]<0.0)
				{
				temp2[j]=-1.0;
				}
			    else
				{
                temp2[j]=0.0;
				}
		   }
 		temp1+=temp2;
		temp2.Mult(x,4.0*3.14159265358979);
		
			for ( j=1;j<=len;j++)
		   {
			   temp2[j]=4.0*sin(temp2[j]);
		   }
		
		signal=temp2-temp1;	
	}
	vector noise(len);//用于存储噪声的向量;
	noise.RandomGauss(mean,variance);
	vector signal_noise;
	signal_noise=signal+noise;
    matrix temp;
			temp.Set(len,1);
			temp.SetColum(signal,1);
            temp.AddColum(noise);
			temp.AddColum(signal_noise);
			data=temp;
			temp.Destroy();
	

	   
}


void Cwavelet::idwt(vector cA, vector cD, vector Lo_R, vector Hi_R, CString mode, int L)
{
    //单尺度一维离散小波逆变换;
	//L为对信号中间长度为L的部分进行重构,默认为0,即对全部信号进行重构;
	vector temp1,temp2;
	temp1=upsaconv(cA, Lo_R, L, mode);
	temp2=upsaconv(cD, Hi_R, L, mode);
	Re_Data.Plus(temp1,temp2);
}

vector Cwavelet::upsaconv(vector x, vector f, int s, CString mode)
{
	int lx=x.Size();
	int lf=f.Size();
	if (mode=="per")
	{
		if (s==0)
		{
			s=2*lx;
		}
		vector temp,temp1;
		dyadup(x,temp,0,1);
		wextend(temp,temp1,"per",lf/2,'b');
		temp.LConvolute(temp1,f);
		temp1.GetMid(temp,lf,2*lx);
		temp.GetMid(temp1,1,s);
		return temp;
	}
	else
	{
        if (s==0)
		{
			s=2*lx-lf+2;
		}
		vector temp,temp1;
		dyadup(x,temp,0,0);
		temp1.LConvolute(temp,f);
		int dd=((temp1.Size()-s)/2)+1;
        temp.GetMid(temp1,dd,s);
		return temp;
	}

}

void Cwavelet::wavedec(vector X, int N, CString wavetype, int n, CString mode)
{
	//多尺度一维小波分解；
	//X表示数据，N表示分解的层数；
	//wavetype ,n,mode可参考dwt()；

	//需要对分解数据进行预处理；
	vector Y;
	preData(X,Y,N);


    vector temp;
	wavedec_L.Set(1);
	 wavedec_L[1]=Y.Size();
	 dwt(Y,wavetype,n,mode);
	 wavedec_C=cD;
	 wavedec_L.Ins(cD.Size(),1);
	 for (int i=2;i<=N;i++)
	 {
		 Y=cA;
		 dwt(Y,wavetype,n,mode);
		 temp=wavedec_C;
		 wavedec_C.Link(cD,temp);
		 wavedec_L.Ins(cD.Size(),1);
	 }
	 temp=wavedec_C;
	 wavedec_C.Link(cA,temp);
	 wavedec_L.Ins(cA.Size(),1);
}




void Cwavelet::appcoef(vector &C, vector &L, int n)
{
    int Len_layer=L.Size()-2;
	int len=0,i=1;
	vector temp;
	Re_Data.GetLeft(C,L[1]);
	for (int j=Len_layer-1;j>=n;j--)
	{
		len+=L[i];
        temp.GetMid(C,len+1,L[i+1]);
		idwt(Re_Data,temp,Lo_R,Hi_R,"per");
		i++;
	}

}

void Cwavelet::waverec(vector &C, vector &L)
{
	appcoef(C,L,0);
	vector temp=Re_Data;
    Re_Data.GetMid(temp,infoData.len_wectend+1,infoData.len_data);

}

double Cwavelet::thselect(vector Coef,vector Coef1)
{
	int len=Coef.Size();//系数的长度；
	double median,s,temp_len=0.0;
	double thr1,thr2;
	vector y(len),y1(len),y2(len),r(len);


//	Coef.Print();
	for (int i=1;i<=len;i++)
	{
		y[i]=fabs(Coef[i]);
	}
	y.SmallToBig();

	if (len%2==0)
	{
		median=(y[len/2]+y[len/2+1])/2.0;
	}
	else
	{
		median=y[len/2+1];
	}
	s=median/0.6745;

	Coef1*=(1.0/s);

	for (i=1;i<=len;i++)
	{
		y[i]=fabs(Coef1[i]);
	}
	y.SmallToBig();
  for (i=1;i<=len;i++)
	{
	  y1[i]=y[i]*y[i];
	  temp_len+=y1[i];
	  y2[i]=temp_len;
	  r[i]=(len-2*i+y2[i]+(len-i)*y1[i])/double(len);
	}
  thr1=sqrt(r.Min());
  thr2=sqrt(2*log(double(len)));
   
  double d=ED(Coef,Coef)*ED(Coef,Coef);
  double eta=(d-len)/double(len);
  double crit=pow(((log(double(len))/log(2.0))),1.5)/(sqrt(len));
  if (eta<crit)
  {
	return thr2*s;
  }
  else
  {
	  return (thr1<thr2?thr1:thr2)*s;
  }
  
}

void Cwavelet::wdencmp(CString mode)
{
    //mode表示选取的滤值方式，soft表示软阈值，hard表示硬阈值；
	int Len_layer=wavedec_L.Size()-2;
	int len=0,i=1,j,k;
	int len_temp;
	vector temp,temp1;
	double th;  //阈值；
	temp1.GetRight(wavedec_C,wavedec_L[Len_layer+1]);
	wavedec_denoisde_C.GetLeft(wavedec_C,wavedec_L[1]);
	for (j=1;j<=Len_layer;j++)
	{
		len+=wavedec_L[i];
        temp.GetMid(wavedec_C,len+1,wavedec_L[i+1]);
		th=thselect(temp,temp1);
	    len_temp=temp.Size();
		for (k=1;k<=len_temp;k++)
		{
			if (mode=="soft")
			{
				if (fabs(temp[k])<th)
				{
					temp[k]=0.0;
				}
				else
				{
					if (temp[k]>=0)
					{
						temp[k]-=th;
					}
					else
					{
						temp[k]+=th;
					}
				}
			}
			else
			{
               if (fabs(temp[k])<th)
					temp[k]=0.0;
			}
			
		}
		wavedec_denoisde_C.Link(temp);
		i++;
	}
	waverec(wavedec_denoisde_C,wavedec_L);
}

void Cwavelet::preData(vector X,vector &Y,int n)
{
  int no_wextend=pow(2,n);
  infoData.len_data=X.Size();
  int len=infoData.len_data;
  if (len%2!=0)
  {   //奇数时应将长度加1，因为"per"延拓在奇数的情况下，会自动延拓最后一位后再进行周期延拓
	  //保证数据长度为偶数；
      len+=1;
  }
  infoData.len_wectend=(no_wextend-len%no_wextend)/2; //表示的是两边延拓的数据长度；
  wextend(X,Y,"per",infoData.len_wectend,'b');
}