function [amps, compspec, sumspec] = spectral_decomp(wl,ap,wlunc,apunc,acs)
% AC-S "un-smoothing" and spectral decomposition method
%
% Ron Zaneveld, WET Labs, Inc., 2005
% Ali Chase, University of Maine, 2014
%
% Vectorized: Guillaume Bourdin, University of Maine, 2020
%
% See the following publication for details on the method:
% Chase, A., et al., Decomposition of in situ particulate absorption
% spectra. Methods in Oceanography (2014), http://dx.doi.org/10.1016/j.mio.2014.02.022
%
% INPUTS:
% wl - wavelengths associated with measured particulate absorption
% spectra
% ap - measured particulate absorption spectra
% wlunc - wavelengths associated with uncertainty of measured
% particulate absorption spectra
% apunc - uncertainty of measured particulate absorption spectra
% acs - 1 (for AC-S data) or 0 (for a different instrument). This code
% was designed for use with particulate absorption spectra measured
% using a WETLabs AC-S instrument, and thus it applies a spectral
% un-smoothing that is specific to that instrument. To use with
% other types of absorption data, input "0" and the correction
% applied for AC-S will be bypassed.
%
% The uncertainty values we use are the standard deviation of one-minute
% binned particulate absorption spectra. If the uncertainty is unknown,
% then all wavelenghts will be treated equally, i.e. no spectral weighting
% will be applied. In this case populate apunc with -999 values.
%
% OUTPUTS:
% amps - the amplitudes of 12 Gaussian function peaks and one non-algal
% particle (NAP) function
% compspec - the component spectra after spectral decomposition; the
% Gaussian and NAP functions multiplied by the 'amps'
% values
% sumspec - the sum of the component spectra. will be simliar to the
% original measured spectrum, but the fit may not be as good in
% parts of the spectrum with higher uncertainty.
%
% The amplitudes of the Gaussian functions represent the relative amounts
% of different phytoplankton pigments and can be compared to HPLC pigment
% concentrations to evaluate the method. See Chase et al. (2014) for
% detail on such evaluation (reference above).
if any(min(apunc,[],2) < 0)
apunc(min(apunc,[],2) < 0, :) = ones(max(length(ap)), size(apunc(min(apunc,[],2) < 0, :), 1));
end
if size(apunc,1) < size(apunc,2)
apunc = apunc';
end
if acs == 1
fprintf('Absorption spectra correction ... ')
% Set up filter factors at every 0.1 nm from 1 to 799 nm, with center
% wavelength at centwavel (i.e. at the data wavelengths)
wavelength = .1:.1:799; % Thus index of 1 nm = 10; 356 nm= 3560;
clear filtfunc
SIG1 = (-9.845*10^-8.*wl.^3 +1.639*10^-4*wl.^2- 7.849*10^-2*wl + 25.24)/2.3547 ;
for i = 1:max(size(wl))
for jkl = 1:max(size(wavelength))
filtfunc(jkl,i) = (1/(sqrt(2*pi)*SIG1(i)))*exp(-0.5*((wavelength(jkl)-wl(i))/SIG1(i)).^2); % First term normalizes area under the curve to 1.
end
end
% Convolve the measurement with the fiter factors add the difference to
% the measured spectrum to get the first corrected spectrum.
% This is the corrected absorption spectrum "abscorr".
minwavel = min(wl);
maxwavel = max(wl);
xixi = minwavel:.1:maxwavel;% The range of centwavel is 0.1 nm.
yiyi = spline(wl, ap, xixi); % Spline the measured data to every 0.1 nm.
% We need data from 0 to 799 nm to multiply by filtfac.
absspec = zeros(size(ap,1), size(wavelength,2));
absspec(:, minwavel*10:maxwavel*10) = yiyi;
absspec(:, 1:minwavel*10-1) = ones(1, size(1:minwavel*10-1,2)) .* absspec(:, minwavel*10);
aspecprime = absspec';
meassignal6 = NaN(size(aspecprime, 2), size(wl, 2));
for i = 1:size(aspecprime, 2)
measur2 = aspecprime(:,i) .* filtfunc; % the measured signal for every filter factor.
meassignal6(i,:) = 0.1 * sum(measur2); % The measured spectrum at a wavelength i is the sum of what a filter measured at
end
abscorr = ap - meassignal6 + ap;
fprintf('Done\n')
else
fprintf('Skip absorption spectra correction\n')
abscorr = ap;
end
% Peak center values ("peak_loc") determined using a interative
% approach that allows the location to vary (uses the matlab
% function LSQNONLIN), and are rounded to nearest integer.
% Sigma values ("lsqsig") are determined similarly. FWHM = sigma*2.355
peak_loc=[406,434,453,470,492,523,550,584,617,638,660,675];
lsqsig=[16,12,12,13,16,14,14,16,13,11,11,10];
onenm = 400:1:720;
fprintf('Gaussian decomposition ... ')
%interpolate the un-smoothed ap spectra to one nm resolution
acorr2onenm = interp1(wl, abscorr', onenm, 'spline');
%define the matrix of component Gaussian functions using the peaks
%and widths (sigma) above
coef2 = exp(-0.5 .* (((onenm .* ones(size(peak_loc,2),1))' - peak_loc .* ...
ones(size(onenm,2),1)) ./ lsqsig) .^ 2);
%define a function for non-algal particles and concatenate this to the Gaussian matrix
coef2nap = exp(-0.01 * (onenm - 400));
coef2 = [coef2nap', coef2];
%normalize both the component functions and the measured ap
%spectrum by the uncertainty (std dev) in the ap spectra
apunc_int = interp1(wlunc, apunc', onenm, 'linear', 'extrap');
acorr2onenm_new = (acorr2onenm ./ apunc_int);
amps = NaN(size(ap,1), size(coef2,2));
sumspec_temp = NaN(size(coef2,1), size(ap,1));
compspec_temp = NaN(size(coef2,1), size(coef2,2), size(ap,1));
for i = 1:size(acorr2onenm_new,2)
coef2_new = coef2 ./ apunc_int(:,i);
%Inversion analysis
amps(i, :) = lsqnonneg(coef2_new, acorr2onenm_new(:, i));
% Using the inverted amplitudes, build a matrix of new component
% spectra (compspec) and the sum of the Gaussian and nap functions (sumspec):
sumspec_temp(:, i) = sum(amps(i, :) .* coef2, 2);
compspec_temp(:, :, i) = amps(i, :) .* coef2;
end
%interpolate back to the original resolution
compspec = interp1(onenm', compspec_temp, wl, 'spline');
sumspec = interp1(onenm, sumspec_temp, wl, 'spline')';
fprintf('Done\n')
end