- Woops, in case anyone is using my attached files: the CSV-file that the .m file is loading is not the one I put up. It shouldn't matter though.
Finding area under peaks for drifting signal
6 views (last 30 days)
Show older comments
Christian Hjuler Christensen
on 23 Jan 2018
Commented: Star Strider
on 25 Jan 2018
EDITED, see bottom, (regarding numerical integration)
Dear all,
How can I find the area under the "larger" peaks, as seen below?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/169345/image.jpeg)
My idea is to:
1) locate where the peaks are (I have tried using the findpeaks() function)
2) add or subtract the minimum value of the particular peak to get a baseline for each peak at y = 0.
3) use some sort of trapezoidal integration (I saw that MATLAB has a built-in function that does this).
When using the findpeaks() function, I'm having a really hard time locating "interesting" peaks. For instance, the first "big" peak has some noise in it and then findpeaks() give me a lot of peaks on it - and then there are other small peaks I'm not interested in.
I hope someone can help!
(If you need it, my data is found in the attached CSV file - I'm only using the first two columns, "Time" and "Channel1V". I have 100 CSV files like the uploaded one, so I really need a smart code ;) )
EDIT:
Q: How to find area under peak using trapz()?
So, using the arguments 'MinPeakDistance' and 'MinPeakProminence' I have been able to figure out which peaks to use. I still need to find the area under the graphs; so far, I've been using the very rough estimate Prominence times Width (where width is the width of the peak taken at the half of the prominence). It should for most purposes be alright, but I'd really appreaciate if anyone has got a good idea of how to numerically integrate this.
11 Comments
Star Strider
on 25 Jan 2018
My pleasure.
I can help with the signal processing. You posted one of your signals, so I will download it, and experiment with it.
Accepted Answer
Star Strider
on 25 Jan 2018
I would use the Signal Processing Toolbox medfilt1() function to estimnate the base line drift, that you can then remove. (None of my other approaches, such as frequency-selective filtering or the Savitzky-Golay filter, gave satisfactory results.)
Try this to remove the baseline drift:
[D,S,R] = xlsread('420B1.1.csv');
T = D(:,1);
Ch1 = D(:,2);
Ch2 = D(:,3);
Ts = mean(diff(T));
L = length(T);
MF = medfilt1(Ch1, 750); % Median Filter
Ch1F = Ch1-MF; % Correct Baseline
figure
plot(T,Ch1F)
grid
I leave the peak detection to you. Here is one possibility:
dCh1F = gradient(Ch1F, Ts); % Derivative
[pkp,locp] = findpeaks(dCh1F, 'MinPeakHeight',0.04, 'MinPeakDistance',500); % Positive Derivative Peaks, ‘locp’ Are Indices
[pkn,locn] = findpeaks(-dCh1F, 'MinPeakHeight',0.03, 'MinPeakDistance',500); % Negative Derivative Peaks, ‘locn’ Are Indices
figure
plot(T,Ch1F)
hold on
plot(T([locp; locn]),Ch1F([locp; locn]), '+r')
hold off
grid
You will have to experiment to get the result you want. This may not be robust to all your data, so you may have to adjust it for each file.
2 Comments
More Answers (1)
AJ Geiger
on 24 Jan 2018
1. instead of using the function maxpeak use
[~,I] = max(f);
to obtain the index of the max peak
2. Nice! Now you know where the max peak is. The next step is to set up your boundaries, so you know where to start and stop your integration. There are several algorithms you can implement. the easiest would be simulated annealing. The algorithm goes as follows.
% Note: E is your data vector
% for the alg I ~= 1
if I == 1 % check condition
I = 2;
end
%
%
% Preform Alg
for t = 0:length(E)-1
% check if next point is
idx = I+t;
Delta_E = E(I+t) - E(I-1+t);
if E > 0, continue, end
%
% Note T: is your cooling function which
% Converges to 0 as t increases. now this
% is basic but T could be as simple as
% T = 1/t or T = 1/(t^2)
% now for tuning part to be complete you
% should say:
% if T == 0, break, end
% however only use it if you develop T
% properly.
%
T = 1/(t^2); % Temperature Function
if exp(Delta_E/T) < Thresh, break, end
end
Now you have found the Right-hand side boundary. For now, you can
3. To be complete, idx_2 = idx;
4. use the function E_New = flip(E);
5. solve for I_New = I - length(E) + 1;
6. then run alg two on E_new and I_New
7. To be complete, idx_1 = idx;
8. get index vector II you wish to integrate over
II = idx_1:idx_2
9. Then solve for the area using II
Area = trapz(E(II));
10. Grab a beer if of legal age and does not effect believes (Optional)
11. Finish beer
12. Explore other integration methods.
2 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!