首页 > 代码库 > 图像处理之基础---一维小波变换,可多次分解

图像处理之基础---一维小波变换,可多次分解

1、题目:一维小波变换,可多次分解

 

2、原理:卷积核变为Daubechies正交小波基h[]和g[]的交替形式。增加了多次分解的功能。

 

3、代码:

 

 

[cpp] view plaincopy
 
  1. #include <stdio.h>  
  2. #include <stdlib.h>  
  3. #include <math.h>  
  4. #define LENGTH 4096//信号长度  
  5. /***************************************************************** 
  6. * 一维卷积函数 
  7. * 说明: 循环卷积,卷积结果的长度与输入信号的长度相同 
  8. * 输入参数: data[],输入信号; h[],Daubechies小波基低通滤波器系数; 
  9. *            g[],Daubechies小波基高通滤波器系数; cov[],卷积结果; 
  10. *            n,输入信号长度; m,卷积核长度. 
  11. * 李承宇, lichengyu2345@126.com 
  12. *  2010-08-22   
  13. *****************************************************************/  
  14. void Covlution(double data[], double h[], double g[], double cov[]  
  15.                , int n, int m)  
  16. {  
  17.     int i = 0;  
  18.     int j = 0;  
  19.     int k = 0;  
  20.   
  21.     //将cov[]清零  
  22.     for(i = 0; i < n; i++)  
  23.     {  
  24.         cov[i] = 0;  
  25.     }  
  26.   
  27.     //****************************************************  
  28.     //奇数行用h[]进行卷积  
  29.     //****************************************************  
  30.     //前m/2+1行  
  31.     i = 0;  
  32.     for(j = 0; j < m/2; j+=2, i+=2)  
  33.     {  
  34.         for(k = m/2-j; k < m; k++ )  
  35.         {  
  36.             cov[i] += data[k-(m/2-j)] * h[k];//k针对core[k]  
  37.         }  
  38.   
  39.         for(k = n-m/2+j; k < n; k++ )  
  40.         {  
  41.             cov[i] += data[k] * h[k-(n-m/2+j)];//k针对data[k]  
  42.         }  
  43.     }  
  44.   
  45.     //中间的n-m行  
  46.     for( ; i <= (n-m)+m/2; i+=2)  
  47.     {  
  48.         for( j = 0; j < m; j++)  
  49.         {  
  50.             cov[i] += data[i-m/2+j] * h[j];  
  51.         }  
  52.     }  
  53.   
  54.     //最后m/2-1行  
  55. //  i = ( (n - m) + m/2 + 1 )/2*2;//**********  
  56.     for(j = 1; j <= m/2; j+=2, i+=2)  
  57.     {  
  58.         for(k = 0; k < j; k++)  
  59.         {  
  60.             cov[i] += data[k] * h[m-j-k];//k针对data[k]  
  61.         }  
  62.   
  63.         for(k = 0; k < m-j; k++)  
  64.         {  
  65.             cov[i] += h[k] * data[n-(m-j)+k];//k针对core[k]  
  66.         }  
  67.     }  
  68.   
  69.     //****************************************************  
  70.     //偶数行用g[]进行卷积  
  71.     //****************************************************  
  72.     //前m/2+1行  
  73.     i = 1;  
  74.     for(j = 0; j < m/2; j+=2, i+=2)  
  75.     {  
  76.         for(k = m/2-j; k < m; k++ )  
  77.         {  
  78.             cov[i] += data[k-(m/2-j)] * g[k];//k针对core[k]  
  79.         }  
  80.   
  81.         for(k = n-m/2+j; k < n; k++ )  
  82.         {  
  83.             cov[i] += data[k] * g[k-(n-m/2+j)];//k针对data[k]  
  84.         }  
  85.     }  
  86.   
  87.     //中间的n-m行  
  88.     for( ; i <= (n-m)+m/2; i+=2)  
  89.     {  
  90.         for( j = 0; j < m; j++)  
  91.         {  
  92.             cov[i] += data[i-m/2+j] * g[j];  
  93.         }  
  94.     }  
  95.   
  96.     //最后m/2-1行  
  97. //  i = ( (n - m) + m/2 + 1 ) ;//*********  
  98.     for(j = 1; j <= m/2; j+=2, i+=2)  
  99.     {  
  100.         for(k = 0; k < j; k++)  
  101.         {  
  102.             cov[i] += data[k] * g[m-j-k];//k针对data[k]  
  103.         }  
  104.   
  105.         for(k = 0; k < m-j; k++)  
  106.         {  
  107.             cov[i] += g[k] * data[n-(m-j)+k];//k针对core[k]  
  108.         }  
  109.     }  
  110. }  
  111.   
  112. /*****************************************************************  
  113. *   排序函数  
  114. *  
  115. *   将卷积后的结果进行排序,使尺度系数和小波系数分开  
  116. *****************************************************************/  
  117. void Sort(double data[], double sort[], int n)  
  118. {  
  119.     for(int i = 0; i < n; i+=2)  
  120.     {  
  121.         sort[i/2] = data[i];  
  122.     }  
  123.   
  124.     for(i = 1; i < n; i+=2)  
  125.     {  
  126.         sort[n/2+i/2] = data[i];  
  127.     }  
  128.   
  129. }  
  130.   
  131. /***************************************************************** 
  132. * 一维小波变换函数 
  133. * 说明: 一维小波变换,可进行多次分解  
  134. * 输入参数: input[],输入信号; output[],小波变换结果,包括尺度系数 
  135. * 和小波系数两部分; temp[],存放中间结果;h[],Daubechies小波基低通滤 
  136. * 波器系数;g[],Daubechies小波基高通滤波器系数;n,输入信号长度; m, 
  137. * Daubechies小波基紧支集长度; nStep,小波变换分解次数 
  138. * 李承宇, lichengyu2345@126.com 
  139. *  2010-08-22   
  140. *****************************************************************/  
  141. void DWT1D(double input[], double output[], double temp[], double h[],   
  142.            double g[], int n, int m, int nStep)  
  143. {  
  144.     int i = 0;  
  145.   
  146.     for(i = 0; i < n; i++)  
  147.     {  
  148.         output[i] = input[i];  
  149.     }  
  150.   
  151.     for(i = 0; i < nStep; i++)  
  152.     {  
  153.         Covlution(output, h, g, temp, n, m);  
  154.         Sort(temp, output, n);  
  155.         n = n/2;  
  156.     }  
  157. }  
  158.   
  159. void main()  
  160. {  
  161.   
  162.     double data[LENGTH];//输入信号  
  163.     double temp[LENGTH];//中间结果  
  164.     double data_output[LENGTH];//一维小波变换后的结果  
  165.     int n = 0;//输入信号长度  
  166.     int m = 6;//Daubechies正交小波基长度  
  167.     int nStep = 6;//分解级数  
  168.     int i = 0;   
  169.     char s[32];//从txt文件中读取一行数据  
  170.   
  171.     static double h[] = {.332670552950, .806891509311, .459877502118,   
  172.         -.135011020010, -.085441273882, .035226291882};  
  173.     static double g[] = {.035226291882, .085441273882, -.135011020010,   
  174.         -.459877502118, .806891509311, -.332670552950};  
  175.   
  176.     //读取输入信号  
  177.     FILE *fp;  
  178.     fp=fopen("data.txt","r");  
  179.     if(fp==NULL) //如果读取失败  
  180.     {  
  181.         printf("错误!找不到要读取的文件/"data.txt/"/n");  
  182.         exit(1);//中止程序  
  183.     }  
  184.   
  185.     while( fgets(s, 32, fp) != NULL )//读取长度n要设置得长一点,要保证读到回车符,这样指针才会定位到下一行?回车符返回的是零值?是,非数字字符经过atoi变换都应该返回零值  
  186.     {  
  187.     //  fscanf(fp,"%d", &data[count]);//一定要有"&"啊!!!最后读了个回车符!适应能力不如atoi啊  
  188.         data[n] = atof(s);  
  189.         n++;  
  190.     }  
  191.   
  192.     //一维小波变换  
  193.     DWT1D(data, data_output, temp, h, g, n, m, nStep);  
  194.   
  195.     //一维小波变换后的结果写入txt文件  
  196.     fp=fopen("test.txt","w");  
  197.   
  198.     //打印一维小波变换后的结果  
  199.     for(i = 0; i < n/pow(2,nStep-1); i++)///pow(2,nStep-1)  
  200.     {  
  201.         printf("%f/n", data_output[i]);  
  202.         fprintf(fp,"%f/n", data_output[i]);  
  203.     }  
  204.   
  205.   
  206.     //关闭文件  
  207.     fclose(fp);  
  208. }  

 

 

 

4、测试结果:

输入信号x(i)为:

取f1 = 5, f2 = 10, f0 = 320, n = 512。x(i)如图1所示:

图1 输入信号

 

各级分解的结果如图2~图7所示,左半部分为尺度系数,右半部分为小波系数:

图2 1级分解结果

 

图3 2级分解结果

图4 3级分解结果

 

图5 4级分解结果

图6 5级分解结果

 

图7 6级分解结果

 

 

图8是各级小波系数和第6级尺度系数的完整结果:

图8 第6级尺度系数和各级小波系数的完整结果

图像处理之基础---一维小波变换,可多次分解