On shifting spectra pixelwise towards lower/higher x-axis values by lsqcurvefit or such - constrain boundaries to integer values

Dear all,

in spectral data evaluation of e. g. a mixture spectrum consisting of peaks of two or more substances, it is a common approach to fit e. g. a pure spectrum of one of the two substances to the mixture spectrum, in order to deconvolute overlapping spectral ranges. Due to molecular interactions however, peaks can shift comparing a pure spectrum and a mixture spectrum, which can be easily corrected by finding maxima of certain peaks in both pure spectrum and mixture spectrum, followed by shifting one of the two. However, this approach is sometimes not accurate, since shifts can occure that lie between the pixels of the detector chip (RamanShiftVector). Therefor the spectral range / xaxis data (RamanShiftVector) can be expanded and the yaxis data / intensities (purespec, mixspec) can be linearly interpolated in order to span more datapoints between the pixels of the detector chip. When trying to solve this problem in lsq-sense however, it is problematic that the shift (pixelstep) must be constrained to be an integer value. When doing so, the solver cannot continue and remains at the initial value. Is there a workaround for that specific problem?

%% assign fit and evaluation boundaries in Raman Shifts
RegionCO2_low = 1320;
RegionCO2_high = 1400;
RegionCH_low = 1350;
RegionCH_high = 1550;
RegionCH_str_low = 2500;
RegionCH_str_high = 3100;

%% define spectral data (purespec, mixspec and x-axis data, which is the Raman Shift)
load('data.mat');
purespec = data.purespec;
RamanShiftVector = data.Raman_shift;
mixspec = data.mixspec;
%% CO2
[~, CO2_low] = min(abs(RamanShiftVector - RegionCO2_low));
[~, CO2_high] = min(abs(RamanShiftVector - RegionCO2_high));
%% comp_a / comp_b: nAlkane or nAlcohol
[~, CH_low] = min(abs(RamanShiftVector - RegionCH_low));
[~, CH_high] = min(abs(RamanShiftVector - RegionCH_high));
%% CH stretching vibration
[~, CH_str_low] = min(abs(RamanShiftVector - RegionCH_str_low));
[~, CH_str_high] = min(abs(RamanShiftVector - RegionCH_str_high));

plot_wish = 1;

%% normalization to CH stretch vibration
mixspec=CH_norming(mixspec, CH_str_low, CH_str_high);
purespec=CH_norming(purespec, CH_str_low, CH_str_high);

%% move purespec by finding the best match with mixspec within assigned fitregion
pixelstep_0 = 10;
moveFun = @(pixelstep, xData)spectral_move_for_fit(purespec,mixspec, RamanShiftVector, xData, pixelstep, CH_str_low, CH_str_high, 0);
pixelstep = lsqcurvefit(moveFun, pixelstep_0, RamanShiftVector(1, CH_str_low:CH_str_high), mixspec(1, CH_str_low:CH_str_high));

function purespec_out = spectral_move_for_fit(purespec,mixspec,RamanShiftVector, xData, pixelstep, CH_str_low, CH_str_high, plot_info)
        %% ensure that pixelstep is an integer
        pixelstep_new = round(pixelstep);

        %% define pseudo RamanShiftVectors that cover more datapoints between original pixels
        old_RamanShiftVector = 1:1:numel(RamanShiftVector);
        new_RamanShiftVector = 1:(1/5):numel(RamanShiftVector);

        %% interpolate according to pseudo RamanShiftVectors between the actual intensities (y-data) for given RamanShiftVector (xdata)
        new_purespec = interp1(old_RamanShiftVector, purespec, new_RamanShiftVector);
        new_mixspec = interp1(old_RamanShiftVector, mixspec, new_RamanShiftVector);

        %% move purespectrum by pixelstep
        moved_purespec = zeros(size(new_purespec));

    if pixelstep_new > 0 % positive step
        for i=(1+pixelstep_new):numel(new_purespec)
        moved_purespec(1,i) = new_purespec(1,i-pixelstep_new);  
        end
    elseif pixelstep_new < 0 % negative step
        for i=1:(numel(new_purespec)+pixelstep_new)
        moved_purespec(1,i) = new_purespec(1, i-pixelstep_new);
        end
    elseif pixelstep_new == 0
        moved_purespec = new_purespec;
    end
    
    %% it must be taken care that the moved pure spectrum is transferred back into the old size of RamanShiftVector (1044 pixels)
    transferred_purespec = interp1(new_RamanShiftVector, moved_purespec, old_RamanShiftVector);
    purespec_out = transferred_purespec(1, CH_str_low:CH_str_high);
    end

    function spec_norm=CH_norming(spec, CH_str_low, CH_str_high)

    %% normalization according to Peak Signal Area of CH-stretch vibration of regarded Fuel
    spec_norm=spec / sum(spec(1,CH_str_low:CH_str_high));

    end

Answer 1

[Harald]

Hi Michael,

least squares solvers are based on gradients and I would therefore not use them for integer problems unless you can treat the variable as float and round later.

The best alternatives I can think of are surrogateopt (I would try this first) and ga from the Global Optimization Toolbox. You would then minimize the norm of the residuals you currently minimize.

Best wishes,

Harald

[michael fechter]

Thank You Harald,

I am going to have a look on it!

With best, Michael

P. S.:

I was missing out that a maximum operation does not work for certain peaks, due to shifting location of maxima from one peak shoulder to another. For which I am trying the described approach.