%
% wavelete smoothing code is obtained at:
% http://bioinformatics.mdanderson.org/cromwell.html
%
%%% spectrum must be a row vector 
function [smoothSpec] = smooth_wavelet(spectrum, t, saturate, L)
% [B S] = waveletSmoothAndBaseleinCorrect(spectrum, threshold, saturation, L)
% computes a denoised spectrum using a wavelet transform, and a baseline corrected
% denoised spectrum using a monotone minimum from a raw spectrum. The wavelet 
% transform is performed using code supplied as part of the Rice Wavelet Toolbox
% (http://www-dsp.rice.edu/software/rwt.shtml). The wavelet transform will only
% work on input vectors whose length is a multiple of a power of 2. To get the
% lengths right, we reverse the end of the spectrum and pad to the next multiple,
% and then truncate the processed spectrum to the original length.
%
% Inputs:
%	spectrum = vector of intensities, assumed to be collected at consecutive time points
%	threshold = multiple of 0.67*MAD at which to threshold the wavelet coefficients
%		we typically use thresholds in the range of 4-10.
%	saturation = index at which saturation ends and baseline correction should start
%	L = power of 2 to use to compute wavelets (default value is L = 10, which says
%		to use blocks of size 1024)
% Outputs:
%	B = vector of denoised and baseline corrected intensities
%	S = vector of denoised intensities
%
% Example:
%	[B S] = waveletSmoothAndBaselineCorrect(spectrum, 6, 2500, 10);
%	noise = spectrum - S;
%	plot(spectrum)
%	hold on
%	plot(S, 'r')

% Copyright (c) 2003, 2004 UT M.D. Anderson Cancer Center. All rights reserved.
% See the accompanying file "license.txt" for details.

% Original version by Kevin Coombes
% $Revision: 2.2 $
% Last modified by $Author: krc $ on $Date: 2004/03/01 21:41:15 $.

if nargin < 4,
   L = 10;				% level for the wavelet denoising
end
if nargin < 3,
   saturate = 2500; % parameter for baseline correction
end

clip = length(spectrum);
spectrum = pad(spectrum, L);	% pad to the next power of 2^l by reflection so the algorithm doesn't barf

% Notes on the wavelet smooth parameters to "denoise"
%	(1) we take q = 8 for the degree of the Daubechies wavelet.
%	(2) the third parameter says to use the undecimated discrete wavelet transform
%		 instead of the regular DWT
%	(3) in the vector
%			options(1) = 0 says to threshold the low pass component. This is irrelevant.
%			options(2) = t is the threshold multiple; this may need to be tweaked
%			options(3) = 0	says to use multiples of the MAD (not the SD). Do not change!
%			options(4) = 1	says to use hard rather than soft thresholding. Do not change!
%			options(5) = L is the level for the wavelet transform
%			options(6) = 0 says to use the threshold as described above; do not change!
smoothSpec = denoise(spectrum, daubcqf(8), 1, [1 t 0 1 L 0]);
smoothSpec = smoothSpec(1:clip);	% restore the original spectrum length

% %%% baseline correction via smart monotone minimum
% baselineCorrected = smoothSpec - monotoneMinimum(smoothSpec, saturate);


