annotate util/Rice Wavelet Toolbox/mdwt.c @ 236:5f4e47b78f2b ver_2.0_alpha1

Added example.
author luisf <luis.figueira@eecs.qmul.ac.uk>
date Thu, 19 Apr 2012 17:59:08 +0100
parents f69ae88b8be5
children
rev   line source
ivan@78 1 /*
ivan@78 2 File Name: mdwt.c
ivan@78 3 Last Modification Date: 06/14/95 12:56:43
ivan@78 4 Current Version: mdwt.c 1.5
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: All software, documentation, and related files in this distribution
ivan@78 9 are Copyright (c) 1994 Rice University
ivan@78 10
ivan@78 11 Permission is granted for use and non-profit distribution providing that this
ivan@78 12 notice be clearly maintained. The right to distribute any portion for profit
ivan@78 13 or as part of any commercial product is specifically reserved for the author.
ivan@78 14
ivan@78 15 Change History: Fixed code such that the result has the same dimension as the
ivan@78 16 input for 1D problems. Also, added some standard error checking.
ivan@78 17 Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
ivan@78 18
ivan@78 19 MATLAB gateway for MDWT.c, discrete wavelet transform
ivan@78 20 */
ivan@78 21 #include <math.h>
ivan@78 22 /*#include <malloc.h>*/
ivan@78 23 #include <stdio.h>
ivan@78 24 #include "mex.h"
ivan@78 25 #include "matrix.h"
ivan@78 26 #if !defined(_WIN32) && !defined(_WIN64)
ivan@78 27 #include <inttypes.h>
ivan@78 28 #endif
ivan@78 29 #define max(A,B) (A > B ? A : B)
ivan@78 30 #define min(A,B) (A < B ? A : B)
ivan@78 31 #define even(x) ((x & 1) ? 0 : 1)
ivan@78 32 #define isint(x) ((x - floor(x)) > 0.0 ? 0 : 1)
ivan@78 33 #define mat(a, i, j) (*(a + (m*(j)+i))) /* macro for matrix indices */
ivan@78 34
ivan@78 35 void mexFunction(const int nlhs,mxArray *plhs[],const int nrhs,const mxArray *prhs[])
ivan@78 36
ivan@78 37 {
ivan@78 38 double *x, *h, *y, *Lf, *Lr;
ivan@78 39 intptr_t m, n, h_col, h_row, lh, L, i, po2, j;
ivan@78 40 double mtest, ntest;
ivan@78 41
ivan@78 42 /* check for correct # of input variables */
ivan@78 43 if (nrhs>3){
ivan@78 44 mexErrMsgTxt("There are at most 3 input parameters allowed!");
ivan@78 45 return;
ivan@78 46 }
ivan@78 47 if (nrhs<2){
ivan@78 48 mexErrMsgTxt("There are at least 2 input parameters required!");
ivan@78 49 return;
ivan@78 50 }
ivan@78 51 x = mxGetPr(prhs[0]);
ivan@78 52 n = mxGetN(prhs[0]);
ivan@78 53 m = mxGetM(prhs[0]);
ivan@78 54 h = mxGetPr(prhs[1]);
ivan@78 55 h_col = mxGetN(prhs[1]);
ivan@78 56 h_row = mxGetM(prhs[1]);
ivan@78 57 if (h_col>h_row)
ivan@78 58 lh = h_col;
ivan@78 59 else
ivan@78 60 lh = h_row;
ivan@78 61 if (nrhs == 3){
ivan@78 62 L = (intptr_t) *mxGetPr(prhs[2]);
ivan@78 63 if (L < 0)
ivan@78 64 mexErrMsgTxt("The number of levels, L, must be a non-negative integer");
ivan@78 65 }
ivan@78 66 else /* Estimate L */ {
ivan@78 67 i=n;j=0;
ivan@78 68 while (even(i)){
ivan@78 69 i=(i>>1);
ivan@78 70 j++;
ivan@78 71 }
ivan@78 72 L=m;i=0;
ivan@78 73 while (even(L)){
ivan@78 74 L=(L>>1);
ivan@78 75 i++;
ivan@78 76 }
ivan@78 77 if(min(m,n) == 1)
ivan@78 78 L = max(i,j);
ivan@78 79 else
ivan@78 80 L = min(i,j);
ivan@78 81 if (L==0){
ivan@78 82 mexErrMsgTxt("Maximum number of levels is zero; no decomposition can be performed!");
ivan@78 83 return;
ivan@78 84 }
ivan@78 85 }
ivan@78 86 /* Check the ROW dimension of input */
ivan@78 87 if(m > 1){
ivan@78 88 mtest = (double) m/pow(2.0, (double) L);
ivan@78 89 if (!isint(mtest))
ivan@78 90 mexErrMsgTxt("The matrix row dimension must be of size m*2^(L)");
ivan@78 91 }
ivan@78 92 /* Check the COLUMN dimension of input */
ivan@78 93 if(n > 1){
ivan@78 94 ntest = (double) n/pow(2.0, (double) L);
ivan@78 95 if (!isint(ntest))
ivan@78 96 mexErrMsgTxt("The matrix column dimension must be of size n*2^(L)");
ivan@78 97 }
ivan@78 98 plhs[0] = mxCreateDoubleMatrix(m,n,mxREAL);
ivan@78 99 y = mxGetPr(plhs[0]);
ivan@78 100 plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
ivan@78 101 Lr = mxGetPr(plhs[1]);
ivan@78 102 *Lr = L;
ivan@78 103 MDWT(x, m, n, h, lh, L, y);
ivan@78 104 }
ivan@78 105 #ifdef __STDC__
ivan@78 106 MDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L, double *y)
ivan@78 107 #else
ivan@78 108 MDWT(x, m, n, h, lh, L, y)
ivan@78 109 double *x, *h, *y;
ivan@78 110 intptr_t m, n, lh, L;
ivan@78 111 #endif
ivan@78 112 {
ivan@78 113 double *h0, *h1, *ydummyl, *ydummyh, *xdummy;
ivan@78 114 int *prob;
ivan@78 115 intptr_t i, j;
ivan@78 116 intptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1;
ivan@78 117
ivan@78 118 xdummy = (double *)mxCalloc(max(m,n)+lh-1,sizeof(double));
ivan@78 119 ydummyl =(double *) (intptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78 120 ydummyh = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double));
ivan@78 121 h0 =(double *)(intptr_t)mxCalloc(lh,sizeof(double));
ivan@78 122 h1 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
ivan@78 123
ivan@78 124
ivan@78 125 /* analysis lowpass and highpass */
ivan@78 126 if (n==1){
ivan@78 127 n = m;
ivan@78 128 m = 1;
ivan@78 129 }
ivan@78 130 for (i=0; i<lh; i++){
ivan@78 131 h0[i] = h[lh-i-1];
ivan@78 132 h1[i] =h[i];
ivan@78 133 }
ivan@78 134 for (i=0; i<lh; i+=2)
ivan@78 135 h1[i] = -h1[i];
ivan@78 136
ivan@78 137 lhm1 = lh - 1;
ivan@78 138 actual_m = 2*m;
ivan@78 139 actual_n = 2*n;
ivan@78 140
ivan@78 141 /* main loop */
ivan@78 142 for (actual_L=1; actual_L <= L; actual_L++){
ivan@78 143 if (m==1)
ivan@78 144 actual_m = 1;
ivan@78 145 else{
ivan@78 146 actual_m = actual_m/2;
ivan@78 147 r_o_a = actual_m/2;
ivan@78 148 }
ivan@78 149 actual_n = actual_n/2;
ivan@78 150 c_o_a = actual_n/2;
ivan@78 151
ivan@78 152 /* go by rows */
ivan@78 153 for (ir=0; ir<actual_m; ir++){ /* loop over rows */
ivan@78 154 /* store in dummy variable */
ivan@78 155 for (i=0; i<actual_n; i++)
ivan@78 156 if (actual_L==1)
ivan@78 157 xdummy[i] = mat(x, ir, i);
ivan@78 158 else
ivan@78 159 xdummy[i] = mat(y, ir, i);
ivan@78 160 /* perform filtering lowpass and highpass*/
ivan@78 161 fpsconv(xdummy, actual_n, h0, h1, lhm1, ydummyl, ydummyh);
ivan@78 162 /* restore dummy variables in matrices */
ivan@78 163 ic = c_o_a;
ivan@78 164 for (i=0; i<c_o_a; i++){
ivan@78 165 mat(y, ir, i) = ydummyl[i];
ivan@78 166 mat(y, ir, ic++) = ydummyh[i];
ivan@78 167 }
ivan@78 168 }
ivan@78 169
ivan@78 170 /* go by columns in case of a 2D signal*/
ivan@78 171 if (m>1){
ivan@78 172 for (ic=0; ic<actual_n; ic++){ /* loop over column */
ivan@78 173 /* store in dummy variables */
ivan@78 174 for (i=0; i<actual_m; i++)
ivan@78 175 xdummy[i] = mat(y, i, ic);
ivan@78 176 /* perform filtering lowpass and highpass*/
ivan@78 177 fpsconv(xdummy, actual_m, h0, h1, lhm1, ydummyl, ydummyh);
ivan@78 178 /* restore dummy variables in matrix */
ivan@78 179 ir = r_o_a;
ivan@78 180 for (i=0; i<r_o_a; i++){
ivan@78 181 mat(y, i, ic) = ydummyl[i];
ivan@78 182 mat(y, ir++, ic) = ydummyh[i];
ivan@78 183 }
ivan@78 184 }
ivan@78 185 }
ivan@78 186 }
ivan@78 187 }
ivan@78 188
ivan@78 189 #ifdef __STDC__
ivan@78 190 fpsconv(double *x_in, intptr_t lx, double *h0, double *h1, intptr_t lhm1,
ivan@78 191 double *x_outl, double *x_outh)
ivan@78 192 #else
ivan@78 193 fpsconv(x_in, lx, h0, h1, lhm1, x_outl, x_outh)
ivan@78 194 double *x_in, *h0, *h1, *x_outl, *x_outh;
ivan@78 195 intptr_t lx, lhm1;
ivan@78 196 #endif
ivan@78 197
ivan@78 198 {
ivan@78 199 intptr_t i, j, ind;
ivan@78 200 double x0, x1;
ivan@78 201
ivan@78 202 for (i=lx; i < lx+lhm1; i++)
ivan@78 203 x_in[i] = *(x_in+(i-lx));
ivan@78 204 ind = 0;
ivan@78 205 for (i=0; i<(lx); i+=2){
ivan@78 206 x0 = 0;
ivan@78 207 x1 = 0;
ivan@78 208 for (j=0; j<=lhm1; j++){
ivan@78 209 x0 = x0 + x_in[i+j]*h0[lhm1-j];
ivan@78 210 x1 = x1 + x_in[i+j]*h1[lhm1-j];
ivan@78 211 }
ivan@78 212 x_outl[ind] = x0;
ivan@78 213 x_outh[ind++] = x1;
ivan@78 214 }
ivan@78 215 }