#include "stdafx.h"
#include "math.h"
#include<fstream>
#include<iostream>
using namespace std;


double Wavelet_Inter_Start(int i,int j)
{ 
	
	double start;
	double length;
         if(i==0){start=0;}
         else{
             double p=pow(2.0,i-1);
               length=1.0/p;
              double temp;
                  if (j%2==0){temp=j/2.0;}
                  else {temp=(j-1)/2.0;}
				  start=temp/p;
              }
 return start;
}
double Wavelet_Inter_End(int i, int j)
{ 
	double p,endpoint;
          if(i==0){endpoint=1;}
          else{
              p=pow(2.0,i-1);
              double length=1.0/p;
              endpoint=Wavelet_Inter_Start(i,j)+length;}
    return endpoint;
}




double Wavelet(int i, int j, double t)
{   
	
	double result=0;
	double temp=0;
	//父小波w(00),w(01)的基函数
	
	if (i==0&&j==0){result=1;}
	if(i==0&&j==1){result=sqrt(3.0)*(2.0*t-1);}
	//母小波w(10),w(11)的基函数 
    if(i==1 && j==0 && t>=0.0 && t<=1.0/2) 
	result=1.0-6.0*t; 
    if(i==1 && j==0 && t>1.0/2 && t<=1.0) 
	result=5.0-6.0*t;
    if(i==1 && j==1 && t>=0.0 && t<=1.0/2.0) 
	result=sqrt(3.0)*(1.0-4.0*t); 
    if(i==1 && j==1 && t>1.0/2.0 && t<=1.0) 
	result=sqrt(3.0)*(4.0*t-3.0);
	//生成的小波空间基函数
	   double t1=Wavelet_Inter_Start(i,j);
	   double t2=(Wavelet_Inter_Start(i,j)+Wavelet_Inter_End(i,j))/2;
	   double t3=Wavelet_Inter_End(i,j);
    int p=j;
    if(i>=2)
	{
	       for(int k=i;k>1;k--)     //记录右移的常数项
		   {

	             if (p>=((int)pow(2.0,k-1)))
				 {
	                  p=p%((int)pow(2.0,k-1));
	                  temp+=pow(2.0,k-2); 
	              }
				 else
				 {
					
					 temp+=0;
				 }

		   }

	   if(t>=t1&&t<=t2)
	   {
		   if(p==0){result=sqrt(pow(2.0,i-1))*(1-6*(pow(2.0,i-1)*t-temp));}
		   if(p==1){result=sqrt(pow(2.0,i-1))*sqrt(3.0)*(1-4*(pow(2.0,i-1)*t-temp));}
	    }
	   else if(t<=t3&&t>t2)
	   {
		   if(p==0){result=sqrt(pow(2.0,i-1))*(5-6*(pow(2.0,i-1)*t-temp));}
		   if(p==1){result=sqrt(pow(2.0,i-1))*sqrt(3.0)*(4*(pow(2.0,i-1)*t-temp)-3);}
	    }
	   else
		   result=0;

	}

	return result;
}


//二重积分高斯五点法
double guass_2D(int l,int m,int p,int q)   //w(a,b),w(c,d)
{   
	    
	    double x[]={-0.9061798459,-0.5384693101,0.0000000000,0.5384693101,0.9061798459};  //高斯点
	    double w[]={0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851 };  //系数
	    double sum=0;
		double a=Wavelet_Inter_Start(l,m);//a,b,c,d为积分上下限
        double b=Wavelet_Inter_End(l,m);
	
	    double c=Wavelet_Inter_Start(p,q);
	    double d=Wavelet_Inter_End(p,q);
	
		double aa=max(a,c);  //积分区间为【aa,bb】
		double bb=min(b,d);
		double cc=(aa+bb)/2;
		
		for(int i=0;i<5;i++)
		{ double s=x[i]*(bb-cc)/2.0+(bb+cc)/2.0, ss=x[i]*(cc-aa)/2.0+(cc+aa)/2.0;
			for (int j=0;j<5;j++)
		     { 
				 double t=x[j]*(bb-cc)/2.0+(bb+cc)/2.0, tt=x[j]*(cc-aa)/2.0+(cc+aa)/2.0;
			     sum+=w[i]*w[j]*(bb-cc)*(cc-aa)/4.0*(Wavelet(l,m,t)*Wavelet(p,q,s)+Wavelet(l,m,t)*Wavelet(p,q,ss)+Wavelet(l,m,tt)*Wavelet(p,q,s)+Wavelet(l,m,tt)*Wavelet(p,q,ss));
			  }
	    }
			return sum;

}

 int _tmain(int argc, _TCHAR* argv[])
{
	int n=4;  //小波的层数
	int nn=(int)pow(2.0,n+1);
long double **A;//A为系数矩阵
	A=new long double *[nn];
	for(int i=0;i<nn;i++)
		A[i]=new long double[nn];
	
	for(int i=0;i<nn;i++)
	{
		for(int j=0;i<nn;i++)
	        A[i]=0;}
	}

//生成矩阵
 
 {
	 for  ( int i=0;i<nn;i++;)
	for(int j=0;j<(i==0?2:pow(2.0,i));j++;)
		{
			int k;
			if(i==0) k=j;else k=(int)pow(2.0,i)+j;
			for (int p=0;p<nn;p++ )
				for(int q=0;q<(p==0?2:pow(2.0,p));q++;)
				{
					int m;
					if(p==0) m=q;else m=(int)pow(2.0,p)+q;
					
					A[m][k]=guass_2D(i,j,p,q);
					
				}
		}
		
 }
ofstream fout; //矩阵的输出
{
for int i, int j;

	fout.open("f:/output_AA.m");
	if(!fout)
		cout<<"不能打开文件"<<endl;
	else
	{
		fout<<"aa=[";
		for(int i=0;i<nn;i++)
		{
			for(int j=0;j<nn;j++)		
				fout<<A[i][j]<<"  ";
			fout<<endl;
		}
		fout<<"];"<<endl;
		
	}
}