標題: 小波包3層分解C語言源程序 [打印本頁]
作者: Larry2018    時間: 2018-1-5 15:53
標題: 小波包3層分解C語言源程序
此程序是用于信號處理分析，突出奇異值的前段處理，對信號進行小波包分解，用C語言實現(xiàn)的，這個程序相比MATLAB程序，沒有進行邊緣的處理，直接根據(jù)卷積原理進行編寫而成的，如有對信號邊緣處理的要求，可以自行改進。程序在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. *　一維卷積函數(shù)
  9. *
  10. *　說明: 循環(huán)卷積,卷積結(jié)果的長度與輸入信號的長度相同
  11. *
  12. *　輸入?yún)?shù): data[],輸入信號; core[],卷積核; cov[],卷積結(jié)果;
  13. *            n,輸入信號長度; m,卷積核長度.
  14. *
  15. *　李承宇, lichengyu2345@126.com
  16. *
  17. *  2010-08-18  
  18. ******************************************************************/
  19. /*void Covlution(double data[], double core[], double cov[], int n, int m)
  20. {
  21.     int i = 0;
  22.     int j = 0;
  23.     int k = 0;

  25.     //將cov[]清零
  26.     for(i = 0; i < n; i++)
  27.     {
  28.         cov[i] = 0;
  29.     }

  31.     //前m/2+1行
  32.     i = 0;
  33.     for(j = 0; j < m/2; j++, i++)
  34.     {
  35.         for(k = m/2-j; k < m; k++ )
  36.         {
  37.             cov[i] += data[k-(m/2-j)] * core[k];//k針對core[k]
  38.         }

  40.         for(k = n-m/2+j; k < n; k++ )
  41.         {
  42.             cov[i] += data[k] * core[k-(n-m/2+j)];//k針對data[k]
  43.         }
  44.     }

  46.     //中間的n-m行
  47.     for( i = m/2; i <= (n-m)+m/2; i++)
  48.     {
  49.         for( j = 0; j < m; j++)
  50.         {
  51.             cov[i] += data[i-m/2+j] * core[j];
  52.         }
  53.     }

  55.     //最后m/2-1行
  56.     i = (n - m) + m/2 + 1;
  57.     for(j = 1; j < m/2; j++, i++)
  58.     {
  59.         for(k = 0; k < j; k++)
  60.         {
  61.             cov[i] += data[k] * core[m-j-k];//k針對data[k]
  62.         }

  64.         for(k = 0; k < m-j; k++)
  65.         {
  66.             cov[i] += core[k] * data[n-(m-j)+k];//k針對core[k]
  67.         }
  68.     }

  70. }
  71. */
  72. //定義一個線性卷積
  73. void Covlution(double data[], double core[], double cov[], int n, int m)
  74. {
  75.          int i = 0;
  76.     int j = 0;
  77.     int t = 0;

  79.     //將cov[]清零
  80.     for(j = 0; j < n+m-1; j++)
  81.     {
  82.         cov[j] = 0;
  83.     }

  85.          for(j=0;j<m+n-1;j++)
  86.          {
  87.                  if(j<=m-1)       //前面m行
  88.                  {
  89.                         for(i=0,t=j;t>=0;i++,t--)
  90.                          cov[j]+=data[i]*core[t];
  91.                  
  92.                  }
  93.                  else if(j<=n-1)    //中間n-m行
  94.                  {
  95.                          for(i=j-m+1,t=m-1;t>=0;i++,t--)
  96.                          cov[j]+=data[i]*core[t];
  97.                  }
  98.                  else     //后面m行
  99.                           {
  100.                                   for(i=j-m+1,t=m-1;i<n;i++,t--)
  101.                                  cov[j]+=data[i]*core[t];        
  102.                           }
  103.          }        

  105. }


  108. /******************************************************************
  109. *　一維小波變換函數(shù)
  110. *
  111. *　說明: 一維小波變換,只變換一次
  112. *
  113. *　輸入?yún)?shù): input[],輸入信號; output[],小波變換結(jié)果，包括尺度系數(shù)和
  114. *　小波系數(shù)兩部分; temp[],存放中間結(jié)果;h[],Daubechies小波基低通濾波器系數(shù);
  115. *　g[],Daubechies小波基高通濾波器系數(shù);n,輸入信號長度; m,Daubechies小波基緊支集長度.
  116. *
  117. *　李承宇, lichengyu2345@126.com
  118. *
  119. *  2010-08-19  
  120. ******************************************************************/
  121. void DWT1D_1(double input[], double output0[], double output1[],double temp[], double h[],
  122.            double g[], int n, int m)
  123. {
  124. //  double temp[LENGTH] = {0};//?????????????
  125.      
  126.     int i = 0;
  127. /*
  128.     //尺度系數(shù)和小波系數(shù)放在一起
  129.     Covlution(input, h, temp, n, m);

  131.     for(i = 0; i < n; i += 2)
  132.     {
  133.         output[i] = te mp[i];
  134.     }

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

  138.     for(i = 1; i < n; i += 2)
  139.     {
  140.         output[i] = temp[i];
  141.     }
  142. */
  143.    //尺度系數(shù)和小波系數(shù)分開
  144.    Covlution(input, h, temp, n, m);

  146.     for(i = 0; i < n+m-1; i += 2)
  147.     {
  148.         output0[i/2] = temp[i];//尺度系數(shù),進行了2抽值 ，即尺度空間，低頻概貌部分
  149.     }

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

  153.     for(i = 1; i < n+m-1; i +=2)
  154.     {
  155.     //   output[i] = temp[i];
  156.            output1[(i-1)/2] = temp[i];//小波系數(shù)，已經(jīng)進行了2抽取，即高頻細節(jié)部分
  157.     }



  161. }

  163. void DWT1D_2(double input[], double AA2[], double DA2[],double AD2[],double DD2[],double temp0[],double temp1[],double temp[], double h[],
  164.            double g[], int n, int m)
  165. {
  166.                 DWT1D_1(input, temp0, temp1,temp, h, g, n, m);
  167.                 int a1=(m+n)/2;
  168.                 DWT1D_1(temp0, AA2, DA2,temp, h, g, a1, m);
  169.                 int d1=(n+m-4)/2;
  170.                 DWT1D_1(temp1, AD2, DD2,temp, h, g, d1, m);
  171.                
  172. }

  174. void DWT1D_3(double input[], double AAA3[], double DAA3[],double ADA3[],double DDA3[],double AAD3[], double DAD3[],double ADD3[],double DDD3[],
  175.                          double temp0[],double temp1[],double temp2[],double temp3[],double temp00[], double temp11[], double temp[], double h[], double g[], int n, int m)
  176. {
  177.                 DWT1D_2(input, temp0, temp1,temp2, temp3,temp00, temp11,temp, h, g, n, m);
  178.                 int aa2=(m+n)/4;
  179.                 DWT1D_1(temp0, AAA3, DAA3,temp, h, g, aa2, m);
  180.                 int da2=(n+3*m-8)/4;
  181.                 DWT1D_1(temp1, ADA3, DDA3,temp, h, g, da2, m);
  182.                 int ad2=(n+3*m-4)/4;
  183.                 DWT1D_1(temp2, AAD3, DAD3,temp, h, g, ad2, m);
  184.                 int dd2=(n+3*m-12)/4;
  185.                 DWT1D_1(temp3, ADD3, DDD3,temp, h, g, dd2, m);
  186.                
  187.                
  188. }


  191. void main()
  192. {

  194.     double data[LENGTH];//輸入信號
  195.     double temp0[(LENGTH+DB_LENGTH)/4]; //先定義了一個中間結(jié)果
  196.     double temp1[(LENGTH+DB_LENGTH*3-8)/4];                        
  197.     double temp2[(LENGTH+DB_LENGTH*3-4)/4];                     
  198.     double temp3[(LENGTH+DB_LENGTH*3-12)/4];                       
  199.     double temp00[(LENGTH+DB_LENGTH)/2];                          
  200.     double temp11[(LENGTH+DB_LENGTH-4)/2];                                                     
  201.     double temp[LENGTH+DB_LENGTH-1];
  202.    
  203.     double AAA3[(LENGTH+DB_LENGTH*5)/8];//一維小波變換后的結(jié)果
  204.     double DAA3[(LENGTH+DB_LENGTH*5-16)/8];
  205.    
  206.         double ADA3[(LENGTH+DB_LENGTH*7-8)/8];   
  207.     double DDA3[(LENGTH+DB_LENGTH*7-24)/8];
  208.    
  209.     double AAD3[(LENGTH+DB_LENGTH*7-4)/8];
  210.         double DAD3[(LENGTH+DB_LENGTH*7-20)/8];
  211.          
  212.         double ADD3[(LENGTH+DB_LENGTH*7-12)/8];
  213.     double DDD3[(LENGTH+DB_LENGTH*7-28)/8];
  214.    
  215.           int aaa3=(LENGTH+DB_LENGTH*5)/8;//一維小波變換后的結(jié)果數(shù)組的長度
  216.     int daa3=(LENGTH+DB_LENGTH*5-16)/8;
  217.    
  218.         int ada3=(LENGTH+DB_LENGTH*7-8)/8;   
  219.     int dda3=(LENGTH+DB_LENGTH*7-24)/8;
  220.    
  221.     int aad3=(LENGTH+DB_LENGTH*7-4)/8;
  222.         int dad3=(LENGTH+DB_LENGTH*7-20)/8;
  223.          
  224.         int add3=(LENGTH+DB_LENGTH*7-12)/8;
  225.     int ddd3=(LENGTH+DB_LENGTH*7-28)/8;
  226.    
  227.     int n = 0;//輸入信號長度
  228.     int m = 8;//Daubechies正交小波基長度
  229.     int i = 0;
  230.     char s[32];//從txt文件中讀取一行數(shù)據(jù)
  231. /*//DB3
  232.     static double h[] = {.332670552950, .806891509311, .459877502118, -.135011020010,
  233.                     -.085441273882, .035226291882};
  234.     static double g[] = {.035226291882, .085441273882, -.135011020010, -.459877502118,
  235.                     .806891509311, -.332670552950};
  236. */
  237. //DB4

  239. static double h[] = {0.2303778133088964,   0.7148465705529154,         0.6308807679398587,-0.0279837694168599, -0.1870348117190931,     0.0308413818355607, 0.0328830116668852, -0.0105974017850690};//h[],Daubechies小波基低通濾波器系數(shù);
  240. static double g[] = {-0.0105974017850690, -0.0328830116668852, 0.0308413818355607, 0.1870348117190931,-0.0279837694168599, -0.6308807679398587, 0.7148465705529154, -0.2303778133088964 };//g[],Daubechies小波基高通濾波器系數(shù)
  241.     //讀取輸入信號
  242.     FILE *fp;
  243.         FILE *fp0,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6,*fp7;
  244.     fp=fopen("heb1.txt","r");
  245.     if(fp==NULL) //如果讀取失敗
  246.     {
  247.         printf("錯誤！找不到要讀取的文件hea1.txt/n");
  248.         exit(1);//中止程序
  249.     }

  251.     while( fgets(s, 32, fp) != NULL )//讀取長度n要設(shè)置得長一點，要保證讀到回車符，這樣指針才會定位到下一行？回車符返回的是零值？是，非數(shù)字字符經(jīng)過atoi變換都應該返回零值
  252.     {
  253.     //  fscanf(fp,"%d", &data[count]);//一定要有"&"啊�。�！最后讀了個回車符！適應能力不如atoi啊
  254.         data[n] = atof(s);
  255.         n++;
  256.     }

  258.     //一維小波變換
  259.     //DWT1D(data, data_output, temp, h, g, n, m);
  260.     DWT1D_3(data, AAA3, DAA3,ADA3,DDA3,AAD3, DAD3,ADD3,DDD3,
  261.                          temp0, temp1, temp2, temp3, temp00, temp11, temp, h, g, n, m);

  263.     //一維小波變換后的結(jié)果寫入txt文件
  264.     fp0=fopen("data_output_db4_aaa3.txt","w");
  265.     fp1=fopen("data_output_db4_daa3.txt","w");
  266.          fp2=fopen("data_output_db4_ada3.txt","w");
  267.     fp3=fopen("data_output_db4_dda3.txt","w");
  268.     fp4=fopen("data_output_db4_aad3.txt","w");
  269.     fp5=fopen("data_output_db4_dad3.txt","w");
  270.          fp6=fopen("data_output_db4_add3.txt","w");
  271.     fp7=fopen("data_output_db4_ddd3.txt","w");
  272.     //打印一維小波變換后的結(jié)果
  273.     for(i = 0; i <aaa3; i++)
  274.     {        //if(i==0)
  275.                 printf("%s\n","AAA3");
  276.         printf("%f\n", AAA3[i]);
  277.         fprintf(fp0,"%f\n", AAA3[i]);
  278. //                }
  279. //        else
  280. //        {
  281. //        printf("%f\n", AAA3[i]);
  282. //        fprintf(fp0,"%f\n", AAA3[i]);
  283. //        }
  284.   
  285.     }

  287.          for(i = 0; i <daa3; i++)
  288.     {//if(i==0)
  289.                 printf("%s\n","DAA3");
  290.         printf("%f\n", DAA3[i]);
  291.         fprintf(fp1,"%f\n", DAA3[i]);
  292. //                }
  293. //                else
  294. //                {
  295. //        printf("%f\n", DAA3[i]);
  296. //        fprintf(fp1,"%f\n", DAA3[i]);
  297. //                }
  298.   
  299.     }
  300.    
  301.         for(i = 0; i <ada3; i++)
  302.     {//if(i==0)
  303.                     printf("%s\n","ADA3");
  304.         printf("%f\n", ADA3[i]);
  305.         fprintf(fp2,"%f\n", ADA3[i]);
  306. //            }
  307. //                else
  308. //                {
  309. //        printf("%f\n", ADA3[i]);
  310. //        fprintf(fp2,"%f\n", ADA3[i]);
  311. //                }
  312.   
  313.     }
  314.    
  315.     for(i = 0; i <dda3; i++)
  316.     {        //if(i==0)
  317.                 printf("%s\n","DDA3");
  318.         printf("%f\n", DDA3[i]);
  319.         fprintf(fp3,"%f\n", DDA3[i]);
  320. //                }
  321. //                else
  322. //                {
  323. //        printf("%f\n", DDA3[i]);
  324. //        fprintf(fp3,"%f\n", DDA3[i]);
  325. //                }
  326.   
  327.     }
  328.    
  329.     for(i = 0; i <aad3; i++)
  330.     {        //if(i==0)
  331.                 printf("%s\n","AAD3");
  332.         printf("%f\n", AAD3[i]);
  333.         fprintf(fp4,"%f\n", AAD3[i]);        
  334. //                }
  335. //                else
  336. //                {
  337. //        printf("%f\n", AAD3[i]);
  338. //        fprintf(fp4,"%f\n", AAD3[i]);
  339. //                }
  340.   
  341.     }
  342.         
  343.         for(i = 0; i <dad3; i++)
  344.     {        //if(i==0)
  345.                 printf("%s\n","DAD3");
  346.         printf("%f\n", DAD3[i]);
  347.         fprintf(fp5,"%f\n", DAD3[i]);
  348. //                }
  349. //                else
  350. //                {        
  351. //        printf("%f\n", DAD3[i]);
  352. //        fprintf(fp5,"%f\n", DAD3[i]);
  353. //                }
  354.   
  355.     }
  356.    
  357.     for(i = 0; i <add3; i++)
  358.     {        //if(i==0)
  359.             printf("%s\n","ADD3");
  360.         printf("%f\n", ADD3[i]);
  361.         fprintf(fp6,"%f\n", ADD3[i]);
  362. //            }
  363. //            else
  364. //            {
  365. //        printf("%f\n", ADD3[i]);
  366. //        fprintf(fp6,"%f\n", ADD3[i]);
  367. //            }
  368. //               
  369.   
  370.     }
  371.    
  372.     for(i = 0; i <ddd3; i++)
  373.     {        //if(i==0)
  374.                 printf("%s\n","DDD3");
  375.         printf("%f\n", DDD3[i]);
  376.         fprintf(fp7,"%f\n", DDD3[i]);
  377. //                }
  378. //                else
  379. //                {
  380. //        printf("%f\n", DDD3[i]);
  381. //        fprintf(fp7,"%f\n", DDD3[i]);
  382. //                }
  383.   
  384.     }
  385.    
  386.     //關(guān)閉文件
  387.     fclose(fp);
  388.     fclose(fp0);
  389.     fclose(fp1);
  390.     fclose(fp2);
  391.     fclose(fp3);
  392.     fclose(fp4);
  393.     fclose(fp5);
  394.     fclose(fp6);
  395.     fclose(fp7);
  396. }


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

  401. /* 嘗試不對
  402.    double A1[(LENGTH+DB_LENGTH)/2];
  403.         double D1[(LENGTH+DB_LENGTH)/2];
  404.         int a1=(LENGTH+DB_LENGTH)/2;
  405.         int d1=(LENGTH+DB_LENGTH)/2;
  406.         double AA2[(LENGTH+DB_LENGTH)/4];
  407.         double DA2[(LENGTH+DB_LENGTH)/4];
  408.         int a2=(LENGTH+DB_LENGTH)/4;
  409.         int d2=(LENGTH+DB_LENGTH)/4;
  410.         double AAA3[(LENGTH+DB_LENGTH)/8];
  411.         double DAA3[(LENGTH+DB_LENGTH)/8];
  412.         int a3=(LENGTH+DB_LENGTH)/8;
  413.         int d3=(LENGTH+DB_LENGTH)/8;
  414.         double AAAA4[(LENGTH+DB_LENGTH)/16];
  415.         double DAAA4[(LENGTH+DB_LENGTH)/16];

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

  419.     for(i = 0; i < n+m-2; i += 2)
  420.     {
  421.         A1[i/2] = temp[i];//一層分解的低頻部分,進行了2抽值 ，即尺度空間，低頻概貌部分
  422.     }
  423.      Covlution(input, g, temp, n, m);

  425.     for(i =0 ; i < n+m-2; i +=2)
  426.     {
  427.            D1[i/2] = temp[i];//一層分解的高頻部分，已經(jīng)進行了2抽取，即高頻細節(jié)部分
  428.     }
  429.    

  431.       
  432.     Covlution(A1, h, temp, a1, m);
  433.     for(i = 0; i < a1+m-2; i += 2)
  434.     {
  435.         AA2[i/2] = temp[i];//一層分解的低頻部分,進行了2抽值 ，即尺度空間，低頻概貌部分
  436.     }
  437.    
  438.     Covlution(A1, g, temp, a1, m);
  439.     for(i = 0; i < a1+m-2; i += 2)
  440.     {
  441.         DA2[i/2] = temp[i];//一層分解的低頻部分,進行了2抽值 ，即尺度空間，低頻概貌部分
  442.     }
  443.    
  444.    
  445.    
  446.      Covlution(AA2, h, temp, a2, m);
  447.     for(i = 0; i < a2+m-2; i += 2)
  448.     {
  449.         AAA3[i/2] = temp[i];//一層分解的低頻部分,進行了2抽值 ，即尺度空間，低頻概貌部分
  450.     }
  451.    
  452.     Covlution(AA2, g, temp, a2, m);
  453.     for(i = 0; i < a2+m-2; i += 2)
  454.     {
  455.         DAA3[i/2] = temp[i];//一層分解的低頻部分,進行了2抽值 ，即尺度空間，低頻概貌部分
  456.     }
  457.    
  458.    
  459.    
  460.    
  461.    
  462.      Covlution(AAA3, h, temp, a3, m);
  463.     for(i = 0; i < a3+m-2; i += 2)
  464.     {
  465.         AAAA4[i/2] = temp[i];//一層分解的低頻部分,進行了2抽值 ，即尺度空間，低頻概貌部分
  466.     }
  467.    
  468.     Covlution(AAA3, g, temp, a3, m);
  469.     for(i = 0; i < a3+m-2; i += 2)
  470.     {
  471.         DAAA4[i/2] = temp[i];//一層分解的低頻部分,進行了2抽值 ，即尺度空間，低頻概貌部分
  472.     }
  473.    
  474.      for(i = 0; i <4116; i++)
  475.     {
  476.             if(i<=259)
  477.       output[i] = AAAA4[i];
  478.       else if(i<=519)
  479.       output[i] = DAAA4[i-260];
  480.       else if(i<=1035)
  481.       output[i] = DAA3[i-520];
  482.       else if(i<=2064)
  483.      output[i] = DA2[i-1036];
  484.      else
  485.      output[i] = DA2[i-2065];
  486.    
  487.     }

  489. */        


復制代碼






歡迎光臨 (http://www.torrancerestoration.com/bbs/) Powered by Discuz! X3.1