Finding peak width at threshold

I have data with two vectors- one vector corresponds to time and other corresponds to voltages. I found the peaks and corresponding time indices of my data using findpeaks function. However, I want to find the peak widths at thresholds (handdrawn by a blackline in the 2nd picture below) rather than half height or halfprominence(appearing in yellow line). Best would be if I can also see the widths on the graph similar to default findpeaks.

While I can find the time indices corresponding to thresholds(which is an array I obtain from another piece of code, not given here) and subtract them from time indices corresponding to maxvol ( my peak value) to obtain the width of rising part of peak, I cannot think of any logic to find the time index in the dropping down part of the peak, where it would again cross the value equal to threshold voltage.

I am rather new to matlab, so have not been able to find a way to do this. Can anyone help me out.

%%
%1. concatenate all the traces
recording = load(sprintf('LSN16-021118-11f%d.mat'));%load data into struct
fields = fieldnames(recording);%get field names of struct
for k = 1:numel(fields); % numel function returns the number of elements in this array
  ind_=  sscanf(fields{k},'Trace_%i_%i_%i_%i');%disentangle name into numbers
  ind(k,:) = ind_';%save indices in long list
end 
relind = ind(:,3);%relevant index is in 3rd column and can change for recordings
[sortind, order] = sort(relind);%sort
CC = [];%initialize a new list
for kk = 1:length(order)
  k=order(kk);
  CC = cat(1,CC,recording.(fields{order(k)}));%concat, CC is current clamp
end
%%  
 % 2. Plot the traces 
spk.time = CC (:,1);
spk.vol = CC (:,2);
spk.sti= CC (:,3);
spk; 
figure(1)
plot (spk.time, spk.vol, 'black')
xlabel 'Time(sec)'
ylabel 'Voltage (V)'  

% 3. find peaks in traces
trace_length=10; % 10 seconds
Fs = length(CC)./trace_length;% Fs= sampling frequency
findpeaks(spk.vol,'minpeakheight',0,'Annotate','extents');
title('Signal Peak Widths')
[maxvol,locs_maxvol,width,p]=findpeaks(spk.vol,'minpeakheight',0);
thresholdvoltage= -0.03;
APheight= maxvol-(thresholdvoltage);
meanAPheight= mean(APheight)*1000
spike_width = mean(width),
figure(2)
plot(spk.time, spk.vol, 'black',(locs_maxvol./Fs),maxvol,'o')

[dpb]

As noted in comment to your comment at <Answer_268729>, need to see what the data are and what you tried...

[Poonam Thakur]

I have edited the question to provide more information. I also included the code and pictures of graphs I obtain. Do you also need to look at the raw data?

Really appreaciate your help here.

[dpb]

Better to use the Save function/button on a figure to write as .jpg format and then attach it...it'll be small enough to be able to see in context.

To be able to really do anything, yeah, we'd need at least a representative sample of the raw data; again, the best way would be to SAVE a section that covers a few of the peaks to be able to see just what that looks like.

Is the above result what you think you want or is it still not actually the desired result?

[Poonam Thakur]

I have provided the graphs as attachments now. So they should be easier to view. I also attached my raw data. I know above result is not what I want because it is giving me peakwidth at halfprominence. I can also get the peak width at half height using default settings in findpeaks. But I am not able to formulate an argument in the code on how to get the width at threshold. In the screenshot where I show a zoomed in picture of a single spike- I desire to find the width of black line (I hand-drew that).

[dpb]

OK, yes, that's much better...thanks! :) For figures, you can use the funny blue-image icon in the 'INSERT' section and insert the .jpg inline instead of as an attachment and that is often handy to have the picture in context...

Let me have a little while and I'll see if can't get the previous to work with your data; it's very clean-looking it appears.

The problem is to return a width at a user-provided level, correct?

[Poonam Thakur]

Yes, I do not have any noise issue, great peak to noise ratio. Just looking to find width at user given level.

Answer 1

[Akira Agata]

One simple way to do your task woudl be utilizing pulsewidth function. By setting 'StateLevels' option to [th - delta, th + delta], where th is your given threshold value and delta is some value.

For example, when setting th = 0.0 and delta = 0.02, the pulse width at y = 0.0 can be calculated as follows (sorry for some Japanese language at the figure due to my MATLAB settings...):

load('LSN16-021118-11f.mat');

t = Trace_1_2_1_1(:,1);
y = Trace_1_2_1_1(:,2);

pulsewidth(y,t,'StateLevels',[-0.02 0.02])

To obtain pulse width value for each pulse, please set the output variable, like:

pWidth = pulsewidth(y,t,'StateLevels',[-0.02 0.02]);

[Poonam Thakur]

Thank you very much Akira. I tried it out and see that I have a much better control of where to get the width. Though I have to find a way to get the appropriate delta. Since my thresholds will be different for each peak, I will need to design a loop to automatically find an appropriate delta value for each peak. Once I figure out a way, I will accept this answer and post the new code.

Thanks a lot.

[dpb]

That's cute function I was unaware of, but will be difficult for the specific case, methinks to find the exact, general application.

Poonam, I ran out of energy last night, unfortunately the previous "trick" doesn't help us here straight as written because the setpoint is below the previous location for many/most -- hence the NaN because was extrapolating.

I'll look at it again this morning after coffee; it appears there are two ways here, to do as findpeaks does and walk along the trace until find the crossings (why TMW didn't implement this with a user-defined threshold is anybody's guess?) or, since your case is pretty-much noise free, you could just independently find the crossings and they'll be where you are looking for.

Patience just a little longer... :)

[dpb]

Actually, I just experimented with pulsewidth and it doesn't do what Poonam wants at all--it finds the equivalent pulse width between the min/max values given, but where creates a square pulse that is centered between the min/max values and the width returned is where that pulse crosses the rise/fall of the peak.

Poonam, otoh, is looking for the base width of the peak at an arbitrary baseline value; to have a square pulse cross the waveform at that location the way pulsewidth derives it would require the lower limit to be at some arbitary location such that it has the crossings at the desired level, not the limits at the level.

I thought I had overlooked it initially and it was, indeed, extremely clever and would work, so I tried the peak and threshold for one...

>> pulsewidth(T(1:4000,2),x(1:4000),'statelevels',[refHt pks(1)])

which results in

from which it's apparent the problem -- to have the width at the location of the lower dotted line one would have to have the lower of the two limits at some arbitrary location around -0.04 or somewhere such that the gray would be the width of that line's intersection points.

One could write an anonymous function and use fsolve to iterate to find that location probably, but it's much better to find directly.

[dpb]

ADDENDUM:

Actually, you can't make it work using the peak at all; it requires there be points on both sides of the peak for the minimum so the lowest value that doesn't return empty is roughly -0.038; to force a wide pulse at base of the the peak one has to move both ends of the limit range arbitrarily -- by trial and error I cheated and set the first point in the wave form to min( T(1,4000,2)) and the upper limit to 0.015. That gets somewhere close but is totally arbitrary and wouldn't be feasible to code as a general solution...

T(1,2)=min(T(1:4000,2));
refHt=T(1,2);
pkHt=-0.015;
pulsewidth(T(1:4000,2),x(1:4000),'statelevels',[refHt pkHt])
ans =
   69.5458
xlim([2400 2800])

results in

[Poonam Thakur]

Thanks a lot for trying. I couldn't try this solution because I ran into problems with part of the code to calculate thresholds. If I am not able to solve that one, maybe another question is coming your way. :)

[dpb]

What's wrong with the solution I posted as the Answer just below?

This was just a sidebar looking at the pulsewidth function I'd not been aware of before; it's cute but really has no bearing on your particular problem.

Answer 2

[dpb]

OK, finally(! :) ) had time to sit down for a few minutes...it's more convoluted than if had just started from scratch, but I thought initially if I picked out the sections from findpeaks that does the locating for the two ways it does that would be quick to modify -- wrong! but, it does work altho could be done quite a lot more simply if want to take the time to streamline it some.

To use as is, use findpeaks to return the peak locations using the array indices as the X values(*), then call the function; it returns a 2D array array of the width low/hi bounds in position units; these can be converted to time by multiplying by the dt or by interpolation for "real" units.

The m-file containing the function and its helpers is attached; ended up long enough owing to using existing structure as base to not paste in its entirety.

load LSN16-021118-11f.mat
x=[1:length(Trace_1_2_1_1(:,1))].';                                 % x position vector
[pks,locs]=findpeaks(Trace_1_2_1_1(:,2),x,'MinPeakHeight',0.025);   % find locations
refHt=-0.030;                                                       % desired crossing level
bnds=getpeakwidth(Trace_1_2_1_1(:,2),x,locs,refHt);                 % find the bounds at level

figure
plot(x,Trace_1_2_1_1(:,2))
for i=1:length(bnds),hL(i)=line(bnds(i,:),[refHt refHt],'color','r');end
xlim([2400 2700])

produced the following figure showing the first peak amplified...

(*) The need for using the index is that one must search the underlying waveform for the crossing so if use the time variable have to convert that internally to be and indes into the array. That's doable with more internal coding, but I took the easy way out and just passed the index vector directly.