home *** CD-ROM | disk | FTP | other *** search
/ ftp.disi.unige.it / 2015-02-11.ftp.disi.unige.it.tar / ftp.disi.unige.it / pub / .person / BarlaA / sw / matlab / Cromwell / waveletSmoothAndBaselineCorrect.m < prev   
Text File  |  2007-12-04  |  3KB  |  64 lines
  1. function [baselineCorrected, smoothSpec] = waveletSmoothAndBaselineCorrect(spectrum, t, saturate, L)
  2. % [B S] = waveletSmoothAndBaseleinCorrect(spectrum, threshold, saturation, L)
  3. % computes a denoised spectrum using a wavelet transform, and a baseline corrected
  4. % denoised spectrum using a monotone minimum from a raw spectrum. The wavelet 
  5. % transform is performed using code supplied as part of the Rice Wavelet Toolbox
  6. % (http://www-dsp.rice.edu/software/rwt.shtml). The wavelet transform will only
  7. % work on input vectors whose length is a multiple of a power of 2. To get the
  8. % lengths right, we reverse the end of the spectrum and pad to the next multiple,
  9. % and then truncate the processed spectrum to the original length.
  10. %
  11. % Inputs:
  12. %    spectrum = vector of intensities, assumed to be collected at consecutive time points
  13. %    threshold = multiple of 0.67*MAD at which to threshold the wavelet coefficients
  14. %        we typically use thresholds in the range of 4-10.
  15. %    saturation = index at which saturation ends and baseline correction should start
  16. %    L = power of 2 to use to compute wavelets (default value is L = 10, which says
  17. %        to use blocks of size 1024)
  18. % Outputs:
  19. %    B = vector of denoised and baseline corrected intensities
  20. %    S = vector of denoised intensities
  21. %
  22. % Example:
  23. %    [B S] = waveletSmoothAndBaselineCorrect(spectrum, 6, 2500, 10);
  24. %    noise = spectrum - S;
  25. %    plot(spectrum)
  26. %    hold on
  27. %    plot(S, 'r')
  28.  
  29. % Copyright (c) 2003, 2004 UT M.D. Anderson Cancer Center. All rights reserved.
  30. % See the accompanying file "license.txt" for details.
  31.  
  32. % Original version by Kevin Coombes
  33. % $Revision: 2.2 $
  34. % Last modified by $Author: krc $ on $Date: 2004/03/01 21:41:15 $.
  35.  
  36. if nargin < 4,
  37.    L = 10;                % level for the wavelet denoising
  38. end
  39. if nargin < 3,
  40.    saturate = 2500; % parameter for baseline correction
  41. end
  42.  
  43. clip = length(spectrum);
  44. spectrum = pad(spectrum, L);    % pad to the next power of 2^l by reflection so the algorithm doesn't barf
  45.  
  46. % Notes on the wavelet smooth parameters to "denoise"
  47. %    (1) we take q = 8 for the degree of the Daubechies wavelet.
  48. %    (2) the third parameter says to use the undecimated discrete wavelet transform
  49. %         instead of the regular DWT
  50. %    (3) in the vector
  51. %            options(1) = 0 says to threshold the low pass component. This is irrelevant.
  52. %            options(2) = t is the threshold multiple; this may need to be tweaked
  53. %            options(3) = 0    says to use multiples of the MAD (not the SD). Do not change!
  54. %            options(4) = 1    says to use hard rather than soft thresholding. Do not change!
  55. %            options(5) = L is the level for the wavelet transform
  56. %            options(6) = 0 says to use the threshold as described above; do not change!
  57. smoothSpec = denoise(spectrum, daubcqf(8), 1, [1 t 0 1 L 0]);
  58. smoothSpec = smoothSpec(1:clip);    % restore the original spectrum length
  59.  
  60. %%% baseline correction via smart monotone minimum
  61. baselineCorrected = smoothSpec - monotoneMinimum(smoothSpec, saturate);
  62.  
  63.  
  64.