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.

convolve(input, bands)
%CONVOLVE Convolves a spectra to a set of spectral response functions
%   C = CONVOLVE(S, B) convolves the specturm or spectra in the structure
%   array B to the set of spectral response functions ("bands") in the
%   array B.
%
%   Known issues:
%   1. Cannot be applied where data are overlapping or not sorted.
%
%   FSF Post Processing Toolbox
%   Field Spectroscopy Facility, Natural Environment Research Council

%   Author: Iain Robinson
%   Contact: fsf@nerc.ac.uk
%   Requirements: MATLAB R2009a or later
%                 the Band class from the FSF Post Processing Toolbox
%   Revision: 1.3.3
%   Date: 2011-09-04

% Convolve a spectrum to a band set.
function output = convolve(input, bands)
for n=1:numel(input)
    % Copy the input spectrum to the output spectrum. This ensures that
    % header data is preserved, as well as other fields in the
    % structure such as the spectrum's name and datetime.
    output(n) = input(n);
    
    % The wavelength scale of the input spectrum.
    inputWavelengths = input(n).wavelength;
    inputValues = input(n).data;
    
    % Check that there are no repeated wavelength values.
    if length(unique(inputWavelengths))~=length(inputWavelengths)
        error('The input spectrum''s wavelength data contains at least one repeated value:\n\t%d\nThis may indicate that there are overlapped data. The convolution cannot be applied to overlapped data.', mode(inputWavelengths))
    end
    
    % Check that input wavelengths are sorted in increasing order.
    if ~issorted(inputWavelengths)
        error('The input spectrum''s wavelength data must be sorted in increasing order of wavelength.');
    end
    
    % The function should also check that the input spectrum is defined at
    % equally spaced wavelength invervals. As yet, it does not do this.
    
    % Interpolate bands to the wavelength scale of the spectrum (if
    % required).
    bandWavelengths = bands(1).wavelength;
    % Check whether interpolation is required...
    if isequal(size(inputWavelengths), size(bandWavelengths)) && isequal(inputWavelengths, bandWavelengths)
        interpolatedBands = bands;
    else
        % Interpolate bands to the wavelength scale of the input spectra.
        % The interpolate function is defined in the Band class.
        interpolatedBands = bands.interpolate(inputWavelengths);
    end
    
    for i=1:length(bands)
        % The output wavelength is the band's centre wavelength.
        if isempty(bands(i).centre)
            % If the band's centre wavelength is not defined the outputs
            % are simply numbered. A better approach may be to estimate the
            % band's centre wavelenght from the band data.
            outputWavelengths(i) = i;
        else
            outputWavelengths(i) = bands(i).centre;
        end
        
        % The output value is the integral of the product
        % of the band with the input spectrum.
        outputValues(i) = trapz(inputWavelengths, inputValues .* interpolatedBands(i).response);
    end
    
    % Copy the input spectrum to the output spectrum, then overwrite
    % wavelength and data with the new values. (Copying the structure
    % ensures that metadata fields such as dsatetime are preserved.)
    output(n).wavelength = outputWavelengths';
    output(n).data = outputValues';
end
end

Contact us