ivan@78: function [xd,xn,option] = denoise(x,h,type,option) ivan@78: % [xd,xn,option] = denoise(x,h,type,option); ivan@78: % ivan@78: % DENOISE is a generic program for wavelet based denoising. ivan@78: % The program will denoise the signal x using the 2-band wavelet ivan@78: % system described by the filter h using either the traditional ivan@78: % discrete wavelet transform (DWT) or the linear shift invariant ivan@78: % discrete wavelet transform (also known as the undecimated DWT ivan@78: % (UDWT)). ivan@78: % ivan@78: % Input: ivan@78: % x : 1D or 2D signal to be denoised ivan@78: % h : Scaling filter to be applied ivan@78: % type : Type of transform (Default: type = 0) ivan@78: % 0 --> Discrete wavelet transform (DWT) ivan@78: % 1 --> Undecimated DWT (UDWT) ivan@78: % option : Default settings is marked with '*': ivan@78: % *type = 0 --> option = [0 3.0 0 0 0 0] ivan@78: % type = 1 --> option = [0 3.6 0 1 0 0] ivan@78: % option(1) : Whether to threshold low-pass part ivan@78: % 0 --> Don't threshold low pass component ivan@78: % 1 --> Threshold low pass component ivan@78: % option(2) : Threshold multiplier, c. The threshold is ivan@78: % computed as: ivan@78: % thld = c*MAD(noise_estimate)). ivan@78: % The default values are: ivan@78: % c = 3.0 for the DWT based denoising ivan@78: % c = 3.6 for the UDWT based denoising ivan@78: % option(3) : Type of variance estimator ivan@78: % 0 --> MAD (mean absolute deviation) ivan@78: % 1 --> STD (classical numerical std estimate) ivan@78: % option(4) : Type of thresholding ivan@78: % 0 --> Soft thresholding ivan@78: % 1 --> Hard thresholding ivan@78: % option(5) : Number of levels, L, in wavelet decomposition. By ivan@78: % setting this to the default value '0' a maximal ivan@78: % decomposition is used. ivan@78: % option(6) : Actual threshold to use (setting this to ivan@78: % anything but 0 will mean that option(3) ivan@78: % is ignored) ivan@78: % ivan@78: % Output: ivan@78: % xd : Estimate of noise free signal ivan@78: % xn : The estimated noise signal (x-xd) ivan@78: % option : A vector of actual parameters used by the ivan@78: % program. The vector is configured the same way as ivan@78: % the input option vector with one added element ivan@78: % option(7) = type. ivan@78: % ivan@78: % HERE'S AN EASY WAY TO RUN THE EXAMPLES: ivan@78: % Cut-and-paste the example you want to run to a new file ivan@78: % called ex.m, for example. Delete out the % at the beginning ivan@78: % of each line in ex.m (Can use search-and-replace in your editor ivan@78: % to replace it with a space). Type 'ex' in matlab and hit return. ivan@78: % ivan@78: % Example 1: ivan@78: % h = daubcqf(6); [s,N] = makesig('Doppler'); n = randn(1,N); ivan@78: % x = s + n/10; % (approximately 10dB SNR) ivan@78: % figure;plot(x);hold on;plot(s,'r'); ivan@78: % ivan@78: % %Denoise x with the default method based on the DWT ivan@78: % [xd,xn,opt1] = denoise(x,h); ivan@78: % figure;plot(xd);hold on;plot(s,'r'); ivan@78: % ivan@78: % %Denoise x using the undecimated (LSI) wavelet transform ivan@78: % [yd,yn,opt2] = denoise(x,h,1); ivan@78: % figure;plot(yd);hold on;plot(s,'r'); ivan@78: % ivan@78: % Example 2: (on an image) ivan@78: % h = daubcqf(6); load lena; ivan@78: % noisyLena = lena + 25 * randn(size(lena)); ivan@78: % figure; colormap(gray); imagesc(lena); title('Original Image'); ivan@78: % figure; colormap(gray); imagesc(noisyLena); title('Noisy Image'); ivan@78: % Denoise lena with the default method based on the DWT ivan@78: % [denoisedLena,xn,opt1] = denoise(noisyLena,h); ivan@78: % figure; colormap(gray); imagesc(denoisedLena); title('denoised Image'); ivan@78: % ivan@78: % ivan@78: % See also: mdwt, midwt, mrdwt, mirdwt, SoftTh, HardTh, setopt ivan@78: % ivan@78: ivan@78: %File Name: denoise.m ivan@78: %Last Modification Date: 04/15/97 10:44:28 ivan@78: %Current Version: denoise.m 2.4 ivan@78: %File Creation Date: Mon Feb 20 08:33:15 1995 ivan@78: %Author: Jan Erik Odegard ivan@78: % ivan@78: %Copyright (c) 2000 RICE UNIVERSITY. All rights reserved. ivan@78: %Created by Jan Erik Odegard, Department of ECE, Rice University. ivan@78: % ivan@78: %This software is distributed and licensed to you on a non-exclusive ivan@78: %basis, free-of-charge. Redistribution and use in source and binary forms, ivan@78: %with or without modification, are permitted provided that the following ivan@78: %conditions are met: ivan@78: % ivan@78: %1. Redistribution of source code must retain the above copyright notice, ivan@78: % this list of conditions and the following disclaimer. ivan@78: %2. Redistribution in binary form must reproduce the above copyright notice, ivan@78: % this list of conditions and the following disclaimer in the ivan@78: % documentation and/or other materials provided with the distribution. ivan@78: %3. All advertising materials mentioning features or use of this software ivan@78: % must display the following acknowledgment: This product includes ivan@78: % software developed by Rice University, Houston, Texas and its contributors. ivan@78: %4. Neither the name of the University nor the names of its contributors ivan@78: % may be used to endorse or promote products derived from this software ivan@78: % without specific prior written permission. ivan@78: % ivan@78: %THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS, ivan@78: %AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, ivan@78: %BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS ivan@78: %FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY ivan@78: %OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, ivan@78: %EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, ivan@78: %PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ivan@78: %OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, ivan@78: %WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR ivan@78: %OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE ivan@78: %USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ivan@78: % ivan@78: %For information on commercial licenses, contact Rice University's Office of ivan@78: %Technology Transfer at techtran@rice.edu or (713) 348-6173 ivan@78: % ivan@78: %Change History: Fixed output of function and an error in the computation ivan@78: % of the threshold for redundant denoising. ivan@78: % ivan@78: % ivan@78: % This code is composed of several of our old codes for ivan@78: % wavelet based denoising. In an effort to make the mess ivan@78: % more manageable we decided to create on code that would ivan@78: % handle all the various wavelet based denoising methods. ivan@78: % However, only time will show (as we discover new and ivan@78: % improved forms of denoising) if we can succeed in our goals. ivan@78: % ivan@78: % ivan@78: ivan@78: if(nargin < 2) ivan@78: error('You need to provide at least 2 inputs: x and h'); ivan@78: end; ivan@78: if(nargin < 3), ivan@78: type = 0; ivan@78: option = []; ivan@78: elseif(nargin < 4) ivan@78: option = []; ivan@78: end; ivan@78: if(isempty(type)), ivan@78: type = 0; ivan@78: end; ivan@78: if(type == 0), ivan@78: default_opt = [0 3.0 0 0 0 0]; ivan@78: elseif(type == 1), ivan@78: default_opt = [0 3.6 0 1 0 0]; ivan@78: else, ivan@78: error(['Unknown denoising method',10,... ivan@78: 'If it is any good we need to have a serious talk :-)']); ivan@78: end; ivan@78: option = setopt(option,default_opt); ivan@78: [mx,nx] = size(x); ivan@78: dim = min(mx,nx); ivan@78: if(dim == 1), ivan@78: n = max(mx,nx); ivan@78: else, ivan@78: n = dim; ivan@78: end; ivan@78: if(option(5) == 0), ivan@78: L = floor(log2(n)); ivan@78: else ivan@78: L = option(5); ivan@78: end; ivan@78: if(type == 0), % Denoising by DWT ivan@78: xd = mdwt(x,h,L); ivan@78: if (option(6) == 0), ivan@78: tmp = xd(floor(mx/2)+1:mx,floor(nx/2)+1:nx); ivan@78: if(option(3) == 0), ivan@78: thld = option(2)*median(abs(tmp(:)))/.67; ivan@78: elseif(option(3) == 1), ivan@78: thld = option(2)*std(tmp(:)); ivan@78: else ivan@78: error('Unknown threshold estimator, Use either MAD or STD'); ivan@78: end; ivan@78: else, ivan@78: thld = option(6); ivan@78: end; ivan@78: if(dim == 1) ivan@78: ix = 1:n/(2^L); ivan@78: ykeep = xd(ix); ivan@78: else ivan@78: ix = 1:mx/(2^L); ivan@78: jx = 1:nx/(2^L); ivan@78: ykeep = xd(ix,jx); ivan@78: end; ivan@78: if(option(4) == 0), ivan@78: xd = SoftTh(xd,thld); ivan@78: elseif(option(4) == 1), ivan@78: xd = HardTh(xd,thld); ivan@78: else, ivan@78: error('Unknown threshold rule. Use either Soft (0) or Hard (1)'); ivan@78: end; ivan@78: if (option(1) == 0), ivan@78: if(dim == 1), ivan@78: xd(ix) = ykeep; ivan@78: else, ivan@78: xd(ix,jx) = ykeep; ivan@78: end; ivan@78: end; ivan@78: xd = midwt(xd,h,L); ivan@78: elseif(type == 1), % Denoising by UDWT ivan@78: [xl,xh] = mrdwt(x,h,L); ivan@78: if(dim == 1), ivan@78: c_offset = 1; ivan@78: else, ivan@78: c_offset = 2*nx + 1; ivan@78: end; ivan@78: if (option(6) == 0), ivan@78: tmp = xh(:,c_offset:c_offset+nx-1); ivan@78: if(option(3) == 0), ivan@78: thld = option(2)*median(abs(tmp(:)))/.67; ivan@78: elseif(option(3) == 1), ivan@78: thld = option(2)*std(tmp(:)); ivan@78: else ivan@78: error('Unknown threshold estimator, Use either MAD or STD'); ivan@78: end; ivan@78: else, ivan@78: thld = option(6); ivan@78: end; ivan@78: if(option(4) == 0), ivan@78: xh = SoftTh(xh,thld); ivan@78: if(option(1) == 1), ivan@78: xl = SoftTh(xl,thld); ivan@78: end; ivan@78: elseif(option(4) == 1), ivan@78: xh = HardTh(xh,thld); ivan@78: if(option(1) == 1), ivan@78: xl = HardTh(xl,thld); ivan@78: end; ivan@78: else, ivan@78: error('Unknown threshold rule. Use either Soft (0) or Hard (1)'); ivan@78: end; ivan@78: xd = mirdwt(xl,xh,h,L); ivan@78: else, % Denoising by unknown method ivan@78: error(['Unknown denoising method',10,... ivan@78: 'If it is any good we need to have a serious talk :-)']); ivan@78: end; ivan@78: option(6) = thld; ivan@78: option(7) = type; ivan@78: xn = x - xd;