mallat program.txt

基于Mallat的二维Haar小波分解与重构算法的C语言实现最大分解级数:5级*/

*用途:基于Mallat的二维Haar小波分解与重构算法的C语言实现最大分解级数:5级*/
#include<iostream.h>      
#include<iomanip.h>
//float data[32][32]; //手工录入时启用此声名
//用作试验的随意数据,在手工录入时注释掉下面的声名和初始化，并启用数据录入函数调用
float data[8][8]={{1,2,3,4,5,6,7,8},{4.22,5.33,3.44,8.66,12,54,6,7},{6,7,1,2,3,2,1,7},{5,2,4,3,5,9,6,3},
{3,2,5,23,5,6,75,3},{2,7,43,54,7,2,6,78},{3.5,4.125,2.625,4,237,6,94,43},{3,6,43,7,84,34,34,33}}; //二维数组的录入函数，手工录入时启用该函数的调用
void Data_Input(int N)
{   cout<<"请输入"<<N<<" X "<<N<<"个象素的值\n";
for(int i=0;i<N;i++)
       {
       for(int j=0;j<N;j++)
              {   cout<<"data["<<i<<"]["<<j<<"]=";  cin>>data[i][j];   }         
       }
}
//Mallat 分解函数
void Mallat_De(int Dim)
{   //分解滤波器系数
float LO_D[]={0.70710678118655,0.70710678118655};
float HI_D[]={-0.70710678118655,0.70710678118655};
float a[32][32],d[32][32];
float a1[32][32],d1[32][32];
float *Atmp=0;
float tmp[34],tmp2[34];
float L[34][34],H[34][34];
float sum=0;
for(int i=0;i<32;i++)   {tmp[i]=0;tmp2[i]=0;}     //暂存区初始化
//进行行变换，求取第一级变换结果L和H
for (i=0;i<Dim;i++)
       {
       for(int j=0;j<Dim;j++)  {       tmp[j]=data[i][j];}   //操作数暂存
       //行变换，求取L
       for(j=0;j<Dim+1;j++)
              {            
              sum=tmp[Dim-j]*LO_D[0]+tmp[Dim-j-1]*LO_D[1];    //卷积
              tmp2[Dim-j]=sum;
              }
       for(j=0;j<Dim;j++)
              {  a[i][j]=tmp2[2*j+1];   }  //抽取(下抽样)
       //行变换，求取H
       for(j=0;j<Dim+1;j++)
              {            
              sum=tmp[Dim-j]*HI_D[0]+tmp[Dim-j-1]*HI_D[1];    //卷积
              tmp2[Dim-j]=sum;
              }
       for(j=0;j<Dim;j++)
              { d[i][j]=tmp2[2*j+1];  }     //抽取(下抽样)       
       }
//将a[][],d[][]结果转置为a1[][],d1[][]。
for(i=0;i<Dim;i++)
       {
       for(int j=0;j<Dim/2;j++)
              {   a1[j][i]=a[i][j];
                  d1[j][i]=d[i][j];
              }
       }
//进行列变换，求得LL
for (i=0;i<Dim/2;i++)
       {
       for(int j=0;j<Dim;j++){       tmp[j]=a1[i][j];}   //操作数暂存
       //列变换，求取LL
       for(j=0;j<Dim+1;j++)
              {            
              sum=tmp[Dim-j]*LO_D[0]+tmp[Dim-j-1]*LO_D[1];    //卷积
              tmp2[Dim-j]=sum;
              }
       for(j=0;j<Dim/2;j++)
              {   data[i][j]=tmp2[2*j+1];  }      //抽取(下抽样),保存分解结果
       }
//进行列变换，求得HL
for (i=0;i<Dim/2;i++)
       {
       for(int j=0;j<Dim;j++){       tmp[j]=a1[i][j];}   //操作数暂存
       //列变换，求取LL
       for(j=0;j<Dim+1;j++)
              {            
              sum=tmp[Dim-j]*HI_D[0]+tmp[Dim-j-1]*HI_D[1];    //卷积
              tmp2[Dim-j]=sum;
              }
       for(j=0;j<Dim/2;j++)
              {   data[i][Dim/2+j]=tmp2[2*j+1];   }    //抽取(下抽样),保存分解结果
       }
//进行列变换，求得LH
for (i=0;i<Dim/2;i++)
       {
       for(int j=0;j<Dim;j++){       tmp[j]=d1[i][j];}   //操作数暂存
       //列变换，求取LL
       for(j=0;j<Dim+1;j++)
              {            
              sum=tmp[Dim-j]*LO_D[0]+tmp[Dim-j-1]*LO_D[1];    //卷积
              tmp2[Dim-j]=sum;
              }
       for(j=0;j<Dim/2;j++)
              {  data[Dim/2+i][j]=tmp2[2*j+1];   }    //抽取(下抽样),保存分解结果
       }
//进行列变换，求得HH
for (i=0;i<Dim/2;i++)
       {
       for(int j=0;j<Dim;j++){       tmp[j]=d1[i][j];}   //操作数暂存
       //列变换，求取LL
       for(j=0;j<Dim+1;j++)
              {            
              sum=tmp[Dim-j]*HI_D[0]+tmp[Dim-j-1]*HI_D[1];    //卷积
              tmp2[Dim-j]=sum;
              }
       for(j=0;j<Dim/2;j++)
              {  data[Dim/2+i][Dim/2+j]=tmp2[2*j+1];  }     //抽取(下抽样),保存分解结果       
       }
}//end Mallat_De
//Mallat  重构函数
void Mallat_Re(int Dim)
{
//重构滤波器系数
float LO_R[]={0.70710678118655,0.70710678118655};
float HI_R[]={0.70710678118655,-0.70710678118655};
float tmp[34];
float sum=0;
float aa[34],ad[34],da[34],dd[34];
float a[32][32],d[32][32];
float a0[32][32],d0[32][32];
//由LL,HL重构L
for (int i=0;i<Dim/2;i++)
       {    for(int j=0;j<Dim;j++)  {aa[j]=0;}            //清空暂存区
            for(j=0;j<Dim/2;j++)   {aa[2*j]=data[i][j];}   //LL插值(上抽样) 
            for(j=0;j<Dim;j++)
                 { if (Dim-2-j<0)    {sum=aa[Dim-1-j]*LO_R[0];}
                   else     {sum=aa[Dim-1-j]*LO_R[0]+aa[Dim-2-j]*LO_R[1]; }   //卷积
                   tmp[Dim-1-j]=sum;       
                 }
       for(j=0;j<Dim;j++)   {  aa[j]=tmp[j];  }               //保存恢复结果
//---
       for(j=0;j<Dim;j++)   {  da[j]=0;      }            //清空暂存区
       for(j=0;j<Dim/2;j++){da[2*j]=data[i][Dim/2+j];}   //HL插值(上抽样)
       for(j=0;j<Dim;j++)
              {     
if (Dim-2-j<0)   {sum=da[Dim-1-j]*HI_R[0];}
                    else   {sum=da[Dim-1-j]*HI_R[0]+da[Dim-j-2]*HI_R[1];}   //卷积
                    tmp[Dim-1-j]=sum;       
}
       for(j=0;j<Dim;j++)   {  da[j]=tmp[j];  }               //保存恢复结果
// 合成第一级低频系数,并转置
       for(j=0;j<Dim;j++)
              {  a[j][i]=aa[j]+da[j];  }
       }
//由LH,HH重构H
for (i=0;i<Dim/2;i++)
       {
       for(int j=0;j<Dim;j++)    {ad[j]=0;}            //清空暂存区
       for(j=0;j<Dim/2;j++)     {ad[2*j]=data[Dim/2+i][j];}   //LL插值(上抽样)
       float sum=0;
       for(j=0;j<Dim;j++)
              {
              if (Dim-2-j<0)  {sum=ad[Dim-1-j]*LO_R[0];}
              else   {sum=ad[Dim-1-j]*LO_R[0]+ad[Dim-2-j]*LO_R[1]; }   //卷积
              tmp[Dim-1-j]=sum;       
              }
       for(j=0;j<Dim;j++) {ad[j]=tmp[j];}               //保存恢复结果
//---
       for(j=0;j<Dim;j++){dd[j]=0;}            //清空暂存区
       for(j=0;j<Dim/2;j++){dd[2*j]=data[Dim/2+i][Dim/2+j];}   //HL插值(上抽样)
       for(j=0;j<Dim;j++)
              {
              if (Dim-2-j<0)   {sum=dd[Dim-1-j]*HI_R[0];}
              else     {sum=dd[Dim-1-j]*HI_R[0]+dd[Dim-j-2]*HI_R[1];}   //卷积
              tmp[Dim-1-j]=sum;       
              }
       for(j=0;j<Dim;j++)    {dd[j]=tmp[j];}               //保存恢复结果
// 合成第一级高频系数,并转置
       for(j=0;j<Dim;j++)
              {   d[j][i]=ad[j]+dd[j];    }
       }
//合成最后信号
for (i=0;i<Dim;i++)
       {
       for (int j=0;j<Dim;j++)
              {     a0[i][j]=0;             //暂存区初始化
                    d0[i][j]=0;
              }     
       for(j=0;j<Dim/2;j++)
              {     a0[i][2*j]=a[i][j];   //L插值(上抽样)
                    d0[i][2*j]=d[i][j];
              }  
       for(j=0;j<Dim;j++)
              {
              if (Dim-2-j<0)   {sum=a0[i][Dim-1-j]*LO_R[0];}
              else    {sum=a0[i][Dim-1-j]*LO_R[0]+a0[i][Dim-2-j]*LO_R[1]; }   //卷积
              tmp[Dim-1-j]=sum;       
              }
       for(j=0;j<Dim;j++){a0[i][j]=tmp[j];}
       for(j=0;j<Dim;j++)
              {
              if (Dim-2-j<0)  {sum=d0[i][Dim-1-j]*HI_R[0];}
              else    {sum=d0[i][Dim-1-j]*HI_R[0]+d0[i][Dim-2-j]*HI_R[1]; }   //卷积
              tmp[Dim-1-j]=sum;       
              }
       for(j=0;j<Dim;j++){d0[i][j]=tmp[j];}
       for(j=0;j<Dim;j++)
               {    data[i][j]=a0[i][j]+d0[i][j];  }      //重构的最终结果       
       }
}//end Mallat_Re
//主函数
void main(void)
{
int Level=3;    //分解级数,注意与试验数据相匹配。手工录入时可不予考虑
int Dim;     //二维矩阵宽度
int n;
char ch;
Dim=1<<Level;   //手工录入数据时注释掉此语句
/*//手工录入数据时，请去掉这一段注释
cout<<"请输入分解的级数(不大于5 的整数)";
cout<<"(矩阵的长度将等于2 的N 次方)"<<endl;
cout<<"N=";
cin>>Level;
Dim=1<<Level;
Data_Input(Dim);  //输入原始数据,手工录入时启用此函数
//  */
//在屏幕上显示原始数据
cout<<"原始数据"<<endl;
for(int i=0;i<Dim;i++)
       {
       for(int j=0;j<Dim;j++)
              {   cout<<setw(7)<<data[i][j]<<"   ";  }
       cout<<endl;
       }
for (i=Level;i>0;i--)
       {    Dim=1<<i;
            Mallat_De(Dim);   //逐次调用Mallat分解算法函数，逐级进行分解
       }
//在屏幕上显示分解结果
Dim=1<<Level;
cout<<"分解结果"<<endl;
for (i=0;i<Dim;i++)
       {
       for (int j=0;j<Dim;j++)
              {   cout<< setw(7) <<data[i][j]<<"   ";   }
       cout<<endl;
       }
for (i=1;i<=Level;i++)
       {    Dim=1<<i;
            Mallat_Re(Dim);   //逐次调用Mallat重构算法函数，逐级进行重构
       }
//在屏幕上显示重构结果
Dim=1<<Level;
cout<<"重构结果"<<endl;
for (i=0;i<Dim;i++)
       {
       for (int j=0;j<Dim;j++)
              {  cout<< setw(7) <<data[i][j]<<"   ";   }
       cout<<endl;
       }
cout<<"显示完毕，请输入一个字符键，并回车退出!"<<endl;
cin>>ch;  //等待键盘输入,以让屏幕显示结果，便于观察屏幕输出
}//end main()