denoise.m
上传用户:speoil
上传日期:2022-06-23
资源大小:224k
文件大小:9k
- 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;