Extracting Zoomed in part of curve and finding rising Edges of data
Show older comments
Hello, I have a plot consisting of "square-ish" TTL pulses and I want to be able to
1: Extract the data thats in the "Zoomed" XLim, and then
2: Find the x-location of each rising edge at a certain y-value (say 3, although I may not have data thats exactly on 3, so the nearest to it)
3: Bin the seperations of each rising edge. As you can see below from the 2nd (zoomed in) plot I am having some drop outs so I want to analyse e.g. if there is any variabity in the drop out seperation as well as the actual rising edge seperation. My actual data consists of thousands of these peaks
(I dont have the signal processing toolbox!)
Heres a sample of my real data
and here's a zoomed in version (this would be what I would want to work with and I obtain it by setting the XLims
So for 1: above, I have this.
ax=app.UIAxes;
newlimits =ax.XLim;
[x,y] = getDataFromGraph(app,ax,1); %See function below
xminidx = find(x >= newlimits(1),1,"first");
xmaxidx = find(x <= newlimits(2),1,"last");
X=x(xminidx:xmaxidx,1); Y=y(xminidx:xmaxidx,1);
D=[X,Y];
head(D,8) %View data
where
function [datax,datay] = getDataFromGraph(app,ax,n)
% ax = axes plot is on
% n is the plot index, 1 is the last plot
b=ax.Children; % b(1) is last object added to ax
datax=(b(n).XData)';
datay=(b(n).YData)';
end
Then for 2, this is where I am struggling, I have tried this:
rising_edges = find(diff(X > 3) > 0)
hold(ax,"on");
xline(ax,X(rising_edges))
But rising_edges is empty.
And for 3 I have no idea how to plot some kind of histogram of the seperations between adjacent pulse
8 Comments
Voss
on 4 Feb 2025
Is rising_edges still empty if you do this?
rising_edges = find(diff(Y > 3) > 0) % use Y not X here
hold(ax,"on");
xline(ax,X(rising_edges))
Jason
on 4 Feb 2025
Thanks Vos, Rising edges now seems correct!
So I now have this:
%Find Rising edges
rising_edges = find(diff(Y > 3) > 0)
hold(ax,"on");
%Plot
plot(ax,X(rising_edges),Y(rising_edges),'r*');
%Calc adjacent pulse seperation
Xsep=X(rising_edges);
format shortG
diff(Xsep) % I think this is the adjacent rising edge spacing?
But why when the threshold is 3 are the r* at lower valuesnd not at 3?
"But why when the threshold is 3 are the r* at lower valuesnd not at 3?"
rising_edges is the indices where Y is about to cross above 3 (i.e., the indices ii where Y(ii) <= 3 and Y(ii+1) > 3). Therefore, Y(rising_edges) is an array of numbers that are <= 3.
A simple example with one rising edge:
Y = linspace(0.5,10,7)
rising_edges = find(diff(Y > 3) > 0)
Y(rising_edges)
If you want to plot the markers with a y-coordinate of exactly 3, then:
plot(ax,X(rising_edges),3*ones(numel(rising_edges),1),'r*');
Jason
on 4 Feb 2025
Voss
on 4 Feb 2025
"Would i need to interp"
You could do that.
"if i wanted the index to instead be the first value >3"
Add 1 to the current definition:
rising_edges = find(diff(Y > 3) > 0) + 1
Steven Lord
on 12 Feb 2025
Just to add one more possible tool to this discussion, you may want to load your data into a live script and use the Find Local Extrema task to interactively experiment with the various options that islocalmax and islocalmin support like min prominence, min separation, and max number of extrema. This will also allow you to see the code that programmatically performs that extrema finding, so you could use the code directly on other data sets after you've interactively tuned the options on one data set.
Jason
on 12 Feb 2025
Steven Lord
on 12 Feb 2025
The "Open the Task" section on the documentation page to which I linked describes two ways to add the task to a live script.
Accepted Answer
T1 = readtable('DS0007_REDUCED.CSV')
t = T1.Col1;
s = T1.Col2;
% figure
% plot(t, s)
% grid
Lvs = islocalmax(s, 'flatselection', 'first', MinProminence=0.5);
Lve = islocalmin(s, 'flatselection', 'last', MinProminence=0.5);
figure
plot(t, s)
hold on
plot(t(Lvs), s(Lvs), '|r')
plot(t(Lve), s(Lve), '|m')
hold off
grid
figure
plot(t, s)
hold on
plot(t(Lvs), s(Lvs), '|r')
plot(t(Lve), s(Lve), '|m')
hold off
grid
xlim([0.011 0.025])
ylim('padded')
Lvt = t >= 0.01;
Ivs = find(Lvs & Lvt);
Ive = find(Lve & Lvt);
for k = 1:numel(Ivs)
idxrng = Ivs(k) - 10 : Ive(k) + 10;
segments{k} = [t(idxrng) s(idxrng)];
end
figure
tiledlayout(4,2)
for k = 1:numel(segments)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
.
29 Comments
Jason
on 5 Feb 2025
Jason
on 5 Feb 2025
Hi, sorry one quick follow on question. Occasionally I get a double count for a single square pulse - how can I ignore the 2nd one
This 2nd one messes up the pulse seperation as the ditance to the next pulse is considered from the 2nd count on this pulse - so I need to get rid of it
Thanks
Star Strider
on 5 Feb 2025
As always, my pleasure!
With respect to the double counts, experiment with the MinProminence value. You may have to nudge it up a bit. I would start by incrementing it from 0.50 to 0.55 by increments of 0.05 or less, until the double counts disappear. That should not significantly affect the pulses without the double counts.
Depending on what you want to do (and if you have the Signal Processing Toolbox), the risetime and falltime functions may also be useful.
.
Thanks. I've noticed that when there are more than 8 pulses the line with idxrng errors and I don't understand why.
for k = 1:numel(Ivs)
idxrng = Ivs(k) - 10 : Ive(k) + 10;
segments{k} = [x(idxrng) y(idxrng)];
end
Is it possible to handle the case when there aremany more pulse and for example just take the 1st 8? (Thats a fantastic graphic btw!!)
BTW, I did try changing the prominence parameter but it didn't have any effect. I also have some occurances like this and was wondering if I should just remove any of the max locations that are within a certain distance 'x1' of another max or a certain distance 'x2' of a min
I crafted the loop and the following tiledlayout plot series simply to show the individual pulses in detail, since the original plot doesn’t permiit that.
To make ‘idxrng’ more robust, change it to:
idxrng = max(1, Ivs(k) - 5) : min(Ive(k) + 5, numel(t));
That should work. (I didn’t do that initially because I didn’t realise you would want to use that part of my code.)
I also changed the offsets from 10 to 5 for clarity. You can change them back (or to other values) if you want to.
If you’re still having problems after making this change, please post the data so that I can work with it.
EDIT —
With respect to the switching transients, I was going to suggest adding MinSeparation, however the documentation indicates that does not work with 'flatselection' because the function considers the flat selection to be a single point. (The problem with any sort of filtering — using smoothdata or an elliptic filter — is that it risks obscuring the sharp transitions of the pulses.) One option would be to use islocalmax and islocalmin without using 'flatselection' and using MinSeparation to return ponly the peaks (these being the initial overshsoots of the switching transients).
Using the current data, that would work liike this —
T1 = readtable('DS0007_REDUCED.CSV')
t = T1.Col1;
s = T1.Col2;
% figure
% plot(t, s)
% grid
% Lvs = islocalmax(s, 'flatselection', 'first', MinProminence=0.5);
% Lve = islocalmin(s, 'flatselection', 'last', MinProminence=0.5);
Lvs = islocalmax(s, MinSeparation=9, MinProminence=0.5); % New Call
Lve = islocalmin(s, MinSeparation=9, MinProminence=0.5); % New Call
figure
plot(t, s)
hold on
plot(t(Lvs), s(Lvs), '|r')
plot(t(Lve), s(Lve), '|m')
hold off
grid
figure
plot(t, s)
hold on
plot(t(Lvs), s(Lvs), '|r')
plot(t(Lve), s(Lve), '|m')
hold off
grid
xlim([0.011 0.025])
ylim('padded')
Lvt = t >= 0.01;
Ivs = find(Lvs & Lvt);
Ive = find(Lve & Lvt);
for k = 1:numel(Ivs)
idxrng = max(1, Ivs(k) - 5) : min(Ive(k) + 5, numel(t)); % New Code
segments{k} = [t(idxrng) s(idxrng)];
end
figure
tiledlayout(4,2)
for k = 1:numel(segments)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
You can take advantage of the switching transients this way.
I originally used 'flatselection' because it seemed the best option, and it usually is for such pulse trains without switching transients, although it is not universally appropriate. The MinProminence name-value pair is still required.
.
Jason
on 5 Feb 2025
As always, my pleasure!
Adapting my code to your new data turned out not to be as straightforward as I thought it would be. There are still problems, most notably when the first oscillation after the switch is higher than the iniitial peak (or valley). This is the situation in pulses 6 and 10 (and perhaps others).
The MinProminence value can only do so much with respect to detecting the ‘correct’ peak, because it then detects (or fails to detect) others that it should not. This causes problems with the indexing. The best solution for this is to increase the sampling frequency (if possible), so that the signal is sampled such that it detects more points and therefore is likely to return the initial peak of the switching transient as always higher than any subsequent peak. Since that information cannot be recovered (or interpolated) from the present data, t+his is likely the best outcome. I encourage you to experiment with the islocalmax MinProminence parameter.
Try this —
T1 = readtable('TESTPULSEDATA.csv')
Ti = mean(diff(T1.Col1)) % Sampling Interval
Fs = 1/Ti % Sampling Frequency
t = T1.Col1;
s = T1.Col2;
% figure
% plot(t, s)
% grid
% Lvs = islocalmax(s, 'flatselection', 'first', MinProminence=0.5);
% Lve = islocalmin(s, 'flatselection', 'last', MinProminence=0.5);
smax = max(s)
Lvm = [0; 0; s > smax*0.9]; % Pulse Peaks
Lvms = strfind(Lvm(:).', [0 0 1]); % Pulse Starts
per = mean(diff(Lvms)) % Pulse Period (Index Units)
Lvs = islocalmax(s, MinSeparation=ceil(per/3), MinProminence=0.275); % New Call
Lve = islocalmin(s, MinSeparation=ceil(per/3), MinProminence=0.275); % New Call
figure
plot(t, s)
hold on
plot(t(Lvs), s(Lvs), '|r')
plot(t(Lve), s(Lve), '|m')
hold off
grid
figure
plot(t, s, '.-') % Plot Showing Sampling Points
hold on
plot(t(Lvs), s(Lvs), 'vr')
plot(t(Lve), s(Lve), '^m')
hold off
grid
xlim([min(t)-0.0001 min(t)+0.005])
ylim('padded')
% Lvt = t >= 0.01;
% Ivs = find(Lvs & Lvt);
% Ive = find(Lve & Lvt);
Ivs = find(Lvs); % Start Indices
Ive = find(Lve); % ENd Indices
maxlen = min([numel(Ive) numel(Ivs)])
CheckIdx = [Ivs(1:maxlen) Ive(1:maxlen) Ive(1:maxlen)-Ivs(1:maxlen)]; % Check Indexing Results
disp(CheckIdx)
Lvgt = CheckIdx(:,3) > 10; % ‘Logical Vector Greater Than’
CheckLen = nnz(Lvgt)
Ivs = Ivs(Lvgt);
Ive = Ive(Lvgt);
for k = 1:numel(Ivs)
idxrng = max(1, Ivs(k) - 5) : min(Ive(k) + 5, numel(t)); % New Code (Changed From Previous)
segments{k} = [t(idxrng) s(idxrng)];
end
NrPlots = numel(Ive)
NrCols = 3;
NrRows = ceil(NrPlots/NrCols)
figure
tiledlayout(NrRows,NrCols)
for k = 1:numel(segments)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
sgtitle('All Pulses')
figure
tiledlayout(3,2)
for k = 1:ceil(numel(Ive)/6):numel(Ive)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
sgtitle('Selected Pulses')
figure
tiledlayout(2,1)
for k = [6 10]
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
This code may still need to be tweaked for other data sets. I made it as robust as I could.
.
Jason
on 5 Feb 2025
Star Strider
on 5 Feb 2025
As always, my pleasure!
This is an atypical application of strfind. It generally works to find specific patterns in strings (string array or character vector), however it turns out to also find specific numerical patterns in numeric vectors. In this instance, I’m looking to see where the logical vector matches the [0 0 1] pattern, indicating that it is detecting a transition. (The results are actually offset by 2 because it returns the index of first instance of the pattern, so here the first 0 in the vector. However since I’m then using it to get the differences in the vectors to estimate the period in index units, the offset doesn’t matter. If I wanted the position of the first 1 in the pattern, I would add 2 to the result.)
I’m using the indices here to get the period of the pulse train to use with the MinSeparation name-value pair.
.
Jason
on 8 Feb 2025
Hi Star Strider, is it possible to suggest a way around this issue that I have found. It occurs when a valley is found before a peak, such as:
I thought the following would handle this:
Ivs = find(Lvs); % Start (Peaks) Indices
Ive = find(Lve); % End (Valleys) Indices
Ive(1,1)
Ivs(1,1)
if Ive(1,1)<=Ivs(1,1) % If first valley is before first peak, delete it
disp('TRUE')
Ive(1,1)=[];
end
But I get the following
ans =
261
ans =
6001
TRUE
A null assignment can have only one non-colon index.
Error in PulseAnalysis/FindPeakMaxMinAUTOButtonPushed (line 332)
Ive(1,1)=[];
Jason
on 8 Feb 2025
I’vee been working on this for a while this morning.
I’m having problems getting an outlier detection approach to work that works in both instances (both data sets). It should be possible.
EDIT — (8 Feb 2025 at 18:14)
I checked my code with respect to both of these data sets, and it turns out that the second version of my code is robust to both of them. I also put the code in a function. At present it returns nothing, however if you want it to return an output, I just need to know what that output is. It should be fairly straightforward to add it.
Try this —
csvfiles = dir('*.csv');
csvfiles = [csvfiles; dir('*.CSV')]
for k = 1:numel(csvfiles)
analysePulseTrain(csvfiles(k).name)
end
function analysePulseTrain(filename)
titlestr = string(strrep(filename,"_","\_"));
T1 = readtable(filename)
t = T1{:,1};
s = T1{:,2};
Ti = mean(diff(t)) % Sampling Interval
Fs = 1/Ti % Sampling Frequency
% figure
% plot(t, s)
% grid
smax = max(s)
Lvm = [0; 0; s > smax*0.9]; % Pulse Peaks
Lvms = strfind(Lvm(:).', [0 0 1]); % Pulse Starts
per = mean(diff(Lvms)) % Pulse Period (Index Units)
Lvs = islocalmax(s, MinSeparation=ceil(per/3), MinProminence=0.275); % New Call
Lve = islocalmin(s, MinSeparation=ceil(per/3), MinProminence=0.275); % New Call
figure
plot(t, s)
hold on
plot(t(Lvs), s(Lvs), '|r')
plot(t(Lve), s(Lve), '|m')
hold off
grid
title(titlestr)
figure
plot(t, s, '.-') % Plot Showing Sampling Points
hold on
plot(t(Lvs), s(Lvs), 'vr')
plot(t(Lve), s(Lve), '^m')
hold off
grid
xlim([min(t)-0.0001 min(t)+0.005])
ylim('padded')
title(titlestr)
Ivs = find(Lvs); % Start Indices
Ive = find(Lve); % ENd Indices
maxlen = min([numel(Ive) numel(Ivs)])
CheckIdx = [Ivs(1:maxlen) Ive(1:maxlen) Ive(1:maxlen)-Ivs(1:maxlen)]; % Check Indexing Results
disp(CheckIdx)
Lvgt = CheckIdx(:,3) > 10; % ‘Logical Vector Greater Than’
CheckLen = nnz(Lvgt)
Ivs = Ivs(Lvgt);
Ive = Ive(Lvgt);
for k = 1:numel(Ivs)
idxrng = max(1, Ivs(k) - 5) : min(Ive(k) + 5, numel(t)); % New Code (Changed From Previous)
segments{k} = [t(idxrng) s(idxrng)];
end
NrPlots = numel(Ive)
NrCols = 3;
NrRows = ceil(NrPlots/NrCols)
figure
tiledlayout(NrRows,NrCols)
for k = 1:numel(segments)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
sgtitle(["All Pulses" titlestr])
figure
tiledlayout(3,2)
for k = 1:ceil(numel(Ive)/6):numel(Ive)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
sgtitle(["Selected Pulses" titlestr])
% figure
% tiledlayout(2,1)
% for k = [6 10]
% nexttile
% plot(segments{k}(:,1), segments{k}(:,2))
% hold on
% plot(t(Ivs(k)), s(Ivs(k)), 'vr')
% plot(t(Ive(k)), s(Ive(k)), '^r')
% hold off
% grid
% title(string(k))
% ylim('padded')
% end
end
.
Thanks so much for all your help!
" At present it returns nothing, however if you want it to return an output, I just need to know what that output is. It should be fairly straightforward to add it. "
As you ask, I want to be able to monitor the pulse seperations and identify any anomalies. The below graphic shows what Im tring to do (so in the pic, I want to draw an xline at the x location on the top plot correspondng to the pulse number picked out in the bottom plot.
So I have added the following to your fantastic code.
% Calc Pulse Seperations, but need to ensure the following
if Ive(1,1)<=Ivs(1,1) %If valley comes before a peak
disp('TRUE')
Ive(1,:)=[]; % Cannot do this: Ive(1,1)=[];
end
f=1000; %y value multiplier (to work in us rather than ms)
sep=diff(x(Ivs));
ax=app.UIAxes2; % Bottom Graph
cl=[0.47 0.25 0.80];
plot(ax,f*sep,'LineWidth',2,'Color',cl); grid(ax,"on");
title(ax,'Pulse Seperations','Color','w','FontSize',14,'FontWeight','normal');
ylim(ax,'padded'); xlabel(ax,'Pulse #'); ylabel(ax,'Period (us)');
ytickformat(ax,'%0.2f');
And the code behind my "Find Pulse Number" is
ax=app.UIAxes; % Top Graph
peaks=app.peaks; % Public property so can access anywhere in GUI
head(peaks,5)
pn=app.PeakNumberEditField.Value;
xloc=peaks(pn)+1
hold(ax,"on");
xline(ax,xloc,'g-'); hold(ax,"off");
(app.peaks) was previously saved in this line of the peak finding routine
Lvs = islocalmax(y, MinSeparation=ceil(per/3), MinProminence=mp); app.peaks=x(Lvs);
Lve = islocalmin(y, MinSeparation=ceil(per/3), MinProminence=mp); app.valleys=x(Lve);
But its not quite showing the xline in the correct location
Jason
on 9 Feb 2025
Star Strider
on 9 Feb 2025
I’m not writing my code in the context of an app, so I’m not sure what changies you’re making to it.
Beyond that, I’m having problems making it robust to everything you want to present to it. (I worked on it for a few hours this morning, then took a break from it.) I’m sure this is possible, just not straightforward.
I’m also running all three (so far) .csv files with it to be certain it works with all of them each time I change it.
Star Strider
on 10 Feb 2025
Edited: Star Strider
on 12 Feb 2025
The significant change here is to get the miidpoint rises and falls of the pulses, and then to use interp1 to find the nearest peak or valley to it (as detected by islocalmax or islocalmin). I have not looked in detaiil at every pulse, however the ones I randomty selected in the plots (that show enough requisite detail) seem to be correct, and the first and last pulse times make sense.
This version actually seems to be robust to every .csv file available to it —
csvfiles = dir('*.csv');
csvfiles = [csvfiles; dir('*.CSV')]
for k = 1:numel(csvfiles)
fn = csvfiles(k).name
analysePulseTrain(csvfiles(k).name)
end
function analysePulseTrain(filename)
titlestr = string(strrep(filename,"_","\_"));
format shortG
T1 = readtable(filename)
t = T1{:,1};
s = T1{:,2};
figure
plot(t,s)
grid
Ti = mean(diff(t)) % Sampling Interval
Fs = 1/Ti % Sampling Frequency
[smin,smax] = bounds(s)
smdn = median([smin smax])
Lvm = [0; 0; s <= smdn]; % Pulse Peaks
Lvms = strfind(Lvm(:).', [0 0 1])-2; % Pulse Starts
per = mean(diff(Lvms)) % Pulse Period (Index Units)
tper = t(ceil(per))
dLvms = [0 diff(Lvms)];
tLvms = t(Lvms).';
Lvme = strfind(Lvm(:).', [1 0 0]); % Pulse Ends
dLvme = [0 diff(Lvme)];
tLvme = t(Lvme).';
Lvs = islocalmax(s, MinSeparation=ceil(per*0.5), MinProminence=0.25); % New Call
Lve = islocalmin(s, MinSeparation=ceil(per*0.5), MinProminence=0.25); % New Call
Ivs = find(Lvs);
Ive = find(Lve);
Ivsi = interp1(Ivs, 1:numel(Ivs), Lvms, 'nearest'); % Interpolate ‘Lvms’ To 'nearest' ‘Ivs’ Index
Ivei = interp1(Ive, 1:numel(Ive), Lvme, 'nearest'); % Interpolate ‘Lvme’ To 'nearest' ‘Ive’ Index
Ivs = Ivs(isfinite(Ivsi));
Ive = Ive(isfinite(Ivei));
figure
plot(t, s)
hold on
% plot(t(Lvs), s(Lvs), '|r')
% plot(t(Lve), s(Lve), '|g')
plot(t(Ivs), s(Ivs), '|r')
plot(t(Ive), s(Ive), '|g')
hold off
grid
title(titlestr)
figure
plot(t, s, '.-') % Plot Showing Sampling Points
hold on
% plot(t(Lvs), s(Lvs), 'vr')
% plot(t(Lve), s(Lve), '^m')
plot(t(Ivs), s(Ivs), 'vr')
plot(t(Ive), s(Ive), '^m')
hold off
grid
xlim([min(t)-0.0001 min(t)+0.015])
ylim('padded')
title(titlestr)
% xline(tLvms, '--r')
% xline(tLvme, '--g')
Ivs = find(Lvs); % Start Indices
Ive = find(Lve); % End Indices
lens = [numel(Ive) numel(Ivs)]
maxlen = min([numel(Ive) numel(Ivs)])
Ivs = Ivs(1:maxlen);
Ive = Ive(1:maxlen);
RelevantTimes = t([Ivs Ive]);
disp("Relevant Times [First Last] —")
disp(RelevantTimes)
disp("Number of Pulse Starts = "+numel(Ivs))
disp("Number of Pulse Ends = "+numel(Ive))
disp("Check Index = ")
CheckIdx = [Ivs(1:maxlen) Ive(1:maxlen) Ive(1:maxlen)-Ivs(1:maxlen)]; % Check Indexing Results
disp(CheckIdx)
Lvgt = CheckIdx(:,3) > 10; % ‘Logical Vector Greater Than’
CheckLen = nnz(Lvgt)
Ivs = Ivs(Lvgt);
Ive = Ive(Lvgt);
for k = 1:numel(Ivs)
idxrng = max(1, Ivs(k) - 5) : min(Ive(k) + 5, numel(t)); % New Code (Changed From Previous)
segments{k} = [t(idxrng) s(idxrng)];
end
NrPlots = numel(Ive)
NrCols = 3;
NrRows = ceil(NrPlots/NrCols)
figure
tiledlayout(NrRows,NrCols)
for k = 1:numel(segments)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
sgtitle(["All Pulses" titlestr])
figure
tiledlayout(3,2)
for k = 1:ceil(numel(Ive)/6):numel(Ive)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
sgtitle(["Selected Pulses" titlestr])
% figure
% tiledlayout(2,1)
% for k = [6 10]
% nexttile
% plot(segments{k}(:,1), segments{k}(:,2))
% hold on
% plot(t(Ivs(k)), s(Ivs(k)), 'vr')
% plot(t(Ive(k)), s(Ive(k)), '^r')
% hold off
% grid
% title(string(k))
% ylim('padded')
% end
end
EDIT — (12 Feb 2025 at 13:22)
Corrected typographical error in code.
.
Star Strider
on 10 Feb 2025
Thank you!
As always, my pleasure!
The ‘m’ is for ‘midpoint’. The ‘s’ and ‘e’ are for ‘start’ and ‘end’ respectively. The initial ‘L’ is for ‘logical’, the ‘I’ for ‘index’ (numerical) and ‘v’ for ‘vector’.
Jason
on 11 Feb 2025
Star Strider
on 11 Feb 2025
It’s appropriate, although I generally don’t do anything wiith GUIs. I’ve simply never needed to use them (although that might change with a project I’m now working on).
When I checked the uitable documentation, I didn’t see any mention of size limits, so that may only be with respect to available memory. I doubt that table arrays have any such limitation, and they do support tall arrays (see Extended Capabilities for that information). The uitable documentation did not mention their excluding tall arrays (or for that matter, anything at all about tall arrays).
Jason
on 12 Feb 2025
Star Strider
on 12 Feb 2025
As always, my pleasure!
It should be an equal sign, and normally something like that would throw an error. (In that instance, I’d have caught it.)
I changed (fixed) it, and re-ran that code, although that change doesn’t seem to make any difference in the results, probably because it didn’t change ‘Ive’, since ‘Ive’ was already defined as (1:maxlen), fortuantely. If it hadn¹t been, the ‘RelevantTimes’ assignment would have thrown a horzcat error, and I‘d have caught it then.
Good pick-up!
.
Im still seeing a few hiccups (as I collect new data)
I fully understand if you have already committed enough time for this.
1: I sometimes get this error
tLvms = t(Lvms).';
Array indices must be positive integers or logical values.
Error in PulseAnalysis/FindPeakMaxMinAUTOButtonPushed (line 303)
tLvms = x(Lvms).';
2: And, i've noticed now that when there is a missing peak (this is exactly what I want to analyse!!)
it captures a false valley.
Perhaps it doesn't matter, as ultimately what Im interested in is the seperations between the peaks which allows me to see how consistent the pulse periods are and highlight any missing peaks. But perhaps these mis-called valleys will show up on the subplots showing a selection of pulses
UPDATE: It does matter, the pulse seperation starts becoming negative
head and tail of CheckIndex:
Check Index =
2 51 49
199 251 52
399 451 52
599 651 52
799 851 52
999002 992054 -6948
999202 992254 -6948
999402 992454 -6948
999602 992654 -6948
999802 992854 -6948
CheckLen =
As a quick fix, changing the Min Prominence value to 2 knocked out those false valleys
Lvs = islocalmax(y, MinSeparation=ceil(per*0.5), MinProminence=0.25);
Lve = islocalmin(y, MinSeparation=ceil(per*0.5), MinProminence=2); % Changed Value HERE
Jason
on 12 Feb 2025
Star Strider
on 12 Feb 2025
Edited: Star Strider
on 12 Feb 2025
The offset should have been +2 insteead of -2, since strfind counts from the first index of the pattern. The first detected instance of the pattern returned should then have been 3 rather than 1, and similarly for the rest of them. (I was probably just tired when I wrote that.)
Try this —
csvfiles = dir('*.csv');
csvfiles = [csvfiles; dir('*.CSV')];
fprintf('\n\nNumber of .csv files in this run: %d\n\n', numel(csvfiles))
for k = 1:numel(csvfiles)
fprintf(['\n\n' repmat('=',1,80) '\n'])
fn = csvfiles(k).name
analysePulseTrain(csvfiles(k).name)
end
function analysePulseTrain(filename)
titlestr = string(strrep(filename,"_","\_"));
format shortG
T1 = readtable(filename)
t = T1{:,1};
s = T1{:,2};
figure
plot(t,s)
grid
title(titlestr)
Ti = mean(diff(t)); % Sampling Interval
Fs = 1/Ti; % Sampling Frequency
[smin,smax] = bounds(s);
smdn = (smax+smin)/2;
Lvm = [0; 0; s <= smdn]; % Pulse Peaks
Lvms = strfind(Lvm(:).', [0 0 1])+2; % Pulse Starts
per = mean(diff(Lvms)); % Pulse Period (Index Units)
tper = t(ceil(per));
fprintf('\n\nSampling Period \t= %10.3E s\nSampling Frequency \t= %10.3E Hz\nAmplitude Minimum \t= %10.3E\nAmplitude Maximum \t= %10.3E\nAmplitude Median \t= %10.3E\n\n', Ti, Fs, smin, smax, smdn)
dLvms = [0 diff(Lvms)];
tLvms = t(Lvms).';
Lvme = strfind(Lvm(:).', [1 0 0]); % Pulse Ends
dLvme = [0 diff(Lvme)];
tLvme = t(Lvme).';
Lvs = islocalmax(s, MinSeparation=ceil(per*0.5), MinProminence=0.25); % New Call
Lve = islocalmin(s, MinSeparation=ceil(per*0.5), MinProminence=0.25); % New Call
Ivs = find(Lvs);
Ive = find(Lve);
Ivsi = interp1(Ivs, 1:numel(Ivs), Lvms, 'nearest'); % Interpolate ‘Lvms’ To 'nearest' ‘Ivs’ Index
Ivei = interp1(Ive, 1:numel(Ive), Lvme, 'nearest'); % Interpolate ‘Lvme’ To 'nearest' ‘Ive’ Index
Ivs = Ivs(isfinite(Ivsi));
Ive = Ive(isfinite(Ivei));
figure
plot(t, s)
hold on
% plot(t(Lvs), s(Lvs), '|r')
% plot(t(Lve), s(Lve), '|g')
plot(t(Ivs), s(Ivs), '|r')
plot(t(Ive), s(Ive), '|g')
hold off
grid
title(titlestr)
figure
plot(t, s, '.-') % Plot Showing Sampling Points
hold on
% plot(t(Lvs), s(Lvs), 'vr')
% plot(t(Lve), s(Lve), '^m')
plot(t(Ivs), s(Ivs), 'vr')
plot(t(Ive), s(Ive), '^m')
hold off
grid
xlim([min(t)-0.0001 min(t)+0.015])
ylim('padded')
title(titlestr)
% xline(tLvms, '--r')
% xline(tLvme, '--g')
Ivs = find(Lvs); % Start Indices
Ive = find(Lve); % End Indices
lens = [numel(Ive) numel(Ivs)]
maxlen = min([numel(Ive) numel(Ivs)])
Ivs = Ivs(1:maxlen);
Ive = Ive(1:maxlen);
RelevantTimes = t([Ivs Ive]);
disp("Relevant Times [First Last] —")
disp(RelevantTimes)
disp("Number of Pulse Starts = "+numel(Ivs))
disp("Number of Pulse Ends = "+numel(Ive))
disp("Check Index = ")
CheckIdx = [Ivs(1:maxlen) Ive(1:maxlen) Ive(1:maxlen)-Ivs(1:maxlen)]; % Check Indexing Results
disp(CheckIdx)
Lvgt = CheckIdx(:,3) > 10; % ‘Logical Vector Greater Than’
CheckLen = nnz(Lvgt)
Ivs = Ivs(Lvgt);
Ive = Ive(Lvgt);
for k = 1:numel(Ivs)
idxrng = max(1, Ivs(k) - 5) : min(Ive(k) + 5, numel(t)); % New Code (Changed From Previous)
segments{k} = [t(idxrng) s(idxrng)];
end
NrPlots = numel(Ive)
NrCols = 3;
NrRows = ceil(NrPlots/NrCols)
figure
tiledlayout(NrRows,NrCols)
for k = 1:numel(segments)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
sgtitle(["All Pulses" titlestr])
figure
tiledlayout(3,2)
for k = 1:ceil(numel(Ive)/6):numel(Ive)
nexttile
plot(segments{k}(:,1), segments{k}(:,2))
hold on
plot(t(Ivs(k)), s(Ivs(k)), 'vr')
plot(t(Ive(k)), s(Ive(k)), '^r')
hold off
grid
title(string(k))
ylim('padded')
end
sgtitle(["Selected Pulses" titlestr])
% figure
% tiledlayout(3,2)
% for k = numel(segments)-5: numel(segments)
% nexttile
% plot(segments{k}(:,1), segments{k}(:,2))
% hold on
% plot(t(Ivs(k)), s(Ivs(k)), 'vr')
% plot(t(Ive(k)), s(Ive(k)), '^r')
% hold off
% grid
% title(string(k))
% ylim('padded')
% end
end
EDIT — Ran all the files again to be certain that my revised code works with all of them. It seems to. Also, some aesthetic changes.
.
Star Strider
on 12 Feb 2025
‘in some instances Im seeing 628 for n1 and only 146 for n2.’
The original data are counted in ‘n1’ however there have to be an equal number of starting and ending indices, annd that is forced in ‘n2’.
My code now seems to be detecting valleys that don’t have a prior peak (or have a very small peak, a small fraction — perhaaps 1% — of the other peak amplitudes in ‘TestTableData5.csv’. The way my code is currently set up, those result in negative pulse lengths for the following pulses and those pulses are eliminated by thresholding the pulse lengths. I’ll have to devise a fix for that, since it didn’t occur in the previous data sets. That will involve either returning those peaks (that might very easily mean returning and filtering a lot of noise, so that might not be practical since those peak amplitudes are similar to the noise amplitude and it may not be possible to distinguish them from the noise) or devising a way of eliminating those valleys from the total. Eliminating them might be easiest.
.
I made this adaptable in order to get it to work. (I’ve actually been working on it for several hours a day over the past few days. These data are not easy to work with, part of that being that they’re noisy. I usually deal with that by filtering or smoothing, however that’s not appropiriate here because filtering would destroy the pulse morphology.) It may be necessary for you to make it further adaptable for data that it may still have problems with (and so it hasn’t yet seen), although it works with all the provided files.
The ‘analysePulseTrain’ function now has an output, that being a table with the starting and ending times and indices for each identified pulse in the file.
Try this —
csvfiles = dir('*.csv');
csvfiles = [csvfiles; dir('*.CSV')];
fprintf('\n\nNumber of .csv files in this run: %d\n\n', numel(csvfiles))
fprintf(['\n\n' repmat('=',1,80)])
for k = 1:numel(csvfiles)
fprintf(['\n\n' repmat('=',1,80) '\n'])
fn{k} = csvfiles(k).name
PulseTable{k,:} = analysePulseTrain(fn{k});
end
for k = 1:numel(PulseTable)
fprintf(['\n\n' repmat('=',1,80) '\n'])
fprintf('\n\nFile: %s\nPulseParameters Size: %5d Rows, %2d Columns\nExcerpt —\n\n',fn{k}, size(PulseTable{k}))
Pulse_Parameters = PulseTable{k}([1:5 end-4:end],:)
end
function PulseParameters = analysePulseTrain(filename)
PulseParameters = [];
titlestr = string(strrep(filename,"_","\_"));
format shortG
T1 = readtable(filename);
fprintf('\n\nFile: %s: %5d Rows, %2d Columns\n\n',filename, size(T1))
t = T1{:,1};
s = T1{:,2};
% figure
% plot(t,s)
% grid
% title(titlestr+" 1")
[Lx,Px] = islocalmax(s);
[Ln,Pn] = islocalmin(s);
% Q0 = [nnz(Lx) nnz(Ln)]
maxkPx = maxk(Px,5);
maxkPn = maxk(Pn,5);
maxkPxm = mean(maxkPx);
maxkPnm = mean(maxkPn);
Pxmax = maxkPxm*0.85;
Pnmax = maxkPnm*0.85;
MinSep = 150;
Lvs = islocalmax(s, MinProminence=Pxmax, MinSeparation=MinSep);
Lve = islocalmin(s, MinProminence=Pnmax, Minseparation=MinSep);
Ivs1 = find(Lvs);
Ive1 = find(Lve);
Ive1 = Ive1(Ive1>Ivs1(1));
% Q3 = [nnz(Lvs) nnz(Lve)]
if numel(Ive1) < numel(Ivs1)
Lve = islocalmin(s, MinProminence=0.05, Minseparation=MinSep);
end
Ivs = find(Lvs);
Ive = find(Lve);
% Q1 = [numel(Ivs) numel(Ive)]
Ive = Ive(Ive>Ivs(1));
% Q2 = [numel(Ivs) numel(Ive)]
% tStart = t(Ivs(1:min(10,numel(Ivs))));
% Ivsv = Ivs(1:min(10,numel(Ivs)))
% Ivev = Ive(1:min(10,numel(Ive)))
figure
plot(t,s)
hold on
plot(t(Ivs), s(Ivs), '|g', MarkerSize=25)
plot(t(Ive), s(Ive), '|r', MarkerSize=25)
hold off
grid
ylim('padded')
title(titlestr+" 1")
PulseParameters = table((1:numel(Ivs)).', t(Ivs), t(Ive), Ivs, Ive, VariableNames=["Pulse","Start Time","Stop Time","Start Index","Stop Index"]);
end
These are the most challenging data sets I’ve encountered here.
.
Jason
on 17 Feb 2025
Star Strider
on 17 Feb 2025
As always, my pleasure!
I wasn’t going to give up on this until I got a reasonable approach, since it seems to be important to what you’re doing. The key difference is that it is now adaptible. For whatever reason, with the MinProminence value set a a higher value, the last ‘valley’ isn’t detected in some files. With a lower value it is, however with some of the other files, that also identifies a lot of noise. So this approach first uses the higher value, and if then are fewer ‘valleys’ that ‘peaks’ (and all the valleys occurring after the first ‘peak’), it uses the lower value and checks those results. (Thus far it hasn’t needed to add any other adaptations.) With this approach, it results in identifying all the peaks and valleys for each file. It then checks to be sure there are an equal number of beginning and end indices, and calculates and creates the ‘PulseParameters’ table.
It would be difficult to have the results table have the same name as the original file (without the extension), however if you want to store the results files in an identifiable way, create a .mat file with that name:
fn = extractBefore(filename,'.');
save([fn ' PulseParameters.mat'], 'PulseParameters')
In the ‘analysePulseTrain’ function, add these lines after creating the table. That should result in a file with a matching name, for example:
DS0007_REDUCED PulseParameters.mat
I thought of that later.
.
More Answers (0)
Categories
Find more on Descriptive Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
