Perform peak fitting, find peaks and label, and integrate certain areas

I have a data set where I need to do a baseline correction, and then integrate certain areas to calculate ratio between them printed as out put.

My baseline correction section seems to be working, but I do not get any peak results nor integration result. I think not finding peaks in the defined regions is the main issue.

I am sharing my code below, along with a sample data set.

% Load Raman data from a text file

data = readmatrix('path to my data.txt');

x = data(:, 1); % Wavenumbers

y = data(:, 2); % Intensity

% Baseline correction using the Asymmetric Least Squares (ALS) algorithm

baseline = als_baseline(y, 0.001, 10^6);

corrected_y = y - baseline;

% Find peaks in the spectrum- This section does not give any output

prominence_threshold = 120;

distance_threshold = 10;

peaks = findpeaks(corrected_y, 'MinPeakProminence', prominence_threshold, 'MinPeakDistance', distance_threshold);

% Perform peak fitting and integration

range1 = [2050, 2200];

range2 = [2800, 3000];

integrated_area1 = 0;

integrated_area2 = 0;

for i = 1:numel(peaks.loc)

position = x(peaks.loc(i));

intensity = corrected_y(peaks.loc(i));

if position >= range1(1) && position <= range1(2)

integrated_area1 = integrated_area1 + intensity;

end

if position >= range2(1) && position <= range2(2)

integrated_area2 = integrated_area2 + intensity;

end

end

result = integrated_area2 / ((integrated_area1 + integrated_area2) + (integrated_area1 + integrated_area2));

fprintf('Integration Result: %f\n', result);

% Define the ALS baseline correction function

function baseline = als_baseline(y, p, lam)

L = numel(y);

D = diff(speye(L), 2);

w = ones(L, 1);

for i = 1:10

W = spdiags(w, 0, L, L);

C = chol(W + lam * (D' * D));

z = C \ (C' \ (w .* y));

w = p * (y > z) + (1 - p) * (y < z);

end

baseline = z;

end

Can you please recommend a better way to do it? Or spotting of any code issues is much appreciated.

Answer 1

[Chiththaka Chaturanga]

Hi,

Thank you so much for the update.

Both areas are having a broad peak (Sample_DataSet). Peak at 2050-2200 is smaller, while 2800-3000 is much bigger.

Yes, I want to measure the total area under the peak of selected areas, then do the following calculation from the results.

Int.1/(Int1+Int2)

result = integrated_area2 / (integrated_area1 + integrated_area2); I have updated this.

To double check, I did a calculation in excel adapting trapazoid rule, and got a very different result with the Sample_DataSet..!!

I also tried to do a peak fitting using the following code.

S = load('Sample_Raman_Data2.mat');

% Perform Lorentzian peak fitting.

S.Fit = PeakFit(S.Data, 'Window', [2050, 2200], 'PeakShape', 'Lorentzian';

It gives an error as

Unrecognized function or variable 'PeakFit'.

I tried to set up corect path and install the required additional packages in the correct place. But nothing was successful unfortunately. Totally lost in this approach here.

[Mathieu NOE]

I need you to clarify a few points before I continue

do you confirm if the goal is to do a peak in both those ranges

range1 = [2050, 2200];

range2 = [2800, 3000];

and then do the integral computation for those ranges, and final result is Int.1/(Int1+Int2) ?

this is not simply picking a few samples like in the original code, it's about saying yes or no we have a bump in the ranges and do the entire integration within this range ?

[Mathieu NOE]

I think I see what you want so I continued the work on the gauss fit ...

but first thing that popped in my mind is that maybe the fit would be better if we denoise the data before the baseline correction.

if we look at what we get :

• if we do the baseline correction on the raw data
• if we do the baseline correction first and smooth afterwards
• or vice versa

it appears that the best solution is indeed to first smooth the data , and then do the baseline correction

the bottom of the curve is then closer to the y = 0 line.

% Load Raman data from a text file
% data = readmatrix('03--Spectrum--006--Spec.Data 1.txt'); 
data = readmatrix('Sample_DataSet.txt'); 
x = data(:, 1);  % Wavenumbers
y = data(:, 2);  % Intensity
% Baseline correction using the Asymmetric Least Squares (ALS) algorithm
baseline = als_baseline(y, 0.001, 10^6);
corrected_y = y - baseline;
% and denoise 
corrected_ys1 = smoothdata(corrected_y,'gaussian',20);

% do the operations in the opposite order 
% denoise the data first
ys = smoothdata(y,'gaussian',20);
% then Baseline correction using the Asymmetric Least Squares (ALS) algorithm
baseline = als_baseline(ys, 0.001, 10^6);
corrected_ys2 = ys - baseline;

figure(1),
plot(x,corrected_y,'b',x,corrected_ys1,'r',x,corrected_ys2,'g');
legend('Raw corrected data','Baseline correction then smoothing','Smoothing then Baseline correction');

% Define the ALS baseline correction function
function baseline = als_baseline(y, p, lam)
    L = numel(y);
    D = diff(speye(L), 2);
    w = ones(L, 1);
    for i = 1:10
        W = spdiags(w, 0, L, L);
        C = chol(W + lam * (D' * D));
        z = C \ (C' \ (w .* y));
        w = p * (y > z) + (1 - p) * (y < z);
    end
    baseline = z;
end

[Mathieu NOE]

based on that , we can now do a better job in term of peak fitting and integrals computation

this is done here

for the data file : 03--Spectrum--006--Spec.Data 1.txt

the result obtained is : Integration Result: 0.844016

for the data file : Sample_DataSet.txt

the result obtained is : Integration Result: 0.814477

code :

% Load Raman data from a text file
% data = readmatrix('03--Spectrum--006--Spec.Data 1.txt'); 
data = readmatrix('Sample_DataSet.txt'); 
x = data(:, 1);  % Wavenumbers
y = data(:, 2);  % Intensity
% Baseline correction using the Asymmetric Least Squares (ALS) algorithm
baseline = als_baseline(y, 0.001, 10^6);
corrected_y = y - baseline;

% it's better to denoise (smooth) the data first and then do the baseline
% correction
% denoise the data
ys = smoothdata(y,'gaussian',20);
% Baseline correction using the Asymmetric Least Squares (ALS) algorithm
baseline = als_baseline(ys, 0.001, 10^6);
corrected_ys = ys - baseline;

figure(1),
plot(x,corrected_y,'b',x,corrected_ys,'r');
legend('Raw corrected data','Smoothing then Baseline correction');


% Find peaks in the spectrum
% put as many lines (rows) as there are individual ranges
ranges =    [2000, 2400;
            2750, 3100]; 
        
for ci = 1:size(ranges,1) % loop on number of ranges
    range = ranges(ci,:);
    ind = (x>=range(1) & x<=range(2));
    xs = x(ind);
    ys = corrected_ys(ind);
    
    % peak fitting (gaussian pulse)
    yfit = my_gauss_fit(xs,ys,xs);
    
    figure(1+ci),
    plot(xs,ys,'b',xs,yfit,'r');
    legend('raw data (smoothed)','gaussian fit');
    
    % compute the integrals (1 and 2)
    integr(ci) = trapz(xs,yfit);
        
end

result = integr(2) / (integr(1) + integr(2)); 

fprintf('Integration Result: %f\n', result);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function yfit = my_gauss_fit(x,y,xfit)

% Define g(x) = a*exp(-(x-mu)^2/(2*sigma^2)):
g = @(A,X) A(1)*exp(-(X-A(2)).^2/(2*A(3)^2));

%% ---Fit Data: Analitical Strategy---
    % Cut Gaussian bell data
    ymax = max(y); xnew=[]; ynew=[];
    ind = y > 0.25*ymax;
    xnew = x(ind);
    ynew = y(ind);

    % Fitting
    ylog=log(ynew); xlog=xnew; B=polyfit(xlog,ylog,2);

    % Compute Parameters
    sigma=sqrt(-1/(2*B(1))); mu=B(2)*sigma^2; a=exp(B(3)+mu^2/(2*sigma^2)); 

    % Plot fitting curve
    A=[a,mu,sigma];
    yfit = g(A,xfit);

end

% Define the ALS baseline correction function
function baseline = als_baseline(y, p, lam)
    L = numel(y);
    D = diff(speye(L), 2);
    w = ones(L, 1);
    for i = 1:10
        W = spdiags(w, 0, L, L);
        C = chol(W + lam * (D' * D));
        z = C \ (C' \ (w .* y));
        w = p * (y > z) + (1 - p) * (y < z);
    end
    baseline = z;
end

[Chiththaka Chaturanga]

Thanks again. This is a great help for me.

Yes, two regions are corresponding to two broad peaks including the boundaries (P1= 2050-2200, P2=2800-3000) which I am interested on analysing the area underneath. Final output is the ratio as shown in your code (P1/(P1+P2)).

Interestingly the code select certain region of the bell shape after curve fitting (% Cut Gaussian bell data), isn't it? So it will be not exactly 2050-2200, but 2055- 2150 for instance? Is it possible to export these x region limits of the each curve analysed? It would be interesting for me to find this region over a number of samples, to determine the actual region of the two peaks that I am interested on.

If this region varies from one sample to another due to the peak distribution, the integration value will also differ from sample to sample, isn't it?

The reason I did not initially use function Trapz was X readings are not having the same distance between two readings. After peak fitting, is it fine to use the function trapz? May be I did not fully understand how this function works.

Thanks a lot for your support again!!!

[Mathieu NOE]

this is a further refinement of the code

now the peak amplitude , x location and width (corresponding to variables a,mu,sigma) are returned from the function my_gauss_fit

[yfit,a,mu,sigma] = my_gauss_fit(x,y,xfit)

you can then do the integral only on the mu+/-3*sigma x range that is typical for a gaussian curve

see below

your initial range is just to find the approximate peak location, then we do the gauss curve fitting and we have from there the +/- 3 sigma x limits

clc
clearvars

% Load Raman data from a text file
% data = readmatrix('03--Spectrum--006--Spec.Data 1.txt'); 
data = readmatrix('Sample_DataSet.txt'); 
x = data(:, 1);  % Wavenumbers
y = data(:, 2);  % Intensity
% Baseline correction using the Asymmetric Least Squares (ALS) algorithm
baseline = als_baseline(y, 0.001, 10^6);
corrected_y = y - baseline;

% it's better to denoise (smooth) the data first and then do the baseline
% correction
% denoise the data
ys = smoothdata(y,'gaussian',20);
% Baseline correction using the Asymmetric Least Squares (ALS) algorithm
baseline = als_baseline(ys, 0.001, 10^6);
corrected_ys = ys - baseline;

figure(1),
plot(x,corrected_y,'b',x,corrected_ys,'r');
legend('Raw corrected data','Smoothing then Baseline correction');


% Find peaks in the spectrum
% put as many lines (rows) as there are individual ranges
% this range must be wider than the typical +/- 3 sigma distribution 
% so take a broad range , the curve fitting below will reduce the
% integration range to x = mu +/- 3 sigma 
initial_range =    [1800, 2500;
                    2600, 3100]; 
        
for ci = 1:size(initial_range,1) % loop on number of ranges
    range = initial_range(ci,:);
    ind = (x>=range(1) & x<=range(2));
    xs = x(ind);
    ys = corrected_ys(ind);
    
    % peak fitting (gaussian pulse)
    [yfit,a,mu,sigma] = my_gauss_fit(xs,ys,xs);
    lower_lim = mu-3*sigma;
    upper_lim = mu+3*sigma;
    
    % refine range within +/- 3 sigma limits (peak is at x = mu ) and do the integration inside
    % these limits
    range = [lower_lim upper_lim];
    ind = (xs>=range(1) & xs<=range(2));
    xs = xs(ind);
    ys = ys(ind);
    yfit = yfit(ind);
    
    figure(1+ci),
    plot(xs,ys,'b',xs,yfit,'r');
    legend('raw data (smoothed)','gaussian fit');
    
    % compute the integrals (1 and 2)
    integr(ci) = trapz(xs,yfit);
        
end

result = integr(2) / (integr(1) + integr(2)); 

fprintf('Integration Result: %f\n', result);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yfit,a,mu,sigma] = my_gauss_fit(x,y,xfit)

% Define g(x) = a*exp(-(x-mu)^2/(2*sigma^2)):
g = @(A,X) A(1)*exp(-(X-A(2)).^2/(2*A(3)^2));

%% ---Fit Data: Analitical Strategy---
    % Cut Gaussian bell data
    ymax = max(y); xnew=[]; ynew=[];
    ind = y > 0.25*ymax;
    xnew = x(ind);
    ynew = y(ind);

    % Fitting
    ylog=log(ynew); xlog=xnew; B=polyfit(xlog,ylog,2);

    % Compute Parameters
    sigma=sqrt(-1/(2*B(1))); mu=B(2)*sigma^2; a=exp(B(3)+mu^2/(2*sigma^2)); 

    % Plot fitting curve
    A=[a,mu,sigma];
    yfit = g(A,xfit);

end

% Define the ALS baseline correction function
function baseline = als_baseline(y, p, lam)
    L = numel(y);
    D = diff(speye(L), 2);
    w = ones(L, 1);
    for i = 1:10
        W = spdiags(w, 0, L, L);
        C = chol(W + lam * (D' * D));
        z = C \ (C' \ (w .* y));
        w = p * (y > z) + (1 - p) * (y < z);
    end
    baseline = z;
end

[Mathieu NOE]

Hello

Problem solved ?

would you mind accepting my answer ? thanks !

[Chiththaka Chaturanga]

Mathieu NOE,

Thank you so much for your contrubution. Code wise, it is soved now now. Somehow peak fitting is not good, which affects my final result. I look into peak fitting, it is likely a higher order polynormial fitting, rather gaussian. So I checked literature to verify this, and most of the publications used gaussian.

Thanks again for your help. I will accept this answer, as now code is fine.

[Mathieu NOE]

ok , good to hear that the problem is solved

i was nevertheless surprised you accepted yourself indeed .. ?

Answer 2

[Mathieu NOE]

hello

fixed some minor bugs

hope it helps

% Load Raman data from a text file
data = readmatrix('03--Spectrum--006--Spec.Data 1.txt'); 
x = data(:, 1);  % Wavenumbers
y = data(:, 2);  % Intensity
% Baseline correction using the Asymmetric Least Squares (ALS) algorithm
baseline = als_baseline(y, 0.001, 10^6);
corrected_y = y - baseline;

% Find peaks in the spectrum- This section does not give any output
prominence_threshold = 10;
distance_threshold = 10;
[peaks,locs] = findpeaks(corrected_y, 'MinPeakProminence', prominence_threshold, 'MinPeakDistance', distance_threshold);
% Perform peak fitting and integration
range1 = [2050, 2200];
range2 = [2800, 3000];
integrated_area1 = 0;
integrated_area2 = 0;
for ci = 1:numel(locs)
    position = x(locs(ci));
    intensity = peaks(ci);
    if position >= range1(1) && position <= range1(2)
        integrated_area1 = integrated_area1 + intensity;
    end
    if position >= range2(1) && position <= range2(2)
        integrated_area2 = integrated_area2 + intensity;
    end
end
result = integrated_area2 / ((integrated_area1 + integrated_area2) + (integrated_area1 + integrated_area2));
fprintf('Integration Result: %f\n', result);

% Define the ALS baseline correction function
function baseline = als_baseline(y, p, lam)
    L = numel(y);
    D = diff(speye(L), 2);
    w = ones(L, 1);
    for i = 1:10
        W = spdiags(w, 0, L, L);
        C = chol(W + lam * (D' * D));
        z = C \ (C' \ (w .* y));
        w = p * (y > z) + (1 - p) * (y < z);
    end
    baseline = z;
end

[Mathieu NOE]

NB that there seems to be no peak in range1 so output is zero

[Chiththaka Chaturanga]

Thanks a lot for your help. Now it is running without any errors, but I always getting output as 0.5, even with different data sets. So I checked the same area with excel, using the same data, and it gives different values there.

I thought this might need proper peak fitting. So I am tring to do Peak fitting with Lorentzian peak fitting before the calculation.

But it gives an error as follows:

Unrecognized function or variable 'PeakFit'.

S.Fit = PeakFit(S.Data, 'Window', [100, 900], 'PeakShape', 'Lorentzian', ...

I have added the PeakFit add-on to path yet it continue to give the error.

What is the right path to add this addon to make it work?

[Mathieu NOE]

I spent some time at this code and figure out what can be the issue - and wonder if the code reflects what the intention is

now, this code with the parameters we gave to findpeaks (prominence_threshold = 10;distance_threshold = 10;) will only return those two peaks

and as the first selected peak does not fullfill the statements vs

range1 = [2050, 2200];

range2 = [2800, 3000];

we are left with only one value used in this computation

if position >= range1(1) && position <= range1(2)

integrated_area1 = integrated_area1 + intensity;

end

if position >= range2(1) && position <= range2(2)

integrated_area2 = integrated_area2 + intensity;

end

so here

integrated_area1 = 0 (no peak valid)

integrated_area2 = 22.2544 (only one single (second) peak valid)

and of course, when you do this

result = integrated_area2 / ((integrated_area1 + integrated_area2) + (integrated_area1 + integrated_area2));

you can only get a result of 0.5 whatever the value of integrated_area2

also , why did you code this denominator with twice the same expression ? this is obviously the same as :

result = integrated_area2 / (2*(integrated_area1 + integrated_area2));

[Mathieu NOE]

I changed the value of to get more peaks being selected so now the result is different from 0.5

in the code I added a few lines to display how many peaks are retained for both ranges

I am still a bit concerned that what you call area_integral is not an area integral, but a simple amplitude sum of a few peaks. Aren't you sippose to do a real integral over the range you define , assuming this range actually contains a bump ?

% Load Raman data from a text file
data = readmatrix('03--Spectrum--006--Spec.Data 1.txt'); 
x = data(:, 1);  % Wavenumbers
y = data(:, 2);  % Intensity
% Baseline correction using the Asymmetric Least Squares (ALS) algorithm
baseline = als_baseline(y, 0.001, 10^6);
corrected_y = y - baseline;
% Find peaks in the spectrum- This section does not give any output
prominence_threshold = 3;
distance_threshold = 10;
[peaks,locs] = findpeaks(corrected_y, 'MinPeakProminence', prominence_threshold, 'MinPeakDistance', distance_threshold);

plot(x,corrected_y,'b',x(locs),peaks,'dr');

% Perform peak fitting and integration
range1 = [2050, 2200];
range2 = [2800, 3000];
integrated_area1 = 0;
integrated_area2 = 0;
n1 = 0; % new code
n2 = 0; % new code

for ci = 1:numel(locs)
    position = x(locs(ci));
    intensity = peaks(ci);
    if position >= range1(1) && position <= range1(2)
        integrated_area1 = integrated_area1 + intensity;
        n1 = n1 + 1;
    end
    if position >= range2(1) && position <= range2(2)
        integrated_area2 = integrated_area2 + intensity;
        n2 = n2 + 1;
    end
end
disp(['number of selected peaks in range 1 = ' num2str(n1)]); % new code
disp(['number of selected peaks in range 2 = ' num2str(n2)]); % new code

% result = integrated_area2 / ((integrated_area1 + integrated_area2) + (integrated_area1 + integrated_area2));
result = integrated_area2 / (2*(integrated_area1 + integrated_area2)); % same thing !
fprintf('Integration Result: %f\n', result);
% Define the ALS baseline correction function
function baseline = als_baseline(y, p, lam)
    L = numel(y);
    D = diff(speye(L), 2);
    w = ones(L, 1);
    for i = 1:10
        W = spdiags(w, 0, L, L);
        C = chol(W + lam * (D' * D));
        z = C \ (C' \ (w .* y));
        w = p * (y > z) + (1 - p) * (y < z);
    end
    baseline = z;
end