LineWavelet.java (8577B)
1 package be.tarsos.dsp.wavelet.lift; 2 3 /** 4 * <p> 5 * Line (with slope) wavelet 6 * </p> 7 * 8 * <p> 9 * The wavelet Lifting Scheme "LineWavelet" wavelet approximates the data set 10 * using a LineWavelet with with slope (in contrast to the HaarWavelet wavelet 11 * where a LineWavelet has zero slope is used to approximate the data). 12 * </p> 13 * 14 * <p> 15 * The predict stage of the LineWavelet wavelet "predicts" that an odd point 16 * will lie midway between its two neighboring even points. That is, that the 17 * odd point will lie on a LineWavelet between the two adjacent even points. The 18 * difference between this "prediction" and the actual odd value replaces the 19 * odd element. 20 * </p> 21 * 22 * <p> 23 * The update stage calculates the average of the odd and even element pairs, 24 * although the method is indirect, since the predict phase has over written the 25 * odd value. 26 * </p> 27 * 28 * <h4> 29 * Copyright and Use</h4> 30 * 31 * <p> 32 * You may use this source code without limitation and without fee as long as 33 * you include: 34 * </p> 35 * <blockquote> This software was written and is copyrighted by Ian Kaplan, Bear 36 * Products International, www.bearcave.com, 2001. </blockquote> 37 * <p> 38 * This software is provided "as is", without any warrenty or claim as to its 39 * usefulness. Anyone who uses this source code uses it at their own risk. Nor 40 * is any support provided by Ian Kaplan and Bear Products International. 41 * <p> 42 * Please send any bug fixes or suggested source changes to: 43 * 44 * <pre> 45 * iank@bearcave.com 46 * </pre> 47 * 48 * @author Ian Kaplan 49 */ 50 public class LineWavelet extends LiftingSchemeBaseWavelet { 51 52 /** 53 * <p> 54 * Calculate an extra "even" value for the LineWavelet wavelet algorithm at 55 * the end of the data series. Here we pretend that the last two values in 56 * the data series are at the x-axis coordinates 0 and 1, respectively. We 57 * then need to calculate the y-axis value at the x-axis coordinate 2. This 58 * point lies on a LineWavelet running through the points at 0 and 1. 59 * </p> 60 * <p> 61 * Given two points, x<sub>1</sub>, y<sub>1</sub> and x<sub>2</sub>, 62 * y<sub>2</sub>, where 63 * </p> 64 * 65 * <pre> 66 * x<sub>1</sub> = 0 67 * x<sub>2</sub> = 1 68 * </pre> 69 * <p> 70 * calculate the point on the LineWavelet at x<sub>3</sub>, y<sub>3</sub>, 71 * where 72 * </p> 73 * 74 * <pre> 75 * x<sub>3</sub> = 2 76 * </pre> 77 * <p> 78 * The "two-point equation" for a LineWavelet given x<sub>1</sub>, 79 * y<sub>1</sub> and x<sub>2</sub>, y<sub>2</sub> is 80 * </p> 81 * 82 * <pre> 83 * . y<sub>2</sub> - y<sub>1</sub> 84 * (y - y<sub>1</sub>) = -------- (x - x<sub>1</sub>) 85 * . x<sub>2</sub> - x<sub>1</sub> 86 * </pre> 87 * <p> 88 * Solving for y 89 * </p> 90 * 91 * <pre> 92 * . y<sub>2</sub> - y<sub>1</sub> 93 * y = -------- (x - x<sub>1</sub>) + y<sub>1</sub> 94 * . x<sub>2</sub> - x<sub>1</sub> 95 * </pre> 96 * <p> 97 * Since x<sub>1</sub> = 0 and x<sub>2</sub> = 1 98 * </p> 99 * 100 * <pre> 101 * . y<sub>2</sub> - y<sub>1</sub> 102 * y = -------- (x - 0) + y<sub>1</sub> 103 * . 1 - 0 104 * </pre> 105 * <p> 106 * or 107 * </p> 108 * 109 * <pre> 110 * y = (y<sub>2</sub> - y<sub>1</sub>)*x + y<sub>1</sub> 111 * </pre> 112 * <p> 113 * We're calculating the value at x<sub>3</sub> = 2, so 114 * </p> 115 * 116 * <pre> 117 * y = 2*y<sub>2</sub> - 2*y<sub>1</sub> + y<sub>1</sub> 118 * </pre> 119 * <p> 120 * or 121 * </p> 122 * 123 * <pre> 124 * y = 2*y<sub>2</sub> - y<sub>1</sub> 125 * </pre> 126 */ 127 private float new_y(float y1, float y2) { 128 float y = 2 * y2 - y1; 129 return y; 130 } 131 132 /** 133 * <p> 134 * Predict phase of LineWavelet Lifting Scheme wavelet 135 * </p> 136 * 137 * <p> 138 * The predict step attempts to "predict" the value of an odd element from 139 * the even elements. The difference between the prediction and the actual 140 * element is stored as a wavelet coefficient. 141 * </p> 142 * <p> 143 * The "predict" step takes place after the split step. The split step will 144 * move the odd elements (b<sub>j</sub>) to the second half of the array, 145 * leaving the even elements (a<sub>i</sub>) in the first half 146 * </p> 147 * 148 * <pre> 149 * a<sub>0</sub>, a<sub>1</sub>, a<sub>1</sub>, a<sub>3</sub>, b<sub>0</sub>, b<sub>1</sub>, b<sub>2</sub>, b<sub>2</sub>, 150 * </pre> 151 * <p> 152 * The predict step of the LineWavelet wavelet "predicts" that the odd 153 * element will be on a LineWavelet between two even elements. 154 * </p> 155 * 156 * <pre> 157 * b<sub>j+1,i</sub> = b<sub>j,i</sub> - (a<sub>j,i</sub> + a<sub>j,i+1</sub>)/2 158 * </pre> 159 * <p> 160 * Note that when we get to the end of the data series the odd element is 161 * the last element in the data series (remember, wavelet algorithms work on 162 * data series with 2<sup>n</sup> elements). Here we "predict" that the odd 163 * element will be on a LineWavelet that runs through the last two even 164 * elements. This can be calculated by assuming that the last two even 165 * elements are located at x-axis coordinates 0 and 1, respectively. The odd 166 * element will be at 2. The <i>new_y()</i> function is called to do this 167 * simple calculation. 168 * </p> 169 */ 170 protected void predict(float[] vec, int N, int direction) { 171 int half = N >> 1; 172 float predictVal; 173 174 for (int i = 0; i < half; i++) { 175 int j = i + half; 176 if (i < half - 1) { 177 predictVal = (vec[i] + vec[i + 1]) / 2; 178 } else if (N == 2) { 179 predictVal = vec[0]; 180 } else { 181 // calculate the last "odd" prediction 182 predictVal = new_y(vec[i - 1], vec[i]); 183 } 184 185 if (direction == forward) { 186 vec[j] = vec[j] - predictVal; 187 } else if (direction == inverse) { 188 vec[j] = vec[j] + predictVal; 189 } else { 190 System.out.println("predictline::predict: bad direction value"); 191 } 192 } 193 } // predict 194 195 /** 196 * <p> 197 * The predict phase works on the odd elements in the second half of the 198 * array. The update phase works on the even elements in the first half of 199 * the array. The update phase attempts to preserve the average. After the 200 * update phase is completed the average of the even elements should be 201 * approximately the same as the average of the input data set from the 202 * previous iteration. The result of the update phase becomes the input for 203 * the next iteration. 204 * </p> 205 * <p> 206 * In a HaarWavelet wavelet the average that replaces the even element is 207 * calculated as the average of the even element and its associated odd 208 * element (e.g., its odd neighbor before the split). This is not possible 209 * in the LineWavelet wavelet since the odd element has been replaced by the 210 * difference between the odd element and the mid-point of its two even 211 * neighbors. As a result, the odd element cannot be recovered. 212 * </p> 213 * <p> 214 * The value that is added to the even element to preserve the average is 215 * calculated by the equation shown below. This equation is given in Wim 216 * Sweldens' journal articles and his tutorial (<i>Building Your Own 217 * Wavelets at Home</i>) and in <i>Ripples in Mathematics</i>. A somewhat 218 * more complete derivation of this equation is provided in <i>Ripples in 219 * Mathematics</i> by A. Jensen and A. la Cour-Harbo, Springer, 2001. 220 * </p> 221 * <p> 222 * The equation used to calculate the average is shown below for a given 223 * iteratin <i>i</i>. Note that the predict phase has already completed, so 224 * the odd values belong to iteration <i>i+1</i>. 225 * </p> 226 * 227 * <pre> 228 * even<sub>i+1,j</sub> = even<sub>i,j</sub> op (odd<sub>i+1,k-1</sub> + odd<sub>i+1,k</sub>)/4 229 * </pre> 230 * <p> 231 * There is an edge problem here, when i = 0 and k = N/2 (e.g., there is no 232 * k-1 element). We assume that the odd<sub>i+1,k-1</sub> is the same as 233 * odd<sub>k</sub>. So for the first element this becomes 234 * 235 * <pre> 236 * (2 * odd<sub>k</sub>)/4 237 * </pre> 238 * <p> 239 * or 240 * </p> 241 * 242 * <pre> 243 * odd<sub>k</sub>/2 244 * </pre> 245 */ 246 protected void update(float[] vec, int N, int direction) { 247 int half = N >> 1; 248 249 for (int i = 0; i < half; i++) { 250 int j = i + half; 251 float val; 252 253 if (i == 0) { 254 val = vec[j] / 2.0f; 255 } else { 256 val = (vec[j - 1] + vec[j]) / 4.0f; 257 } 258 if (direction == forward) { 259 vec[i] = vec[i] + val; 260 } else if (direction == inverse) { 261 vec[i] = vec[i] - val; 262 } else { 263 System.out.println("update: bad direction value"); 264 } 265 } // for 266 } 267 268 } // LineWavelet