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

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

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

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

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

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

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

图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 }

内容版权声明:除非注明,否则皆为本站原创文章。

转载注明出处:https://www.heiqu.com/zwgjgf.html