ivan@78: /*
ivan@78: File Name: MRDWT.c
ivan@78: Last Modification Date:	09/21/95	15:42:59
ivan@78: Current Version: MRDWT.c	2.4
ivan@78: File Creation Date: Wed Oct 12 08:44:43 1994
ivan@78: Author: Markus Lang  <lang@jazz.rice.edu>
ivan@78: 
ivan@78: Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
ivan@78: Created by Markus Lang, Department of ECE, Rice University. 
ivan@78: 
ivan@78: This software is distributed and licensed to you on a non-exclusive 
ivan@78: basis, free-of-charge. Redistribution and use in source and binary forms, 
ivan@78: with or without modification, are permitted provided that the following 
ivan@78: conditions are met:
ivan@78: 
ivan@78: 1. Redistribution of source code must retain the above copyright notice, 
ivan@78:    this list of conditions and the following disclaimer.
ivan@78: 2. Redistribution in binary form must reproduce the above copyright notice, 
ivan@78:    this list of conditions and the following disclaimer in the 
ivan@78:    documentation and/or other materials provided with the distribution.
ivan@78: 3. All advertising materials mentioning features or use of this software 
ivan@78:    must display the following acknowledgment: This product includes 
ivan@78:    software developed by Rice University, Houston, Texas and its contributors.
ivan@78: 4. Neither the name of the University nor the names of its contributors 
ivan@78:    may be used to endorse or promote products derived from this software 
ivan@78:    without specific prior written permission.
ivan@78: 
ivan@78: THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS, 
ivan@78: AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 
ivan@78: BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 
ivan@78: FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY 
ivan@78: OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
ivan@78: EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
ivan@78: PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
ivan@78: OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
ivan@78: WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
ivan@78: OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE 
ivan@78: USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
ivan@78: 
ivan@78: For information on commercial licenses, contact Rice University's Office of 
ivan@78: Technology Transfer at techtran@rice.edu or (713) 348-6173
ivan@78: 
ivan@78: Change History: Fixed the code such that 1D vectors passed to it can be in
ivan@78:                 either passed as a row or column vector. Also took care of 
ivan@78: 		the code such that it will compile with both under standard
ivan@78: 		C compilers as well as for ANSI C compilers
ivan@78: 		Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
ivan@78: 
ivan@78: MATLAB description:
ivan@78: %[yl,yh] = mrdwt(x,h,L);
ivan@78: % 
ivan@78: % function computes the redundant discrete wavelet transform y for a 1D or
ivan@78: % 2D input signal . redundant means here that the subsampling after each
ivan@78: % stage is omitted. yl contains the lowpass and yl the highpass
ivan@78: % components. In case of a 2D signal the ordering in yh is [lh hl hh lh hl
ivan@78: % ... ] (first letter refers to row, second to column filtering). 
ivan@78: %
ivan@78: %    Input:
ivan@78: %	x    : finite length 1D or 2D signal (implicitely periodized)
ivan@78: %       h    : scaling filter
ivan@78: %       L    : number of levels. in case of a 1D signal length(x) must be
ivan@78: %              divisible by 2^L; in case of a 2D signal the row and the
ivan@78: %              column dimension must be divisible by 2^L.
ivan@78: %   
ivan@78: %    Output:
ivan@78: %       yl   : lowpass component
ivan@78: %       yh   : highpass components
ivan@78: %
ivan@78: % see also: mdwt, midwt, mirdwt
ivan@78: 
ivan@78: 
ivan@78: */
ivan@78: 
ivan@78: #include <math.h>
ivan@78: #include <stdio.h>
ivan@78: #include <inttypes.h>
ivan@78: 
ivan@78: /*#define mat(a, i, j) (a[m*(j)+i]) */
ivan@78: #define max(a, b) ((a) > (b) ? (a) : (b))
ivan@78: #define mat(a, i, j) (*(a + (m*(j)+i))) 
ivan@78: 
ivan@78: 
ivan@78: #ifdef __STDC__
ivan@78: MRDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L,
ivan@78:       double *yl, double *yh)
ivan@78: #else
ivan@78: MRDWT(x, m, n, h, lh, L, yl, yh)
ivan@78: double *x, *h, *yl, *yh;
ivan@78: uintptr_t m, n, lh, L;
ivan@78: #endif
ivan@78: {
ivan@78:   double *tmp;
ivan@78:   double  *h0, *h1, *ydummyll, *ydummylh, *ydummyhl;
ivan@78:   double *ydummyhh, *xdummyl , *xdummyh;
ivan@78:   long i, j;
ivan@78:   uintptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o;
ivan@78:   uintptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
ivan@78:   xdummyl = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
ivan@78:   xdummyh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
ivan@78:   ydummyll = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78:   ydummylh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78:   ydummyhl = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78:   ydummyhh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78:   h0 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
ivan@78:   h1 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
ivan@78: 
ivan@78:   if (n==1){
ivan@78:     n = m;
ivan@78:     m = 1;
ivan@78:   }  
ivan@78:   /* analysis lowpass and highpass */
ivan@78:   for (i=0; i<lh; i++){
ivan@78:     h0[i] = h[lh-i-1];
ivan@78:     h1[i] =h[i];
ivan@78:   }
ivan@78:   for (i=0; i<lh; i+=2)
ivan@78:     h1[i] = -h1[i];
ivan@78:   
ivan@78:   actual_m = 2*m;
ivan@78:   actual_n = 2*n;
ivan@78:   for (i=0; i<m*n; i++)
ivan@78:     yl[i] = x[i];
ivan@78:   
ivan@78:   /* main loop */
ivan@78:   sample_f = 1;
ivan@78:   for (actual_L=1; actual_L <= L; actual_L++){
ivan@78:     actual_m = actual_m/2;
ivan@78:     actual_n = actual_n/2;
ivan@78:     /* actual (level dependent) column offset */
ivan@78:     if (m==1)
ivan@78:       c_o_a = n*(actual_L-1);
ivan@78:     else
ivan@78:       c_o_a = 3*n*(actual_L-1);
ivan@78:     c_o_a_p2n = c_o_a + 2*n;
ivan@78:     
ivan@78:     /* go by rows */
ivan@78:     n_cb = n/actual_n;                 /* # of column blocks per row */
ivan@78:     for (ir=0; ir<m; ir++){            /* loop over rows */
ivan@78:       for (n_c=0; n_c<n_cb; n_c++){    /* loop within one row */      
ivan@78: 	/* store in dummy variable */
ivan@78: 	ic = -sample_f + n_c;
ivan@78: 	for (i=0; i<actual_n; i++){    
ivan@78: 	  ic = ic + sample_f;
ivan@78: 	  xdummyl[i] = mat(yl, ir, ic);  
ivan@78: 	}
ivan@78: 	/* perform filtering lowpass/highpass */
ivan@78: 	fpconv(xdummyl, actual_n, h0, h1, lh, ydummyll, ydummyhh); 
ivan@78: 	/* restore dummy variables in matrices */
ivan@78: 	ic = -sample_f + n_c;
ivan@78: 	for  (i=0; i<actual_n; i++){    
ivan@78: 	  ic = ic + sample_f;
ivan@78: 	  mat(yl, ir, ic) = ydummyll[i];  
ivan@78: 	  mat(yh, ir, c_o_a+ic) = ydummyhh[i];  
ivan@78: 	} 
ivan@78:       }
ivan@78:     }
ivan@78:       
ivan@78:     /* go by columns in case of a 2D signal*/
ivan@78:     if (m>1){
ivan@78:       n_rb = m/actual_m;                 /* # of row blocks per column */
ivan@78:       for (ic=0; ic<n; ic++){            /* loop over column */
ivan@78: 	for (n_r=0; n_r<n_rb; n_r++){    /* loop within one column */
ivan@78: 	  /* store in dummy variables */
ivan@78: 	  ir = -sample_f + n_r;
ivan@78: 	  for (i=0; i<actual_m; i++){    
ivan@78: 	    ir = ir + sample_f;
ivan@78: 	    xdummyl[i] = mat(yl, ir, ic);  
ivan@78: 	    xdummyh[i] = mat(yh, ir,c_o_a+ic);  
ivan@78: 	  }
ivan@78: 	  /* perform filtering: first LL/LH, then HL/HH */
ivan@78: 	  fpconv(xdummyl, actual_m, h0, h1, lh, ydummyll, ydummylh); 
ivan@78: 	  fpconv(xdummyh, actual_m, h0, h1, lh, ydummyhl, ydummyhh); 
ivan@78: 	  /* restore dummy variables in matrices */
ivan@78: 	  ir = -sample_f + n_r;
ivan@78: 	  for (i=0; i<actual_m; i++){    
ivan@78: 	    ir = ir + sample_f;
ivan@78: 	    mat(yl, ir, ic) = ydummyll[i];  
ivan@78: 	    mat(yh, ir, c_o_a+ic) = ydummylh[i];  
ivan@78: 	    mat(yh, ir,c_o_a+n+ic) = ydummyhl[i];  
ivan@78: 	    mat(yh, ir, c_o_a_p2n+ic) = ydummyhh[i];  
ivan@78: 	  }
ivan@78: 	}
ivan@78:       }
ivan@78:     }
ivan@78:     sample_f = sample_f*2;
ivan@78:   }
ivan@78: }
ivan@78: 
ivan@78: #ifdef __STDC__
ivan@78: fpconv(double *x_in, uintptr_t lx, double *h0, double *h1, uintptr_t lh,
ivan@78:        double *x_outl, double *x_outh)
ivan@78: #else
ivan@78: fpconv(x_in, lx, h0, h1, lh, x_outl, x_outh)
ivan@78: double *x_in, *h0, *h1, *x_outl, *x_outh;
ivan@78: uintptr_t lx, lh;
ivan@78: #endif
ivan@78: {
ivan@78:   uintptr_t i, j;
ivan@78:   double x0, x1;
ivan@78: 
ivan@78:   for (i=lx; i < lx+lh-1; i++)
ivan@78:     x_in[i] = x_in[i-lx];
ivan@78:   for (i=0; i<lx; i++){
ivan@78:     x0 = 0;
ivan@78:     x1 = 0;
ivan@78:     for (j=0; j<lh; j++){
ivan@78:       x0 = x0 + x_in[j+i]*h0[lh-1-j];
ivan@78:       x1 = x1 + x_in[j+i]*h1[lh-1-j];
ivan@78:     }
ivan@78:     x_outl[i] = x0;
ivan@78:     x_outh[i] = x1;
ivan@78:   }
ivan@78: }