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

Question

Chiththaka Chaturanga 2023-6-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1982739-perform-peak-fitting-find-peaks-and-label-and-integrate-certain-areas

评论： Mathieu NOE 2023-6-28

03--Spectrum--006--Spec.Data 1.txt

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

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Chiththaka Chaturanga 2023-6-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1982739-perform-peak-fitting-find-peaks-and-label-and-integrate-certain-areas#answer_1257399

Sample_DataSet.txt

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Mathieu NOE 2023-6-16

在 MATLAB Online 中打开

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 2023-6-16

编辑：Mathieu NOE 2023-6-16

在 MATLAB Online 中打开

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 2023-6-16

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 2023-6-20

在 MATLAB Online 中打开

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

Chiththaka Chaturanga 2023-6-28

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 2023-6-28

ok , good to hear that the problem is solved

i was nevertheless surprised you accepted yourself indeed .. ?

请先登录，再进行评论。

Answer 2

Mathieu NOE 2023-6-14

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1982739-perform-peak-fitting-find-peaks-and-label-and-integrate-certain-areas#answer_1255654

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mathieu NOE 2023-6-15

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 2023-6-15

在 MATLAB Online 中打开

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

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论