annotate util/Rice Wavelet Toolbox/mdwt_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: MDWT.c
ivan@78 3 Last Modification Date: 06/14/95 13:15:44
ivan@78 4 Current Version: MDWT.c 2.4
ivan@78 5 File Creation Date: Wed Oct 19 10:51:58 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 %y = mdwt(x,h,L);
ivan@78 50 %
ivan@78 51 % function computes the discrete wavelet transform y for a 1D or 2D input
ivan@78 52 % signal x.
ivan@78 53 %
ivan@78 54 % Input:
ivan@78 55 % x : finite length 1D or 2D signal (implicitely periodized)
ivan@78 56 % h : scaling filter
ivan@78 57 % L : number of levels. in case of a 1D signal length(x) must be
ivan@78 58 % divisible by 2^L; in case of a 2D signal the row and the
ivan@78 59 % column dimension must be divisible by 2^L.
ivan@78 60 %
ivan@78 61 % see also: midwt, mrdwt, mirdwt
ivan@78 62 */
ivan@78 63
ivan@78 64 #include <math.h>
ivan@78 65 #include <stdio.h>
ivan@78 66 #include <inttypes.h>
ivan@78 67
ivan@78 68 #define max(A,B) (A > B ? A : B)
ivan@78 69 #define mat(a, i, j) (*(a + (m*(j)+i))) /* macro for matrix indices */
ivan@78 70
ivan@78 71 #ifdef __STDC__
ivan@78 72 MDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L, double *y)
ivan@78 73 #else
ivan@78 74 MDWT(x, m, n, h, lh, L, y)
ivan@78 75 double *x, *h, *y;
ivan@78 76 uintptr_t m, n, lh, L;
ivan@78 77 #endif
ivan@78 78 {
ivan@78 79 double *h0, *h1, *ydummyl, *ydummyh, *xdummy;
ivan@78 80 int *prob;
ivan@78 81 uintptr_t i, j;
ivan@78 82 uintptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1;
ivan@78 83
ivan@78 84 xdummy = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
ivan@78 85 ydummyl =(double *) (uintptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78 86 ydummyh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78 87 h0 =(double *)(uintptr_t)mxCalloc(lh,sizeof(double));
ivan@78 88 h1 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
ivan@78 89
ivan@78 90
ivan@78 91 /* analysis lowpass and highpass */
ivan@78 92 if (n==1){
ivan@78 93 n = m;
ivan@78 94 m = 1;
ivan@78 95 }
ivan@78 96 for (i=0; i<lh; i++){
ivan@78 97 h0[i] = h[lh-i-1];
ivan@78 98 h1[i] =h[i];
ivan@78 99 }
ivan@78 100 for (i=0; i<lh; i+=2)
ivan@78 101 h1[i] = -h1[i];
ivan@78 102
ivan@78 103 lhm1 = lh - 1;
ivan@78 104 actual_m = 2*m;
ivan@78 105 actual_n = 2*n;
ivan@78 106
ivan@78 107 /* main loop */
ivan@78 108 for (actual_L=1; actual_L <= L; actual_L++){
ivan@78 109 if (m==1)
ivan@78 110 actual_m = 1;
ivan@78 111 else{
ivan@78 112 actual_m = actual_m/2;
ivan@78 113 r_o_a = actual_m/2;
ivan@78 114 }
ivan@78 115 actual_n = actual_n/2;
ivan@78 116 c_o_a = actual_n/2;
ivan@78 117
ivan@78 118 /* go by rows */
ivan@78 119 for (ir=0; ir<actual_m; ir++){ /* loop over rows */
ivan@78 120 /* store in dummy variable */
ivan@78 121 for (i=0; i<actual_n; i++)
ivan@78 122 if (actual_L==1)
ivan@78 123 xdummy[i] = mat(x, ir, i);
ivan@78 124 else
ivan@78 125 xdummy[i] = mat(y, ir, i);
ivan@78 126 /* perform filtering lowpass and highpass*/
ivan@78 127 fpsconv(xdummy, actual_n, h0, h1, lhm1, ydummyl, ydummyh);
ivan@78 128 /* restore dummy variables in matrices */
ivan@78 129 ic = c_o_a;
ivan@78 130 for (i=0; i<c_o_a; i++){
ivan@78 131 mat(y, ir, i) = ydummyl[i];
ivan@78 132 mat(y, ir, ic++) = ydummyh[i];
ivan@78 133 }
ivan@78 134 }
ivan@78 135
ivan@78 136 /* go by columns in case of a 2D signal*/
ivan@78 137 if (m>1){
ivan@78 138 for (ic=0; ic<actual_n; ic++){ /* loop over column */
ivan@78 139 /* store in dummy variables */
ivan@78 140 for (i=0; i<actual_m; i++)
ivan@78 141 xdummy[i] = mat(y, i, ic);
ivan@78 142 /* perform filtering lowpass and highpass*/
ivan@78 143 fpsconv(xdummy, actual_m, h0, h1, lhm1, ydummyl, ydummyh);
ivan@78 144 /* restore dummy variables in matrix */
ivan@78 145 ir = r_o_a;
ivan@78 146 for (i=0; i<r_o_a; i++){
ivan@78 147 mat(y, i, ic) = ydummyl[i];
ivan@78 148 mat(y, ir++, ic) = ydummyh[i];
ivan@78 149 }
ivan@78 150 }
ivan@78 151 }
ivan@78 152 }
ivan@78 153 }
ivan@78 154
ivan@78 155 #ifdef __STDC__
ivan@78 156 fpsconv(double *x_in, uintptr_t lx, double *h0, double *h1, uintptr_t lhm1,
ivan@78 157 double *x_outl, double *x_outh)
ivan@78 158 #else
ivan@78 159 fpsconv(x_in, lx, h0, h1, lhm1, x_outl, x_outh)
ivan@78 160 double *x_in, *h0, *h1, *x_outl, *x_outh;
ivan@78 161 uintptr_t lx, lhm1;
ivan@78 162 #endif
ivan@78 163
ivan@78 164 {
ivan@78 165 uintptr_t i, j, ind;
ivan@78 166 double x0, x1;
ivan@78 167
ivan@78 168 for (i=lx; i < lx+lhm1; i++)
ivan@78 169 x_in[i] = *(x_in+(i-lx));
ivan@78 170 ind = 0;
ivan@78 171 for (i=0; i<(lx); i+=2){
ivan@78 172 x0 = 0;
ivan@78 173 x1 = 0;
ivan@78 174 for (j=0; j<=lhm1; j++){
ivan@78 175 x0 = x0 + x_in[i+j]*h0[lhm1-j];
ivan@78 176 x1 = x1 + x_in[i+j]*h1[lhm1-j];
ivan@78 177 }
ivan@78 178 x_outl[ind] = x0;
ivan@78 179 x_outh[ind++] = x1;
ivan@78 180 }
ivan@78 181 }