Area between baseline and data peak

Hi,

I am trying to determine the area under a peak curve, see image. The baseline illustrated in image is not zero, and I would like to subtract this from the area under the peak. The data is from a picture and been summed to obtain this profile. I want to determine the intensity by summing the area under the peak.

1. Is there a way to determine the baseline (y- value) automatically by code.

2. How do I determine the area minus the baseline? I use sum to determine the area under the peak curve.

Answer 1

[Star Strider]

Here is how I would do it:

% Create Data
x = 0:480;
y = 1E4*(0.5+rand(size(x)));
y(400:450) = 1.3E+5*exp(-(x(400:450)-425).^2/50)+y(400:450);

% Generate Statistics & Identify Peak
ysts = [mean(y) median(y);  mean(y)-1.96*std(y)/sqrt(length(y))   mean(y)+1.96*std(y)/sqrt(length(y))];
[ypk, pki] = max(y);
pkiv = find(y >= ysts(2,2));                % Find Peak Indices

Iy = cumtrapz(x, y-ysts(1,1));              % Integrate Entire Record
Iydb = [ones(pkiv(1),1)  x(1:pkiv(1))']\Iy(1:pkiv(1))';
Iyd = Iy' - [ones(size(x))'  x']*Iydb;      % Detrend Using Baseline

Ipk = Iy(max(pkiv)) - Iy(min(pkiv));        % Find Peak Area

% Plot Data & Detrended Integral & Selected Statistics
figure(1)
plot(x, y)
hold on
plot(x, Iyd, '-g')
hold off
legend('Data', 'Detrended Integral', 'Location', 'NW')
text(100, 1E+6, sprintf('Peak Area = %23.15E', Ipk))

The idea is reasonably straightforward:

1. Identify the peak as y-values greater than the upper 95% confidence interval for the mean of all the data;
2. Integrate the entire record using cumtrapz;
3. Do a linear regression on all the integrated data up to the first index of the peak;
4. Use that regression to detrend the integrated data;
5. Use the identified indices of the peak to calculate the approximate integral of the peak.

This approach may not generalise well if most of the baseline information is not at x-values less than the peak, because it depends on that information to detrend the integral.

The plot produces:

[Muthu Annamalai]

I was mentally thinking do a cumsum and then diff the result to find the peak. Nice solution there!

[Star Strider]

Thank you!

[Lizan]

Thank you for your reply. =)

I don't understand the "detrended integral" line. What is this line showing?

I am confused about what this code does:

Iydb = [ones(pkiv(1),1)  x(1:pkiv(1))']\Iy(1:pkiv(1))';
Iyd = Iy' - [ones(size(x))'  x']*Iydb;      % Detrend Using Baseline

Also, what is the reason for using ystat(1,1) for substring from y for the integration?

Iy = cumtrapz(x, y-ysts(1,1));              % Integrate Entire Record

[Star Strider]

Those two lines calculate the trend line ( ‘Iydb’ ) and detrend ( ‘Iyd’ ) the integral ‘Iy’. (The integral of a constant creates a linear trend. Subtracting that trend straightens out the integral line. It may not strictly be necessary, but it improves the accuracy of the peak integral calculation.)

I calculated a few statistics in ‘ysts’ and then referred to them rather than recalculating them. The value of ‘ysts(1,1)’ is the mean. I subtracted the mean from the signal before integrating it because it essentially zeros the baseline, improving the accuracy of the integral and peak area calculation. This does not affect the integral of the peak, since that is measured from the chosen threshold. (The other values in ‘ysts’ are the median and the upper and lower 95% confidence intervals on the mean. I later used the upper limit as the threshold to define the peak. I experimented with the other values in zeroing the initial baseline and defining the integral before deciding on the best option, that I posted.)

[Lizan]

I see, thank you!

[Star Strider]

My pleasure!

[Image Analyst]

When you do \Iy you're doing matrix division. Did you really want element-by-element division? If so, use dot slash instead of slash.

Also ones() is returning a column vector, so make sure x is a row vector, since you transpose x to be a column vector and want to stitch it to the right of ones.

Also be away of automatic expansion of your ly vector since you're dividing a pkiv(1)-by-2 vector by a vector array.

Answer 2

[Andrew Reibold]

The area under a curve is the definition of an integral. You will want to take the integral and subtract the baseline * baseline_width to find the area.

I'm not sure if its really scientific to determine the baseline just from raw data. Maybe this entire unique data set has some factor that shifts every point in some direction which makes you perceive the baseline somewhere that it isn't. I'm sure you could take the mean of non-peaked data or something, but normally the baseline is determined from the original function, historical data, or whatever you are stepping off from.

[Lizan]

Can you determinate the integral from a set of data points (x,y) in matlab? I don't have a function of this curve.

Also how do you establish the baseline? Mean-value? Can you do that without the peak?

[Image Analyst]

Lizan, you could speed along a solution to your problem by just posting your data.

[Lizan]

Some data examples:

Answer 3

[Image Analyst]

You can find the peak (at least for this plot) very easily with max(). Then you can just "fall down" the peak on each side (with a for loop) and find the place (index) where the data is no longer decreasing. Then you can just fit a line between those two valley locations, say with polyfit() or whatever, to get the baseline.

[Lizan]

How would I find out that the curve is no longer decreasing?

[Image Analyst]

signal(i) will be greater than or equal to signal(i-1).