Data analysis - How to recognize local peaks/different areas?
Show older comments
Hello MatLab Community,
i'm completely new to MatLab, but i think it is the right tool for my current project. I have a CNC milling machine do multiple drillings in a test piece. With MatLab i want to read my data and export just mean values for each drilling.
I am done with a simple linear compensation of temperatur drift during the process. Data looks like this:
Now i would like to get 5 average values out of this, without having to enter start and end lines/times. The goal is to get these values in a short amount of time for multiple datasets, so it can be done directly after the machining.
My first idea was to find the beginnings and endings by looking at the derivation, but i`m not sure how to implement this and if there is maybe a simpler solution like a function that could be used.
Thank you for your help!
The code looks like this:
%% Import
clear; clc;
%Tabelle2 = readtable("d12_z4_ap7_ae3_vc110_n2920_f0,12_vf1400-1.txt");
Tabelle = readtable("_N_DREHOPT_FRAESEN_D6_AP_12_MPF_PL_2023_06_29-15_00_43.txt");
name = "Fraesen_D6_AP12_2023_06_29"
%% Auswertung
Startline = 1;
Endline = 155;
Data = Tabelle(:,["Var2" "Var23"]);
Data = renamevars(Data,["Var2","Var23"],["Time","Torsion"])
Data([end-Endline:end],:) = [];
Data([1:Startline],:) = [];
%Data.Torsion = movmean(Data.Torsion,5);
%Data.Time = movmean(Data.Time,5);
%linear compensation of temperature drift
meanStart = mean(Data.Torsion(1:20));
meanEnd = mean(Data.Torsion(end-10:end));
OffsetLine = linspace(meanStart,meanEnd,length(Data.Torsion));
Offset = Data.Torsion-OffsetLine';
%% Output
%original plot
fig1 = figure(1);
axes1 = axes;
plot(Data.Time, Data.Torsion)
grid on
%plot with compensation
fig2 = figure(2)
plot(Data.Time, Offset)
grid on
saveas(fig2, name, "png")
Accepted Answer
Star Strider
on 4 Sep 2023
Edited: Star Strider
on 5 Sep 2023
It would help to have the data.
It is relatively straightforward to get the start and end times and then store the relevant data for later analysis.
Try something like this —
t = datetime('15:00:00','InputFormat','HH:mm:ss') + seconds(1:370).';
tv = second(t);
y = rectpuls(tv-50,50)*1.7+randn(size(t))*0.0051; % Signal Processing Toolbox Function
Lv = y>1.5;
pkStart = strfind(Lv.', [0 1])+1; % Start Indices
pkEnd = strfind(Lv.', [1 0]); % End Indices
for k = 1:numel(pkStart) % Store Data For Each Peak
idxrng = pkStart(k) : pkEnd(k);
pkData{k,1} = t(idxrng);
pkData{k,2} = y(idxrng);
end
for k = 1:size(pkData,1) % Process Data
pkMean(k,:) = mean(pkData{k,2});
pkDur(k,:) = second(pkData{k,1}(end)) - second(pkData{k,1}(1));
pkArea(k,:) = trapz(second(pkData{k,1}), pkData{k,2});
end
Results = table(pkMean,pkDur,pkArea)
PeakMean = pkMean;
PeakDuration = pkDur;
figure
plot(t, y)
hold on
plot(t(pkStart), y(pkStart), 'sg')
plot(t(pkEnd), y(pkEnd), 'sr')
hold off
grid
EDIT — (5 Sep 2023 at 14:55)
An alternative method to find the start and end indices (again assuming column vectors) is:
pkStart = find(diff(+[Lv(1); Lv]) > 0);
pkEnd = find(diff(+[Lv(1); Lv]) < 0);
That works just as well (although the strfind ‘hack’ can be useful for finding patterns). The ‘+[Lv(1); Lv]’ argument to diff first ‘pads’ it with the first element of ‘Lv’ (to compensate for how diff works in order to get the indices correct), then converts the logical vector ‘Lv’ to a numerical vector by preceeding it with a unary ‘+’ (uplus, +), and then finds the indices that satisfy each condition.
Added ‘pkArea’ calculation to the processing loop, since ‘areas’ was mentioned in the title.
.
7 Comments
William Rose
on 4 Sep 2023
He made a simulated data set, with 370 points, quite similar to yours.
He noticed that the level y=1.5 would be a good value for a threhold, so he used it:
Lv = y>1.5;
The line above creates a vector of logical values (1's and 0's), with length=370. It is 1 for each element of y>1.5.
He detects upward threshold crossings by finding the indices in Lv where a 0 is followed by a 1.
pkStart = strfind(Lv.', [0 1])+1;
He detects downward threshold crossings by finding the indices in Lv where a 1 is followed by a 0.
I didn't remember that you could search a vector of logical values as if it were a string of 1's and 0's. I think I have seen @Star Strider use a similar idea for other problems, and it is good to be reminded.
After that it is pretty straightforward.
There are two places in your figure where the trace goes backward in time. I don;t know how that happens, but that's what it looks like. One of those points is above the threshold. Depending on how the data is represented in the file, that could throw off the algorithm a bit, so I would check.
Star Strider
on 12 Sep 2023
Moved: Star Strider
on 12 Sep 2023
Thresholding is easier, however if you want the exact starting and ending points of the peaks, that is certainly possible.
Try this —
T1 = readtable('_N_DREHOPT_FRA...-15_00_43.txt')
T2 = readtable('_N_DREHOPT_BOH...-13_49_42.txt')
% t = datetime('15:00:00','InputFormat','HH:mm:ss') + seconds(1:370).';
% tv = second(t);
% y = rectpuls(tv-50,50)*1.7+randn(size(t))*0.0051; % Signal Processing Toolbox Function
xy{1} = {T1{:,2} T1{:,23}};
xy{2} = {T2{:,2} T2{:,23}};
for k1 = 1:numel(xy)
t = xy{k1}{:,1};
y = xy{k1}{:,2};
ymax = max(y);
ymin = y(1);
y(y<=0.9*y(1)) = y(1);
[h,m,s] = hms(t);
[pks,locs] = findpeaks(y, 'MinPeakHeight',mean([ymax ymin]));
thrshld = 25;
dps = diff([0; locs]);
pkStart = locs(dps>thrshld);
dpe = diff([locs; numel(y)]);
pkEnd = locs(dpe>thrshld);
figure
plot(t, y)
hold on
hp1 = plot(t(pkStart), y(pkStart), 'sg', 'DisplayName','Pulse Start');
hp2 = plot(t(pkEnd), y(pkEnd), 'sr', 'DisplayName','Pulse End');
hold off
grid
xlabel('Time')
ylabel('Signal')
legend([hp1(1) hp2(1)], 'Location','best')
ylim([y(1) max(y)].*[0.9 1.1])
for k = 1:numel(pkStart) % Store Data For Each Peak
idxrng = pkStart(k) : pkEnd(k);
pkData{k,1} = t(idxrng);
pkData{k,2} = y(idxrng);
end
for k = 1:size(pkData,1) % Process Data
[h2,m2,s2] = hms(pkData{k,1});
pkMean(k,:) = mean(pkData{k,2});
pkDur(k,:) = s2(end) - s2(1);
pkArea(k,:) = trapz(s2, pkData{k,2});
end
Results = table(pkMean,pkDur,pkArea)
end
EDIT — (12 Sep 2023 at 20:42)
Eliminated thresholding approach.
.
Jonas Chlechowitz
on 25 Sep 2023
Hi,
thank you again for your answer and especially your patience. Im sorry it took so long for me to try out your input. It worked perfectly until i included my drift compensation.
In line 32 "%table to cell format" i convert the table to a 1x1 cell format, since the [h,m,s] was not working otherwise. I'm not sure why we have to do this. If you could tell me, this would be great, but for now it is working :)
When just using the original table "Tabelle" for the "xy{1}" pkStart & Ends are recognized like shown in your awnser. But when i use my "Data.Time" and "Data.Torsion" for the "xy{1}", pkEnds start drifting downwards...
%% Import
clear; clc;
%Tabelle = readtable("d12_z4_ap7_ae3_vc110_n2920_f0,12_vf1400-1.txt");
Tabelle = readtable("_N_DREHOPT_BOHREN 424_MPF_PL_2023_06_29-13_49_42.txt");
name = "Fraesen D6 AP12 2023-06-29"
%% Auswertung
%identify and rename important columns
Data = Tabelle(:,["Var2" "Var23"]);
%Data = {Tabelle{:,2} Tabelle{:,23}};
Data = renamevars(Data,["Var2","Var23"],["Time","Torsion"])
%manual start & end distances
Startline = 1; %forcing manual control of start & end to avoid mistakes in temp compensation (e.g. wrong "meanEnd" caused by spindle shutdown -> negative peak)
Endline = 1;
%delete manually set start and end lines
Data([end-Endline:end],:) = [];
Data([1:Startline],:) = [];
%linear compensation of temperature drift
meanStart = mean(Data.Torsion(1:20));
meanEnd = mean(Data.Torsion(end-10:end));
OffsetLine = linspace(meanStart,meanEnd,length(Data.Torsion));
Data.Torsion = Data.Torsion-OffsetLine';
%table to cell format
xy{1} = {Data{:,"Time"} Data{:,"Torsion"}};
%xy{1} = {Tabelle{:,2} Tabelle{:,23}}; <-- working with uncorrected Data only :/
t = xy{1}{:,1};
y = xy{1}{:,2};
ymax = max(y);
ymin = y(1) %mean(Data.Torsion(1:20));
y(y<=0.9*ymin) = y(1); % = ymin;
[h,m,s] = hms(t);
[pks,locs] = findpeaks(y, 'MinPeakHeight',mean([ymax ymin]));
thrshld = 25;
dps = diff([0; locs]);
pkStart = locs(dps>thrshld);
dpe = diff([locs; numel(y)]);
pkEnd = locs(dpe>thrshld);
%Store Data for each peak
for k = 1:numel(pkStart)
idxrng = pkStart(k) : pkEnd(k);
pkData{k,1} = t(idxrng);
pkData{k,2} = y(idxrng);
end
%Calculate duration, mean & area
for k = 1:size(pkData,1)
[h2,m2,s2] = hms(pkData{k,1});
pkMean(k,:) = mean(pkData{k,2});
pkDur(k,:) = s2(end) - s2(1);
pkArea(k,:) = trapz(s2, pkData{k,2});
end
%% Output
%plot
figure
plot(t, y)
hold on
hp1 = plot(t(pkStart), y(pkStart), 'sg', 'DisplayName','Pulse Start');
hp2 = plot(t(pkEnd), y(pkEnd), 'sr', 'DisplayName','Pulse End');
hold off
grid
xlabel('Time')
ylabel('Signal')
legend([hp1(1) hp2(1)], 'Location','best')
ylim([y(1)-1 max(y)].*[0.9 1.1])
%duration, mean & area
Results = table(pkMean,pkDur,pkArea)
After about two weeks, I had to go back and remember what Iw as doing here.
We are back to a version of thresholding again with those changes, because there is simply no way around it with these data. (Thresholding augments the peak finding. I am using it in conjunction with the ismember function to more closely define the peaks, especially at the ends.)
I also experimented with findchangepts and ischange. They came up with the same results my original code (revised here) does.
This is the best I can do with these data.
Files = dir('*.txt');
Files.name
for kf = 1:numel(Files)
Tabelle = readtable(Files(kf).name)
%Tabelle = readtable("d12_z4_ap7_ae3_vc110_n2920_f0,12_vf1400-1.txt");
% Tabelle = readtable("_N_DREHOPT_BOHREN 424_MPF_PL_2023_06_29-13_49_42.txt");
% name = "Fraesen D6 AP12 2023-06-29"
%% Auswertung
%identify and rename important columns
Data = Tabelle(:,["Var2" "Var23"]);
%Data = {Tabelle{:,2} Tabelle{:,23}};
Data = renamevars(Data,["Var2","Var23"],["Time","Torsion"])
%manual start & end distances
Startline = 1; %forcing manual control of start & end to avoid mistakes in temp compensation (e.g. wrong "meanEnd" caused by spindle shutdown -> negative peak)
Endline = 1;
%delete manually set start and end lines
Data([end-Endline:end],:) = [];
Data([1:Startline],:) = [];
%linear compensation of temperature drift
meanStart = mean(Data.Torsion(1:20));
meanEnd = mean(Data.Torsion(end-10:end));
OffsetLine = linspace(meanStart,meanEnd,length(Data.Torsion));
Data.Torsion = Data.Torsion-OffsetLine';
%table to cell format
xy{1} = {Data{:,"Time"} Data{:,"Torsion"}};
%xy{1} = {Tabelle{:,2} Tabelle{:,23}}; <-- working with uncorrected Data only :/
t = xy{1}{:,1};
y = xy{1}{:,2};
ymax = max(y);
ymin = y(1); %mean(Data.Torsion(1:20));
y(y<=0.9*ymin) = y(1); % = ymin;
[h,m,s] = hms(t);
[pks,locs] = findpeaks(y, 'MinPeakHeight',mean([ymax ymin]));
% LocsStats = [mean(diff(locs)); std(diff(locs)); median(diff(locs)); min(diff(locs)); max(diff(locs))]
% PksStats = [mean(pks); std(pks); median(pks); std(pks)/sqrt(numel(pks)); 1.96*std(pks)]
thrshld = 20;
puls = y > mean(pks)-1.96*std(pks);
Lvlocs = ismember(locs, find(puls)); % Threshold The Peak Definitions (Logical Vector)
locs = locs(Lvlocs);
tx = 1:numel(t);
%
% figure
% plot(locs, pks, '.-')
% hold on
% plot(tx, +puls*max(pks))
% hold off
% grid
% title("kf = "+string(kf))
dps = diff([0; locs]);
pkStart = locs(dps>thrshld);
dpe = diff([locs; numel(y)]);
pkEnd = locs(dpe>thrshld);
%Store Data for each peak
for k = 1:numel(pkStart)
idxrng = pkStart(k) : pkEnd(k);
pkData{k,1} = t(idxrng);
pkData{k,2} = y(idxrng);
end
%Calculate duration, mean & area
for k = 1:size(pkData,1)
[h2,m2,s2] = hms(pkData{k,1});
pkMean(k,:) = mean(pkData{k,2});
pkDur(k,:) = s2(end) - s2(1);
pkArea(k,:) = trapz(s2, pkData{k,2});
end
%% Output
%plot
figure
plot(t, y)
hold on
hp1 = plot(t(pkStart), y(pkStart), 'sg', 'DisplayName','Pulse Start');
hp2 = plot(t(pkEnd), y(pkEnd), 'sr', 'DisplayName','Pulse End');
% plot(t, +puls*max(pks)) % Thresholding Diagnostic
hold off
grid
xlabel('Time')
ylabel('Signal')
legend([hp1(1) hp2(1)], 'Location','best')
ylim([y(1)-1 max(y)].*[0.9 1.1])
%duration, mean & area
Results = table(pkMean,pkDur,pkArea)
end
.
Jonas Chlechowitz
on 26 Sep 2023
Star Strider
on 26 Sep 2023
As always, my pleasure!
The correction essentially appeared to be a baseline shift (I experimented with different approaches — these data are not easy to work with), and that could distort the peaks themselves by changing their locations with respect to criteria used to determine the peaks and peak locations. Intuitively, it would not seem that a simple subtraction, such as the baseline correction, would do that, however it seems to. I have no other explanation for it.
I thought a bit more about this, especially the detrending step, and came up with some improvements.
This is robust to these two data sets, so I hope it works as well with others.
The newly revised code —
Files = dir('*.txt');
Files.name
for kf = 1:numel(Files)
Tabelle = readtable(Files(kf).name)
%Tabelle = readtable("d12_z4_ap7_ae3_vc110_n2920_f0,12_vf1400-1.txt");
% Tabelle = readtable("_N_DREHOPT_BOHREN 424_MPF_PL_2023_06_29-13_49_42.txt");
% name = "Fraesen D6 AP12 2023-06-29"
%% Auswertung
%identify and rename important columns
Data = Tabelle(:,["Var2" "Var23"]);
%Data = {Tabelle{:,2} Tabelle{:,23}};
Data = renamevars(Data,["Var2","Var23"],["Time","Torsion"])
%manual start & end distances
Startline = 1; %forcing manual control of start & end to avoid mistakes in temp compensation (e.g. wrong "meanEnd" caused by spindle shutdown -> negative peak)
Endline = 1;
%delete manually set start and end lines
Data([end-Endline:end],:) = [];
Data([1:Startline],:) = [];
%linear compensation of temperature drift
meanStart = mean(Data.Torsion(1:20));
meanEnd = mean(Data.Torsion(end-10:end));
OffsetLine = linspace(meanStart,meanEnd,length(Data.Torsion));
Data.Torsion = Data.Torsion-OffsetLine';
%table to cell format
xy{1} = {Data{:,"Time"} Data{:,"Torsion"}};
%xy{1} = {Tabelle{:,2} Tabelle{:,23}}; <-- working with uncorrected Data only :/
t = xy{1}{:,1};
y = xy{1}{:,2};
[h,m,s] = hms(t);
tv = s+cumsum(diff([m(1);m])>0)*60; % Cumulative Seconds (Some Functions Do Not Work With 'datetime' Or 'duration' Arrays)
% figure
% plot(t, tv)
% grid
y = max([zeros(size(y)) y],[],2); % Delete Values < 0
Lv1 = islocalmin(y, 'MinProminence',0.005); % Used In The Detrending Polynomial
Lv1 = Lv1 & y<1; % Thresholding
% p = polyfit(tv(Lv1), y(Lv1), 3);
% blc = polyval(p, tv)
% figure
% plot(t, y)
% hold on
% plot(t(Lv1), y(Lv1), '.r')
% hold off
% grid
%
p = polyfit(tv(Lv1), y(Lv1), 3); % Calculate Detrending Polynomial
blc = polyval(p, tv); % Evaluate Detrending Polynomial
y = y-blc; % Detrend 'y'
Lv2 = islocalmax(y, 'MinProminence',0.01); % Determine Peaks
Lv2 = Lv2 & y>1; % Threshold Peaks
% figure
% plot(t, y)
% hold on
% plot(t(Lv2), y(Lv2), '.r')
% hold off
% grid
% xline(Lv2)
% return
ymax = max(y);
ymin = y(1); %mean(Data.Torsion(1:20));
y(y<=0.9*ymin) = y(1); % = ymin;
[h,m,s] = hms(t);
pks = y(Lv2);
locs = find(Lv2);
% [pks,locs] = findpeaks(y, 'MinPeakHeight',mean([ymax ymin]));
% LocsStats = [mean(diff(locs)); std(diff(locs)); median(diff(locs)); min(diff(locs)); max(diff(locs))]
% PksStats = [mean(pks); std(pks); median(pks); std(pks)/sqrt(numel(pks)); 1.96*std(pks)]
thrshld = 20;
puls = y > mean(pks)-1.96*std(pks);
Lvlocs = ismember(locs, find(puls)); % Threshold The Peak Definitions (Logical Vector)
locs = locs(Lvlocs);
% tx = 1:numel(t);
% %
% % figure
% % plot(locs, pks, '.-')
% % hold on
% % plot(tx, +puls*max(pks))
% % hold off
% % grid
% % title("kf = "+string(kf))
dps = diff([0; locs]);
pkStart = locs(dps>thrshld);
dpe = diff([locs; numel(y)]);
pkEnd = locs(dpe>thrshld);
%Store Data for each peak
for k = 1:numel(pkStart)
idxrng = pkStart(k) : pkEnd(k);
pkData{k,1} = t(idxrng);
pkData{k,2} = y(idxrng);
end
%Calculate duration, mean & area
for k = 1:size(pkData,1)
[h2,m2,s2] = hms(pkData{k,1});
pkMean(k,:) = mean(pkData{k,2});
pkDur(k,:) = s2(end) - s2(1);
pkArea(k,:) = trapz(s2, pkData{k,2});
end
%% Output
%plot
figure
plot(t, y)
hold on
hp1 = plot(t(pkStart), y(pkStart), 'sg', 'DisplayName','Pulse Start');
hp2 = plot(t(pkEnd), y(pkEnd), 'sr', 'DisplayName','Pulse End');
% plot(t, +puls*max(pks)) % Thresholding Diagnostic
hold off
grid
xlabel('Time')
ylabel('Signal')
legend([hp1(1) hp2(1)], 'Location','best')
ylim([y(1)-1 max(y)].*[0.9 1.1])
%duration, mean & area
Results = table(pkMean,pkDur,pkArea)
end
.
More Answers (3)
William Rose
on 12 Sep 2023
Edited: William Rose
on 12 Sep 2023
[edit: I attached each data file twice, by mistake. Removing the extra files.]
You should fix the points that have incorrect time stamps, i.e. time stamps that are smaller than the preceding point, because they can result in incorrect peak detections.
T1=readtable('Jonas1.txt');
t1=T1.(2); x1=T1.(23);
% plot the time step size versus time
figure; subplot(211); plot(t1(2:end),diff(t1),'-r.')
You can see the bad time steps in the red plot below. Fix as follows.
t1fix=t1;
for i=2:length(t1)-1
if t1(i)<t1(i-1)
t1fix(i)=(t1(i-1)+t1(i+1))/2;
end
end
Plot the time step size versus time, for the fixed t1 (blue plot below):
subplot(212); plot(t1fix(2:end),diff(t1fix),'-b.')
That looks better. The two unusually large positive time steps in the blue plot above should not be a problem.
Likewise for the Jonas2 data:
T2=readtable('Jonas2.txt');
t2=T2.(2); x2=T2.(23);
% plot the time step size versus time
figure; subplot(211); plot(t2(2:end),diff(t2),'-r.')
t2fix=t2;
for i=2:length(t2)-1
if t2(i)<t2(i-1)
t2fix(i)=(t2(i-1)+t2(i+1))/2;
end
end
subplot(212); plot(t2fix(2:end),diff(t2fix),'-b.')
The plots above demonstrate that we have fixed the time vectors.
You said you don't like the threshold approach which @Star Strider proposed, because it can find the start points too early, and it can find the end points too late. I like the method of @Star Strider. You could define the start as one point after the upward threshold crossing, and define the end as one point before the downward threshold crossing, which would adress your concern, at least to some extent.
The attached script fixes the time vector. The script estimates the mean value of each peak, using @Star Strider's peak-finding method, but with a 1 point delay at the start and a 1-point advance at the end of each peak. The script plots the data and the fitted peak values. The results for files Jonas1 and Jonas2 are below. The code finds the peak locaitons and peak mean values quite well.
2 Comments
Jonas Chlechowitz
on 25 Sep 2023
William Rose
on 25 Sep 2023
Regarding the time stamp solution, you are welcome.
You wrote: "I will try to stick with the approach of not having a threshold, because the frequency of mesaurements may depend on the machine." This implies that a threshold-based approach is undesirable when the frequency of measurement varies. I do not understand why that is the case. I am also not sure what method you intend to use, since the method of @Star Strider, and my similar script, both use thresholds.
You said in your original post: "i would like to get 5 average values out of this, without having to enter start and end lines/times. The goal is to get these values in a short amount of time for multiple datasets, so it can be done directly after the machining". My script and the script of @Star Strider do provide the mean values in the variable pkMean.
My script uses hard-coded thresholds of 1.7 and 15 for the two datasets. The script of @Star Strider uses a more sophisticated automatic method. Here's what his code does, I think: Find all peaks in the original signal that are above the midpoint between the min and the max. Then find the mean and SD of those peaks, and exclude any "peaks" that are >1.96 SDs below the mean. With this new set of peaks, identify "start" points as peaks that occur more than thrshld=20 samples after the immediately preceding peak. Identify "end" points as peaks that occur more than 20 samples before the next peak. This automatic method gives nice results for the two files.
Jonas Chlechowitz
on 12 Sep 2023
Hello again,
thank you very much for your help. I was not able to continue with this until now. But my employer agreed, that it worth investing some time, so i am able to start daily work on this next week :) ...I tried the first method by including:
Lv = Offset>1.7;
pkStart = strfind(Lv.', [0 1])+1;
pkEnd = strfind(Lv.', [1 0]);
for k = 1:numel(pkStart) % Store Data For Each Peak
idxrng = pkStart(k) : pkEnd(k);
pkData{k,1} = Data.Time(idxrng);
pkData{k,2} = Offset(idxrng);
end
for k = 1:size(pkData,1) % Process Data
pkMean(k,:) = mean(pkData{k,2});
end
pkMean
I have attached this exact data and another dataset to this comment. It represents a lot of channels the machine exports, the only interesting ones being Var2 and Var23, when imported as a table. I dont think that a threshold is ideal for this. In the following picture we can see, that it detects pkStart & End a little too early/late. Besides other datasets have different pkMeans and maybe dont start just as sharp.(2nd dataset) Sorry for not including the data right away!
The values going back in time are a problem of the machine. If they are a problem for finding start & end indices, i would just delete all rows going back in time.
You can use the ischange function to identify where the mean of the data changes. Using the same approach for generating sample data that Star Strider's accepted answer does:
t = datetime('15:00:00','InputFormat','HH:mm:ss') + seconds(1:370).';
tv = second(t);
y = rectpuls(tv-50,50)*1.7+randn(size(t))*0.0051;
plot(t, y, '-');
the ischange function will indicate the points where the mean of the data changes (using the default change indication method.)
c = ischange(y, 'SamplePoints', t);
figure
plot(t, y, '-o', 'MarkerFaceColor', 'r', 'MarkerIndices', find(c))
Note that ischange detects the start of the peak and the start of the valley, so if you want the two edges of each peak you'd need to determine if the values of the selected points are high or low and "back up" one step if they're low (to catch the top of the "cliff" rather than the bottom.)
functionValues = y(c) % Note that these alternate high / low / high / low
lowerPoints = functionValues < mean(functionValues); % Split the high points from the low points
loc = find(c); % All the change points
loc(lowerPoints) = loc(lowerPoints)-1; % Move the lower ones to the top of the cliff
figure
plot(t, y, 'o-', 'MarkerFaceColor', 'r', 'MarkerIndices', loc)
Categories
Find more on Time-Frequency Analysis in Help Center and File Exchange
on 4 Sep 2023
on 28 Sep 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
