Mercurial > hg > smallbox
view util/Rice Wavelet Toolbox/mdwt.c @ 240:04fb7249011e
Added tag ver_2.0_beta for changeset 71128ec3e532
author | luisf <luis.figueira@eecs.qmul.ac.uk> |
---|---|
date | Wed, 25 Apr 2012 13:06:47 +0100 |
parents | f69ae88b8be5 |
children |
line wrap: on
line source
/* File Name: mdwt.c Last Modification Date: 06/14/95 12:56:43 Current Version: mdwt.c 1.5 File Creation Date: Wed Oct 12 08:44:43 1994 Author: Markus Lang <lang@jazz.rice.edu> Copyright: All software, documentation, and related files in this distribution are Copyright (c) 1994 Rice University Permission is granted for use and non-profit distribution providing that this notice be clearly maintained. The right to distribute any portion for profit or as part of any commercial product is specifically reserved for the author. Change History: Fixed code such that the result has the same dimension as the input for 1D problems. Also, added some standard error checking. Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995 MATLAB gateway for MDWT.c, discrete wavelet transform */ #include <math.h> /*#include <malloc.h>*/ #include <stdio.h> #include "mex.h" #include "matrix.h" #if !defined(_WIN32) && !defined(_WIN64) #include <inttypes.h> #endif #define max(A,B) (A > B ? A : B) #define min(A,B) (A < B ? A : B) #define even(x) ((x & 1) ? 0 : 1) #define isint(x) ((x - floor(x)) > 0.0 ? 0 : 1) #define mat(a, i, j) (*(a + (m*(j)+i))) /* macro for matrix indices */ void mexFunction(const int nlhs,mxArray *plhs[],const int nrhs,const mxArray *prhs[]) { double *x, *h, *y, *Lf, *Lr; intptr_t m, n, h_col, h_row, lh, L, i, po2, j; double mtest, ntest; /* check for correct # of input variables */ if (nrhs>3){ mexErrMsgTxt("There are at most 3 input parameters allowed!"); return; } if (nrhs<2){ mexErrMsgTxt("There are at least 2 input parameters required!"); return; } x = mxGetPr(prhs[0]); n = mxGetN(prhs[0]); m = mxGetM(prhs[0]); h = mxGetPr(prhs[1]); h_col = mxGetN(prhs[1]); h_row = mxGetM(prhs[1]); if (h_col>h_row) lh = h_col; else lh = h_row; if (nrhs == 3){ L = (intptr_t) *mxGetPr(prhs[2]); if (L < 0) mexErrMsgTxt("The number of levels, L, must be a non-negative integer"); } else /* Estimate L */ { i=n;j=0; while (even(i)){ i=(i>>1); j++; } L=m;i=0; while (even(L)){ L=(L>>1); i++; } if(min(m,n) == 1) L = max(i,j); else L = min(i,j); if (L==0){ mexErrMsgTxt("Maximum number of levels is zero; no decomposition can be performed!"); return; } } /* Check the ROW dimension of input */ if(m > 1){ mtest = (double) m/pow(2.0, (double) L); if (!isint(mtest)) mexErrMsgTxt("The matrix row dimension must be of size m*2^(L)"); } /* Check the COLUMN dimension of input */ if(n > 1){ ntest = (double) n/pow(2.0, (double) L); if (!isint(ntest)) mexErrMsgTxt("The matrix column dimension must be of size n*2^(L)"); } plhs[0] = mxCreateDoubleMatrix(m,n,mxREAL); y = mxGetPr(plhs[0]); plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL); Lr = mxGetPr(plhs[1]); *Lr = L; MDWT(x, m, n, h, lh, L, y); } #ifdef __STDC__ MDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L, double *y) #else MDWT(x, m, n, h, lh, L, y) double *x, *h, *y; intptr_t m, n, lh, L; #endif { double *h0, *h1, *ydummyl, *ydummyh, *xdummy; int *prob; intptr_t i, j; intptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1; xdummy = (double *)mxCalloc(max(m,n)+lh-1,sizeof(double)); ydummyl =(double *) (intptr_t)mxCalloc(max(m,n),sizeof(double)); ydummyh = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double)); h0 =(double *)(intptr_t)mxCalloc(lh,sizeof(double)); h1 = (double *)(intptr_t)mxCalloc(lh,sizeof(double)); /* analysis lowpass and highpass */ if (n==1){ n = m; m = 1; } for (i=0; i<lh; i++){ h0[i] = h[lh-i-1]; h1[i] =h[i]; } for (i=0; i<lh; i+=2) h1[i] = -h1[i]; lhm1 = lh - 1; actual_m = 2*m; actual_n = 2*n; /* main loop */ for (actual_L=1; actual_L <= L; actual_L++){ if (m==1) actual_m = 1; else{ actual_m = actual_m/2; r_o_a = actual_m/2; } actual_n = actual_n/2; c_o_a = actual_n/2; /* go by rows */ for (ir=0; ir<actual_m; ir++){ /* loop over rows */ /* store in dummy variable */ for (i=0; i<actual_n; i++) if (actual_L==1) xdummy[i] = mat(x, ir, i); else xdummy[i] = mat(y, ir, i); /* perform filtering lowpass and highpass*/ fpsconv(xdummy, actual_n, h0, h1, lhm1, ydummyl, ydummyh); /* restore dummy variables in matrices */ ic = c_o_a; for (i=0; i<c_o_a; i++){ mat(y, ir, i) = ydummyl[i]; mat(y, ir, ic++) = ydummyh[i]; } } /* go by columns in case of a 2D signal*/ if (m>1){ for (ic=0; ic<actual_n; ic++){ /* loop over column */ /* store in dummy variables */ for (i=0; i<actual_m; i++) xdummy[i] = mat(y, i, ic); /* perform filtering lowpass and highpass*/ fpsconv(xdummy, actual_m, h0, h1, lhm1, ydummyl, ydummyh); /* restore dummy variables in matrix */ ir = r_o_a; for (i=0; i<r_o_a; i++){ mat(y, i, ic) = ydummyl[i]; mat(y, ir++, ic) = ydummyh[i]; } } } } } #ifdef __STDC__ fpsconv(double *x_in, intptr_t lx, double *h0, double *h1, intptr_t lhm1, double *x_outl, double *x_outh) #else fpsconv(x_in, lx, h0, h1, lhm1, x_outl, x_outh) double *x_in, *h0, *h1, *x_outl, *x_outh; intptr_t lx, lhm1; #endif { intptr_t i, j, ind; double x0, x1; for (i=lx; i < lx+lhm1; i++) x_in[i] = *(x_in+(i-lx)); ind = 0; for (i=0; i<(lx); i+=2){ x0 = 0; x1 = 0; for (j=0; j<=lhm1; j++){ x0 = x0 + x_in[i+j]*h0[lhm1-j]; x1 = x1 + x_in[i+j]*h1[lhm1-j]; } x_outl[ind] = x0; x_outh[ind++] = x1; } }