How can I modify the findpeaks command to give me a width at 80% of peak height (instead of at the build in 'halfheight')?

I am trying to analyze peaks but and need to determine the width of signal at 80% of the peak height. Matlab has a built in criterion when you specify findpeaks(...,'WidthReference",'halfheight') and this gives you 50%.

I have tried to modify the code but have so far been unsuccessful. If anybody can help I would really appreciate it.

Thank you

Answer 1

[dpb]

I'd suggest not modifying findpeaks but to postprocess the results instead--a sample from one of the examples of the basic process I'd undertake:

load mtlb                       % a sample dataset
select = mtlb(1001:1200);       % a subsection that's got a few peaks...
x=1:length(select);             % corollary variable of bin locations
[pks,locs,wdths]=findpeaks(select,x,'MinPeakHeight',2);       % the result for just a few >2

% The "engine" to postprocess...

fnInterpL=@(p,l,w) interp1(select(l-fix(w):l),l-fix(w):l,0.8*p);  % interpolate lower side
fnInterpH=@(p,l,w) interp1(select(l:l+fix(w)),l:l+fix(w),0.8*p);  % and upper...

This can be used as:

>> [arrayfun(fnInterpL,pks,locs.',wdths.') arrayfun(fnInterpH,pks,locs.',wdths.')]
ans =
 51.3944   53.5555
 82.5567   84.8145
145.4656  146.3182
>>

NB: the transpose on the locations/widths vectors; for some reason findpeaks doesn't return all three in same orientation as arrayfun needs.

One could generalize the above a little more by passing in the limit percent as well. Also note the above has the anonymous functions with the builtin signal; probably in a general case would be best to just go ahead and do an m-file where can pass it as well...

ADDENDUM

The above are the boundaries; the width is the difference of course:

>> w80=[arrayfun(fnInterpH,pks,locs.',wdths.')-arrayfun(fnInterpL,pks,locs.',wdths.') ]
w80 =
  2.1611
  2.2578
  0.8526
>>

[dpb]

Had plotted just to see what the signal was and to ensure had signs, etc., right above...

The upper shows the peaks>2 and the bounds returned by findpeaks for the first in red. The lower is the composite of the lower and upper section for interp1 plotted where the index is the dependent variable and the 80% magnitude is the vertical red line for the first peak. Inspection shows the crossing points correlate with the returned min/max values.

Have to split the interpolation into the two pieces since the independent variable being the amplitude is double-valued; have to work on the two sides separately.

One could also do a find for locations that are above/below the threshold and bracket them; this "gets 'er done!" in basically one step.

[dpb]

I'll note it does seem a worthwhile enhancement request to allow an additional input parameter (as if there weren't enough already!) to findpeaks to set the FW computation to a given value (FWUM maybe for FW "User" Max would be a reasonable acronym).

[Poonam Thakur]

Hi,

I understand this is an old thread but taking my chance here. I am intersted in finding peak width at threshold. I tried to use the above code using my threshold value instead of 0.8*p. However, it returned me arrays full of NaN. My data consists of very long arrays and has multiple peaks.

Thank you.

[dpb]

Would probably be best to post new Q? with the specific data (representative subset, of course).

Would need to see the specific code you tried, obviously, to see if could spot any problems there; NaN is a typical result from interp1 when it tries to extrapolate so one could infer that perhaps the ranges to it in the above anonymous function aren't bounding the input points being looked for...

Answer 2

[David Zimmerman]

Hello thread,

First I want to apologize, I did not see the response posted by dpb and I have not personally tryed it but it may be a much neater code than my own. However, this is what I ended up writing. It begins by using findpeaks to get the location and value of each peak. It also set a minimum peak distance and a minimum peak height of 2. Then it proceeds to test if terminal points exist, higher than peaks and then goes through multiple peaks, interpolating right and left and eliminating double peaks that overlap (keeping only the larger peak value). My analysis only looks at peaks with a value greater than 2 so the code begins with an if statement but I suppressed it for you (however, you should change the MinPeakDistance and MinPeakHeight to your own values) [written with matlab 2015]

David Z.

select = inpaint_nans(rate_matrix_left(i_plot,:));
% if (max_y_L >= 2) == 1
            [pks,locs]=findpeaks(rate_matrix_left(i_plot,:),'MinPeakDistance',10,'MinPeakHeight',2);
            
            num = length(pks);
            cell = ones(1,num)*cell_count;
            Low = zeros(1,num);
            High = zeros(1,num);
            L = zeros(1,num);
            R = zeros(1,num);
            w80_L = zeros(1,num);
            
            % checks if there exists a terminal point greater than the peaks
            if isempty(locs) == 0 && max_y_L > max(pks)
                if max_y_L == select(72)
                    locs(num) = 72;                    % come back to bc i think should be num+1
                    pks(num) = select(72);
                elseif max_y_L == select(1)
                    locs(num) = 1;
                    pks(num) = select(1);
                end
            end
            
            % if there is a value greater than 2, but there are not peaks,
            % if statement checks if maximum value occurs at
            % first or last bin
            if isempty(locs) == 1
                if max_y_L == select(72)
                    locs = 72;
                    pks = select(72);
                    cell = cell_count;
                elseif max_y_L == select(1)
                    locs = 1;
                    pks = select(1);
                    cell = cell_count;
                end
            end

            % INTERPOLATION
            for i = 1:length(locs)
                % interpolate left
                L(i) = locs(i);
                while select(L(i)) >= .19*pks(i) && L(i) >= 2
                    L(i) = L(i) - 1;
                end
                if L(i) == 1
                    Low(i) = 1;
                else
                    Low(i) = interp1(select(L(i):L(i)+2),L(i):L(i)+2,0.2*pks(i));
                end
                
                % interpolate right
                R(i) = locs(i);
                while select(R(i)) >= .19*pks(i) && R(i) <= 71
                    R(i) = R(i) + 1;
                end
                if R(i) == 72
                    High(i) = 72;
                else
 %                  High(i) = interp1(select(locs(i):R(i)),locs(i):R(i),0.2*pks(i));
                    High(i) = interp1(select(R(i)-2:R(i)),R(i)-2:R(i),0.2*pks(i));
                end
            end
            
            for i = 1:length(locs)
                if isnan(High(i)) == 1
                    High(i) = 72;
                end
                if isnan(Low(i)) == 1
                    Low(i) = 1;
                end
                w80_L(i)= High(i) - Low(i);
            end
            
            % to get rid of double peaks
            if length(locs)>1
            i = 0;            
            while i <= length(locs)-1               
                if High(i)>locs(i+1) || Low(i+1)<locs(i)
                    if pks(i+1)<pks(i)
                        pks(i+1) = [];
                        locs(i+1) = [];
                        Low(i+1) = [];
                        High(i+1) = [];
                        w80_L(i+1) = [];
                        cell(i+1) = [];
                    else
                        pks(i) = [];
                        locs(i) = [];
                        Low(i) = [];
                        High(i) = [];
                        w80_L(i) = [];
                        cell(i) = [];
                    end
                end
            end
            end
                
            for i = 1:length(locs) 
                hold on
                scatter(Low(i),pks(i)*.2)
                scatter(High(i),pks(i)*.2)
            end
% else
%    [pks,locs, Low, High, w80_L] = deal(NaN);
%    cell = cell_count;
% end

[MARWA NOOR]

hi,

i tried this code ,it is correct but i do not understand how zeros function work here

..i know that is to generate matrix of zeros ,

[David Zimmerman]

hi Marwa,

speaking of the code I posted above, the zeros are being used to preallocate space for 5 values. The values L and R are used during the interpolation, such that after you find he peak you move left ('L') and right ('R') while the value Low and High are created to store the 80% peak hight value, left and right of the peak respectively. Lastly, the value w80_L is the width, so how far away (in x) the lower 80% peak hight value is from the higher 80% peak hight value.