一维小波变换的C++实现


  将小波展开系数当成离散信号,尺度函数和小波函数的MRA方程系数看成数字滤波器组,根据Mallat快速算法的原理,小波变换对数据的处理方法可简化成对信号逐级采样和滤波的过程。

图1 小波变换的滤波器实现

(a)分解算法 (b)重构算法

  一层小波分解算法流程如图2所示,信号将先经过小波分解低通滤波器和高通滤波器,随后被降采样,实现数据重构。而滤波算法可简化为待处理信号与滤波器数组卷积的过程,为了保证卷积前和卷积后数组的长度相同,结合小波变换中数组延拓的思想,在实际编程过程中,可以将超过信号长度的那段数据以前端对齐的方式与前面一段数据相加。将卷积后的数组每2个点采样一次,即可获得小波分解后的尺度系数和小波系数。

图2 一层小波分解算法

(X:待分解数组;H,G:小波分解滤波器;C,D:小波重构后数组)

  小波重构算法是小波分解算法的逆运算,其流程为升采样和滤波,最后数据相加实现重构。小波重构算法中滤波可视为系数与小波重构滤波器的卷积,与小波正变换类似,在重构算法中,需要将卷积后的数组末位对齐相加,获得与原数组长度相同的卷积结果。将小波系数和尺度系数以2为步长进行升采样,将获得的新数组分别经过小波重构低通滤波器和高通滤波器,再将滤波后的两组数据相加,即实现了一层小波重构。

  1 #define LENGTH 512
  2 #define LEVEL 4
  3 #define L_core 6
  4 
  5 static void Covlution(double data[], double core[], double cov[], int LEN)
  6 {
  7         double temp[LENGTH + L_core - 1] = {0};
  8         int i = 0;
  9         int j = 0;
 10 
 11     for(i = 0; i < LEN; i++)
 12         {
 13                 for(j = 0; j  < L_core; j++)
 14                 {
 15             temp[i + j] += data[i] * core[j];
 16         }
 17     }
 18 
 19         for(i = 0; i < LEN; i++)
 20         {
 21                 if(i < L_core - 1)
 22                         cov[i] = temp[i] + temp[LEN + i];
 23                 else
 24                         cov[i] = temp[i];
 25         }
 26 
 27 }
 28 
 29 static void Covlution2(double data[], double core[], double cov[], int LEN)
 30 {
 31         double temp[LENGTH + L_core - 1] = {0};
 32         int i = 0;
 33         int j = 0;
 34 
 35     for(i = 0; i < LEN; i++)
 36         {
 37                 for(j = 0; j  < L_core; j++)
 38                 {
 39             temp[i + j] += data[i] * core[j];
 40         }
 41     }
 42 
 43         for(i = 0; i < LEN; i++)
 44         {
 45                 if(i < L_core - 1)
 46                         cov[i + LEN - L_core + 1] = temp[i] + temp[LEN + i];
 47                 else
 48                         cov[i - L_core + 1] = temp[i];
 49         }
 50 
 51 }
 52 
 53 static void DWT1D(double input[], double output[], double LF[], double HF[], int l)
 54 {
 55     int i = 0;
 56         double temp[LENGTH] = {0};
 57         int LEN = LENGTH / pow(2, l - 1);
 58 
 59     Covlution(input, LF, temp, LEN);
 60     for(i = 1; i < LEN; i += 2)
 61     {
 62         output[i/2] = temp[i];
 63     }
 64 
 65     Covlution(input, HF, temp, LEN);
 66     for(i = 1; i < LEN; i += 2)
 67     {
 68         output[LEN/2 + i/2] = temp[i];
 69     }
 70 }
 71 
 72 static void DWT(double input[], double output[], double LF[], double HF[], int len[])
 73 {
 74         int i;
 75         int j;
 76 
 77         len[0] = len[1] = LENGTH / pow(2, LEVEL);
 78         for(i = 2; i <= LEVEL; i++) len[i] = len[i - 1] * 2;
 79 
 80         DWT1D(input, output, LF, HF, 1);
 81         for(i = 2; i <= LEVEL; i++)
 82     {
 83                 for(j = 0; j < len[LEVEL + 2 - i]; j++) input[j] = output[j];
 84         DWT1D(input, output, LF, HF, i);
 85     }
 86 }
 87 
 88 static void IDWT1D(double input[], double output[], double LF[], double HF[], int l, int flag)
 89 {
 90     int i = 0;
 91         double temp[LENGTH] = {0};
 92         int LEN = l * 2;
 93 
 94         if(flag) Covlution2(input, HF, temp, LEN);
 95         else Covlution2(input, LF, temp, LEN);
 96 
 97         for(i = 0; i < LEN; i++)
 98     {
 99         output[i] = temp[i];
100     }
101 }
102 
103 static void IDWT(double input[], double output[], double LF[], double HF[], int len[], int level)
104 {
105         int i;
106         int j;
107         for(j = 0; j < len[LEVEL + 1 - level]; j++)
108         {
109                 output[2 * j] = 0;
110                 output[2 * j + 1] = input[j];
111         }
112         for(j = 0; j < 2 * len[LEVEL + 1 - level]; j++)
113         {
114                 input[j] = output[j];
115         }
116     IDWT1D(input, output, LF, HF, len[LEVEL + 1 - level], 1);
117 
118         for(i = level - 1; i > 0; i--)
119     {
120                 for(j = 0; j < len[LEVEL + 1 - i]; j++)
121                 {
122                         input[2 * j] = 0;
123                         input[2 * j + 1] = output[j];
124                 }
125         IDWT1D(input, output, LF, HF, len[LEVEL + 1 - i], 0);
126     }
127 }

 用C++算法实现的小波变换结果与MATLAB实现的小波变换结果对比(心电信号,db5小波,5层分解)。计算的结果几乎相差无异,大多数位于中间的数据计算结果完全相同,两端的数据则存在一定误差,这可能与对信号采取的延拓方式不同有关。

 

 

 

(a)小波变换结果(红线:MATLAB计算结果;蓝线:C++算法计算结果)

(b)小波变换误差(两个计算结果相减)

参考文献:网上各类资料,时间有点久,当时忘记记下了。

 


免责声明!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系本站邮箱yoyou2525@163.com删除。



 
粤ICP备18138465号  © 2018-2025 CODEPRJ.COM