Thread Subject: Area of hysteresis loops

Please could I have help determining the area of series (30+) of hysteresis loops from experimental data? My data is plotted as Load (L) on the y axis against Displacement (D) on the x axis. I used the peakdet function (http://billauer.co.il/peakdet.html) to generate new variables which provide the indices for the minima and maxima of my experimental data. My specific questions are:

a. How do I measure the area of each hystereis loop?
b. How do I plot these values against time?

I am a novice MATLAB user and could not follow previous posts on hysteresis loop area calculations.

Thank you.

"Etienne O'Brien" <etienneobrien@yahoo.com> wrote in message <ic3ofg$jed$1@fred.mathworks.com>...
> Please could I have help determining the area of series (30+) of hysteresis loops from experimental data? My data is plotted as Load (L) on the y axis against Displacement (D) on the x axis. I used the peakdet function (http://billauer.co.il/peakdet.html) to generate new variables which provide the indices for the minima and maxima of my experimental data. My specific questions are:
>
> a. How do I measure the area of each hystereis loop?
> b. How do I plot these values against time?
>
> I am a novice MATLAB user and could not follow previous posts on hysteresis loop area calculations.
>
> Thank you.

Use sum(trapz(x,y)) on each section defined by a change in sign in diff(y) (same as peakdet I presume).

Then depending on your hysteresis loops, which we haven't seen, subtract the bottom part from the top part to get the enclosed area.

If you post sample data it'll be much easier for us to be helpful and we'll be able to provide more accurate feedback/results.

Thank you Sean,

What's the best way for me to post some data as I can't see a file attachment option here.

Etienne

"Etienne O'Brien" <etienneobrien@yahoo.com> wrote in message <ic6hhm$nqe$1@fred.mathworks.com>...
> Thank you Sean,
>
> What's the best way for me to post some data as I can't see a file attachment option here.
>
> Etienne

The one we always used was drop.io but it's been shut down.

Maybe this one?
http://www.4shared.com/

Any other free file hosting site as long as you don't have to be a member to download.

Thank you Sean,

I've uploaded some files. The 'data' file is fairly large but the others very small.

Although I have managed to use 'peakdet', as I have been using MATLAB for a very short time I don't yet understand the inner workings of this function. Therefore I couldn't follow your comment 'Use sum(trapz(x,y)) on each section defined by a change in sign in diff(y) (same as peakdet I presume)'

I will be grateful your your guidance.

http://www.4shared.com/file/lGZ6cltQ/data.html
http://www.4shared.com/file/TD3Nay-5/Hysteresis_plot_with_Peaks_and.html
http://www.4shared.com/file/P35dMyvd/index_Peaks__red_.html
http://www.4shared.com/file/2w3Z-Br5/index_Valleys__green_.html
http://www.4shared.com/file/CtX8DQdH/Key_to_columns_in_data_variabl.html

Etienne

"Etienne O'Brien" <etienneobrien@yahoo.com> wrote in message <ic9rff$kn1$1@fred.mathworks.com>...
> Thank you Sean,
>
> I've uploaded some files. The 'data' file is fairly large but the others very small.
>
> Although I have managed to use 'peakdet', as I have been using MATLAB for a very short time I don't yet understand the inner workings of this function. Therefore I couldn't follow your comment 'Use sum(trapz(x,y)) on each section defined by a change in sign in diff(y) (same as peakdet I presume)'
>
> I will be grateful your your guidance.
>
> http://www.4shared.com/file/lGZ6cltQ/data.html
> http://www.4shared.com/file/TD3Nay-5/Hysteresis_plot_with_Peaks_and.html
> http://www.4shared.com/file/P35dMyvd/index_Peaks__red_.html
> http://www.4shared.com/file/2w3Z-Br5/index_Valleys__green_.html
> http://www.4shared.com/file/CtX8DQdH/Key_to_columns_in_data_variabl.html
>
> Etienne

Etienne, this looks a lot like the load-displacement data we get!

First I would recommend smoothing it. If you zoom in, you'll see that there are lots of little overlaps and bubbles, noise. To account I'll use a moving average filter. You also need to remove the NaN values you have. The code is included, this is fully functional for your data. I've commented what I've done to try and make it easy to understand. You should be able to run this as long as your data and index variables are in your workspace. It should also be fully functional if you wish to make it a function to use on other data (i.e. nothing is hardwired).
Also, when copying from here make sure to reconnect lines broken by the newsreader.

%%%%
%% Extract and Filter
theDisp = filter(ones(1,5),5,data(:,4).'); %5-element moving average
theLoad = filter(ones(1,5),5,data(:,5).');
theLoad(isnan(theLoad)) = []; %Remove NaNs.
theDisp(isnan(theDisp)) = [];

if ~isequal(length(theLoad),length(theDisp))
    %Insure equal length
    %You now want to recalculate the peaks and valleys since the data has
    %changed. I'm not going to do this, since your method seems to be working.
    %I'll use the previous ones, for demonstration etc.
    error('After removing NaNs the data is no longer the same length');
end

%% Engine & Visualization
theAreas = zeros(length(index_Valleys)-1,1); %preallocate
figure;
hold off
for ii = 1:length(theAreas)
    %Engine
    iispan = index_Valleys(ii):index_Peaks(ii); %span of data between current valley and peak
    ip1span = index_Peaks(ii):index_Valleys(ii+1); %span of data between current peak and next valley
    disp_ip1span = theDisp(ip1span); %extracted next part (need to extract to analyze)
    load_ip1span = theLoad(ip1span);
        %Remove parts of the next span that fall outside the range of the current span
        %I.e. remove data that's not below the "ceiling" of the hysteresis loop
        %You may not actually wish to do this. If you don't comment or delete these two lines.
    load_ip1span(disp_ip1span<theDisp(index_Valleys(ii))|disp_ip1span>theDisp(index_Peaks(ii))) = [];
    disp_ip1span(disp_ip1span<theDisp(index_Valleys(ii))|disp_ip1span>theDisp(index_Peaks(ii))) = [];
    Ceiling = trapz(theDisp(iispan),theLoad(iispan)); %Top integral
    Basement = trapz(disp_ip1span,load_ip1span); %Bottom integral
    theAreas(ii) = Ceiling-Basement; %Area between.
    
    %Visualization
    plot(theDisp,theLoad,'b-')
    hold on
    plot(theDisp(index_Valleys(ii)),theLoad(index_Valleys(ii)),'gp');
    plot(theDisp(index_Peaks(ii)),theLoad(index_Peaks(ii)),'rp');
    fill(sort([theDisp(iispan) disp_ip1span]),sort([theLoad(iispan) load_ip1span]),'m');
    hold off
    pause(.5)
end
%SCd 11/22/2010
%%%%%

Thank you very much Sean for this code and your comments. I got it to work. I am lucky to receive advice from someone doing similar tests.

May I ask some related questions?
a. I hadn’t considered filtering these data before your recommended it (I am fairly new to biomech testing as well as MATLAB). Is it always good practise to filter and remove NaNs before further analysis of this type (i.e. data collecated at a similar sampling rate and with similar noise)? For example, I also test my samples to failure and measure max force and area under their respective force-displacement plots.
b. I see that the variable ‘Areas’ contains 28 values which would seem to be the correct number given that I have 29 peaks and 29 valleys. However, the plotting of the pink filled in areas seems to stop before it reaches the end, i.e some loops do not appear to be fully ‘filled in’. Is this as it should be? I attach a file here: http://www.4shared.com/file/FkLBDv2q/26_Nov_2010_Final_fig_in_hyste.html
c. Many authors, it seems to me, express hysteresis as the ratio of the area within the loop to the area below the loading (ascending) curve only. Could you guide me on how to calculate the area under the loading curve too please? I assume that this calculation will also generate a list of 28 values which I can use with the 28 Area values which I already have.

Less important questions
d. Although the code seems to work fine, I can’t open ‘theDisp’ or ‘theLoad’ variables – MATLAB repeatedly gives a black window and then crashes. Is this unimportant?
e. If I zoom in a lot on the plotted filtered data, at each peak and valley I can still see some ‘bubbles’ however, is it the case that by filtering much more I may lose important data?

Thank you.

Etienne

"Etienne O'Brien" <etienneobrien@yahoo.com> wrote in message <icp7fr$g2e$1@fred.mathworks.com>...
> Thank you very much Sean for this code and your comments. I got it to work. I am lucky to receive advice from someone doing similar tests.
>
> May I ask some related questions?
> a. I hadn’t considered filtering these data before your recommended it (I am fairly new to biomech testing as well as MATLAB). Is it always good practise to filter and remove NaNs before further analysis of this type (i.e. data collecated at a similar sampling rate and with similar noise)? For example, I also test my samples to failure and measure max force and area under their respective force-displacement plots.
> b. I see that the variable ‘Areas’ contains 28 values which would seem to be the correct number given that I have 29 peaks and 29 valleys. However, the plotting of the pink filled in areas seems to stop before it reaches the end, i.e some loops do not appear to be fully ‘filled in’. Is this as it should be? I attach a file here: http://www.4shared.com/file/FkLBDv2q/26_Nov_2010_Final_fig_in_hyste.html
> c. Many authors, it seems to me, express hysteresis as the ratio of the area within the loop to the area below the loading (ascending) curve only. Could you guide me on how to calculate the area under the loading curve too please? I assume that this calculation will also generate a list of 28 values which I can use with the 28 Area values which I already have.
>
> Less important questions
> d. Although the code seems to work fine, I can’t open ‘theDisp’ or ‘theLoad’ variables – MATLAB repeatedly gives a black window and then crashes. Is this unimportant?
> e. If I zoom in a lot on the plotted filtered data, at each peak and valley I can still see some ‘bubbles’ however, is it the case that by filtering much more I may lose important data?
>
> Thank you.
>
> Etienne

Hi Etienne,

Sorry it's taken so long to get back to you, Thanksgiving vacation here.

a. To decide if you need to filter perhaps you should compare unfiltered results to filtered results and see if there is any significant difference? You could in many different ways, and that's its own topic, one in which I'm not particularly educated. However, having too long a span will round your peaks which is bad.
Instead of removing the nan values directly like I did, you could interpolate over them using interp1 or similar. I don't know where they actually were, the code just threw an error the first time and for demonstration ease I got rid of them. As far as noise removal, I'm sure that you could measure the signal to noise ratio of your sensors and design a filter around that. There are others on here who would offer better advice.

b. There is a bug in the area/filling code. The first valley is a greater index than the first peak, and thus the first peak is unnecessary and useless. That is corrected in the the code below. If you're going to generalize this then you should check to see which is greater index_Peaks(1) or index_Valleys(1). If index_Peaks(1)>index_Valleys(1) then do what we had before, else do what's below.

c. Area under the loading curve is the variable 'Ceiling' in my code. I first calculated the entire area under the loading curve (trapezoidal integral), and then subtracted the area under the unloading curve. You could store the ceilings if you wanted.
>> Ceiling(ii) = ...

d. This is expected, the variables are vectors over a million elements long! If you want to see what's inside them, extract pieces and view those, it would take a VERY long time to inspect a million element long vector.

e. Yes, as I alluded to earlier you don't want to lose peaks/valleys. I would recommend posting that as it's own thread with a title: "How to filter/smooth this" or similar.

Good Luck!

%%%
%% Engine & Visualization
theAreas = zeros(length(index_Valleys)-1,1); %preallocate
figure;
hold off
for ii = 1:length(theAreas)
    %Engine
    iispan = index_Valleys(ii):index_Peaks(ii+1); %span of data between current valley and peak
    ip1span = index_Peaks(ii+1):index_Valleys(ii+1); %span of data between current peak and next valley
    disp_ip1span = theDisp(ip1span); %extracted next part (need to extract to analyze)
    load_ip1span = theLoad(ip1span);
        %Remove parts of the next span that fall outside the range of the current span
        %I.e. remove data that's not below the "ceiling" of the hysteresis loop
        %You may not actually wish to do this. If you don't comment or delete these two lines.
    load_ip1span(disp_ip1span<theDisp(index_Valleys(ii))|disp_ip1span>theDisp(index_Peaks(ii))) = [];
    disp_ip1span(disp_ip1span<theDisp(index_Valleys(ii))|disp_ip1span>theDisp(index_Peaks(ii))) = [];
    Ceiling = trapz(theDisp(iispan),theLoad(iispan)); %Top integral
    Basement = trapz(disp_ip1span,load_ip1span); %Bottom integral
    theAreas(ii) = Ceiling-Basement; %Area between.
    
    %Visualization
    plot(theDisp,theLoad,'b-')
    hold on
    plot(theDisp(index_Valleys(ii)),theLoad(index_Valleys(ii)),'gp');
    plot(theDisp(index_Peaks(ii)),theLoad(index_Peaks(ii)),'rp');
    fill([theDisp(iispan) disp_ip1span],[theLoad(iispan) load_ip1span],'m');
    hold off
    pause(.2)
end
%SCd
%%%

Thanks Sean,

I'm a little stuck with my experimental conditions (quite separate from MATLAB) so haven't been able to try your solution properly. I hope to sort this out and then say thanks finally or ask again.

I would like to say I found this code helpful.
Though after four months, I found a bug that's been messing around with my data.
The line that reads:

 theAreas(ii) = Ceiling-Basement; %Area between.

Should read Ceiling + Basement;

Computing the riemann sum for the bottom half of the loop will yield a negative value as the trapz function is going backwards in indices.

I believe this is the case, as I was trying to divide the average area by the area of the rectangle of maxs and mins, and got a value greater than 1.

Andrew

"Andrew Ho" <andrew_ho@hmc.edu> wrote in message <iv0082$c3l$1@newscl01ah.mathworks.com>...
> I would like to say I found this code helpful.
> Though after four months, I found a bug that's been messing around with my data.
> The line that reads:
>
> theAreas(ii) = Ceiling-Basement; %Area between.
>
> Should read Ceiling + Basement;
>
> Computing the riemann sum for the bottom half of the loop will yield a negative value as the trapz function is going backwards in indices.
>
> I believe this is the case, as I was trying to divide the average area by the area of the rectangle of maxs and mins, and got a value greater than 1.
>
> Andrew


Good catch Andrew; you are correct. On that note: you could do the ascending and descending trapezoidal integral of each loop all at once. I.e. index_valleys(ii):index_valleys(ii+1)
or peaks; your choice!

Hi Sean,

May I ask for your help again? It's only now (18 months later) that I've been able to try out your code on experimental data (as opposed to data from sham samples) and I must admit that I still don't understand its inner workings.

My specific queries are:
(1) Am I using the code correctly/are my results sensible?
(2) I added correction from Andrew Ho but that gave me negative areas so I removed it. I take it that this correction wasn't appropriate for my data?
(3) When I try to divide the area of the loops by area under the loading (ascending) curve I don't get a meaningful result.

I've uploaded some example images/data here:
http://www.4shared.com/photo/FY40da2R/Etienne_area_per_cycle.html
http://www.4shared.com/file/kBglDNBI/Etienne_data.html
http://www.4shared.com/photo/fOsGI6Gt/Etienne_filled_areas.html

The code as I'm using it now (EOB refers to me) is below.

Thank you.

Etienne

%% Engine & Visualization
theAreas = zeros(length(index_Valleys)-1,1); %preallocate
figure;
hold off
for ii = 1:length(theAreas)
    %Engine
    iispan = index_Valleys(ii+1):index_Peaks(ii+1); %span of data between current valley and peak. EOB added a '+1' (inverted commas for illustration only) to the statement iispan = index_Valleys(ii'+1')
    ip1span = index_Peaks(ii+1):index_Valleys(ii+1); %span of data between current peak and next valley
    disp_ip1span = theDisp2(ip1span); %extracted next part (need to extract to analyze)
    load_ip1span = theLoad2(ip1span);
        %Remove parts of the next span that fall outside the range of the current span
        %I.e. remove data that's not below the "ceiling" of the hysteresis loop
        %You may not actually wish to do this. If you don't comment or delete these two lines.
    load_ip1span(disp_ip1span<theDisp2(index_Valleys(ii))|disp_ip1span>theDisp2(index_Peaks(ii))) = [];
    disp_ip1span(disp_ip1span<theDisp2(index_Valleys(ii))|disp_ip1span>theDisp2(index_Peaks(ii))) = [];
    Ceiling = trapz(theDisp2(iispan),theLoad2(iispan)); %Top integral
    Basement = trapz(disp_ip1span,load_ip1span); %Bottom integral
    theAreas(ii) = Ceiling-Basement; %Area between.
   
    %Visualization
    plot(theDisp2,theLoad2,'b-')
    axis([1.8 2.05 0 5])% EOB addition
    hold on
    plot(theDisp2(index_Valleys(ii)),theLoad2(index_Valleys(ii)),'gp');
    plot(theDisp2(index_Peaks(ii)),theLoad2(index_Peaks(ii)),'rp');
    fill([theDisp2(iispan) disp_ip1span],[theLoad2(iispan) load_ip1span],'m');
    hold off
    pause(.2)
end
%SCd
%%%

%% Plotting Area of hysteresis loops against cycle number (EOB addition)
figure;
plot(theAreas,'-*')
xlabel('Loop number','FontSize',16)
ylabel('Loop area (Nmm)','FontSize',16)
title('Area of hysteresis loops','FontSize',16)
grid on

%% Determining and plotting Hysteresis Index (EOB addition)
HysteresisIndex=(theAreas/Ceiling)*100;
figure;
plot(HysteresisIndex,'-*')
xlabel('Loop number','FontSize',16)
ylabel('Hysteresis index(%)','FontSize',16)
title('Hysteresis index','FontSize',16)
grid on












"Sean de " <sean.dewolski@nospamplease.umit.maine.edu> wrote in message <iv02df$hia$1@newscl01ah.mathworks.com>...
> "Andrew Ho" <andrew_ho@hmc.edu> wrote in message <iv0082$c3l$1@newscl01ah.mathworks.com>...
> > I would like to say I found this code helpful.
> > Though after four months, I found a bug that's been messing around with my data.
> > The line that reads:
> >
> > theAreas(ii) = Ceiling-Basement; %Area between.
> >
> > Should read Ceiling + Basement;
> >
> > Computing the riemann sum for the bottom half of the loop will yield a negative value as the trapz function is going backwards in indices.
> >
> > I believe this is the case, as I was trying to divide the average area by the area of the rectangle of maxs and mins, and got a value greater than 1.
> >
> > Andrew
>
>
> Good catch Andrew; you are correct. On that note: you could do the ascending and descending trapezoidal integral of each loop all at once. I.e. index_valleys(ii):index_valleys(ii+1)
> or peaks; your choice!