To find the FWHM and area of pulses in a signal waveform

A = xlsread('z_data.xlsx');
z1 = A(:,1);
z2 = A(:,2);

%% finding beginning of pulse and peak 
z11 = z1;
z21 = smooth(z2,10);        
z3 = diff(z2);
z4 = diff(z1);
z5 = (diff(z21)./diff(z11)) > 0;
z6 = diff(z5);
z7 = logical(diff(z5));
z8 = diff(z5) < 0 ;
z9 = logical(diff(z8));

%% plots
figure
subplot(6,1,1);
plot(z11,z21);
hold on
plot(z11(z7),z21(z7),'.g','markersize',20);
hold off
subplot(6,1,2);
plot(z5);
subplot(6,1,3);
plot(z6);
subplot(6,1,4);
plot(z7);
subplot(6,1,5);
plot(z8);
subplot(6,1,6);
plot(z9);

I have made use the derivative and logical operations one after the other and hence found out the start point and peak of the pulse.

There after not able to find half maximums of each pulse, and the point where the pulses end.

Also i do not have a formula through which the area of such varing waveform could be found out.

[SHANTANU KSHIRSAGAR]

this the file containig the data which has been used

Answer 1

[Star Strider]

Try this:

M = readmatrix('z_data.xlsx');
t = M(:,1);
A = M(:,2);

As = smoothdata(A, 'movmean',15);                                       % Remove Noise
zci = @(v) find(v(:).*circshift(v(:), 1, 1) <= 0);                      % Returns Approximate Zero-Crossing Indices Of Argument Vector (>= R2016b)
dAdt = gradient(As)./gradient(t);                                       % Time Derivative
zx = zci(dAdt);                                                         % Zero-Crossing Indices
[pks,locs,fwhm,prmns] = findpeaks(A, 'MinPeakDistance',50);             % Return Full-Width-Half-Maximum Of Peaks

zx = [1; zx];                                                           % Modify ‘zx’ To Include First Index
Idx = zx(1:2:end);                                                      % Beginnings Of Each Peak
TotArea = cumtrapz(t,A);                                                % Integral
PkAreas = diff(TotArea(Idx));                                           % Peak Areas
PkLocs = locs(1:end-1);                                                 % Eliminate Last (Incomplete) Peak
PkTime = t(PkLocs);

Results = table(t(PkLocs), fwhm(1:end-1), PkAreas, 'VariableNames',{'PeakTime', 'FWHM', 'Area'});

Note: This code is reasonably robust, although fragile in that it may need to be modified slightly for different data sets. If there are incomplete peaks (that do not end in zero-crossings), or it begins with a complete peak (so that the initial index appears twice), those will need to be considered and corrected for in the calculation of the indices (‘zx’ and ‘locs’). It works as written for the data presented.

If you want to see graphically what the code does, plot this example:

figure                                                                  % Information Plot (Not Needed In Code)
plot(t, A)
hold on
plot(t, dAdt, '-r')
plot(t(zx), dAdt(zx), '+g')
plot(t(locs), pks, '^r')
hold off
grid
legend('Data', 'Detivative', 'Zero-Crossings', 'Peaks', 'Location','N')
xlim([ 0 750])
ylim([-5  30])

[SHANTANU KSHIRSAGAR]

Thank you Star Strider, the answer it was really helpful, i will see that from the method you followed, i will make use also to find the area under each pulse .

[Star Strider]

As always, my pleasure!

The ‘Results’ table has the time of the peaks, the full-width-half-maximum values computed by findpeaks, and the areas of the pulses.

[SHANTANU KSHIRSAGAR]

yes thats true but i do not have the signal processing tollbox, so i will need to modify it to a simple code , or i would have used it earlier itself. Let's see..

[Star Strider]

I did not realise that.

I will see if I can revise my code to work without the Signal Processing Toolbox.

[Star Strider]

Here you go:

M = readmatrix('z_data.xlsx');
t = M(:,1);
A = M(:,2);

As = smoothdata(A, 'movmean',15);                                       % Remove Noise
zci = @(v) find(v(:).*circshift(v(:), 1, 1) <= 0);                      % Returns Approximate Zero-Crossing Indices Of Argument Vector (>= R2016b)
dAdt = gradient(As)./gradient(t);                                       % Time Derivative
zx = zci(dAdt); 
locs = find(islocalmax(A, 'MinSeparation', 50));

zx = [1; zx];                                                           % Modify ‘zx’ To Include First Index
Idx = zx(1:2:end);                                                      % Beginnings Of Each Peak
TotArea = cumtrapz(t,A);                                                % Integral
PkAreas = diff(TotArea(Idx));                                           % Peak Areas
PkLocs = locs(1:end-1);                                                 % Eliminate Last (Incomplete) Peak
PkTime = t(PkLocs);
for k = 1:numel(Idx)-1
    thm(1) = interp1(A(Idx(k):locs(k)), t(Idx(k):locs(k)), A(locs(k))/2);
    thm(2) = interp1(A(locs(k):Idx(k+1)), t(locs(k):Idx(k+1)), A(locs(k))/2);
    fwhm(k,:) = diff(thm);
end

Results = table(t(PkLocs), fwhm, PkAreas, 'VariableNames',{'PeakTime', 'FWHM', 'Area'});

It took a bit to figure out how best to calculate the ‘fwhm’ vector in the absence of findpeaks. (You should have all the other functions I use here.) The loop uses interp1 to calculate the ‘t’ values for the full-width-half-maximum in the rising and descending parts of each pulse.

My previous comments are unchanged.

The slightly revised plot code:

figure                                                                  % Information Plot (Not Needed In Code)
plot(t, A)
hold on
plot(t, dAdt, '-r')
plot(t(zx), dAdt(zx), '+g')
plot(t(locs), A(locs), '^r')
hold off
grid
legend('Data', 'Derivative', 'Zero-Crossings', 'Peaks', 'Location','N')
xlim([ 0 750])
ylim([-5  30])

[SHANTANU KSHIRSAGAR]

Thank you Star Strider this is simplified code , it works.

The use of zci = @(v) find(v(:).*circshift(v(:), 1, 1) <= 0) is little difficult to understand for me.

The area calculated is area under each indivisual pulse right, i just needed to confirm.

[Star Strider]

As always, my pleasure.

The ‘zci’ function shifts the argument vector by +1 in the first (row) dimension, then multiplies it with the original vector. Where the signs of the two are the same, the result is positive, where a zero-crossing occurs, it is negative, so detecting the negative values detects the zero-crossing. The original result is a logical vector, the find call returns a numeric vector of the zero-crossing indices. The only problem that can arise with it is if the first and last values of a vector are of different signs, in that event it will return the last value as a zero-crossing when it is not actually one. A simple test for that is to multiply the first and last values. If that is positive, there is no problem. If it is negative, remove the last value returned by ‘zci’.

[SHANTANU KSHIRSAGAR]

I am getting many number of zero crossings on using ourf method. I have attached .jpg file of the same, so not able to calculate exact FWHM and area. One of the reason it is happening is the change in positive pulse width.getting error:

Index exceeds the number of array elements (0), extptids = Xq < X(1) | Xq > X(end); in the thm(2) execution.

[Star Strider]

The ‘comp_1.fig’. signal is not even remotely close to the sort of signal my code was originally designed to analyse. Even when I filter it I cannot do anything with it. For the others, they do not resemble closely enough the signal in your original Question, and my code requires that for it to work correctly. My code is not designed to work with them, and I am not certain that I could write code that could work with them. They have transients that cannot be eliminated by filtering, and are interfering with the detection methods my code uses.

I am sorry, however I cannot do anything with those data.

One thing that you will need to do is to force the data to be column vectors, with:

t = t(:);
A = A(:);

My code will not work with row vectors.

Also, I cannot find this line:

extptids = Xq < X(1) | Xq > X(end);

anywhere in my code. My code executes without error (in R2020a) on those data, so I do not know what is throwing it.

I doubt if it is possible to write code that would be robust to data in these figure files. I will see if I can devise something, however I can make no promises.

[SHANTANU KSHIRSAGAR]

Will the use of fourier transforms help in this case?

[Star Strider]

I tried that in order to design a frequency-selective filter. The Fourier transform did not produce any useful information.

[SHANTANU KSHIRSAGAR]

ok Star, So the fourier transform method too, will not be useful here. I will look for another alternative closer to what you had suggested. I also tried taking derivative twice but it did not help.

[Star Strider]

I worked on those signals for a while and was unable to do anything with them. The findpeaks results are not a problem (except with the first set that is completely unusable), however the areas can not be reliably calculated because it is not possible to accurately detect the beginnings and ends of the pulses.

[SHANTANU KSHIRSAGAR]

yes finding the biginnings and the ending of pulses is the main problem. Thanks Star Strider.

I also got following error:

The grid vectors must contain unique points.

Error in interp1 (line 151)

F = griddedInterpolant(X,V,method);

[SHANTANU KSHIRSAGAR]

we are getting output of approximate zero cressings by using the function handle "zci".

can we filter out the zci's we need by using find and some minseperation in any way.

[Star Strider]

I have been thinking about this, and came up with an alternative approach:

F = openfig('comp_2.fig');                          % Choose One
F = openfig('comp_3.fig');                          % Choose One
Lines = findobj(F, 'Type','Line');
M = Lines(4);
t = M.XData;
A = M.YData;
Idx = islocalmin(A, 'FlatSelection','last');        % Works With These Two
NIdx = find(Idx);

figure(2)
plot(t, A)
hold on
plot(t(Idx), A(Idx), '+r')
hold off
grid
xlim([0 1000])
ylim([-2  2])

This accurately locates the start of each pulse, at least for these two data sets. (The first data set is a lost cause, however you may be able to get this to work with it as well.) That will work for the integration. The findpeaks function will return the other parameters-of-interest. (I did not test it on the signal from the orginal post. It may also work for that one.) You will need to integrate this into the rest of the code.

The islocalmin function was introduced in R2017b, so you have it and will be able to use it.

[SHANTANU KSHIRSAGAR]

I am getting some repeated data points in the arrays and so getting an error using interp1 , that there should be unique points only in the array for using it. Can we do something for that. I have attached the .mat file for the same.

[Star Strider]

I don’t remember what I was doing with the interpolations two months ago, however they aren’t appropriate for this data set.

Use this (slightly revised) version of the original code:

D = load('non_uniqueness.mat');
A1 = D.A1;
t1 = D.t1;

figure
plot(t1, A1)
grid

t = t1;
A = A1;
As = smoothdata(A, 'movmean',50);                                       % Remove Noise
zci = @(v) find(v(:).*circshift(v(:), 1, 1) <= 0);                      % Returns Approximate Zero-Crossing Indices Of Argument Vector (>= R2016b)
dAdt = gradient(As)./gradient(t);                                       % Time Derivative
zx = zci(dAdt);                                                         % Zero-Crossing Indices
[pks,locs,fwhm,prmns] = findpeaks(As, 'MinPeakDistance',50);             % Return Full-Width-Half-Maximum Of Peaks
% zx = [1; zx];                                                           % Modify ‘zx’ To Include First Index
Idx = zx(1:2:end);                                                      % Beginnings Of Each Peak
TotArea = cumtrapz(t,A);                                                % Integral
PkAreas = diff(TotArea(Idx));                                           % Peak Areas
PkAreas2 = PkAreas(PkAreas>0);
% PkLocs = locs(1:end-1);                                                 % Eliminate Last (Incomplete) Peak
% PkTime = t(PkLocs);
PkTime = t(locs);
Results = table(PkTime, fwhm, PkAreas, 'VariableNames',{'PeakTime', 'FWHM', 'Area'})

to get:

Results =

  24×3 table

    PeakTime     FWHM      Area 
    ________    ______    ______

      6.875     26.256    237.76
     75.313     50.043    273.97
     142.19     50.038    474.01
     208.75     50.055     502.8
     275.31     50.046    477.57
     342.19     50.038    484.04
     408.75     50.047    505.02
     475.31     50.049    465.72
     542.19     50.065    506.83
     608.75      50.05    487.37
     675.31     50.048    482.41
     742.19     50.047    520.94
     808.75     50.053    473.42
     875.31     50.051    489.55
     941.88     50.054     513.8
     1008.8     50.041    469.69
     1075.3     50.052    509.54
     1142.2      50.06    505.52
     1208.8     50.047    488.93
     1275.3      50.05    517.39
     1341.9     50.057    495.93
     1408.8     50.029     490.4
     1475.3     50.057    534.83
     1541.9     50.043    485.69

That works, and appears to produce consistent results.

.

[SHANTANU KSHIRSAGAR]

Yeah thats really simple and robust nice technique, you posted it earlier too, but in the absence of signal processing toolbox thats what is my main problem.

[Star Strider]

I’d forgotten about the Signal Processing Toolbox problem.

Try this:

D = load('non_uniqueness.mat');
A1 = D.A1;
t1 = D.t1;

figure
plot(t1, A1)
grid

t = t1;
A = A1;
As = smoothdata(A, 'movmean',50);                                       % Remove Noise
zci = @(v) find(v(:).*circshift(v(:), 1, 1) <= 0);                      % Returns Approximate Zero-Crossing Indices Of Argument Vector (>= R2016b)
dAdt = gradient(As)./gradient(t);                                       % Time Derivative
zx = zci(dAdt); 
zx = zx([0; diff(zx)]> 50);
locs = find(islocalmax(As, 'MinSeparation', 50));
locix = locs > min(zx);                                                 % Restrict ‘locs’ To After First Zero-Crossing
locs = locs(locix);
Idx = zx(1:2:end);                                                      % Beginnings Of Each Peak
TotArea = cumtrapz(t,A);                                                % Integral
PkAreas = diff(TotArea(Idx));                                           % Peak Areas
PkLocs = locs(1:end-1);                                                 % Eliminate Last (Incomplete) Peak
PkTime = t(PkLocs);
for k = 1:numel(Idx)-1
    Q1 = [k Idx(k) locs(k) Idx(k+1)  As(locs(k))]
    thm(1) = interp1(As(Idx(k):locs(k)), t(Idx(k):locs(k)), As(locs(k))/2);
    thm(2) = interp1(As(locs(k):Idx(k+1)), t(locs(k):Idx(k+1)), As(locs(k))/2);
    fwhm(k,:) = diff(thm);
end
Results = table(t(PkLocs), fwhm, PkAreas, 'VariableNames',{'PeakTime', 'FWHM', 'Area'})

There is no robust way to deal with the noise in every data set. This works for this one. The results are different from what findpeaks returns, however they appear to be reasonable when looking at the data.

[SHANTANU KSHIRSAGAR]

Thank it worked for the same data efficiently. I am Trying to work it out for number of signals, its a difficult thing but.

[Star Strider]

As always, my pleasure!

Answer 2

[SHANTANU KSHIRSAGAR]

Thanks Star Strider for quick and such an explainatory answer. It also cleared some new concepts to me other than the proposed problem in matlab.

[Star Strider]

As alwasy, my pleasure!