Code covered by the BSD License  

Highlights from
Field Spectroscopy Facility Post Processing Toolbox

from Field Spectroscopy Facility Post Processing Toolbox by Iain Robinson
A toolbox for importing and processing optical spectra acquired with portable spectroradiometers.

featuredepth(spectrum, leftWavelength, rightWavelength)
%FEATUREDEPTH Plots spectrum in region of an absorption feature and estimates its depth.
%   [w, D] = featuredepth(S, w1, w2) finds the local minimum within the
%   wavelength region w1 to w2 estimates the depth of minimum relative to
%   a straight line (continuum) plotted between w1 and w2.
%
%   SEE ALSO
%   The definition of feature depth is given in Clark et al.
%
function [featureWavelength, depth] = featuredepth(spectrum, leftWavelength, rightWavelength)
    plotWavelengthRange = [leftWavelength - 200, rightWavelength + 200];
    plotReflectanceRange = [ 0, 1 ];
    
    plot(spectrum.wavelength, spectrum.data, 'Color', 'b');
    axis([plotWavelengthRange plotReflectanceRange]);
    xlabel('Wavelength in nm');
    ylabel('Reflectance');
    if isfield(spectrum, 'name')
        title(spectrum.name);
    else
        title('Specturm of unknown sample');
    end
    
    % Mark on the continuum interval points.
    leftReflectance = spline(spectrum.wavelength, spectrum.data, leftWavelength);
    rightReflectance = spline(spectrum.wavelength, spectrum.data, rightWavelength);    
    line(leftWavelength, leftReflectance, 'Marker', 's', 'Color', 'r');
    line(rightWavelength, rightReflectance, 'Marker', 's', 'Color', 'r');
    
    % Calculate the continuum (which is a straight line).
    gradient = (rightReflectance - leftReflectance) / (rightWavelength - leftWavelength);
    intercept = (leftReflectance*rightWavelength - rightReflectance*leftWavelength) / (rightWavelength - leftWavelength);
    continuum = gradient * spectrum.wavelength + intercept;
    
    % Plot the continuum.
    hold on;
    plot(spectrum.wavelength, continuum, 'Color', 'r');
    
    % Find the feature wavelength. This is the minimum value of the reflectance
    % within the continuum interval.
    % a is the start index of the continuum interval.
    % b is the end index of the continuum interval.
    intervalWavelength = spectrum.wavelength(spectrum.wavelength > leftWavelength & spectrum.wavelength < rightWavelength);
    intervalReflectance = spectrum.data(spectrum.wavelength > leftWavelength & spectrum.wavelength < rightWavelength);
    featureWavelength = intervalWavelength(intervalReflectance == min(intervalReflectance));
    featureReflectance = spectrum.data(spectrum.wavelength == featureWavelength);
    if numel(featureWavelength) > 1
        error('Feature does not have a single minimum value.');
    end
    
    % Find the continuum reflectance. This is the reflectance at the
    % feature wavelength.
    continuumReflectance = gradient*featureWavelength + intercept;
    
    % Plot a line to show the feature wavelength.
    line([featureWavelength, featureWavelength], [featureReflectance, continuumReflectance], 'Color', 'c');
    
    % Calculate the depth of the feature.
    depth = 1 - (featureReflectance / continuumReflectance);
    
    % Switch off graphics holding.
    hold off;
end

Contact us