Mercurial > hg > smallbox

--- a/SMALLboxSetup.m	Fri Mar 25 14:01:50 2011 +0000
+++ b/SMALLboxSetup.m	Mon Mar 28 11:19:41 2011 +0100
@@ -103,18 +103,9 @@
             'will not work.\n\n']);
     end

-    cd([Sparco_path, FS, 'sparco-1.2', FS, 'tools' ,FS, 'rwt'])
-    fprintf('Compiling the Rice Wavelet Toolbox MEX interfaces...');
-    try
-        if exist('mdwt'  ,'file')~=3, mex mdwt.c   mdwt_r.c;   end
-        if exist('midwt' ,'file')~=3, mex midwt.c  midwt_r.c;  end
-        if exist('mrdwt' ,'file')~=3, mex mrdwt.c  mrdwt_r.c;  end
-        if exist('mirdwt','file')~=3, mex mirdwt.c mirdwt_r.c; end
-        fprintf('SPARCO Installation Successful!\n');
-    catch
-        warning('Could not compile Rice Wavelet Toolbox MEX interfaces.');
-    end
+
     cd(SMALL_path);
+    fprintf('SPARCO Installation Successful!\n');
 else
     fprintf('\n ******************************************************************');
     fprintf('\n\n SPARCO and Rice Wavelet Toolbox are already installed');
@@ -354,28 +345,39 @@
     fprintf('\n\n matlab_midi (http://www.kenschutte.com/midi/) is already installed');
 end

-cd([pwd,FS,'util',FS,'ksvd utils']);
+cd([SMALL_path,FS,'util',FS,'ksvd utils']);
 make
 cd(SMALL_path);
 %%
-
+cd([SMALL_path, FS, 'util', FS, 'Rice Wavelet Toolbox'])
+    fprintf('Compiling the Rice Wavelet Toolbox MEX interfaces...');
+    try
+        if exist('mdwt'  ,'file')~=3, mex mdwt.c;   end
+        if exist('midwt' ,'file')~=3, mex midwt.c;  end
+        if exist('mrdwt' ,'file')~=3, mex mrdwt.c;  end
+        if exist('mirdwt','file')~=3, mex mirdwt.c; end
+        fprintf('Rice Wavelet Toolbox Installation Successful!\n\n');
+    catch
+        warning('Could not compile Rice Wavelet Toolbox MEX interfaces.\n');
+    end
+cd(SMALL_path);

 fprintf('\n ******************************************************************');
 fprintf('\n\n Initialising SMALLbox Examples Setup');

-    % Need to do a bit of temporary housekeeping first.
-    cd(SMALL_path);
-  try
-    cd(['Problems',FS,'private']);
-    if exist('addtocols'  ,'file')~=3,
-    fprintf('\n Compiling MEX interfaces for SMALL examples \n');
-    make;
-    end
-    fprintf('\n SMALLbox Problems Installation Successful! \n');
-  catch
-    fprintf('\n SMALLbox Problems Installation Failed \n');
-  end
-    cd(SMALL_path);
+%     % Need to do a bit of temporary housekeeping first.
+%     cd(SMALL_path);
+%   try
+%     cd(['Problems',FS,'private']);
+%     if exist('addtocols'  ,'file')~=3,
+%     fprintf('\n Compiling MEX interfaces for SMALL examples \n');
+%     make;
+%     end
+%     fprintf('\n SMALLbox Problems Installation Successful! \n');
+%   catch
+%     fprintf('\n SMALLbox Problems Installation Failed \n');
+%   end
+%     cd(SMALL_path);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/Contents.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,18 @@
+% Rice Wavelet Toolbox
+% Version 2.4   Dec 2002
+%
+% Wavelet Transforms.
+%     mdwt - Discrete orthogonal wavelet transform using Mallat alg. (1D and 2D)
+%     midwt - Inverse discrete orthogonal wavelet transform
+%     mrdwt - Undecimated (redundant) discrete wavelet transform (1D and 2D)
+%     mirdwt - Inverse undecimated discrete wavelet transform
+%     daubcqf - Daubecheis filter coefficients
+%
+% Wavelet Domain Processing.
+%     denoise - Denoise signal and images by thresholding wavelet coefficients
+%     HardTh - Hard thresholding
+%     SoftTh - Soft thresholding
+%
+% Other.
+%     makesig - Create Donoho-Johnstone test signals
+%     compile - compile the Rice Wavelet toolbox
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/HardTh.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,71 @@
+function  x = HardTh(y,thld)
+%    x = HardTh(y,thld);
+%
+%    HARDTH hard thresholds the input signal y with the threshold value
+%    thld.
+%
+%    Input:
+%       y    : 1D or 2D signal to be thresholded
+%       thld : threshold value
+%
+%    Output:
+%       x : Hard thresholded output (x = (abs(y)>thld).*y)
+%
+%  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
+%  Cut-and-paste the example you want to run to a new file
+%  called ex.m, for example. Delete out the % at the beginning
+%  of each line in ex.m (Can use search-and-replace in your editor
+%  to replace it with a space). Type 'ex' in matlab and hit return.
+%
+%
+%    Example:
+%       y = makesig('WernerSorrows',8);
+%       thld = 1;
+%       x = HardTh(y,thld)
+%       x = 1.5545 5.3175 0 1.6956  -1.2678 0 1.7332 0
+%
+%    See also: SoftTh
+%
+
+%File Name: HardTh.m
+%Last Modification Date: 8/15/95	17:49:37
+%Current Version: HardTh.m	2.4
+%File Creation Date: Mon Jan 31 09:42:50 1994
+%Author: Haitao Guo  <harry@jazz.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Haitao Guo, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+
+x = (abs(y) > thld).*y;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/INSTALL	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,26 @@
+In order to install this distribution of Rice Wavelet Tools version 2.4
+released - <Dec 1 2002>
+
+1. Properly set up your system to create MEX-files. Please refer to the
+   "Matlab Application Program Guide" to properly set up of your matlab
+   and C-compiler to be able to compile C-mex files on your system.
+   All reference documentations are available on the MathWorks web page:
+   www.mathworks.com
+
+2. Make a toolbox directory and uncompress/extract all the files.
+   For example, in the unix environment,
+
+        gunzip rwt.tar.gz
+        tar xvf rwt.tar
+
+3. Run MATLAB and change to the temporary directory containing the files.
+
+4. Compile the toolbox by executing the Matlab command
+
+        compile
+
+5. Add the toolbox directory to your Matlab path.
+
+6. For further instructions, please refer to the README file.
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/LICENSE	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,38 @@
+This "rice-wlet-tools", version 2.4
+Released - <Dec 1 2002>
+
+CONDITIONS FOR USE:
+Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+
+This software is distributed and licensed to you on a non-exclusive
+basis, free-of-charge. Redistribution and use in source and binary forms,
+with or without modification, are permitted provided that the following
+conditions are met:
+
+1. Redistribution of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistribution in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. All advertising materials mentioning features or use of this software
+   must display the following acknowledgment: This product includes
+   software developed by Rice University, Houston, Texas and its contributors.
+4. Neither the name of the University nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+For information on commercial licenses, contact Rice University's Office of
+Technology Transfer at techtran@rice.edu or (713) 348-6173
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/README	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,38 @@
+Rice Wavelet Tools version 2.4
+Released - <Dec 1 2002>
+
+INSTALLATION:
+To install this distribution of Rice Wavelet Tools see the INSTALL file.
+
+SOURCE:
+     www.dsp.rice.edu/software/RWT
+
+EMAIL:
+For bug reports and questions, send email to wlet-tools@rice.edu
+
+CONDITIONS FOR USE:
+See the LICENSE file
+
+TOOLBOX FUNCTIONS:
+
+ Wavelet Transforms
+     mdwt - Discrete orthogonal wavelet transform using the Mallat algorithm (1D and 2D)
+     midwt - Inverse discrete orthogonal wavelet transform
+     mrdwt - Undecimated (redundant) discrete wavelet transform (1D and 2D)
+     mirdwt - Inverse undecimated discrete wavelet transform
+     daubcqf - Daubechies filter coefficients
+
+ Wavelet Domain Processing
+     denoise - Denoise signals and images by thresholding wavelet coefficients
+     HardTh - Hard thresholding
+     SoftTh - Soft thresholding
+
+ Other
+     makesig - Create Donoho-Johnstone test signals
+     compile - Compile the Rice Wavelet Toolbox
+
+Functions omitted in this version of toolbox can be found in
+version 2.01 at www.dsp.rice.edu/software/RWT/RWT-2.01.tar.Z
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/README.sparco	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,21 @@
+What are these funny ".mhelp" files?!
+-------------------------------------
+
+Some users of Matlab 7.0.4 have reported that the Rice Wavelet Toolbox
+MEX interfaces aren't being found.  Instead, the .m files (which only
+contain the help text) are taking precedence, and function calls to
+these routines cause an error.
+
+A simple workaround is to move these .m help files so that Matlab
+doesn't have to choose between them and the MEX versions.  We've
+renamed these files:
+
+mdwt.m    to  mdwt.mhelp
+midwt.m   to  midwt.mhelpm
+mirdwt.m  to  mirdwt.mhelpm
+mrdwt.dll to  mrdwt.mhelpm
+
+-- The Sparco development team
+   (Michael P. Friedlander and Ewout van den Berg)
+$Id: README.sparco 830 2008-03-25 22:28:04Z mpf $
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/SoftTh.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,79 @@
+function  x = SoftTh(y,thld)
+%    x = SoftTh(y,thld);
+%
+%    SOFTTH soft thresholds the input signal y with the threshold value
+%    thld.
+%
+%    Input:
+%       y    : 1D or 2D signal to be thresholded
+%       thld : Threshold value
+%
+%    Output:
+%       x : Soft thresholded output (x = sign(y)(|y|-thld)_+)
+%
+%  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
+%  Cut-and-paste the example you want to run to a new file
+%  called ex.m, for example. Delete out the % at the beginning
+%  of each line in ex.m (Can use search-and-replace in your editor
+%  to replace it with a space). Type 'ex' in matlab and hit return.
+%
+%
+%    Example:
+%       y = makesig('Doppler',8);
+%       thld = 0.2;
+%       x = SoftTh(y,thld)
+%       x = 0 0 0 -0.0703 0 0.2001 0.0483 0
+%
+%    See also: HardTh
+%
+%    Reference:
+%       "De-noising via Soft-Thresholding" Tech. Rept. Statistics,
+%       Stanford, 1992. D.L. Donoho.
+%
+
+%File Name: SoftTh.m
+%Last Modification Date: 8/15/95	17:49:48
+%Current Version: SoftTh.m	2.4
+%File Creation Date: Mon Mar  7 10:38:45 1994
+%Author: Haitao Guo  <harry@jazz.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Haitao Guo, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+%
+%Change History:
+%
+
+x = abs(y);
+x = sign(y).*(x >= thld).*(x - thld);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/compile.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,42 @@
+%    COMPILE compiles the c files and generates mex files.
+%
+
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Haitao Guo, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+
+mex mdwt.c
+mex midwt.c
+mex mrdwt.c
+mex mirdwt.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/daubcqf.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,106 @@
+function [h_0,h_1] = daubcqf(N,TYPE)
+%    [h_0,h_1] = daubcqf(N,TYPE);
+%
+%    Function computes the Daubechies' scaling and wavelet filters
+%    (normalized to sqrt(2)).
+%
+%    Input:
+%       N    : Length of filter (must be even)
+%       TYPE : Optional parameter that distinguishes the minimum phase,
+%              maximum phase and mid-phase solutions ('min', 'max', or
+%              'mid'). If no argument is specified, the minimum phase
+%              solution is used.
+%
+%    Output:
+%       h_0 : Minimal phase Daubechies' scaling filter
+%       h_1 : Minimal phase Daubechies' wavelet filter
+%
+%    Example:
+%       N = 4;
+%       TYPE = 'min';
+%       [h_0,h_1] = daubcqf(N,TYPE)
+%       h_0 = 0.4830 0.8365 0.2241 -0.1294
+%       h_1 = 0.1294 0.2241 -0.8365 0.4830
+%
+%    Reference: "Orthonormal Bases of Compactly Supported Wavelets",
+%                CPAM, Oct.89
+%
+
+%File Name: daubcqf.m
+%Last Modification Date: 01/02/96	15:12:57
+%Current Version: daubcqf.m	2.4
+%File Creation Date: 10/10/88
+%Author: Ramesh Gopinath  <ramesh@dsp.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Ramesh Gopinath, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+
+if(nargin < 2),
+  TYPE = 'min';
+end;
+if(rem(N,2) ~= 0),
+  error('No Daubechies filter exists for ODD length');
+end;
+K = N/2;
+a = 1;
+p = 1;
+q = 1;
+h_0 = [1 1];
+for j  = 1:K-1,
+  a = -a * 0.25 * (j + K - 1)/j;
+  h_0 = [0 h_0] + [h_0 0];
+  p = [0 -p] + [p 0];
+  p = [0 -p] + [p 0];
+  q = [0 q 0] + a*p;
+end;
+q = sort(roots(q));
+qt = q(1:K-1);
+if TYPE=='mid',
+  if rem(K,2)==1,
+    qt = q([1:4:N-2 2:4:N-2]);
+  else
+    qt = q([1 4:4:K-1 5:4:K-1 N-3:-4:K N-4:-4:K]);
+  end;
+end;
+h_0 = conv(h_0,real(poly(qt)));
+h_0 = sqrt(2)*h_0/sum(h_0); 	%Normalize to sqrt(2);
+if(TYPE=='max'),
+  h_0 = fliplr(h_0);
+end;
+if(abs(sum(h_0 .^ 2))-1 > 1e-4)
+  error('Numerically unstable for this value of "N".');
+end;
+h_1 = rot90(h_0,2);
+h_1(1:2:N)=-h_1(1:2:N);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/denoise.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,245 @@
+function [xd,xn,option] = denoise(x,h,type,option)
+%    [xd,xn,option] = denoise(x,h,type,option);
+%
+%    DENOISE is a generic program for wavelet based denoising.
+%    The program will denoise the signal x using the 2-band wavelet
+%    system described by the filter h using either the traditional
+%    discrete wavelet transform (DWT) or the linear shift invariant
+%    discrete wavelet transform (also known as the undecimated DWT
+%    (UDWT)).
+%
+%    Input:
+%       x         : 1D or 2D signal to be denoised
+%       h         : Scaling filter to be applied
+%       type      : Type of transform (Default: type = 0)
+%                   0 --> Discrete wavelet transform (DWT)
+%                   1 --> Undecimated DWT (UDWT)
+%       option    : Default settings is marked with '*':
+%                   *type = 0 --> option = [0 3.0 0 0 0 0]
+%                   type = 1 --> option = [0 3.6 0 1 0 0]
+%       option(1) : Whether to threshold low-pass part
+%                   0 --> Don't threshold low pass component
+%                   1 --> Threshold low pass component
+%       option(2) : Threshold multiplier, c. The threshold is
+%                   computed as:
+%                     thld = c*MAD(noise_estimate)).
+%                   The default values are:
+%                     c = 3.0 for the DWT based denoising
+%                     c = 3.6 for the UDWT based denoising
+%       option(3) : Type of variance estimator
+%                   0 --> MAD (mean absolute deviation)
+%                   1 --> STD (classical numerical std estimate)
+%       option(4) : Type of thresholding
+%                   0 --> Soft thresholding
+%                   1 --> Hard thresholding
+%       option(5) : Number of levels, L, in wavelet decomposition. By
+%                   setting this to the default value '0' a maximal
+%                   decomposition is used.
+%       option(6) : Actual threshold to use (setting this to
+%                   anything but 0 will mean that option(3)
+%                   is ignored)
+%
+%    Output:
+%       xd     : Estimate of noise free signal
+%       xn     : The estimated noise signal (x-xd)
+%       option : A vector of actual parameters used by the
+%                program. The vector is configured the same way as
+%                the input option vector with one added element
+%                option(7) = type.
+%
+%  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
+%  Cut-and-paste the example you want to run to a new file
+%  called ex.m, for example. Delete out the % at the beginning
+%  of each line in ex.m (Can use search-and-replace in your editor
+%  to replace it with a space). Type 'ex' in matlab and hit return.
+%
+%    Example 1:
+%       h = daubcqf(6); [s,N] = makesig('Doppler'); n = randn(1,N);
+%       x = s + n/10;     % (approximately 10dB SNR)
+%       figure;plot(x);hold on;plot(s,'r');
+%
+%       %Denoise x with the default method based on the DWT
+%       [xd,xn,opt1] = denoise(x,h);
+%       figure;plot(xd);hold on;plot(s,'r');
+%
+%       %Denoise x using the undecimated (LSI) wavelet transform
+%       [yd,yn,opt2] = denoise(x,h,1);
+%       figure;plot(yd);hold on;plot(s,'r');
+%
+% Example 2: (on an image)
+%      h = daubcqf(6);  load lena;
+%      noisyLena = lena + 25 * randn(size(lena));
+%      figure; colormap(gray); imagesc(lena); title('Original Image');
+%       figure; colormap(gray); imagesc(noisyLena); title('Noisy Image');
+%       Denoise lena with the default method based on the DWT
+%      [denoisedLena,xn,opt1] = denoise(noisyLena,h);
+%      figure; colormap(gray); imagesc(denoisedLena); title('denoised Image');
+%
+%
+%    See also: mdwt, midwt, mrdwt, mirdwt, SoftTh, HardTh, setopt
+%
+
+%File Name: denoise.m
+%Last Modification Date: 04/15/97	10:44:28
+%Current Version: denoise.m	2.4
+%File Creation Date: Mon Feb 20 08:33:15 1995
+%Author: Jan Erik Odegard  <odegard@ece.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Jan Erik Odegard, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+%
+%Change History: Fixed output of function and an error in the computation
+%                of the threshold for redundant denoising.
+%                <Jan Erik Odegard> <Mon Jul 31, 1995>
+%
+%                This code is composed of several of our old codes for
+%                wavelet based denoising. In an effort to make the mess
+%                more manageable we decided to create on code that would
+%                handle all the various wavelet based denoising methods.
+%                However, only time will show (as we discover new and
+%                improved forms of denoising) if we can succeed in our goals.
+%                <Jan Erik Odegard> <Thu May 11, 1995>
+%
+
+if(nargin < 2)
+  error('You need to provide at least 2 inputs: x and h');
+end;
+if(nargin < 3),
+  type = 0;
+  option = [];
+elseif(nargin < 4)
+  option = [];
+end;
+if(isempty(type)),
+  type = 0;
+end;
+if(type == 0),
+  default_opt = [0 3.0 0 0 0 0];
+elseif(type == 1),
+  default_opt = [0 3.6 0 1 0 0];
+else,
+  error(['Unknown denoising method',10,...
+	  'If it is any good we need to have a serious talk :-)']);
+end;
+option = setopt(option,default_opt);
+[mx,nx] = size(x);
+dim = min(mx,nx);
+if(dim == 1),
+  n = max(mx,nx);
+else,
+  n = dim;
+end;
+if(option(5) == 0),
+  L = floor(log2(n));
+else
+  L = option(5);
+end;
+if(type == 0), 			% Denoising by DWT
+  xd = mdwt(x,h,L);
+  if (option(6) == 0),
+    tmp = xd(floor(mx/2)+1:mx,floor(nx/2)+1:nx);
+    if(option(3) == 0),
+      thld = option(2)*median(abs(tmp(:)))/.67;
+    elseif(option(3) == 1),
+      thld = option(2)*std(tmp(:));
+    else
+      error('Unknown threshold estimator, Use either MAD or STD');
+    end;
+  else,
+    thld = option(6);
+  end;
+  if(dim == 1)
+    ix = 1:n/(2^L);
+    ykeep = xd(ix);
+  else
+    ix = 1:mx/(2^L);
+    jx = 1:nx/(2^L);
+    ykeep = xd(ix,jx);
+  end;
+  if(option(4) == 0),
+    xd = SoftTh(xd,thld);
+  elseif(option(4) == 1),
+    xd = HardTh(xd,thld);
+  else,
+    error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
+  end;
+  if (option(1) == 0),
+    if(dim == 1),
+      xd(ix) = ykeep;
+    else,
+      xd(ix,jx) = ykeep;
+    end;
+  end;
+  xd = midwt(xd,h,L);
+elseif(type == 1), 			% Denoising by UDWT
+  [xl,xh] = mrdwt(x,h,L);
+  if(dim == 1),
+    c_offset = 1;
+  else,
+    c_offset = 2*nx + 1;
+  end;
+  if (option(6) == 0),
+    tmp = xh(:,c_offset:c_offset+nx-1);
+    if(option(3) == 0),
+      thld = option(2)*median(abs(tmp(:)))/.67;
+    elseif(option(3) == 1),
+      thld = option(2)*std(tmp(:));
+    else
+      error('Unknown threshold estimator, Use either MAD or STD');
+    end;
+  else,
+    thld = option(6);
+  end;
+  if(option(4) == 0),
+    xh = SoftTh(xh,thld);
+    if(option(1) == 1),
+      xl = SoftTh(xl,thld);
+    end;
+  elseif(option(4) == 1),
+    xh = HardTh(xh,thld);
+    if(option(1) == 1),
+      xl = HardTh(xl,thld);
+    end;
+  else,
+    error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
+  end;
+  xd = mirdwt(xl,xh,h,L);
+else, 					% Denoising by unknown method
+  error(['Unknown denoising method',10,...
+         'If it is any good we need to have a serious talk :-)']);
+end;
+option(6) = thld;
+option(7) = type;
+xn = x - xd;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/makesig.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,211 @@
+function [x,N] = makesig(SigName,N)
+% [x,N] = makesig(SigName,N) Creates artificial test signal identical to the
+%     standard test signals proposed and used by D. Donoho and I. Johnstone
+%     in WaveLab (- a matlab toolbox developed by Donoho et al. the statistics
+%     department at Stanford University).
+%
+%    Input:  SigName - Name of the desired signal (Default 'all')
+%                        'AllSig' (Returns a matrix with all the signals)
+%                        'HeaviSine'
+%                        'Bumps'
+%                        'Blocks'
+%                        'Doppler'
+%                        'Ramp'
+%                        'Cusp'
+%                        'Sing'
+%                        'HiSine'
+%                        'LoSine'
+%                        'LinChirp'
+%                        'TwoChirp'
+%                        'QuadChirp'
+%                        'MishMash'
+%                        'Werner Sorrows' (Heisenberg)
+%                        'Leopold' (Kronecker)
+%            N       - Length in samples of the desired signal (Default 512)
+%
+%    Output: x   - vector/matrix of test signals
+%            N   - length of signal returned
+%
+%    See also:
+%
+%    References:
+%            WaveLab can be accessed at
+%            www_url: http://playfair.stanford.edu/~wavelab/
+%            Also see various articles by D.L. Donoho et al. at
+%            web_url: http://playfair.stanford.edu/
+
+%File Name: makesig.m
+%Last Modification Date: 08/30/95	15:52:03
+%Current Version: makesig.m	2.4
+%File Creation Date: Thu Jun  8 10:31:11 1995
+%Author: Jan Erik Odegard  <odegard@ece.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Jan Erik Odegard, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+%
+%Change History: This m-file is a copy of the  code provided with WaveLab
+%                customized to be consistent with RWT.
+%                Jan Erik Odegard <odegard@ece.rice.edu> Thu Jun  8 1995
+%
+
+if(nargin < 1)
+  SigName = 'AllSig';
+  N = 512;
+elseif(nargin == 1)
+  N = 512;
+end;
+t = (1:N) ./N;
+x = [];
+y = [];
+if(strcmp(SigName,'HeaviSine') | strcmp(SigName,'AllSig')),
+  y = 4.*sin(4*pi.*t);
+  y = y - sign(t - .3) - sign(.72 - t);
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'Bumps') | strcmp(SigName,'AllSig')),
+  pos = [ .1 .13 .15 .23 .25 .40 .44 .65  .76 .78 .81];
+  hgt = [ 4  5   3   4  5  4.2 2.1 4.3  3.1 5.1 4.2];
+  wth = [.005 .005 .006 .01 .01 .03 .01 .01  .005 .008 .005];
+  y = zeros(size(t));
+  for j =1:length(pos)
+    y = y + hgt(j)./( 1 + abs((t - pos(j))./wth(j))).^4;
+  end
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'Blocks') | strcmp(SigName,'AllSig')),
+  pos = [ .1 .13 .15 .23 .25 .40 .44 .65  .76 .78 .81];
+  hgt = [4 (-5) 3 (-4) 5 (-4.2) 2.1 4.3  (-3.1) 2.1 (-4.2)];
+  y = zeros(size(t));
+  for j=1:length(pos)
+    y = y + (1 + sign(t-pos(j))).*(hgt(j)/2) ;
+  end
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'Doppler') | strcmp(SigName,'AllSig')),
+  y = sqrt(t.*(1-t)).*sin((2*pi*1.05) ./(t+.05));
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'Ramp') | strcmp(SigName,'AllSig')),
+  y = t - (t >= .37);
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'Cusp') | strcmp(SigName,'AllSig')),
+  y = sqrt(abs(t - .37));
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'Sing') | strcmp(SigName,'AllSig')),
+  k = floor(N * .37);
+  y = 1 ./abs(t - (k+.5)/N);
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'HiSine') | strcmp(SigName,'AllSig')),
+  y = sin( pi * (N * .6902) .* t);
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'LoSine') | strcmp(SigName,'AllSig')),
+  y = sin( pi * (N * .3333) .* t);
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'LinChirp') | strcmp(SigName,'AllSig')),
+  y = sin(pi .* t .* ((N .* .125) .* t));
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'TwoChirp') | strcmp(SigName,'AllSig')),
+  y = sin(pi .* t .* (N .* t)) + sin((pi/3) .* t .* (N .* t));
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'QuadChirp') | strcmp(SigName,'AllSig')),
+  y = sin( (pi/3) .* t .* (N .* t.^2));
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'MishMash') | strcmp(SigName,'AllSig')),
+  % QuadChirp + LinChirp + HiSine
+  y = sin( (pi/3) .* t .* (N .* t.^2)) ;
+  y = y +  sin( pi * (N * .6902) .* t);
+  y = y +  sin(pi .* t .* (N .* .125 .* t));
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'WernerSorrows') | strcmp(SigName,'AllSig')),
+  y = sin( pi .* t .* (N/2 .* t.^2)) ;
+  y = y +  sin( pi * (N * .6902) .* t);
+  y = y +  sin(pi .* t .* (N .* t));
+  pos = [ .1 .13 .15 .23 .25 .40 .44 .65  .76 .78 .81];
+  hgt = [ 4  5   3   4  5  4.2 2.1 4.3  3.1 5.1 4.2];
+  wth = [.005 .005 .006 .01 .01 .03 .01 .01  .005 .008 .005];
+  for j =1:length(pos)
+    y = y + hgt(j)./( 1 + abs((t - pos(j))./wth(j))).^4;
+  end
+end;
+x = [x;y];
+y = [];
+if(strcmp(SigName,'Leopold') | strcmp(SigName,'AllSig')),
+  y = (t == floor(.37 * N)/N); 		% Kronecker
+end;
+x = [x;y];
+y = [];
+
+%  disp(sprintf('MakeSignal: I don*t recognize << %s>>',SigName))
+%  disp('Allowable SigNames are:')
+%  disp('AllSig'),
+%  disp('HeaviSine'),
+%  disp('Bumps'),
+%  disp('Blocks'),
+%  disp('Doppler'),
+%  disp('Ramp'),
+%  disp('Cusp'),
+%  disp('Crease'),
+%  disp('Sing'),
+%  disp('HiSine'),
+%  disp('LoSine'),
+%  disp('LinChirp'),
+%  disp('TwoChirp'),
+%  disp('QuadChirp'),
+%  disp('MishMash'),
+%  disp('WernerSorrows'),
+%  disp('Leopold'),
+%end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mdwt.c	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,215 @@
+/*
+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;
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mdwt.mhelp	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,132 @@
+function [y,L] = mdwt(x,h,L);
+%    [y,L] = mdwt(x,h,L);
+%
+%    Function computes the discrete wavelet transform y for a 1D or 2D input
+%    signal x using the scaling filter h.
+%
+%    Input:
+%	x : finite length 1D or 2D signal (implicitly periodized)
+%       h : scaling filter
+%       L : number of levels. In the case of a 1D signal, length(x) must be
+%           divisible by 2^L; in the case of a 2D signal, the row and the
+%           column dimension must be divisible by 2^L. If no argument is
+%           specified, a full DWT is returned for maximal possible L.
+%
+%    Output:
+%       y : the wavelet transform of the signal
+%           (see example to understand the coefficients)
+%       L : number of decomposition levels
+%
+%    1D Example:
+%       x = makesig('LinChirp',8);
+%       h = daubcqf(4,'min');
+%       L = 2;
+%       [y,L] = mdwt(x,h,L)
+%
+%    1D Example's  output and explanation:
+%
+%       y = [1.1097 0.8767 0.8204 -0.5201 -0.0339 0.1001 0.2201 -0.1401]
+%       L = 2
+%
+%    The coefficients in output y are arranged as follows
+%
+%       y(1) and y(2) : Scaling coefficients (lowest frequency)
+%       y(3) and y(4) : Band pass wavelet coefficients
+%       y(5) to y(8)  : Finest scale wavelet coefficients (highest frequency)
+%
+%    2D Example:
+%
+%       load test_image
+%       h = daubcqf(4,'min');
+%       L = 1;
+%       [y,L] = mdwt(test_image,h,L);
+%
+%    2D Example's  output and explanation:
+%
+%       The coefficients in y are arranged as follows.
+%
+%              .------------------.
+%              |         |        |
+%              |    4    |   2    |
+%              |         |        |
+%              |   L,L   |   H,L  |
+%              |         |        |
+%              --------------------
+%              |         |        |
+%              |    3    |   1    |
+%              |         |        |
+%              |   L,H   |  H,H   |
+%              |         |        |
+%              `------------------'
+%
+%       where
+%            1 : High pass vertically and high pass horizontally
+%            2 : Low pass vertically and high pass horizontally
+%            3 : High pass vertically and low  pass horizontally
+%            4 : Low pass vertically and Low pass horizontally
+%                (scaling coefficients)
+%
+%
+%
+%
+%    See also: midwt, mrdwt, mirdwt
+%
+
+%File Name: mdwt.m
+%Last Modification Date: 08/07/95	15:13:25
+%Current Version: mdwt.m	2.4
+%File Creation Date: Wed Oct 19 10:51:58 1994
+%Author: Markus Lang  <lang@jazz.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Markus Lang, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+%
+%Change History:
+%
+%Modification #1
+%Mon Aug  7 11:42:11 CDT 1995
+%Rebecca Hindman <hindman@ece.rice.edu>
+%Added L to function line so that it can be displayed as an output
+%
+%Change History:
+%
+%Modification #1
+%Thu Mar  2 13:07:11 CDT 2000
+%Ramesh Neelamani<neelsh@ece.rice.edu>
+%Revamped the help file
+%
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mdwt_r.c	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,181 @@
+/*
+File Name: MDWT.c
+Last Modification Date:	06/14/95	13:15:44
+Current Version: MDWT.c	2.4
+File Creation Date: Wed Oct 19 10:51:58 1994
+Author: Markus Lang  <lang@jazz.rice.edu>
+
+Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+Created by Markus Lang, Department of ECE, Rice University.
+
+This software is distributed and licensed to you on a non-exclusive
+basis, free-of-charge. Redistribution and use in source and binary forms,
+with or without modification, are permitted provided that the following
+conditions are met:
+
+1. Redistribution of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistribution in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. All advertising materials mentioning features or use of this software
+   must display the following acknowledgment: This product includes
+   software developed by Rice University, Houston, Texas and its contributors.
+4. Neither the name of the University nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+For information on commercial licenses, contact Rice University's Office of
+Technology Transfer at techtran@rice.edu or (713) 348-6173
+
+Change History: Fixed the code such that 1D vectors passed to it can be in
+                either passed as a row or column vector. Also took care of
+		the code such that it will compile with both under standard
+		C compilers as well as for ANSI C compilers
+		Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
+
+%y = mdwt(x,h,L);
+%
+% function computes the discrete wavelet transform y for a 1D or 2D input
+% signal x.
+%
+%    Input:
+%	x    : finite length 1D or 2D signal (implicitely periodized)
+%       h    : scaling filter
+%       L    : number of levels. in case of a 1D signal length(x) must be
+%              divisible by 2^L; in case of a 2D signal the row and the
+%              column dimension must be divisible by 2^L.
+%
+% see also: midwt, mrdwt, mirdwt
+*/
+
+#include <math.h>
+#include <stdio.h>
+#include <inttypes.h>
+
+#define max(A,B) (A > B ? A : B)
+#define mat(a, i, j) (*(a + (m*(j)+i)))  /* macro for matrix indices */
+
+#ifdef __STDC__
+MDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L, double *y)
+#else
+MDWT(x, m, n, h, lh, L, y)
+double *x, *h, *y;
+uintptr_t m, n, lh, L;
+#endif
+{
+  double  *h0, *h1, *ydummyl, *ydummyh, *xdummy;
+  int *prob;
+  uintptr_t i, j;
+  uintptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1;
+
+  xdummy = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyl =(double *) (uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  h0 =(double *)(uintptr_t)mxCalloc(lh,sizeof(double));
+  h1 = (double *)(uintptr_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, uintptr_t lx, double *h0, double *h1, uintptr_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;
+uintptr_t lx, lhm1;
+#endif
+
+{
+ uintptr_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;
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/midwt.c	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,226 @@
+/*
+File Name: midwt.c
+Last Modification Date:	06/14/95	12:55:58
+Current Version: midwt.c	1.4
+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
+
+*/
+#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)
+
+
+
+void mexFunction(int nlhs,mxArray *plhs[],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;
+  }
+  y = 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);
+  x = mxGetPr(plhs[0]);
+  plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
+  Lr = mxGetPr(plhs[1]);
+  *Lr = L;
+  MIDWT(x, m, n, h, lh, L, y);
+}
+
+#define mat(a, i, j) (*(a + (m*(j)+i)))  /* macro for matrix indices */
+
+
+#ifdef __STDC__
+MIDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L, double *y)
+#else
+MIDWT(x, m, n, h, lh, L, y)
+double *x, *h, *y;
+intptr_t m, n, lh, L;
+#endif
+{
+  double  *g0, *g1, *ydummyl, *ydummyh, *xdummy;
+  intptr_t i, j;
+  intptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1, lhhm1, sample_f;
+  xdummy = (double *)mxCalloc(max(m,n),sizeof(double));
+  ydummyl = (double *)mxCalloc(max(m,n)+lh/2-1,sizeof(double));
+  ydummyh = (double *)(intptr_t)mxCalloc(max(m,n)+lh/2-1,sizeof(double));
+  g0 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
+  g1 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
+
+  if (n==1){
+    n = m;
+    m = 1;
+  }
+  /* synthesis lowpass and highpass */
+  for (i=0; i<lh; i++){
+    g0[i] = h[i];
+    g1[i] = h[lh-i-1];
+  }
+  for (i=1; i<=lh; i+=2)
+    g1[i] = -g1[i];
+
+  lhm1 = lh - 1;
+  lhhm1 = lh/2 - 1;
+  /* 2^L */
+  sample_f = 1;
+  for (i=1; i<L; i++)
+    sample_f = sample_f*2;
+
+  if (m>1)
+    actual_m = m/sample_f;
+  else
+    actual_m = 1;
+  actual_n = n/sample_f;
+
+  for (i=0; i<(m*n); i++)
+    x[i] = y[i];
+
+  /* main loop */
+  for (actual_L=L; actual_L >= 1; actual_L--){
+    r_o_a = actual_m/2;
+    c_o_a = actual_n/2;
+
+    /* 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 */
+	ir = r_o_a;
+	for (i=0; i<r_o_a; i++){
+	  ydummyl[i+lhhm1] = mat(x, i, ic);
+	  ydummyh[i+lhhm1] = mat(x, ir++, ic);
+	}
+	/* perform filtering lowpass and highpass*/
+	bpsconv(xdummy, r_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh);
+	/* restore dummy variables in matrix */
+	for (i=0; i<actual_m; i++)
+	  mat(x, i, ic) = xdummy[i];
+      }
+    }
+    /* go by rows */
+    for (ir=0; ir<actual_m; ir++){            /* loop over rows */
+      /* store in dummy variable */
+      ic = c_o_a;
+      for  (i=0; i<c_o_a; i++){
+	ydummyl[i+lhhm1] = mat(x, ir, i);
+	ydummyh[i+lhhm1] = mat(x, ir, ic++);
+      }
+      /* perform filtering lowpass and highpass*/
+      bpsconv(xdummy, c_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh);
+      /* restore dummy variables in matrices */
+      for (i=0; i<actual_n; i++)
+        mat(x, ir, i) = xdummy[i];
+    }
+    if (m==1)
+      actual_m = 1;
+    else
+      actual_m = actual_m*2;
+    actual_n = actual_n*2;
+  }
+}
+
+#ifdef __STDC__
+bpsconv(double *x_out, intptr_t lx, double *g0, double *g1, intptr_t lhm1,
+	intptr_t lhhm1, double *x_inl, double *x_inh)
+#else
+bpsconv(x_out, lx, g0, g1, lhm1, lhhm1, x_inl, x_inh)
+double *x_inl, *x_inh, *g0, *g1, *x_out;
+intptr_t lx, lhm1, lhhm1;
+#endif
+{
+  intptr_t i, j, ind, tj;
+  double x0, x1;
+
+  for (i=lhhm1-1; i > -1; i--){
+    x_inl[i] = x_inl[lx+i];
+    x_inh[i] = x_inh[lx+i];
+  }
+  ind = 0;
+  for (i=0; i<(lx); i++){
+    x0 = 0;
+    x1 = 0;
+    tj = -2;
+    for (j=0; j<=lhhm1; j++){
+      tj+=2;
+      x0 = x0 + x_inl[i+j]*g0[lhm1-1-tj] + x_inh[i+j]*g1[lhm1-1-tj] ;
+      x1 = x1 + x_inl[i+j]*g0[lhm1-tj] + x_inh[i+j]*g1[lhm1-tj] ;
+    }
+    x_out[ind++] = x0;
+    x_out[ind++] = x1;
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/midwt.mhelp	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,93 @@
+function [y,L] = midwt(x,h,L);
+%    [x,L] = midwt(y,h,L);
+%
+%    Function computes the inverse discrete wavelet transform x for a 1D or
+%    2D input signal y using the scaling filter h.
+%
+%    Input:
+%	y : finite length 1D or 2D input signal (implicitly periodized)
+%           (see function mdwt to find the structure of y)
+%       h : scaling filter
+%       L : number of levels. In the case of a 1D signal, length(x) must be
+%           divisible by 2^L; in the case of a 2D signal, the row and the
+%           column dimension must be divisible by 2^L.  If no argument is
+%           specified, a full inverse DWT is returned for maximal possible
+%           L.
+%
+%    Output:
+%       x : periodic reconstructed signal
+%       L : number of decomposition levels
+%
+%    1D Example:
+%       xin = makesig('LinChirp',8);
+%       h = daubcqf(4,'min');
+%       L = 1;
+%       [y,L] = mdwt(xin,h,L);
+%       [x,L] = midwt(y,h,L)
+%
+%    1D Example's  output:
+%
+%       x = 0.0491 0.1951 0.4276 0.7071 0.9415 0.9808 0.6716 0.0000
+%       L = 1
+%
+%    See also: mdwt, mrdwt, mirdwt
+%
+
+%
+%
+%File Name: midwt.m
+%Last Modification Date: 08/07/95	15:13:52
+%Current Version: midwt.m	2.4
+%File Creation Date: Wed Oct 19 10:51:58 1994
+%Author: Markus Lang  <lang@jazz.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Markus Lang, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+%
+%Change History:
+%
+%Modification #1
+%Mon Aug  7 11:52:33 CDT 1995
+%Rebecca Hindman <hindman@ece.rice.edu>
+%Added L to function line so that it can be displayed as an output
+%
+%Thu Mar  2 13:07:11 CDT 2000
+%Ramesh Neelamani<neelsh@ece.rice.edu>
+%Revamped the help file
+%
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/midwt_r.c	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,193 @@
+/*
+File Name: MIDWT.c
+Last Modification Date:	06/14/95	13:01:15
+Current Version: MIDWT.c	2.4
+File Creation Date: Wed Oct 12 08:44:43 1994
+Author: Markus Lang  <lang@jazz.rice.edu>
+
+Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+Created by Markus Lang, Department of ECE, Rice University.
+
+This software is distributed and licensed to you on a non-exclusive
+basis, free-of-charge. Redistribution and use in source and binary forms,
+with or without modification, are permitted provided that the following
+conditions are met:
+
+1. Redistribution of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistribution in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. All advertising materials mentioning features or use of this software
+   must display the following acknowledgment: This product includes
+   software developed by Rice University, Houston, Texas and its contributors.
+4. Neither the name of the University nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+For information on commercial licenses, contact Rice University's Office of
+Technology Transfer at techtran@rice.edu or (713) 348-6173
+
+Change History: Fixed the code such that 1D vectors passed to it can be in
+                either passed as a row or column vector. Also took care of
+		the code such that it will compile with both under standard
+		C compilers as well as for ANSI C compilers
+		Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
+
+		Fix minor bug to allow maximum number of levels
+
+decription of the matlab call:
+%y = midwt(x,h,L);
+%
+% function computes the inverse discrete wavelet transform y for a 1D or 2D
+% input signal x.
+%
+%    Input:
+%	x    : finite length 1D or 2D input signal (implicitely periodized)
+%       h    : scaling filter
+%       L    : number of levels. in case of a 1D signal length(x) must be
+%              divisible by 2^L; in case of a 2D signal the row and the
+%              column dimension must be divisible by 2^L.
+%
+% see also: mdwt, mrdwt, mirdwt
+
+*/
+
+#include <math.h>
+#include <stdio.h>
+#include <inttypes.h>
+
+#define max(A,B) (A > B ? A : B)
+#define mat(a, i, j) (*(a + (m*(j)+i)))  /* macro for matrix indices */
+
+
+#ifdef __STDC__
+MIDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L, double *y)
+#else
+MIDWT(x, m, n, h, lh, L, y)
+double *x, *h, *y;
+uintptr_t m, n, lh, L;
+#endif
+{
+  double  *g0, *g1, *ydummyl, *ydummyh, *xdummy;
+  uintptr_t i, j;
+  uintptr_t actual_L, actual_m, actual_n, r_o_a, c_o_a, ir, ic, lhm1, lhhm1, sample_f;
+  xdummy = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyl = (double *)(uintptr_t)mxCalloc(max(m,n)+lh/2-1,sizeof(double));
+  ydummyh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh/2-1,sizeof(double));
+  g0 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
+  g1 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
+
+  if (n==1){
+    n = m;
+    m = 1;
+  }
+  /* synthesis lowpass and highpass */
+  for (i=0; i<lh; i++){
+    g0[i] = h[i];
+    g1[i] = h[lh-i-1];
+  }
+  for (i=1; i<=lh; i+=2)
+    g1[i] = -g1[i];
+
+  lhm1 = lh - 1;
+  lhhm1 = lh/2 - 1;
+  /* 2^L */
+  sample_f = 1;
+  for (i=1; i<L; i++)
+    sample_f = sample_f*2;
+
+  if (m>1)
+    actual_m = m/sample_f;
+  else
+    actual_m = 1;
+  actual_n = n/sample_f;
+
+  for (i=0; i<(m*n); i++)
+    x[i] = y[i];
+
+  /* main loop */
+  for (actual_L=L; actual_L >= 1; actual_L--){
+    r_o_a = actual_m/2;
+    c_o_a = actual_n/2;
+
+    /* 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 */
+	ir = r_o_a;
+	for (i=0; i<r_o_a; i++){
+	  ydummyl[i+lhhm1] = mat(x, i, ic);
+	  ydummyh[i+lhhm1] = mat(x, ir++, ic);
+	}
+	/* perform filtering lowpass and highpass*/
+	bpsconv(xdummy, r_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh);
+	/* restore dummy variables in matrix */
+	for (i=0; i<actual_m; i++)
+	  mat(x, i, ic) = xdummy[i];
+      }
+    }
+    /* go by rows */
+    for (ir=0; ir<actual_m; ir++){            /* loop over rows */
+      /* store in dummy variable */
+      ic = c_o_a;
+      for  (i=0; i<c_o_a; i++){
+	ydummyl[i+lhhm1] = mat(x, ir, i);
+	ydummyh[i+lhhm1] = mat(x, ir, ic++);
+      }
+      /* perform filtering lowpass and highpass*/
+      bpsconv(xdummy, c_o_a, g0, g1, lhm1, lhhm1, ydummyl, ydummyh);
+      /* restore dummy variables in matrices */
+      for (i=0; i<actual_n; i++)
+        mat(x, ir, i) = xdummy[i];
+    }
+    if (m==1)
+      actual_m = 1;
+    else
+      actual_m = actual_m*2;
+    actual_n = actual_n*2;
+  }
+}
+
+#ifdef __STDC__
+bpsconv(double *x_out, uintptr_t lx, double *g0, double *g1, uintptr_t lhm1,
+	uintptr_t lhhm1, double *x_inl, double *x_inh)
+#else
+bpsconv(x_out, lx, g0, g1, lhm1, lhhm1, x_inl, x_inh)
+double *x_inl, *x_inh, *g0, *g1, *x_out;
+uintptr_t lx, lhm1, lhhm1;
+#endif
+{
+  uintptr_t i, j, ind, tj;
+  double x0, x1;
+
+  for (i=lhhm1-1; i > -1; i--){
+    x_inl[i] = x_inl[lx+i];
+    x_inh[i] = x_inh[lx+i];
+  }
+  ind = 0;
+  for (i=0; i<(lx); i++){
+    x0 = 0;
+    x1 = 0;
+    tj = -2;
+    for (j=0; j<=lhhm1; j++){
+      tj+=2;
+      x0 = x0 + x_inl[i+j]*g0[lhm1-1-tj] + x_inh[i+j]*g1[lhm1-1-tj] ;
+      x1 = x1 + x_inl[i+j]*g0[lhm1-tj] + x_inh[i+j]*g1[lhm1-tj] ;
+    }
+    x_out[ind++] = x0;
+    x_out[ind++] = x1;
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mirdwt.c	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,261 @@
+/*
+File Name: mirdwt.c
+Last Modification Date:	%G%	%U%
+Current Version: %M%	%I%
+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
+
+*/
+
+#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)
+
+
+void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
+
+{
+  double *x, *h,  *yl, *yh, *Lf, *Lr;
+  intptr_t m, n, mh, nh, h_col, h_row, lh, L, i, po2, j;
+  double mtest, ntest;
+
+  /* check for correct # of input variables */
+  if (nrhs>4){
+    mexErrMsgTxt("There are at most 4 input parameters allowed!");
+    return;
+  }
+  if (nrhs<3){
+    mexErrMsgTxt("There are at least 3 input parameters required!");
+    return;
+  }
+  yl = mxGetPr(prhs[0]);
+  n = mxGetN(prhs[0]);
+  m = mxGetM(prhs[0]);
+  yh = mxGetPr(prhs[1]);
+  nh = mxGetN(prhs[1]);
+  mh = mxGetM(prhs[1]);
+  h = mxGetPr(prhs[2]);
+  h_col = mxGetN(prhs[2]);
+  h_row = mxGetM(prhs[2]);
+  if (h_col>h_row)
+    lh = h_col;
+  else
+    lh = h_row;
+  if (nrhs == 4){
+    L = (intptr_t) *mxGetPr(prhs[3]);
+    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 for consistency of rows and columns of yl, yh */
+  if (min(m,n) > 1){
+    if((m != mh) | (3*n*L != nh)){
+      mexErrMsgTxt("Dimensions of first two input matrices not consistent!");
+      return;
+    }
+  }
+  else{
+    if((m != mh) | (n*L != nh)){
+      mexErrMsgTxt("Dimensions of first two input vectors not consistent!");{
+	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);
+  x = mxGetPr(plhs[0]);
+  plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);
+  Lr = mxGetPr(plhs[1]);
+  *Lr = L;
+  MIRDWT(x, m, n, h, lh, L, yl, yh);
+}
+#define mat(a, i, j) (*(a + (m*(j)+i)))  /* macro for matrix indices */
+
+#ifdef __STDC__
+MIRDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L,
+       double *yl, double *yh)
+#else
+MIRDWT(x, m, n, h, lh, L, yl, yh)
+double *x, *h, *yl, *yh;
+intptr_t m, n, lh, L;
+#endif
+{
+  double  *g0, *g1, *ydummyll, *ydummylh, *ydummyhl;
+  double *ydummyhh, *xdummyl , *xdummyh, *xh;
+  intptr_t i, j;
+  intptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o, lhm1;
+  intptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
+  xh = (double *)(intptr_t)mxCalloc(m*n,sizeof(double));
+  xdummyl = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double));
+  xdummyh = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyll = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummylh = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyhl = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyhh = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  g0 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
+  g1 = (double *)(intptr_t)mxCalloc(lh,sizeof(double));
+
+  if (n==1){
+    n = m;
+    m = 1;
+  }
+  /* analysis lowpass and highpass */
+  for (i=0; i<lh; i++){
+    g0[i] = h[i]/2;
+    g1[i] = h[lh-i-1]/2;
+  }
+  for (i=1; i<=lh; i+=2)
+    g1[i] = -g1[i];
+
+  lhm1 = lh - 1;
+  /* 2^L */
+  sample_f = 1;
+  for (i=1; i<L; i++)
+    sample_f = sample_f*2;
+  actual_m = m/sample_f;
+  actual_n = n/sample_f;
+  /* restore yl in x */
+  for (i=0;i<m*n;i++)
+    x[i] = yl[i];
+
+  /* main loop */
+  for (actual_L=L; actual_L >= 1; actual_L--){
+    /* actual (level dependent) column offset */
+    if (m==1)
+      c_o_a = n*(actual_L-1);
+    else
+      c_o_a = 3*n*(actual_L-1);
+    c_o_a_p2n = c_o_a + 2*n;
+
+    /* go by columns in case of a 2D signal*/
+    if (m>1){
+      n_rb = m/actual_m;                 /* # of row blocks per column */
+      for (ic=0; ic<n; ic++){            /* loop over column */
+	for (n_r=0; n_r<n_rb; n_r++){    /* loop within one column */
+	  /* store in dummy variables */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){
+	    ir = ir + sample_f;
+	    ydummyll[i+lhm1] = mat(x, ir, ic);
+	    ydummylh[i+lhm1] = mat(yh, ir, c_o_a+ic);
+	    ydummyhl[i+lhm1] = mat(yh, ir,c_o_a+n+ic);
+	    ydummyhh[i+lhm1] = mat(yh, ir, c_o_a_p2n+ic);
+	  }
+	  /* perform filtering and adding: first LL/LH, then HL/HH */
+	  bpconv(xdummyl, actual_m, g0, g1, lh, ydummyll, ydummylh);
+	  bpconv(xdummyh, actual_m, g0, g1, lh, ydummyhl, ydummyhh);
+	  /* store dummy variables in matrices */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){
+	    ir = ir + sample_f;
+	    mat(x, ir, ic) = xdummyl[i];
+	    mat(xh, ir, ic) = xdummyh[i];
+	  }
+	}
+      }
+    }
+
+    /* go by rows */
+    n_cb = n/actual_n;                 /* # of column blocks per row */
+    for (ir=0; ir<m; ir++){            /* loop over rows */
+      for (n_c=0; n_c<n_cb; n_c++){    /* loop within one row */
+	/* store in dummy variable */
+	ic = -sample_f + n_c;
+	for  (i=0; i<actual_n; i++){
+	  ic = ic + sample_f;
+	  ydummyll[i+lhm1] = mat(x, ir, ic);
+	  if (m>1)
+	    ydummyhh[i+lhm1] = mat(xh, ir, ic);
+	  else
+	    ydummyhh[i+lhm1] = mat(yh, ir, c_o_a+ic);
+	}
+	/* perform filtering lowpass/highpass */
+	bpconv(xdummyl, actual_n, g0, g1, lh, ydummyll, ydummyhh);
+	/* restore dummy variables in matrices */
+	ic = -sample_f + n_c;
+	for (i=0; i<actual_n; i++){
+	  ic = ic + sample_f;
+	  mat(x, ir, ic) = xdummyl[i];
+	}
+      }
+    }
+    sample_f = sample_f/2;
+    actual_m = actual_m*2;
+    actual_n = actual_n*2;
+  }
+}
+
+#ifdef __STDC__
+bpconv(double *x_out, intptr_t lx, double *g0, double *g1, intptr_t lh,
+       double *x_inl, double *x_inh)
+#else
+bpconv(x_out, lx, g0, g1, lh, x_inl, x_inh)
+double *x_inl, *x_inh, *g0, *g1, *x_out;
+intptr_t lx, lh;
+#endif
+{
+  intptr_t i, j;
+  double x0;
+
+  for (i=lh-2; i > -1; i--){
+    x_inl[i] = x_inl[lx+i];
+    x_inh[i] = x_inh[lx+i];
+  }
+  for (i=0; i<lx; i++){
+    x0 = 0;
+    for (j=0; j<lh; j++)
+      x0 = x0 + x_inl[j+i]*g0[lh-1-j] +
+	x_inh[j+i]*g1[lh-1-j];
+    x_out[i] = x0;
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mirdwt.mhelp	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,138 @@
+function [x,L] = mirdwt(yl,yh,h,L);
+%    function [x,L] = mirdwt(yl,yh,h,L);
+%
+%    Function computes the inverse redundant discrete wavelet
+%    transform x  for a 1D or 2D input signal. (Redundant means here
+%    that the sub-sampling after each stage of the forward transform
+%    has been omitted.) yl contains the lowpass and yl the highpass
+%    components as computed, e.g., by mrdwt. In the case of a 2D
+%    signal, the ordering in
+%    yh is [lh hl hh lh hl ... ] (first letter refers to row, second
+%    to column filtering).
+%
+%    Input:
+%       yl : lowpass component
+%       yh : highpass components
+%       h  : scaling filter
+%       L  : number of levels. In the case of a 1D signal,
+%            length(yl) must  be divisible by 2^L;
+%            in the case of a 2D signal, the row and
+%            the column dimension must be divisible by 2^L.
+%
+%    Output:
+%	     x : finite length 1D or 2D signal
+%	     L : number of levels
+%
+%  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
+%  Cut-and-paste the example you want to run to a new file
+%  called ex.m, for example. Delete out the % at the beginning
+%  of each line in ex.m (Can use search-and-replace in your editor
+%  to replace it with a space). Type 'ex' in matlab and hit return.
+%
+%
+%    Example 1:
+%    xin = makesig('Leopold',8);
+%    h = daubcqf(4,'min');
+%    L = 1;
+%    [yl,yh,L] = mrdwt(xin,h,L);
+%    [x,L] = mirdwt(yl,yh,h,L)
+%    x = 0.0000 1.0000 0.0000 -0.0000 0 0 0 -0.0000
+%    L = 1
+%
+%    Example 2:
+%    load lena;
+%    h = daubcqf(4,'min');
+%    L = 2;
+%    [ll_lev2,yh,L] = mrdwt(lena,h,L); % lena is a 256x256 matrix
+%    N = 256;
+%    lh_lev1 = yh(:,1:N);
+%    hl_lev1 = yh(:,N+1:2*N);
+%    hh_lev1 = yh(:,2*N+1:3*N);
+%    lh_lev2 = yh(:,3*N+1:4*N);
+%    hl_lev2 = yh(:,4*N+1:5*N);
+%    hh_lev2 = yh(:,5*N+1:6*N);
+%    figure; colormap(gray); imagesc(lena); title('Original Image');
+%    figure; colormap(gray); imagesc(ll_lev2); title('LL Level 2');
+%    figure; colormap(gray); imagesc(hh_lev2); title('HH Level 2');
+%    figure; colormap(gray); imagesc(hl_lev2); title('HL Level 2');
+%    figure; colormap(gray); imagesc(lh_lev2); title('LH Level 2');
+%    figure; colormap(gray); imagesc(hh_lev1); title('HH Level 1');
+%    figure; colormap(gray); imagesc(hl_lev2); title('HL Level 1');
+%    figure; colormap(gray); imagesc(lh_lev2); title('LH Level 1');
+%    [lena_Hat,L] = mirdwt(ll_lev2,yh,h,L);
+%    figure; colormap(gray); imagesc(lena_Hat);
+%                            title('Reconstructed Image');
+%
+%    See also: mdwt, midwt, mrdwt
+%
+%    Warning! min(size(yl))/2^L should be greater than length(h)
+%
+
+%File Name: mirdwt.m
+%Last Modification Date: 08/07/95	15:14:21
+%Current Version: mirdwt.m	2.4
+%File Creation Date: Wed Oct 19 10:51:58 1994
+%Author: Markus Lang  <lang@jazz.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Markus Lang, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+%
+%Change History:
+%
+%Modification #1
+%Mon Aug  7 15:09:51 CDT 1995
+%Rebecca Hindman <hindman@ece.rice.edu>
+%Added L to function line so that it can be displayed as an output
+%
+%Modification #2
+%Thursday Mar 2 2000
+% Added Example 2
+% Felix Fernandes <felixf@rice.edu>
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mirdwt_r.c	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,219 @@
+/*
+File Name: MIRDWT.c
+Last Modification Date:	06/14/95	16:22:45
+Current Version: MIRDWT.c	2.4
+File Creation Date: Wed Oct 12 08:44:43 1994
+Author: Markus Lang  <lang@jazz.rice.edu>
+
+Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+Created by Markus Lang, Department of ECE, Rice University.
+
+This software is distributed and licensed to you on a non-exclusive
+basis, free-of-charge. Redistribution and use in source and binary forms,
+with or without modification, are permitted provided that the following
+conditions are met:
+
+1. Redistribution of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistribution in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. All advertising materials mentioning features or use of this software
+   must display the following acknowledgment: This product includes
+   software developed by Rice University, Houston, Texas and its contributors.
+4. Neither the name of the University nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+For information on commercial licenses, contact Rice University's Office of
+Technology Transfer at techtran@rice.edu or (713) 348-6173
+
+Change History: Fixed the code such that 1D vectors passed to it can be in
+                either passed as a row or column vector. Also took care of
+		the code such that it will compile with both under standard
+		C compilers as well as for ANSI C compilers
+		Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
+
+                Fix minor bug to allow maximum number of levels
+
+MATLAB description:
+%function x = mirdwt(yl,yh,h,L);
+%
+% function computes the inverse redundant discrete wavelet transform y for a
+% 1D or  2D input signal. redundant means here that the subsampling after
+% each stage of the forward transform has been omitted. yl contains the
+% lowpass and yl the highpass components as computed, e.g., by mrdwt. In
+% case of a 2D signal the ordering in yh is [lh hl hh lh hl ... ] (first
+% letter refers to row, second to column filtering).
+%
+%    Input:
+%       yl   : lowpass component
+%       yh   : highpass components
+%       h    : scaling filter
+%       L    : number of levels. in case of a 1D signal length(yl) must be
+%              divisible by 2^L; in case of a 2D signal the row and the
+%              column dimension must be divisible by 2^L.
+%
+%    Output:
+%	x    : finite length 1D or 2D signal
+%
+% see also: mdwt, midwt, mrdwt
+
+*/
+#include <math.h>
+#include <stdio.h>
+#include <inttypes.h>
+
+#define max(a, b) ((a) > (b) ? (a) : (b))
+#define mat(a, i, j) (*(a + (m*(j)+i)))  /* macro for matrix indices */
+
+#ifdef __STDC__
+MIRDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L,
+       double *yl, double *yh)
+#else
+MIRDWT(x, m, n, h, lh, L, yl, yh)
+double *x, *h, *yl, *yh;
+uintptr_t m, n, lh, L;
+#endif
+{
+  double  *g0, *g1, *ydummyll, *ydummylh, *ydummyhl;
+  double *ydummyhh, *xdummyl , *xdummyh, *xh;
+  long i, j;
+  uintptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o, lhm1;
+  uintptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
+  xh = (double *)(uintptr_t)mxCalloc(m*n,sizeof(double));
+  xdummyl = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  xdummyh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyll = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummylh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyhl = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyhh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  g0 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
+  g1 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
+
+  if (n==1){
+    n = m;
+    m = 1;
+  }
+  /* analysis lowpass and highpass */
+  for (i=0; i<lh; i++){
+    g0[i] = h[i]/2;
+    g1[i] = h[lh-i-1]/2;
+  }
+  for (i=1; i<=lh; i+=2)
+    g1[i] = -g1[i];
+
+  lhm1 = lh - 1;
+  /* 2^L */
+  sample_f = 1;
+  for (i=1; i<L; i++)
+    sample_f = sample_f*2;
+  actual_m = m/sample_f;
+  actual_n = n/sample_f;
+  /* restore yl in x */
+  for (i=0;i<m*n;i++)
+    x[i] = yl[i];
+
+  /* main loop */
+  for (actual_L=L; actual_L >= 1; actual_L--){
+    /* actual (level dependent) column offset */
+    if (m==1)
+      c_o_a = n*(actual_L-1);
+    else
+      c_o_a = 3*n*(actual_L-1);
+    c_o_a_p2n = c_o_a + 2*n;
+
+    /* go by columns in case of a 2D signal*/
+    if (m>1){
+      n_rb = m/actual_m;                 /* # of row blocks per column */
+      for (ic=0; ic<n; ic++){            /* loop over column */
+	for (n_r=0; n_r<n_rb; n_r++){    /* loop within one column */
+	  /* store in dummy variables */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){
+	    ir = ir + sample_f;
+	    ydummyll[i+lhm1] = mat(x, ir, ic);
+	    ydummylh[i+lhm1] = mat(yh, ir, c_o_a+ic);
+	    ydummyhl[i+lhm1] = mat(yh, ir,c_o_a+n+ic);
+	    ydummyhh[i+lhm1] = mat(yh, ir, c_o_a_p2n+ic);
+	  }
+	  /* perform filtering and adding: first LL/LH, then HL/HH */
+	  bpconv(xdummyl, actual_m, g0, g1, lh, ydummyll, ydummylh);
+	  bpconv(xdummyh, actual_m, g0, g1, lh, ydummyhl, ydummyhh);
+	  /* store dummy variables in matrices */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){
+	    ir = ir + sample_f;
+	    mat(x, ir, ic) = xdummyl[i];
+	    mat(xh, ir, ic) = xdummyh[i];
+	  }
+	}
+      }
+    }
+
+    /* go by rows */
+    n_cb = n/actual_n;                 /* # of column blocks per row */
+    for (ir=0; ir<m; ir++){            /* loop over rows */
+      for (n_c=0; n_c<n_cb; n_c++){    /* loop within one row */
+	/* store in dummy variable */
+	ic = -sample_f + n_c;
+	for  (i=0; i<actual_n; i++){
+	  ic = ic + sample_f;
+	  ydummyll[i+lhm1] = mat(x, ir, ic);
+	  if (m>1)
+	    ydummyhh[i+lhm1] = mat(xh, ir, ic);
+	  else
+	    ydummyhh[i+lhm1] = mat(yh, ir, c_o_a+ic);
+	}
+	/* perform filtering lowpass/highpass */
+	bpconv(xdummyl, actual_n, g0, g1, lh, ydummyll, ydummyhh);
+	/* restore dummy variables in matrices */
+	ic = -sample_f + n_c;
+	for (i=0; i<actual_n; i++){
+	  ic = ic + sample_f;
+	  mat(x, ir, ic) = xdummyl[i];
+	}
+      }
+    }
+    sample_f = sample_f/2;
+    actual_m = actual_m*2;
+    actual_n = actual_n*2;
+  }
+}
+
+#ifdef __STDC__
+bpconv(double *x_out, uintptr_t lx, double *g0, double *g1, uintptr_t lh,
+       double *x_inl, double *x_inh)
+#else
+bpconv(x_out, lx, g0, g1, lh, x_inl, x_inh)
+double *x_inl, *x_inh, *g0, *g1, *x_out;
+uintptr_t lx, lh;
+#endif
+{
+  uintptr_t i, j;
+  double x0;
+
+  for (i=lh-2; i > -1; i--){
+    x_inl[i] = x_inl[lx+i];
+    x_inh[i] = x_inh[lx+i];
+  }
+  for (i=0; i<lx; i++){
+    x0 = 0;
+    for (j=0; j<lh; j++)
+      x0 = x0 + x_inl[j+i]*g0[lh-1-j] +
+	x_inh[j+i]*g1[lh-1-j];
+    x_out[i] = x0;
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mrdwt.c	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,242 @@
+/*
+File Name: mrdwt.c
+Last Modification Date:	%G%	%U%
+Current Version: %M%	%I%
+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
+*/
+
+#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)
+
+
+void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
+
+{
+  double *x, *h,  *yl, *yh, *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);
+  yl = mxGetPr(plhs[0]);
+  if (min(m,n) == 1)
+    plhs[1] = mxCreateDoubleMatrix(m,L*n,mxREAL);
+  else
+    plhs[1] = mxCreateDoubleMatrix(m,3*L*n,mxREAL);
+  yh = mxGetPr(plhs[1]);
+  plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL);
+  Lr = mxGetPr(plhs[2]);
+  *Lr = L;
+  MRDWT(x, m, n, h, lh, L, yl, yh);
+}
+#define mat(a, i, j) (*(a + (m*(j)+i)))
+
+
+#ifdef __STDC__
+MRDWT(double *x, intptr_t m, intptr_t n, double *h, intptr_t lh, intptr_t L,
+      double *yl, double *yh)
+#else
+MRDWT(x, m, n, h, lh, L, yl, yh)
+double *x, *h, *yl, *yh;
+intptr_t m, n, lh, L;
+#endif
+{
+  double *tmp;
+  double  *h0, *h1, *ydummyll, *ydummylh, *ydummyhl;
+  double *ydummyhh, *xdummyl , *xdummyh;
+  long i, j;
+  intptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o;
+  intptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
+  xdummyl = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  xdummyh = (double *)(intptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyll = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummylh = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyhl = (double *)(intptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyhh = (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));
+
+  if (n==1){
+    n = m;
+    m = 1;
+  }
+  /* analysis lowpass and highpass */
+  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];
+
+  actual_m = 2*m;
+  actual_n = 2*n;
+  for (i=0; i<m*n; i++)
+    yl[i] = x[i];
+
+  /* main loop */
+  sample_f = 1;
+  for (actual_L=1; actual_L <= L; actual_L++){
+    actual_m = actual_m/2;
+    actual_n = actual_n/2;
+    /* actual (level dependent) column offset */
+    if (m==1)
+      c_o_a = n*(actual_L-1);
+    else
+      c_o_a = 3*n*(actual_L-1);
+    c_o_a_p2n = c_o_a + 2*n;
+
+    /* go by rows */
+    n_cb = n/actual_n;                 /* # of column blocks per row */
+    for (ir=0; ir<m; ir++){            /* loop over rows */
+      for (n_c=0; n_c<n_cb; n_c++){    /* loop within one row */
+	/* store in dummy variable */
+	ic = -sample_f + n_c;
+	for (i=0; i<actual_n; i++){
+	  ic = ic + sample_f;
+	  xdummyl[i] = mat(yl, ir, ic);
+	}
+	/* perform filtering lowpass/highpass */
+	fpconv(xdummyl, actual_n, h0, h1, lh, ydummyll, ydummyhh);
+	/* restore dummy variables in matrices */
+	ic = -sample_f + n_c;
+	for  (i=0; i<actual_n; i++){
+	  ic = ic + sample_f;
+	  mat(yl, ir, ic) = ydummyll[i];
+	  mat(yh, ir, c_o_a+ic) = ydummyhh[i];
+	}
+      }
+    }
+
+    /* go by columns in case of a 2D signal*/
+    if (m>1){
+      n_rb = m/actual_m;                 /* # of row blocks per column */
+      for (ic=0; ic<n; ic++){            /* loop over column */
+	for (n_r=0; n_r<n_rb; n_r++){    /* loop within one column */
+	  /* store in dummy variables */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){
+	    ir = ir + sample_f;
+	    xdummyl[i] = mat(yl, ir, ic);
+	    xdummyh[i] = mat(yh, ir,c_o_a+ic);
+	  }
+	  /* perform filtering: first LL/LH, then HL/HH */
+	  fpconv(xdummyl, actual_m, h0, h1, lh, ydummyll, ydummylh);
+	  fpconv(xdummyh, actual_m, h0, h1, lh, ydummyhl, ydummyhh);
+	  /* restore dummy variables in matrices */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){
+	    ir = ir + sample_f;
+	    mat(yl, ir, ic) = ydummyll[i];
+	    mat(yh, ir, c_o_a+ic) = ydummylh[i];
+	    mat(yh, ir,c_o_a+n+ic) = ydummyhl[i];
+	    mat(yh, ir, c_o_a_p2n+ic) = ydummyhh[i];
+	  }
+	}
+      }
+    }
+    sample_f = sample_f*2;
+  }
+}
+
+#ifdef __STDC__
+fpconv(double *x_in, intptr_t lx, double *h0, double *h1, intptr_t lh,
+       double *x_outl, double *x_outh)
+#else
+fpconv(x_in, lx, h0, h1, lh, x_outl, x_outh)
+double *x_in, *h0, *h1, *x_outl, *x_outh;
+intptr_t lx, lh;
+#endif
+{
+  intptr_t i, j;
+  double x0, x1;
+
+  for (i=lx; i < lx+lh-1; i++)
+    x_in[i] = x_in[i-lx];
+  for (i=0; i<lx; i++){
+    x0 = 0;
+    x1 = 0;
+    for (j=0; j<lh; j++){
+      x0 = x0 + x_in[j+i]*h0[lh-1-j];
+      x1 = x1 + x_in[j+i]*h1[lh-1-j];
+    }
+    x_outl[i] = x0;
+    x_outh[i] = x1;
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mrdwt.mhelp	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,120 @@
+function [yl,yh,L] = mrdwt(x,h,L);
+%    [yl,yh,L] = mrdwt(x,h,L);
+%
+%    Function computes the redundant discrete wavelet transform y
+%    for a 1D  or 2D input signal. (Redundant means here that the
+%    sub-sampling after each stage is omitted.) yl contains the
+%    lowpass and yh the highpass components. In the case of a 2D
+%    signal, the ordering in yh is
+%    [lh hl hh lh hl ... ] (first letter refers to row, second to
+%    column filtering).
+%
+%    Input:
+%	     x : finite length 1D or 2D signal (implicitly periodized)
+%       h : scaling filter
+%       L : number of levels. In the case of a 1D
+%           length(x) must be  divisible by 2^L;
+%           in the case of a 2D signal, the row and the
+%           column dimension must be divisible by 2^L.
+%           If no argument is
+%           specified, a full DWT is returned for maximal possible L.
+%
+%    Output:
+%       yl : lowpass component
+%       yh : highpass components
+%       L  : number of levels
+%
+%  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
+%  Cut-and-paste the example you want to run to a new file
+%  called ex.m, for example. Delete out the % at the beginning
+%  of each line in ex.m (Can use search-and-replace in your editor
+%  to replace it with a space). Type 'ex' in matlab and hit return.
+%
+%
+%    Example 1::
+%    x = makesig('Leopold',8);
+%    h = daubcqf(4,'min');
+%    L = 1;
+%    [yl,yh,L] = mrdwt(x,h,L)
+%    yl =  0.8365  0.4830 0 0 0 0 -0.1294 0.2241
+%    yh = -0.2241 -0.1294 0 0 0 0 -0.4830 0.8365
+%    L = 1
+%    Example 2:
+%    load lena;
+%    h = daubcqf(4,'min');
+%    L = 2;
+%    [ll_lev2,yh,L] = mrdwt(lena,h,L); % lena is a 256x256 matrix
+%    N = 256;
+%    lh_lev1 = yh(:,1:N);
+%    hl_lev1 = yh(:,N+1:2*N);
+%    hh_lev1 = yh(:,2*N+1:3*N);
+%    lh_lev2 = yh(:,3*N+1:4*N);
+%    hl_lev2 = yh(:,4*N+1:5*N);
+%    hh_lev2 = yh(:,5*N+1:6*N);
+%    figure; colormap(gray); imagesc(lena); title('Original Image');
+%    figure; colormap(gray); imagesc(ll_lev2); title('LL Level 2');
+%    figure; colormap(gray); imagesc(hh_lev2); title('HH Level 2');
+%    figure; colormap(gray); imagesc(hl_lev2); title('HL Level 2');
+%    figure; colormap(gray); imagesc(lh_lev2); title('LH Level 2');
+%    figure; colormap(gray); imagesc(hh_lev1); title('HH Level 1');
+%    figure; colormap(gray); imagesc(hl_lev2); title('HL Level 1');
+%    figure; colormap(gray); imagesc(lh_lev2); title('LH Level 1');
+%
+%    See also: mdwt, midwt, mirdwt
+%
+%    Warning! min(size(x))/2^L should be greater than length(h)
+%
+
+%File Name: mrdwt.m
+%Last Modification Date: 08/07/95	15:14:47
+%Current Version: mrdwt.m	2.4
+%File Creation Date: Wed Oct 19 10:51:58 1994
+%Author: Markus Lang  <lang@jazz.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Markus Lang, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+%
+%Change History:
+%
+%Modification #1
+%Mon Aug  7 14:26:26 CDT 1995
+%Rebecca Hindman <hindman@ece.rice.edu>
+%Added L to function line so that it can be displayed as an output
+%
+%Thursday Mar 2 2000
+% Added Example 2
+% Felix Fernandes <felixf@rice.edu>
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/mrdwt_r.c	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,213 @@
+/*
+File Name: MRDWT.c
+Last Modification Date:	09/21/95	15:42:59
+Current Version: MRDWT.c	2.4
+File Creation Date: Wed Oct 12 08:44:43 1994
+Author: Markus Lang  <lang@jazz.rice.edu>
+
+Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+Created by Markus Lang, Department of ECE, Rice University.
+
+This software is distributed and licensed to you on a non-exclusive
+basis, free-of-charge. Redistribution and use in source and binary forms,
+with or without modification, are permitted provided that the following
+conditions are met:
+
+1. Redistribution of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistribution in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. All advertising materials mentioning features or use of this software
+   must display the following acknowledgment: This product includes
+   software developed by Rice University, Houston, Texas and its contributors.
+4. Neither the name of the University nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+For information on commercial licenses, contact Rice University's Office of
+Technology Transfer at techtran@rice.edu or (713) 348-6173
+
+Change History: Fixed the code such that 1D vectors passed to it can be in
+                either passed as a row or column vector. Also took care of
+		the code such that it will compile with both under standard
+		C compilers as well as for ANSI C compilers
+		Jan Erik Odegard <odegard@ece.rice.edu> Wed Jun 14 1995
+
+MATLAB description:
+%[yl,yh] = mrdwt(x,h,L);
+%
+% function computes the redundant discrete wavelet transform y for a 1D or
+% 2D input signal . redundant means here that the subsampling after each
+% stage is omitted. yl contains the lowpass and yl the highpass
+% components. In case of a 2D signal the ordering in yh is [lh hl hh lh hl
+% ... ] (first letter refers to row, second to column filtering).
+%
+%    Input:
+%	x    : finite length 1D or 2D signal (implicitely periodized)
+%       h    : scaling filter
+%       L    : number of levels. in case of a 1D signal length(x) must be
+%              divisible by 2^L; in case of a 2D signal the row and the
+%              column dimension must be divisible by 2^L.
+%
+%    Output:
+%       yl   : lowpass component
+%       yh   : highpass components
+%
+% see also: mdwt, midwt, mirdwt
+
+
+*/
+
+#include <math.h>
+#include <stdio.h>
+#include <inttypes.h>
+
+/*#define mat(a, i, j) (a[m*(j)+i]) */
+#define max(a, b) ((a) > (b) ? (a) : (b))
+#define mat(a, i, j) (*(a + (m*(j)+i)))
+
+
+#ifdef __STDC__
+MRDWT(double *x, uintptr_t m, uintptr_t n, double *h, uintptr_t lh, uintptr_t L,
+      double *yl, double *yh)
+#else
+MRDWT(x, m, n, h, lh, L, yl, yh)
+double *x, *h, *yl, *yh;
+uintptr_t m, n, lh, L;
+#endif
+{
+  double *tmp;
+  double  *h0, *h1, *ydummyll, *ydummylh, *ydummyhl;
+  double *ydummyhh, *xdummyl , *xdummyh;
+  long i, j;
+  uintptr_t actual_L, actual_m, actual_n, c_o_a, ir, n_c, n_cb, n_c_o;
+  uintptr_t ic, n_r, n_rb, n_r_o, c_o_a_p2n, sample_f;
+  xdummyl = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  xdummyh = (double *)(uintptr_t)mxCalloc(max(m,n)+lh-1,sizeof(double));
+  ydummyll = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummylh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyhl = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  ydummyhh = (double *)(uintptr_t)mxCalloc(max(m,n),sizeof(double));
+  h0 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
+  h1 = (double *)(uintptr_t)mxCalloc(lh,sizeof(double));
+
+  if (n==1){
+    n = m;
+    m = 1;
+  }
+  /* analysis lowpass and highpass */
+  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];
+
+  actual_m = 2*m;
+  actual_n = 2*n;
+  for (i=0; i<m*n; i++)
+    yl[i] = x[i];
+
+  /* main loop */
+  sample_f = 1;
+  for (actual_L=1; actual_L <= L; actual_L++){
+    actual_m = actual_m/2;
+    actual_n = actual_n/2;
+    /* actual (level dependent) column offset */
+    if (m==1)
+      c_o_a = n*(actual_L-1);
+    else
+      c_o_a = 3*n*(actual_L-1);
+    c_o_a_p2n = c_o_a + 2*n;
+
+    /* go by rows */
+    n_cb = n/actual_n;                 /* # of column blocks per row */
+    for (ir=0; ir<m; ir++){            /* loop over rows */
+      for (n_c=0; n_c<n_cb; n_c++){    /* loop within one row */
+	/* store in dummy variable */
+	ic = -sample_f + n_c;
+	for (i=0; i<actual_n; i++){
+	  ic = ic + sample_f;
+	  xdummyl[i] = mat(yl, ir, ic);
+	}
+	/* perform filtering lowpass/highpass */
+	fpconv(xdummyl, actual_n, h0, h1, lh, ydummyll, ydummyhh);
+	/* restore dummy variables in matrices */
+	ic = -sample_f + n_c;
+	for  (i=0; i<actual_n; i++){
+	  ic = ic + sample_f;
+	  mat(yl, ir, ic) = ydummyll[i];
+	  mat(yh, ir, c_o_a+ic) = ydummyhh[i];
+	}
+      }
+    }
+
+    /* go by columns in case of a 2D signal*/
+    if (m>1){
+      n_rb = m/actual_m;                 /* # of row blocks per column */
+      for (ic=0; ic<n; ic++){            /* loop over column */
+	for (n_r=0; n_r<n_rb; n_r++){    /* loop within one column */
+	  /* store in dummy variables */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){
+	    ir = ir + sample_f;
+	    xdummyl[i] = mat(yl, ir, ic);
+	    xdummyh[i] = mat(yh, ir,c_o_a+ic);
+	  }
+	  /* perform filtering: first LL/LH, then HL/HH */
+	  fpconv(xdummyl, actual_m, h0, h1, lh, ydummyll, ydummylh);
+	  fpconv(xdummyh, actual_m, h0, h1, lh, ydummyhl, ydummyhh);
+	  /* restore dummy variables in matrices */
+	  ir = -sample_f + n_r;
+	  for (i=0; i<actual_m; i++){
+	    ir = ir + sample_f;
+	    mat(yl, ir, ic) = ydummyll[i];
+	    mat(yh, ir, c_o_a+ic) = ydummylh[i];
+	    mat(yh, ir,c_o_a+n+ic) = ydummyhl[i];
+	    mat(yh, ir, c_o_a_p2n+ic) = ydummyhh[i];
+	  }
+	}
+      }
+    }
+    sample_f = sample_f*2;
+  }
+}
+
+#ifdef __STDC__
+fpconv(double *x_in, uintptr_t lx, double *h0, double *h1, uintptr_t lh,
+       double *x_outl, double *x_outh)
+#else
+fpconv(x_in, lx, h0, h1, lh, x_outl, x_outh)
+double *x_in, *h0, *h1, *x_outl, *x_outh;
+uintptr_t lx, lh;
+#endif
+{
+  uintptr_t i, j;
+  double x0, x1;
+
+  for (i=lx; i < lx+lh-1; i++)
+    x_in[i] = x_in[i-lx];
+  for (i=0; i<lx; i++){
+    x0 = 0;
+    x1 = 0;
+    for (j=0; j<lh; j++){
+      x0 = x0 + x_in[j+i]*h0[lh-1-j];
+      x1 = x1 + x_in[j+i]*h1[lh-1-j];
+    }
+    x_outl[i] = x0;
+    x_outh[i] = x1;
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/Rice Wavelet Toolbox/setopt.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,70 @@
+function option = setopt(opt_par,default);
+%    option = setopt(opt_par,default);
+%
+%    SETOPT can modify a default option vector with user specified options.
+%
+%    Input:
+%       opt_par : Users desired option vector
+%       default : Program default option vector
+%
+%    Output:
+%       option : New option vector
+%
+%    Example:
+%       opt_par = [1 2 3 4];
+%       default = [1 1 1 1];
+%       option = setopt(opt_par,default)
+%       option = 1     2     3     4
+%
+
+%File Name: setopt.m
+%Last Modification Date: 10/26/95	13:06:18
+%Current Version: setopt.m	2.4
+%File Creation Date: Thu Nov 11 10:40:22 1993
+%Author: Jan Erik Odegard  <odegard@ece.rice.edu>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Jan Erik Odegard, Department of ECE, Rice University.
+%
+%This software is distributed and licensed to you on a non-exclusive
+%basis, free-of-charge. Redistribution and use in source and binary forms,
+%with or without modification, are permitted provided that the following
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice,
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software
+%   must display the following acknowledgment: This product includes
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors
+%   may be used to endorse or promote products derived from this software
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of
+%Technology Transfer at techtran@rice.edu or (713) 348-6173
+%
+%Change History:
+%
+
+if (nargin < 2)
+  error('You need two inputs');
+end;
+len = length(opt_par);
+option = zeros(size(default));
+option(1:len) = opt_par(1:len);
+option = option + (option == 0).*default;