annotate util/Rice Wavelet Toolbox/midwt.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: midwt.c
ivan@78 3 Last Modification Date: 06/14/95 12:55:58
ivan@78 4 Current Version: midwt.c 1.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: 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 */
ivan@78 20 #include <math.h>
ivan@78 21 /*#include <malloc.h>*/
ivan@78 22 #include <stdio.h>
ivan@78 23 #include "mex.h"
ivan@78 24 #include "matrix.h"
ivan@78 25 #if !defined(_WIN32) && !defined(_WIN64)
ivan@78 26 #include <inttypes.h>
ivan@78 27 #endif
ivan@78 28 #define max(A,B) (A > B ? A : B)
ivan@78 29 #define min(A,B) (A < B ? A : B)
ivan@78 30 #define even(x) ((x & 1) ? 0 : 1)
ivan@78 31 #define isint(x) ((x - floor(x)) > 0.0 ? 0 : 1)
ivan@78 32
ivan@78 33
ivan@78 34
ivan@78 35 void mexFunction(int nlhs,mxArray *plhs[],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 y = 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 x = 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 MIDWT(x, m, n, h, lh, L, y);
ivan@78 104 }
ivan@78 105
ivan@78 106 #define mat(a, i, j) (*(a + (m*(j)+i))) /* macro for matrix indices */
ivan@78 107
ivan@78 108
ivan@78 109 #ifdef __STDC__
ivan@78 110 MIDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L, double *y)
ivan@78 111 #else
ivan@78 112 MIDWT(x, m, n, h, lh, L, y)
ivan@78 113 double *x, *h, *y;
ivan@78 114 intptr_t m, n, lh, L;
ivan@78 115 #endif
ivan@78 116 {
ivan@78 117 double *g0, *g1, *ydummyl, *ydummyh, *xdummy;
ivan@78 118 intptr_t i, j;
ivan@78 119 intptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1, lhhm1, sample_f;
ivan@78 120 xdummy = (double *)mxCalloc(max(m,n),sizeof(double));
ivan@78 121 ydummyl = (double *)mxCalloc(max(m,n)+lh/2-1,sizeof(double));
ivan@78 122 ydummyh = (double *)(intptr_t)mxCalloc(max(m,n)+lh/2-1,sizeof(double));
ivan@78 123 g0 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
ivan@78 124 g1 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
ivan@78 125
ivan@78 126 if (n==1){
ivan@78 127 n = m;
ivan@78 128 m = 1;
ivan@78 129 }
ivan@78 130 /* synthesis lowpass and highpass */
ivan@78 131 for (i=0; i<lh; i++){
ivan@78 132 g0[i] = h[i];
ivan@78 133 g1[i] = h[lh-i-1];
ivan@78 134 }
ivan@78 135 for (i=1; i<=lh; i+=2)
ivan@78 136 g1[i] = -g1[i];
ivan@78 137
ivan@78 138 lhm1 = lh - 1;
ivan@78 139 lhhm1 = lh/2 - 1;
ivan@78 140 /* 2^L */
ivan@78 141 sample_f = 1;
ivan@78 142 for (i=1; i<L; i++)
ivan@78 143 sample_f = sample_f*2;
ivan@78 144
ivan@78 145 if (m>1)
ivan@78 146 actual_m = m/sample_f;
ivan@78 147 else
ivan@78 148 actual_m = 1;
ivan@78 149 actual_n = n/sample_f;
ivan@78 150
ivan@78 151 for (i=0; i<(m*n); i++)
ivan@78 152 x[i] = y[i];
ivan@78 153
ivan@78 154 /* main loop */
ivan@78 155 for (actual_L=L; actual_L >= 1; actual_L--){
ivan@78 156 r_o_a = actual_m/2;
ivan@78 157 c_o_a = actual_n/2;
ivan@78 158
ivan@78 159 /* go by columns in case of a 2D signal*/
ivan@78 160 if (m>1){
ivan@78 161 for (ic=0; ic<actual_n; ic++){ /* loop over column */
ivan@78 162 /* store in dummy variables */
ivan@78 163 ir = r_o_a;
ivan@78 164 for (i=0; i<r_o_a; i++){
ivan@78 165 ydummyl[i+lhhm1] = mat(x, i, ic);
ivan@78 166 ydummyh[i+lhhm1] = mat(x, ir++, ic);
ivan@78 167 }
ivan@78 168 /* perform filtering lowpass and highpass*/
ivan@78 169 bpsconv(xdummy, r_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh);
ivan@78 170 /* restore dummy variables in matrix */
ivan@78 171 for (i=0; i<actual_m; i++)
ivan@78 172 mat(x, i, ic) = xdummy[i];
ivan@78 173 }
ivan@78 174 }
ivan@78 175 /* go by rows */
ivan@78 176 for (ir=0; ir<actual_m; ir++){ /* loop over rows */
ivan@78 177 /* store in dummy variable */
ivan@78 178 ic = c_o_a;
ivan@78 179 for (i=0; i<c_o_a; i++){
ivan@78 180 ydummyl[i+lhhm1] = mat(x, ir, i);
ivan@78 181 ydummyh[i+lhhm1] = mat(x, ir, ic++);
ivan@78 182 }
ivan@78 183 /* perform filtering lowpass and highpass*/
ivan@78 184 bpsconv(xdummy, c_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh);
ivan@78 185 /* restore dummy variables in matrices */
ivan@78 186 for (i=0; i<actual_n; i++)
ivan@78 187 mat(x, ir, i) = xdummy[i];
ivan@78 188 }
ivan@78 189 if (m==1)
ivan@78 190 actual_m = 1;
ivan@78 191 else
ivan@78 192 actual_m = actual_m*2;
ivan@78 193 actual_n = actual_n*2;
ivan@78 194 }
ivan@78 195 }
ivan@78 196
ivan@78 197 #ifdef __STDC__
ivan@78 198 bpsconv(double *x_out, intptr_t lx, double *g0, double *g1, intptr_t lhm1,
ivan@78 199 intptr_t lhhm1, double *x_inl, double *x_inh)
ivan@78 200 #else
ivan@78 201 bpsconv(x_out, lx, g0, g1, lhm1, lhhm1, x_inl, x_inh)
ivan@78 202 double *x_inl, *x_inh, *g0, *g1, *x_out;
ivan@78 203 intptr_t lx, lhm1, lhhm1;
ivan@78 204 #endif
ivan@78 205 {
ivan@78 206 intptr_t i, j, ind, tj;
ivan@78 207 double x0, x1;
ivan@78 208
ivan@78 209 for (i=lhhm1-1; i > -1; i--){
ivan@78 210 x_inl[i] = x_inl[lx+i];
ivan@78 211 x_inh[i] = x_inh[lx+i];
ivan@78 212 }
ivan@78 213 ind = 0;
ivan@78 214 for (i=0; i<(lx); i++){
ivan@78 215 x0 = 0;
ivan@78 216 x1 = 0;
ivan@78 217 tj = -2;
ivan@78 218 for (j=0; j<=lhhm1; j++){
ivan@78 219 tj+=2;
ivan@78 220 x0 = x0 + x_inl[i+j]*g0[lhm1-1-tj] + x_inh[i+j]*g1[lhm1-1-tj] ;
ivan@78 221 x1 = x1 + x_inl[i+j]*g0[lhm1-tj] + x_inh[i+j]*g1[lhm1-tj] ;
ivan@78 222 }
ivan@78 223 x_out[ind++] = x0;
ivan@78 224 x_out[ind++] = x1;
ivan@78 225 }
ivan@78 226 }