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