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 / basecorr2.m < prev    next >
Text File  |  2007-12-04  |  3KB  |  94 lines
  1. function bc = basecorr2(raw, rawMass, M, L, thld)
  2. %
  3. % bc = basecorr2(raw, rawMass, M, L, thld)
  4. % wrapper to perform wavelet denoising, baseline correction, and
  5. % peak finding on a collection of spectra.
  6. %
  7. % inputs are:
  8. %    raw        % matrix of intensities
  9. %    rawMass    % vector of masses
  10. %    M            % number of data points to ignore at beginning
  11. %    L            % wavelet level
  12. %    thld        % multiple of MAD at which to set the threshold
  13. % output is:
  14. %    bc            % matrix of baseline corrected intensities
  15. %
  16. % side effects:
  17. % creates the collowing binary MATLAB (*.mat) files in the current directory
  18. %    baselineCorrected    % the matrix bc
  19. %    noises                % matrix of noise estimates
  20. %    normalized            % matrix of normalized baselin corrected intensities
  21. %    meanBase                % vector containing mean baseline corrected spectrum
  22. %    meanNorm                % vector containing mean normalized spectrum
  23. %    peakLocations        % matrix of peaks (uses trivialPeakFinder)
  24. %    numpeaks                % vector of the number of peaks pe spectrum
  25.  
  26. % Copyright (c) 2003, 2004 UT M.D. Anderson Cancer Center. All rights reserved.
  27. % See the accompanying file "license.txt" for details.
  28.  
  29. w = size(raw);                % how much data are we handling?
  30. n = floor(w(2)/2^L);        % number of length 2^L pieces available
  31.  
  32. echo on
  33. status = 'start baseline correction'
  34. tic
  35. echo off
  36. %bc = repmat(0, [w(1) n*2^L]);
  37. %noises = repmat(0, [w(1) n*2^L]);
  38. bc = repmat(0, [w(1) w(2)]);
  39. noises = repmat(0, [w(1) w(2)]);
  40. for i = 1:w(1),
  41.    [b, s] = waveletSmoothAndBaselineCorrect(raw(i,:), thld, M, L);
  42.    bc(i,:) = b;
  43. %   noises(i, :) = raw(i, 1:(n*2^L)) - s;
  44.    noises(i, :) = raw(i, :) - s;
  45. end
  46. echo on
  47. toc
  48. status = 'end baseline correction'
  49. echo off
  50.  
  51. lm = length(rawMass);
  52. bc = bc(:, 1:lm);
  53. noises = noises(:, 1:lm);
  54. save baselineCorrected bc -V4
  55. save noises noises -V4
  56. sc = sum(bc(:, M:end)')';
  57.  
  58. nc = bc ./ repmat(sc, [1 lm]) * 10^4;
  59. echo on
  60. status = 'finished normalizing'
  61. echo off
  62.  
  63. save normalized nc -V4
  64.  
  65. meanBase = mean(bc);
  66. meanNorm = mean(nc);
  67.  
  68. save meanBase meanBase -V4
  69. save meanNorm meanNorm -V4
  70.  
  71. echo on
  72. status = 'start peak finding'
  73. tic
  74. echo off
  75. peaks = repmat(0, [w(1) 1000]);
  76. numpeaks = repmat(0, [w(1) 1]);
  77. for i = 1:w(1),
  78.    i
  79.    to = trivialPeakFinder(nc(i,:));
  80.    numpeaks(i) = length(to);
  81.    if length(to) > 500,
  82.       to = to(1:500);
  83.    end
  84.    peaks(i, 1:length(to)) = to;
  85. end
  86. echo on
  87. toc
  88. peaks = peaks(:, 1:max(numpeaks));
  89. status = 'end peak finding'
  90. echo off
  91.  
  92. save peakLocations peaks -V4
  93. save numpeaks numpeaks -V4
  94.