annotate util/Rice Wavelet Toolbox/denoise.m @ 78:f69ae88b8be5

added Rice Wavelet Toolbox with my modification, so it can be compiled on newer systems.
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Fri, 25 Mar 2011 15:27:33 +0000
parents
children
rev   line source
ivan@78 1 function [xd,xn,option] = denoise(x,h,type,option)
ivan@78 2 % [xd,xn,option] = denoise(x,h,type,option);
ivan@78 3 %
ivan@78 4 % DENOISE is a generic program for wavelet based denoising.
ivan@78 5 % The program will denoise the signal x using the 2-band wavelet
ivan@78 6 % system described by the filter h using either the traditional
ivan@78 7 % discrete wavelet transform (DWT) or the linear shift invariant
ivan@78 8 % discrete wavelet transform (also known as the undecimated DWT
ivan@78 9 % (UDWT)).
ivan@78 10 %
ivan@78 11 % Input:
ivan@78 12 % x : 1D or 2D signal to be denoised
ivan@78 13 % h : Scaling filter to be applied
ivan@78 14 % type : Type of transform (Default: type = 0)
ivan@78 15 % 0 --> Discrete wavelet transform (DWT)
ivan@78 16 % 1 --> Undecimated DWT (UDWT)
ivan@78 17 % option : Default settings is marked with '*':
ivan@78 18 % *type = 0 --> option = [0 3.0 0 0 0 0]
ivan@78 19 % type = 1 --> option = [0 3.6 0 1 0 0]
ivan@78 20 % option(1) : Whether to threshold low-pass part
ivan@78 21 % 0 --> Don't threshold low pass component
ivan@78 22 % 1 --> Threshold low pass component
ivan@78 23 % option(2) : Threshold multiplier, c. The threshold is
ivan@78 24 % computed as:
ivan@78 25 % thld = c*MAD(noise_estimate)).
ivan@78 26 % The default values are:
ivan@78 27 % c = 3.0 for the DWT based denoising
ivan@78 28 % c = 3.6 for the UDWT based denoising
ivan@78 29 % option(3) : Type of variance estimator
ivan@78 30 % 0 --> MAD (mean absolute deviation)
ivan@78 31 % 1 --> STD (classical numerical std estimate)
ivan@78 32 % option(4) : Type of thresholding
ivan@78 33 % 0 --> Soft thresholding
ivan@78 34 % 1 --> Hard thresholding
ivan@78 35 % option(5) : Number of levels, L, in wavelet decomposition. By
ivan@78 36 % setting this to the default value '0' a maximal
ivan@78 37 % decomposition is used.
ivan@78 38 % option(6) : Actual threshold to use (setting this to
ivan@78 39 % anything but 0 will mean that option(3)
ivan@78 40 % is ignored)
ivan@78 41 %
ivan@78 42 % Output:
ivan@78 43 % xd : Estimate of noise free signal
ivan@78 44 % xn : The estimated noise signal (x-xd)
ivan@78 45 % option : A vector of actual parameters used by the
ivan@78 46 % program. The vector is configured the same way as
ivan@78 47 % the input option vector with one added element
ivan@78 48 % option(7) = type.
ivan@78 49 %
ivan@78 50 % HERE'S AN EASY WAY TO RUN THE EXAMPLES:
ivan@78 51 % Cut-and-paste the example you want to run to a new file
ivan@78 52 % called ex.m, for example. Delete out the % at the beginning
ivan@78 53 % of each line in ex.m (Can use search-and-replace in your editor
ivan@78 54 % to replace it with a space). Type 'ex' in matlab and hit return.
ivan@78 55 %
ivan@78 56 % Example 1:
ivan@78 57 % h = daubcqf(6); [s,N] = makesig('Doppler'); n = randn(1,N);
ivan@78 58 % x = s + n/10; % (approximately 10dB SNR)
ivan@78 59 % figure;plot(x);hold on;plot(s,'r');
ivan@78 60 %
ivan@78 61 % %Denoise x with the default method based on the DWT
ivan@78 62 % [xd,xn,opt1] = denoise(x,h);
ivan@78 63 % figure;plot(xd);hold on;plot(s,'r');
ivan@78 64 %
ivan@78 65 % %Denoise x using the undecimated (LSI) wavelet transform
ivan@78 66 % [yd,yn,opt2] = denoise(x,h,1);
ivan@78 67 % figure;plot(yd);hold on;plot(s,'r');
ivan@78 68 %
ivan@78 69 % Example 2: (on an image)
ivan@78 70 % h = daubcqf(6); load lena;
ivan@78 71 % noisyLena = lena + 25 * randn(size(lena));
ivan@78 72 % figure; colormap(gray); imagesc(lena); title('Original Image');
ivan@78 73 % figure; colormap(gray); imagesc(noisyLena); title('Noisy Image');
ivan@78 74 % Denoise lena with the default method based on the DWT
ivan@78 75 % [denoisedLena,xn,opt1] = denoise(noisyLena,h);
ivan@78 76 % figure; colormap(gray); imagesc(denoisedLena); title('denoised Image');
ivan@78 77 %
ivan@78 78 %
ivan@78 79 % See also: mdwt, midwt, mrdwt, mirdwt, SoftTh, HardTh, setopt
ivan@78 80 %
ivan@78 81
ivan@78 82 %File Name: denoise.m
ivan@78 83 %Last Modification Date: 04/15/97 10:44:28
ivan@78 84 %Current Version: denoise.m 2.4
ivan@78 85 %File Creation Date: Mon Feb 20 08:33:15 1995
ivan@78 86 %Author: Jan Erik Odegard <odegard@ece.rice.edu>
ivan@78 87 %
ivan@78 88 %Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
ivan@78 89 %Created by Jan Erik Odegard, Department of ECE, Rice University.
ivan@78 90 %
ivan@78 91 %This software is distributed and licensed to you on a non-exclusive
ivan@78 92 %basis, free-of-charge. Redistribution and use in source and binary forms,
ivan@78 93 %with or without modification, are permitted provided that the following
ivan@78 94 %conditions are met:
ivan@78 95 %
ivan@78 96 %1. Redistribution of source code must retain the above copyright notice,
ivan@78 97 % this list of conditions and the following disclaimer.
ivan@78 98 %2. Redistribution in binary form must reproduce the above copyright notice,
ivan@78 99 % this list of conditions and the following disclaimer in the
ivan@78 100 % documentation and/or other materials provided with the distribution.
ivan@78 101 %3. All advertising materials mentioning features or use of this software
ivan@78 102 % must display the following acknowledgment: This product includes
ivan@78 103 % software developed by Rice University, Houston, Texas and its contributors.
ivan@78 104 %4. Neither the name of the University nor the names of its contributors
ivan@78 105 % may be used to endorse or promote products derived from this software
ivan@78 106 % without specific prior written permission.
ivan@78 107 %
ivan@78 108 %THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS,
ivan@78 109 %AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
ivan@78 110 %BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
ivan@78 111 %FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY
ivan@78 112 %OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
ivan@78 113 %EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
ivan@78 114 %PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
ivan@78 115 %OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
ivan@78 116 %WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
ivan@78 117 %OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE
ivan@78 118 %USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
ivan@78 119 %
ivan@78 120 %For information on commercial licenses, contact Rice University's Office of
ivan@78 121 %Technology Transfer at techtran@rice.edu or (713) 348-6173
ivan@78 122 %
ivan@78 123 %Change History: Fixed output of function and an error in the computation
ivan@78 124 % of the threshold for redundant denoising.
ivan@78 125 % <Jan Erik Odegard> <Mon Jul 31, 1995>
ivan@78 126 %
ivan@78 127 % This code is composed of several of our old codes for
ivan@78 128 % wavelet based denoising. In an effort to make the mess
ivan@78 129 % more manageable we decided to create on code that would
ivan@78 130 % handle all the various wavelet based denoising methods.
ivan@78 131 % However, only time will show (as we discover new and
ivan@78 132 % improved forms of denoising) if we can succeed in our goals.
ivan@78 133 % <Jan Erik Odegard> <Thu May 11, 1995>
ivan@78 134 %
ivan@78 135
ivan@78 136 if(nargin < 2)
ivan@78 137 error('You need to provide at least 2 inputs: x and h');
ivan@78 138 end;
ivan@78 139 if(nargin < 3),
ivan@78 140 type = 0;
ivan@78 141 option = [];
ivan@78 142 elseif(nargin < 4)
ivan@78 143 option = [];
ivan@78 144 end;
ivan@78 145 if(isempty(type)),
ivan@78 146 type = 0;
ivan@78 147 end;
ivan@78 148 if(type == 0),
ivan@78 149 default_opt = [0 3.0 0 0 0 0];
ivan@78 150 elseif(type == 1),
ivan@78 151 default_opt = [0 3.6 0 1 0 0];
ivan@78 152 else,
ivan@78 153 error(['Unknown denoising method',10,...
ivan@78 154 'If it is any good we need to have a serious talk :-)']);
ivan@78 155 end;
ivan@78 156 option = setopt(option,default_opt);
ivan@78 157 [mx,nx] = size(x);
ivan@78 158 dim = min(mx,nx);
ivan@78 159 if(dim == 1),
ivan@78 160 n = max(mx,nx);
ivan@78 161 else,
ivan@78 162 n = dim;
ivan@78 163 end;
ivan@78 164 if(option(5) == 0),
ivan@78 165 L = floor(log2(n));
ivan@78 166 else
ivan@78 167 L = option(5);
ivan@78 168 end;
ivan@78 169 if(type == 0), % Denoising by DWT
ivan@78 170 xd = mdwt(x,h,L);
ivan@78 171 if (option(6) == 0),
ivan@78 172 tmp = xd(floor(mx/2)+1:mx,floor(nx/2)+1:nx);
ivan@78 173 if(option(3) == 0),
ivan@78 174 thld = option(2)*median(abs(tmp(:)))/.67;
ivan@78 175 elseif(option(3) == 1),
ivan@78 176 thld = option(2)*std(tmp(:));
ivan@78 177 else
ivan@78 178 error('Unknown threshold estimator, Use either MAD or STD');
ivan@78 179 end;
ivan@78 180 else,
ivan@78 181 thld = option(6);
ivan@78 182 end;
ivan@78 183 if(dim == 1)
ivan@78 184 ix = 1:n/(2^L);
ivan@78 185 ykeep = xd(ix);
ivan@78 186 else
ivan@78 187 ix = 1:mx/(2^L);
ivan@78 188 jx = 1:nx/(2^L);
ivan@78 189 ykeep = xd(ix,jx);
ivan@78 190 end;
ivan@78 191 if(option(4) == 0),
ivan@78 192 xd = SoftTh(xd,thld);
ivan@78 193 elseif(option(4) == 1),
ivan@78 194 xd = HardTh(xd,thld);
ivan@78 195 else,
ivan@78 196 error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
ivan@78 197 end;
ivan@78 198 if (option(1) == 0),
ivan@78 199 if(dim == 1),
ivan@78 200 xd(ix) = ykeep;
ivan@78 201 else,
ivan@78 202 xd(ix,jx) = ykeep;
ivan@78 203 end;
ivan@78 204 end;
ivan@78 205 xd = midwt(xd,h,L);
ivan@78 206 elseif(type == 1), % Denoising by UDWT
ivan@78 207 [xl,xh] = mrdwt(x,h,L);
ivan@78 208 if(dim == 1),
ivan@78 209 c_offset = 1;
ivan@78 210 else,
ivan@78 211 c_offset = 2*nx + 1;
ivan@78 212 end;
ivan@78 213 if (option(6) == 0),
ivan@78 214 tmp = xh(:,c_offset:c_offset+nx-1);
ivan@78 215 if(option(3) == 0),
ivan@78 216 thld = option(2)*median(abs(tmp(:)))/.67;
ivan@78 217 elseif(option(3) == 1),
ivan@78 218 thld = option(2)*std(tmp(:));
ivan@78 219 else
ivan@78 220 error('Unknown threshold estimator, Use either MAD or STD');
ivan@78 221 end;
ivan@78 222 else,
ivan@78 223 thld = option(6);
ivan@78 224 end;
ivan@78 225 if(option(4) == 0),
ivan@78 226 xh = SoftTh(xh,thld);
ivan@78 227 if(option(1) == 1),
ivan@78 228 xl = SoftTh(xl,thld);
ivan@78 229 end;
ivan@78 230 elseif(option(4) == 1),
ivan@78 231 xh = HardTh(xh,thld);
ivan@78 232 if(option(1) == 1),
ivan@78 233 xl = HardTh(xl,thld);
ivan@78 234 end;
ivan@78 235 else,
ivan@78 236 error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
ivan@78 237 end;
ivan@78 238 xd = mirdwt(xl,xh,h,L);
ivan@78 239 else, % Denoising by unknown method
ivan@78 240 error(['Unknown denoising method',10,...
ivan@78 241 'If it is any good we need to have a serious talk :-)']);
ivan@78 242 end;
ivan@78 243 option(6) = thld;
ivan@78 244 option(7) = type;
ivan@78 245 xn = x - xd;