|
@ -22,19 +22,28 @@ namespace wvlt{ |
|
|
return x[0] * coef[0] + x[stride] * coef[1] + x[2*stride] * coef[2] + x[3*stride] * coef[3]; |
|
|
return x[0] * coef[0] + x[stride] * coef[1] + x[2*stride] * coef[2] + x[3*stride] * coef[3]; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
// will overwrite x, x1 and x2 are next elements, or wrap around
|
|
|
// will calculate part of wavelete transform (in place!)
|
|
|
// size is size of vector x (so x[size-1] is valid)
|
|
|
// size is size of vector x (so x[size-1] is valid)
|
|
|
inline void wavelet_mul(double* x, double x1, double x2, unsigned int size, unsigned int stride){ |
|
|
// does not calculate "last two" elements (it does not assume periodicity)
|
|
|
assert(is_pow_of_two(size) && is_pow_of_two(stride) && 4*stride <= size); |
|
|
// calculates size/stride - 2 elements of the output
|
|
|
|
|
|
inline void wavelet_mul_base(double* x, unsigned int size, unsigned int stride){ |
|
|
|
|
|
assert(x && is_even(size) && is_pow_of_two(stride) && 4*stride <= size); |
|
|
|
|
|
|
|
|
for(int i = 0; i < size - 2*stride; i += 2*stride){ |
|
|
for(unsigned int i = 0; i < size - 2*stride; i += 2*stride){ |
|
|
double y1 = inner_product(x + i, evn_coef, stride); |
|
|
double y1 = inner_product(x + i, evn_coef, stride); |
|
|
double y2 = inner_product(x + i, odd_coef, stride); |
|
|
double y2 = inner_product(x + i, odd_coef, stride); |
|
|
x[i] = y1; |
|
|
x[i] = y1; |
|
|
x[i+stride] = y2; |
|
|
x[i+stride] = y2; |
|
|
} |
|
|
} |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
// x1 and x2 are next elements, or wrap around
|
|
|
|
|
|
// calculates size/stride elements of the output
|
|
|
|
|
|
inline void wavelet_mul(double* x, double x1, double x2, unsigned int size, unsigned int stride){ |
|
|
|
|
|
assert(x && is_even(size) && is_pow_of_two(stride) && 4*stride <= size); |
|
|
|
|
|
wavelet_mul_base(x, size, stride); |
|
|
|
|
|
|
|
|
int i = size - 2*stride; |
|
|
unsigned int i = size - 2*stride; |
|
|
double y1 = x[i] * evn_coef[0] + x[i+stride] * evn_coef[1] + x1 * evn_coef[2] + x2 * evn_coef[3]; |
|
|
double y1 = x[i] * evn_coef[0] + x[i+stride] * evn_coef[1] + x1 * evn_coef[2] + x2 * evn_coef[3]; |
|
|
double y2 = x[i] * odd_coef[0] + x[i+stride] * odd_coef[1] + x1 * odd_coef[2] + x2 * odd_coef[3]; |
|
|
double y2 = x[i] * odd_coef[0] + x[i+stride] * odd_coef[1] + x1 * odd_coef[2] + x2 * odd_coef[3]; |
|
|
x[i] = y1; |
|
|
x[i] = y1; |
|
@ -44,16 +53,16 @@ namespace wvlt{ |
|
|
// will overwrite x, x2 and x1 are previous elements, or wrap around
|
|
|
// will overwrite x, x2 and x1 are previous elements, or wrap around
|
|
|
// size is size of vector x (so x[size-1] is valid)
|
|
|
// size is size of vector x (so x[size-1] is valid)
|
|
|
inline void wavelet_inv(double* x, double x1, double x2, unsigned int size, unsigned int stride){ |
|
|
inline void wavelet_inv(double* x, double x1, double x2, unsigned int size, unsigned int stride){ |
|
|
assert(is_pow_of_two(size) && is_pow_of_two(stride) && 4*stride <= size); |
|
|
assert(x && is_even(size) && is_pow_of_two(stride) && 4*stride <= size); |
|
|
|
|
|
|
|
|
for(int i = size - 2*stride; i >= 2*stride; i -= 2*stride){ |
|
|
for(unsigned int i = size - 2*stride; i >= 2*stride; i -= 2*stride){ |
|
|
double y1 = inner_product(x + i - 2*stride, evn_coef_inv, stride); |
|
|
double y1 = inner_product(x + i - 2*stride, evn_coef_inv, stride); |
|
|
double y2 = inner_product(x + i - 2*stride, odd_coef_inv, stride); |
|
|
double y2 = inner_product(x + i - 2*stride, odd_coef_inv, stride); |
|
|
x[i] = y1; |
|
|
x[i] = y1; |
|
|
x[i+stride] = y2; |
|
|
x[i+stride] = y2; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
int i = 0; |
|
|
unsigned int i = 0; |
|
|
double y1 = x2 * evn_coef_inv[0] + x1 * evn_coef_inv[1] + x[i] * evn_coef_inv[2] + x[i+stride] * evn_coef_inv[3]; |
|
|
double y1 = x2 * evn_coef_inv[0] + x1 * evn_coef_inv[1] + x[i] * evn_coef_inv[2] + x[i+stride] * evn_coef_inv[3]; |
|
|
double y2 = x2 * odd_coef_inv[0] + x1 * odd_coef_inv[1] + x[i] * odd_coef_inv[2] + x[i+stride] * odd_coef_inv[3]; |
|
|
double y2 = x2 * odd_coef_inv[0] + x1 * odd_coef_inv[1] + x[i] * odd_coef_inv[2] + x[i+stride] * odd_coef_inv[3]; |
|
|
x[i] = y1; |
|
|
x[i] = y1; |
|
@ -61,15 +70,15 @@ namespace wvlt{ |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
inline void wavelet(double* x, unsigned int size){ |
|
|
inline void wavelet(double* x, unsigned int size){ |
|
|
assert(is_pow_of_two(size) && size >= 4); |
|
|
assert(x && is_pow_of_two(size) && size >= 4); |
|
|
for(int i = 1; i <= size / 4; i <<= 1){ |
|
|
for(unsigned int i = 1; i <= size / 4; i <<= 1){ |
|
|
wavelet_mul(x, x[0], x[i], size, i); |
|
|
wavelet_mul(x, x[0], x[i], size, i); |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
inline void unwavelet(double* x, unsigned int size){ |
|
|
inline void unwavelet(double* x, unsigned int size){ |
|
|
assert(is_pow_of_two(size) && size >= 4); |
|
|
assert(x && is_pow_of_two(size) && size >= 4); |
|
|
for(int i = size / 4; i >= 1; i >>= 1){ |
|
|
for(unsigned int i = size / 4; i >= 1; i >>= 1){ |
|
|
wavelet_inv(x, x[size-i], x[size-2*i], size, i); |
|
|
wavelet_inv(x, x[size-i], x[size-2*i], size, i); |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|