比特138

 找回密码
 立即注册
搜索
查看: 52|回复: 0

小波包三层分解C语言源程序

[复制链接]

485

主题

485

帖子

105

积分

初级会员

Rank: 2

金币
105
发表于 2020-5-2 00:46:09 | 显示全部楼层 |阅读模式
此程序是用于信号处理分析,突出奇异值的前段处理,对信号进行小波包分解,用C语言实现的,这个程序相比MATLAB程序,没有进行边缘的处理,直接根据卷积原理进行编写而成的,如有对信号边缘处理的要求,可以自行改进。程序在DEV C++软件中运行成功,能够进行,程序中有注释部分,欢迎大家改进提高。
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #define LENGTH 4096//信号长度
  6. #define DB_LENGTH 8  //Daubechies小波基紧支集长度
  7. /******************************************************************
  8. * 一维卷积函数
  9. *
  10. * 说明: 循环卷积,卷积结果的长度与输入信号的长度相同
  11. *
  12. * 输入参数: data[],输入信号; core[],卷积核; cov[],卷积结果;
  13. *            n,输入信号长度; m,卷积核长度.
  14. *
  15. ******************************************************************/
  16. /*void Covlution(double data[], double core[], double cov[], int n, int m)
  17. {
  18.     int i = 0;
  19.     int j = 0;
  20.     int k = 0;

  21.     //将cov[]清零
  22.     for(i = 0; i < n; i++)
  23.     {
  24.         cov[i] = 0;
  25.     }

  26.     //前m/2+1行
  27.     i = 0;
  28.     for(j = 0; j < m/2; j++, i++)
  29.     {
  30.         for(k = m/2-j; k < m; k++ )
  31.         {
  32.             cov[i] += data[k-(m/2-j)] * core[k];//k针对core[k]
  33.         }

  34.         for(k = n-m/2+j; k < n; k++ )
  35.         {
  36.             cov[i] += data[k] * core[k-(n-m/2+j)];//k针对data[k]
  37.         }
  38.     }

  39.     //中间的n-m行
  40.     for( i = m/2; i <= (n-m)+m/2; i++)
  41.     {
  42.         for( j = 0; j < m; j++)
  43.         {
  44.             cov[i] += data[i-m/2+j] * core[j];
  45.         }
  46.     }

  47.     //最后m/2-1行
  48.     i = (n - m) + m/2 + 1;
  49.     for(j = 1; j < m/2; j++, i++)
  50.     {
  51.         for(k = 0; k < j; k++)
  52.         {
  53.             cov[i] += data[k] * core[m-j-k];//k针对data[k]
  54.         }

  55.         for(k = 0; k < m-j; k++)
  56.         {
  57.             cov[i] += core[k] * data[n-(m-j)+k];//k针对core[k]
  58.         }
  59.     }

  60. }
  61. */
  62. //定义一个线性卷积
  63. void Covlution(double data[], double core[], double cov[], int n, int m)
  64. {
  65.          int i = 0;
  66.     int j = 0;
  67.     int t = 0;

  68.     //将cov[]清零
  69.     for(j = 0; j < n+m-1; j++)
  70.     {
  71.         cov[j] = 0;
  72.     }

  73.          for(j=0;j<m+n-1;j++)
  74.          {
  75.                  if(j<=m-1)       //前面m行
  76.                  {
  77.                         for(i=0,t=j;t>=0;i++,t--)
  78.                          cov[j]+=data[i]*core[t];
  79.                  
  80.                  }
  81.                  else if(j<=n-1)    //中间n-m行
  82.                  {
  83.                          for(i=j-m+1,t=m-1;t>=0;i++,t--)
  84.                          cov[j]+=data[i]*core[t];
  85.                  }
  86.                  else     //后面m行
  87.                           {
  88.                                   for(i=j-m+1,t=m-1;i<n;i++,t--)
  89.                                  cov[j]+=data[i]*core[t];        
  90.                           }
  91.          }        

  92. }


  93. /******************************************************************
  94. * 一维小波变换函数
  95. *
  96. * 说明: 一维小波变换,只变换一次
  97. *
  98. * 输入参数: input[],输入信号; output[],小波变换结果,包括尺度系数和
  99. * 小波系数两部分; temp[],存放中间结果;h[],Daubechies小波基低通滤波器系数;
  100. * g[],Daubechies小波基高通滤波器系数;n,输入信号长度; m,Daubechies小波基紧支集长度.
  101. *
  102. * 李承宇, lichengyu2345@126.com
  103. *
  104. *  2010-08-19  
  105. ******************************************************************/
  106. void DWT1D_1(double input[], double output0[], double output1[],double temp[], double h[],
  107.            double g[], int n, int m)
  108. {
  109. //  double temp[LENGTH] = {0};//?????????????
  110.      
  111.     int i = 0;
  112. /*
  113.     //尺度系数和小波系数放在一起
  114.     Covlution(input, h, temp, n, m);

  115.     for(i = 0; i < n; i += 2)
  116.     {
  117.         output[i] = te mp[i];
  118.     }

  119.     Covlution(input, g, temp, n, m);

  120.     for(i = 1; i < n; i += 2)
  121.     {
  122.         output[i] = temp[i];
  123.     }
  124. */
  125.    //尺度系数和小波系数分开
  126.    Covlution(input, h, temp, n, m);

  127.     for(i = 0; i < n+m-1; i += 2)
  128.     {
  129.         output0[i/2] = temp[i];//尺度系数,进行了2抽值 ,即尺度空间,低频概貌部分
  130.     }

  131. Covlution(input, g, temp, n, m);

  132.     for(i = 1; i < n+m-1; i +=2)
  133.     {
  134.     //   output[i] = temp[i];
  135.            output1[(i-1)/2] = temp[i];//小波系数,已经进行了2抽取,即高频细节部分
  136.     }



  137. }

  138. void DWT1D_2(double input[], double AA2[], double DA2[],double AD2[],double DD2[],double temp0[],double temp1[],double temp[], double h[],
  139.            double g[], int n, int m)
  140. {
  141.                 DWT1D_1(input, temp0, temp1,temp, h, g, n, m);
  142.                 int a1=(m+n)/2;
  143.                 DWT1D_1(temp0, AA2, DA2,temp, h, g, a1, m);
  144.                 int d1=(n+m-4)/2;
  145.                 DWT1D_1(temp1, AD2, DD2,temp, h, g, d1, m);
  146.                
  147. }

  148. void DWT1D_3(double input[], double AAA3[], double DAA3[],double ADA3[],double DDA3[],double AAD3[], double DAD3[],double ADD3[],double DDD3[],
  149.                          double temp0[],double temp1[],double temp2[],double temp3[],double temp00[], double temp11[], double temp[], double h[], double g[], int n, int m)
  150. {
  151.                 DWT1D_2(input, temp0, temp1,temp2, temp3,temp00, temp11,temp, h, g, n, m);
  152.                 int aa2=(m+n)/4;
  153.                 DWT1D_1(temp0, AAA3, DAA3,temp, h, g, aa2, m);
  154.                 int da2=(n+3*m-8)/4;
  155.                 DWT1D_1(temp1, ADA3, DDA3,temp, h, g, da2, m);
  156.                 int ad2=(n+3*m-4)/4;
  157.                 DWT1D_1(temp2, AAD3, DAD3,temp, h, g, ad2, m);
  158.                 int dd2=(n+3*m-12)/4;
  159.                 DWT1D_1(temp3, ADD3, DDD3,temp, h, g, dd2, m);
  160.                
  161.                
  162. }


  163. void main()
  164. {

  165.     double data[LENGTH];//输入信号
  166.     double temp0[(LENGTH+DB_LENGTH)/4]; //先定义了一个中间结果
  167.     double temp1[(LENGTH+DB_LENGTH*3-8)/4];                        
  168.     double temp2[(LENGTH+DB_LENGTH*3-4)/4];                     
  169.     double temp3[(LENGTH+DB_LENGTH*3-12)/4];                       
  170.     double temp00[(LENGTH+DB_LENGTH)/2];                          
  171.     double temp11[(LENGTH+DB_LENGTH-4)/2];                                                     
  172.     double temp[LENGTH+DB_LENGTH-1];
  173.    
  174.     double AAA3[(LENGTH+DB_LENGTH*5)/8];//一维小波变换后的结果
  175.     double DAA3[(LENGTH+DB_LENGTH*5-16)/8];
  176.    
  177.         double ADA3[(LENGTH+DB_LENGTH*7-8)/8];   
  178.     double DDA3[(LENGTH+DB_LENGTH*7-24)/8];
  179.    
  180.     double AAD3[(LENGTH+DB_LENGTH*7-4)/8];
  181.         double DAD3[(LENGTH+DB_LENGTH*7-20)/8];
  182.          
  183.         double ADD3[(LENGTH+DB_LENGTH*7-12)/8];
  184.     double DDD3[(LENGTH+DB_LENGTH*7-28)/8];
  185.    
  186.           int aaa3=(LENGTH+DB_LENGTH*5)/8;//一维小波变换后的结果数组的长度
  187.     int daa3=(LENGTH+DB_LENGTH*5-16)/8;
  188.    
  189.         int ada3=(LENGTH+DB_LENGTH*7-8)/8;   
  190.     int dda3=(LENGTH+DB_LENGTH*7-24)/8;
  191.    
  192.     int aad3=(LENGTH+DB_LENGTH*7-4)/8;
  193.         int dad3=(LENGTH+DB_LENGTH*7-20)/8;
  194.          
  195.         int add3=(LENGTH+DB_LENGTH*7-12)/8;
  196.     int ddd3=(LENGTH+DB_LENGTH*7-28)/8;
  197.    
  198.     int n = 0;//输入信号长度
  199.     int m = 8;//Daubechies正交小波基长度
  200.     int i = 0;
  201.     char s[32];//从txt文件中读取一行数据
  202. /*//DB3
  203.     static double h[] = {.332670552950, .806891509311, .459877502118, -.135011020010,
  204.                     -.085441273882, .035226291882};
  205.     static double g[] = {.035226291882, .085441273882, -.135011020010, -.459877502118,
  206.                     .806891509311, -.332670552950};
  207. */
  208. //DB4

  209. static double h[] = {0.2303778133088964,   0.7148465705529154,         0.6308807679398587,-0.0279837694168599, -0.1870348117190931,     0.0308413818355607, 0.0328830116668852, -0.0105974017850690};//h[],Daubechies小波基低通滤波器系数;
  210. static double g[] = {-0.0105974017850690, -0.0328830116668852, 0.0308413818355607, 0.1870348117190931,-0.0279837694168599, -0.6308807679398587, 0.7148465705529154, -0.2303778133088964 };//g[],Daubechies小波基高通滤波器系数
  211.     //读取输入信号
  212.     FILE *fp;
  213.         FILE *fp0,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6,*fp7;
  214.     fp=fopen("heb1.txt","r");
  215.     if(fp==NULL) //如果读取失败
  216.     {
  217.         printf("错误!找不到要读取的文件hea1.txt/n");
  218.         exit(1);//中止程序
  219.     }

  220.     while( fgets(s, 32, fp) != NULL )//读取长度n要设置得长一点,要保证读到回车符,这样指针才会定位到下一行?回车符返回的是零值?是,非数字字符经过atoi变换都应该返回零值
  221.     {
  222.     //  fscanf(fp,"%d", &data[count]);//一定要有"&"啊!!!最后读了个回车符!适应能力不如atoi啊
  223.         data[n] = atof(s);
  224.         n++;
  225.     }

  226.     //一维小波变换
  227.     //DWT1D(data, data_output, temp, h, g, n, m);
  228.     DWT1D_3(data, AAA3, DAA3,ADA3,DDA3,AAD3, DAD3,ADD3,DDD3,
  229.                          temp0, temp1, temp2, temp3, temp00, temp11, temp, h, g, n, m);

  230.     //一维小波变换后的结果写入txt文件
  231.     fp0=fopen("data_output_db4_aaa3.txt","w");
  232.     fp1=fopen("data_output_db4_daa3.txt","w");
  233.          fp2=fopen("data_output_db4_ada3.txt","w");
  234.     fp3=fopen("data_output_db4_dda3.txt","w");
  235.     fp4=fopen("data_output_db4_aad3.txt","w");
  236.     fp5=fopen("data_output_db4_dad3.txt","w");
  237.          fp6=fopen("data_output_db4_add3.txt","w");
  238.     fp7=fopen("data_output_db4_ddd3.txt","w");
  239.     //打印一维小波变换后的结果
  240.     for(i = 0; i <aaa3; i++)
  241.     {        //if(i==0)
  242.                 printf("%s\n","AAA3");
  243.         printf("%f\n", AAA3[i]);
  244.         fprintf(fp0,"%f\n", AAA3[i]);
  245. //                }
  246. //        else
  247. //        {
  248. //        printf("%f\n", AAA3[i]);
  249. //        fprintf(fp0,"%f\n", AAA3[i]);
  250. //        }
  251.   
  252.     }

  253.          for(i = 0; i <daa3; i++)
  254.     {//if(i==0)
  255.                 printf("%s\n","DAA3");
  256.         printf("%f\n", DAA3[i]);
  257.         fprintf(fp1,"%f\n", DAA3[i]);
  258. //                }
  259. //                else
  260. //                {
  261. //        printf("%f\n", DAA3[i]);
  262. //        fprintf(fp1,"%f\n", DAA3[i]);
  263. //                }
  264.   
  265.     }
  266.    
  267.         for(i = 0; i <ada3; i++)
  268.     {//if(i==0)
  269.                     printf("%s\n","ADA3");
  270.         printf("%f\n", ADA3[i]);
  271.         fprintf(fp2,"%f\n", ADA3[i]);
  272. //            }
  273. //                else
  274. //                {
  275. //        printf("%f\n", ADA3[i]);
  276. //        fprintf(fp2,"%f\n", ADA3[i]);
  277. //                }
  278.   
  279.     }
  280.    
  281.     for(i = 0; i <dda3; i++)
  282.     {        //if(i==0)
  283.                 printf("%s\n","DDA3");
  284.         printf("%f\n", DDA3[i]);
  285.         fprintf(fp3,"%f\n", DDA3[i]);
  286. //                }
  287. //                else
  288. //                {
  289. //        printf("%f\n", DDA3[i]);
  290. //        fprintf(fp3,"%f\n", DDA3[i]);
  291. //                }
  292.   
  293.     }
  294.    
  295.     for(i = 0; i <aad3; i++)
  296.     {        //if(i==0)
  297.                 printf("%s\n","AAD3");
  298.         printf("%f\n", AAD3[i]);
  299.         fprintf(fp4,"%f\n", AAD3[i]);        
  300. //                }
  301. //                else
  302. //                {
  303. //        printf("%f\n", AAD3[i]);
  304. //        fprintf(fp4,"%f\n", AAD3[i]);
  305. //                }
  306.   
  307.     }
  308.         
  309.         for(i = 0; i <dad3; i++)
  310.     {        //if(i==0)
  311.                 printf("%s\n","DAD3");
  312.         printf("%f\n", DAD3[i]);
  313.         fprintf(fp5,"%f\n", DAD3[i]);
  314. //                }
  315. //                else
  316. //                {        
  317. //        printf("%f\n", DAD3[i]);
  318. //        fprintf(fp5,"%f\n", DAD3[i]);
  319. //                }
  320.   
  321.     }
  322.    
  323.     for(i = 0; i <add3; i++)
  324.     {        //if(i==0)
  325.             printf("%s\n","ADD3");
  326.         printf("%f\n", ADD3[i]);
  327.         fprintf(fp6,"%f\n", ADD3[i]);
  328. //            }
  329. //            else
  330. //            {
  331. //        printf("%f\n", ADD3[i]);
  332. //        fprintf(fp6,"%f\n", ADD3[i]);
  333. //            }
  334. //               
  335.   
  336.     }
  337.    
  338.     for(i = 0; i <ddd3; i++)
  339.     {        //if(i==0)
  340.                 printf("%s\n","DDD3");
  341.         printf("%f\n", DDD3[i]);
  342.         fprintf(fp7,"%f\n", DDD3[i]);
  343. //                }
  344. //                else
  345. //                {
  346. //        printf("%f\n", DDD3[i]);
  347. //        fprintf(fp7,"%f\n", DDD3[i]);
  348. //                }
  349.   
  350.     }
  351.    
  352.     //关闭文件
  353.     fclose(fp);
  354.     fclose(fp0);
  355.     fclose(fp1);
  356.     fclose(fp2);
  357.     fclose(fp3);
  358.     fclose(fp4);
  359.     fclose(fp5);
  360.     fclose(fp6);
  361.     fclose(fp7);
  362. }


  363. /* run this program using the console pauser or add your own getch, system("pause") or input loop */

  364. /* 尝试不对
  365.    double A1[(LENGTH+DB_LENGTH)/2];
  366.         double D1[(LENGTH+DB_LENGTH)/2];
  367.         int a1=(LENGTH+DB_LENGTH)/2;
  368.         int d1=(LENGTH+DB_LENGTH)/2;
  369.         double AA2[(LENGTH+DB_LENGTH)/4];
  370.         double DA2[(LENGTH+DB_LENGTH)/4];
  371.         int a2=(LENGTH+DB_LENGTH)/4;
  372.         int d2=(LENGTH+DB_LENGTH)/4;
  373.         double AAA3[(LENGTH+DB_LENGTH)/8];
  374.         double DAA3[(LENGTH+DB_LENGTH)/8];
  375.         int a3=(LENGTH+DB_LENGTH)/8;
  376.         int d3=(LENGTH+DB_LENGTH)/8;
  377.         double AAAA4[(LENGTH+DB_LENGTH)/16];
  378.         double DAAA4[(LENGTH+DB_LENGTH)/16];

  379.     Covlution(input, h, temp, n, m);

  380.     for(i = 0; i < n+m-2; i += 2)
  381.     {
  382.         A1[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ,即尺度空间,低频概貌部分
  383.     }
  384.      Covlution(input, g, temp, n, m);

  385.     for(i =0 ; i < n+m-2; i +=2)
  386.     {
  387.            D1[i/2] = temp[i];//一层分解的高频部分,已经进行了2抽取,即高频细节部分
  388.     }
  389.    

  390.       
  391.     Covlution(A1, h, temp, a1, m);
  392.     for(i = 0; i < a1+m-2; i += 2)
  393.     {
  394.         AA2[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ,即尺度空间,低频概貌部分
  395.     }
  396.    
  397.     Covlution(A1, g, temp, a1, m);
  398.     for(i = 0; i < a1+m-2; i += 2)
  399.     {
  400.         DA2[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ,即尺度空间,低频概貌部分
  401.     }
  402.    
  403.    
  404.    
  405.      Covlution(AA2, h, temp, a2, m);
  406.     for(i = 0; i < a2+m-2; i += 2)
  407.     {
  408.         AAA3[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ,即尺度空间,低频概貌部分
  409.     }
  410.    
  411.     Covlution(AA2, g, temp, a2, m);
  412.     for(i = 0; i < a2+m-2; i += 2)
  413.     {
  414.         DAA3[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ,即尺度空间,低频概貌部分
  415.     }
  416.    
  417.    
  418.    
  419.    
  420.    
  421.      Covlution(AAA3, h, temp, a3, m);
  422.     for(i = 0; i < a3+m-2; i += 2)
  423.     {
  424.         AAAA4[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ,即尺度空间,低频概貌部分
  425.     }
  426.    
  427.     Covlution(AAA3, g, temp, a3, m);
  428.     for(i = 0; i < a3+m-2; i += 2)
  429.     {
  430.         DAAA4[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ,即尺度空间,低频概貌部分
  431.     }
  432.    
  433.      for(i = 0; i <4116; i++)
  434.     {
  435.             if(i<=259)
  436.       output[i] = AAAA4[i];
  437.       else if(i<=519)
  438.       output[i] = DAAA4[i-260];
  439.       else if(i<=1035)
  440.       output[i] = DAA3[i-520];
  441.       else if(i<=2064)
  442.      output[i] = DA2[i-1036];
  443.      else
  444.      output[i] = DA2[i-2065];
  445.    
  446.     }

  447. */        


复制代码


回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

手机版|小黑屋|比特138 |网站地图

GMT+8, 2020-10-1 04:25 , Processed in 0.084072 second(s), 7 queries , Redis On.

Powered by Discuz! X3.4

Copyright © 2001-2020, Tencent Cloud.

快速回复 返回顶部 返回列表