dwt/
transform.rs

1use num::Float;
2
3use Operation;
4use wavelet::Wavelet;
5
6/// The transform.
7pub trait Transform<T> {
8    /// Perform the transform.
9    ///
10    /// The number of points should be divisible by `2^level`. If the operation
11    /// is forward, the data are replaced by the approximation and detail
12    /// coefficients stored in the first and second halves of `data`,
13    /// respectively. If the operation is inverse, the data are assumed to be
14    /// stored according to the above convention.
15    fn transform(&mut self, operation: Operation, wavelet: &Wavelet<T>, level: usize);
16}
17
18impl<T> Transform<T> for [T] where T: Float {
19    fn transform(&mut self, operation: Operation, wavelet: &Wavelet<T>, level: usize) {
20        if level == 0 {
21            return;
22        }
23        let n = self.len();
24        assert!(n % (1 << level) == 0);
25        let mut work = Vec::with_capacity(n);
26        unsafe { work.set_len(n) };
27        match operation {
28            Operation::Forward => {
29                for i in 0..level {
30                    forward_step(self, wavelet, n >> i, &mut work);
31                }
32            },
33            Operation::Inverse => {
34                for i in 0..level {
35                    inverse_step(self, wavelet, n >> (level - i - 1), &mut work);
36                }
37            },
38        }
39    }
40}
41
42macro_rules! copy(
43    ($source:ident, $destination:ident, $n:expr) => ({
44        use std::ptr::copy_nonoverlapping as copy;
45        unsafe { copy($source.as_ptr(), $destination.as_mut_ptr(), $n) };
46    });
47);
48
49macro_rules! zero(
50    ($buffer:expr) => ({
51        use std::ptr::write_bytes as write;
52        unsafe { write($buffer.as_mut_ptr(), 0, $buffer.len()) };
53    });
54);
55
56#[inline(always)]
57pub fn forward_step<T>(data: &mut [T], wavelet: &Wavelet<T>, n: usize, work: &mut [T])
58    where T: Float
59{
60    zero!(work);
61    let nm = wavelet.length * n - wavelet.offset;
62    let nh = n >> 1;
63    for i in 0..nh {
64        let (mut h, mut g) = (T::zero(), T::zero());
65        let k = 2 * i + nm;
66        for j in 0..wavelet.length {
67            let k = (k + j) % n;
68            h = h + wavelet.dec_lo[j] * data[k];
69            g = g + wavelet.dec_hi[j] * data[k];
70        }
71        work[i] = work[i] + h;
72        work[i + nh] = work[i + nh] + g;
73    }
74    copy!(work, data, n);
75}
76
77#[inline(always)]
78pub fn inverse_step<T>(data: &mut [T], wavelet: &Wavelet<T>, n: usize, work: &mut [T])
79    where T: Float
80{
81    zero!(work);
82    let nm = wavelet.length * n - wavelet.offset;
83    let nh = n >> 1;
84    for i in 0..nh {
85        let (h, g) = (data[i], data[i + nh]);
86        let k = 2 * i + nm;
87        for j in 0..wavelet.length {
88            let k = (k + j) % n;
89            work[k] = work[k] + wavelet.rec_lo[j] * h + wavelet.rec_hi[j] * g;
90        }
91    }
92    copy!(work, data, n);
93}
dwt/transform.rs

dwt/
transform.rs
