annotate util/Rice Wavelet Toolbox/mirdwt_r.c @ 78:f69ae88b8be5

added Rice Wavelet Toolbox with my modification, so it can be compiled on newer systems.
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Fri, 25 Mar 2011 15:27:33 +0000
parents
children
rev   line source
ivan@78 1 /*
ivan@78 2 File Name: MIRDWT.c
ivan@78 3 Last Modification Date: 06/14/95 16:22:45
ivan@78 4 Current Version: MIRDWT.c 2.4
ivan@78 5 File Creation Date: Wed Oct 12 08:44:43 1994
ivan@78 6 Author: Markus Lang <lang@jazz.rice.edu>
ivan@78 7
ivan@78 8 Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
ivan@78 9 Created by Markus Lang, Department of ECE, Rice University.
ivan@78 10
ivan@78 11 This software is distributed and licensed to you on a non-exclusive
ivan@78 12 basis, free-of-charge. Redistribution and use in source and binary forms,
ivan@78 13 with or without modification, are permitted provided that the following
ivan@78 14 conditions are met:
ivan@78 15
ivan@78 16 1. Redistribution of source code must retain the above copyright notice,
ivan@78 17 this list of conditions and the following disclaimer.
ivan@78 18 2. Redistribution in binary form must reproduce the above copyright notice,
ivan@78 19 this list of conditions and the following disclaimer in the
ivan@78 20 documentation and/or other materials provided with the distribution.
ivan@78 21 3. All advertising materials mentioning features or use of this software
ivan@78 22 must display the following acknowledgment: This product includes
ivan@78 23 software developed by Rice University, Houston, Texas and its contributors.
ivan@78 24 4. Neither the name of the University nor the names of its contributors
ivan@78 25 may be used to endorse or promote products derived from this software
ivan@78 26 without specific prior written permission.
ivan@78 27
ivan@78 28 THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
ivan@78 29 AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
ivan@78 30 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
ivan@78 31 FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
ivan@78 32 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
ivan@78 33 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
ivan@78 34 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
ivan@78 35 OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
ivan@78 36 WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
ivan@78 37 OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
ivan@78 38 USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
ivan@78 39
ivan@78 40 For information on commercial licenses, contact Rice University's Office of
ivan@78 41 Technology Transfer at techtran@rice.edu or (713) 348-6173
ivan@78 42
ivan@78 43 Change History: Fixed the code such that 1D vectors passed to it can be in
ivan@78 44 either passed as a row or column vector. Also took care of
ivan@78 45 the code such that it will compile with both under standard
ivan@78 46 C compilers as well as for ANSI C compilers
ivan@78 47 Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
ivan@78 48
ivan@78 49 Fix minor bug to allow maximum number of levels
ivan@78 50
ivan@78 51 MATLAB description:
ivan@78 52 %function x = mirdwt(yl,yh,h,L);
ivan@78 53 %
ivan@78 54 % function computes the inverse redundant discrete wavelet transform y for a
ivan@78 55 % 1D or 2D input signal. redundant means here that the subsampling after
ivan@78 56 % each stage of the forward transform has been omitted. yl contains the
ivan@78 57 % lowpass and yl the highpass components as computed, e.g., by mrdwt. In
ivan@78 58 % case of a 2D signal the ordering in yh is [lh hl hh lh hl ... ] (first
ivan@78 59 % letter refers to row, second to column filtering).
ivan@78 60 %
ivan@78 61 % Input:
ivan@78 62 % yl : lowpass component
ivan@78 63 % yh : highpass components
ivan@78 64 % h : scaling filter
ivan@78 65 % L : number of levels. in case of a 1D signal length(yl) must be
ivan@78 66 % divisible by 2^L; in case of a 2D signal the row and the
ivan@78 67 % column dimension must be divisible by 2^L.
ivan@78 68 %
ivan@78 69 % Output:
ivan@78 70 % x : finite length 1D or 2D signal
ivan@78 71 %
ivan@78 72 % see also: mdwt, midwt, mrdwt
ivan@78 73
ivan@78 74 */
ivan@78 75 #include <math.h>
ivan@78 76 #include <stdio.h>
ivan@78 77 #include <inttypes.h>
ivan@78 78
ivan@78 79 #define max(a, b) ((a) > (b) ? (a) : (b))
ivan@78 80 #define mat(a, i, j) (*(a + (m*(j)+i))) /* macro for matrix indices */
ivan@78 81
ivan@78 82 #ifdef __STDC__
ivan@78 83 MIRDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L,
ivan@78 84 double *yl, double *yh)
ivan@78 85 #else
ivan@78 86 MIRDWT(x, m, n, h, lh, L, yl, yh)
ivan@78 87 double *x, *h, *yl, *yh;
ivan@78 88 uintptr_t m, n, lh, L;
ivan@78 89 #endif
ivan@78 90 {
ivan@78 91 double *g0, *g1, *ydummyll, *ydummylh, *ydummyhl;
ivan@78 92 double *ydummyhh, *xdummyl , *xdummyh, *xh;
ivan@78 93 long i, j;
ivan@78 94 uintptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o, lhm1;
ivan@78 95 uintptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
ivan@78 96 xh = (double *)(uintptr_t)mxCalloc(m*n,sizeof(double));
ivan@78 97 xdummyl = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78 98 xdummyh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78 99 ydummyll = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
ivan@78 100 ydummylh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
ivan@78 101 ydummyhl = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
ivan@78 102 ydummyhh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
ivan@78 103 g0 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
ivan@78 104 g1 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
ivan@78 105
ivan@78 106 if (n==1){
ivan@78 107 n = m;
ivan@78 108 m = 1;
ivan@78 109 }
ivan@78 110 /* analysis lowpass and highpass */
ivan@78 111 for (i=0; i<lh; i++){
ivan@78 112 g0[i] = h[i]/2;
ivan@78 113 g1[i] = h[lh-i-1]/2;
ivan@78 114 }
ivan@78 115 for (i=1; i<=lh; i+=2)
ivan@78 116 g1[i] = -g1[i];
ivan@78 117
ivan@78 118 lhm1 = lh - 1;
ivan@78 119 /* 2^L */
ivan@78 120 sample_f = 1;
ivan@78 121 for (i=1; i<L; i++)
ivan@78 122 sample_f = sample_f*2;
ivan@78 123 actual_m = m/sample_f;
ivan@78 124 actual_n = n/sample_f;
ivan@78 125 /* restore yl in x */
ivan@78 126 for (i=0;i<m*n;i++)
ivan@78 127 x[i] = yl[i];
ivan@78 128
ivan@78 129 /* main loop */
ivan@78 130 for (actual_L=L; actual_L >= 1; actual_L--){
ivan@78 131 /* actual (level dependent) column offset */
ivan@78 132 if (m==1)
ivan@78 133 c_o_a = n*(actual_L-1);
ivan@78 134 else
ivan@78 135 c_o_a = 3*n*(actual_L-1);
ivan@78 136 c_o_a_p2n = c_o_a + 2*n;
ivan@78 137
ivan@78 138 /* go by columns in case of a 2D signal*/
ivan@78 139 if (m>1){
ivan@78 140 n_rb = m/actual_m; /* # of row blocks per column */
ivan@78 141 for (ic=0; ic<n; ic++){ /* loop over column */
ivan@78 142 for (n_r=0; n_r<n_rb; n_r++){ /* loop within one column */
ivan@78 143 /* store in dummy variables */
ivan@78 144 ir = -sample_f + n_r;
ivan@78 145 for (i=0; i<actual_m; i++){
ivan@78 146 ir = ir + sample_f;
ivan@78 147 ydummyll[i+lhm1] = mat(x, ir, ic);
ivan@78 148 ydummylh[i+lhm1] = mat(yh, ir, c_o_a+ic);
ivan@78 149 ydummyhl[i+lhm1] = mat(yh, ir,c_o_a+n+ic);
ivan@78 150 ydummyhh[i+lhm1] = mat(yh, ir, c_o_a_p2n+ic);
ivan@78 151 }
ivan@78 152 /* perform filtering and adding: first LL/LH, then HL/HH */
ivan@78 153 bpconv(xdummyl, actual_m, g0, g1, lh, ydummyll, ydummylh);
ivan@78 154 bpconv(xdummyh, actual_m, g0, g1, lh, ydummyhl, ydummyhh);
ivan@78 155 /* store dummy variables in matrices */
ivan@78 156 ir = -sample_f + n_r;
ivan@78 157 for (i=0; i<actual_m; i++){
ivan@78 158 ir = ir + sample_f;
ivan@78 159 mat(x, ir, ic) = xdummyl[i];
ivan@78 160 mat(xh, ir, ic) = xdummyh[i];
ivan@78 161 }
ivan@78 162 }
ivan@78 163 }
ivan@78 164 }
ivan@78 165
ivan@78 166 /* go by rows */
ivan@78 167 n_cb = n/actual_n; /* # of column blocks per row */
ivan@78 168 for (ir=0; ir<m; ir++){ /* loop over rows */
ivan@78 169 for (n_c=0; n_c<n_cb; n_c++){ /* loop within one row */
ivan@78 170 /* store in dummy variable */
ivan@78 171 ic = -sample_f + n_c;
ivan@78 172 for (i=0; i<actual_n; i++){
ivan@78 173 ic = ic + sample_f;
ivan@78 174 ydummyll[i+lhm1] = mat(x, ir, ic);
ivan@78 175 if (m>1)
ivan@78 176 ydummyhh[i+lhm1] = mat(xh, ir, ic);
ivan@78 177 else
ivan@78 178 ydummyhh[i+lhm1] = mat(yh, ir, c_o_a+ic);
ivan@78 179 }
ivan@78 180 /* perform filtering lowpass/highpass */
ivan@78 181 bpconv(xdummyl, actual_n, g0, g1, lh, ydummyll, ydummyhh);
ivan@78 182 /* restore dummy variables in matrices */
ivan@78 183 ic = -sample_f + n_c;
ivan@78 184 for (i=0; i<actual_n; i++){
ivan@78 185 ic = ic + sample_f;
ivan@78 186 mat(x, ir, ic) = xdummyl[i];
ivan@78 187 }
ivan@78 188 }
ivan@78 189 }
ivan@78 190 sample_f = sample_f/2;
ivan@78 191 actual_m = actual_m*2;
ivan@78 192 actual_n = actual_n*2;
ivan@78 193 }
ivan@78 194 }
ivan@78 195
ivan@78 196 #ifdef __STDC__
ivan@78 197 bpconv(double *x_out, uintptr_t lx, double *g0, double *g1, uintptr_t lh,
ivan@78 198 double *x_inl, double *x_inh)
ivan@78 199 #else
ivan@78 200 bpconv(x_out, lx, g0, g1, lh, x_inl, x_inh)
ivan@78 201 double *x_inl, *x_inh, *g0, *g1, *x_out;
ivan@78 202 uintptr_t lx, lh;
ivan@78 203 #endif
ivan@78 204 {
ivan@78 205 uintptr_t i, j;
ivan@78 206 double x0;
ivan@78 207
ivan@78 208 for (i=lh-2; i > -1; i--){
ivan@78 209 x_inl[i] = x_inl[lx+i];
ivan@78 210 x_inh[i] = x_inh[lx+i];
ivan@78 211 }
ivan@78 212 for (i=0; i<lx; i++){
ivan@78 213 x0 = 0;
ivan@78 214 for (j=0; j<lh; j++)
ivan@78 215 x0 = x0 + x_inl[j+i]*g0[lh-1-j] +
ivan@78 216 x_inh[j+i]*g1[lh-1-j];
ivan@78 217 x_out[i] = x0;
ivan@78 218 }
ivan@78 219 }