Removing noise from Hysteresis loops

Hi, I need a help in removing noise from Hysteresis loops, I have generated code which loads multiple files(containing Force and Displacement data) one by one and plots their Hysteresis loops but i need to remove noise and the smoothening curve should fit the orginal shape, I have used a filter in following code but it distorts the orginal shape of loops which effects the area under the curve for calculation of damping.This filter works well for one data file but disfigures the other file, so I want a solution that would work for my all files. Any help would be highly appreciated.

Following is my code:

clc
clear all
projectdir = 'F:\NUST\SAMPLE 9\Measured'
dinfo = dir(fullfile(projectdir))
dinfo([dinfo.isdir]) = []    %get rid of all directories including . and ..
%nfiles = length(dinfo)
folder = 'F:\NUST\IMAGE\hysteresis';
for j = 5 :6
    
  filename = fullfile(projectdir, dinfo(j).name)
  delimiter = ';';
startRow = 1000;
endRow = 2000;
formatSpec = '%*s%*s%f%f%f%f%f%f%f%f%f%[^\n\r]';
  f1 = fopen(filename, 'r')
  dataArray = textscan(f1, formatSpec, endRow-startRow+1, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines', startRow-1, 'ReturnOnError', false);
  fclose(f1)
  %% Allocate imported array to column variable names
  Displacement = dataArray{:, 1}
Force = (dataArray{:, 2})*1000

%%%%%%%%%%%%%Plots%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure (j)
%plot(Testtime,Force)
%b=fir1(20,0.1,'low');
%[b,a]=butter(10,0.2,'low');
[b,a]=butter(10,0.2,'low');
%y1=filter(b,1,Force);%Fir
y2=filter(b,a,Force);%butter
plot(Displacement,y2)
%plot(Displacement,Force)
title(dinfo(j).name)
%title(['S1  R',sprintf('%d',j)])
xlabel('Displacement(mm)') 
ylabel('Force(N)')
%saveas(h,sprintf('FIG%d.png',j));
hold on
%saveas(plot, fullfile(folder, sprintf('fig%02d.jpg', j)));

pngFileName = sprintf('Fig_%d.jpg', j);
  fullFileName = fullfile(folder, pngFileName);
    
  % Then save it
  %saveas(gcf, fullFileName)
  


end

Answer 1

[Star Strider]

The easiest way to remove the noise is likely to fit your data to a function that describes the hysteresis.

See: Interpolation when y data is not strictly a function of x for an illustration of how to do that.

[Sadia]

Thank you Sir for your valuable replay. I have modified my code using your illustration but ascending and descending curves are very narrow, seems like they are not following the data accurately. i am attaching figures of both usng my old code and new one. and following is my modified code. kindly assesist me if I am doing something wrong.

clc
clear all
projectdir = 'F:\NUST\SAMPLE 9\Measured'
dinfo = dir(fullfile(projectdir))
dinfo([dinfo.isdir]) = []    %get rid of all directories including . and ..
%nfiles = length(dinfo)
folder = 'F:\NUST\IMAGE\hysteresis';
for j = 5 :6
    
  filename = fullfile(projectdir, dinfo(j).name)
  delimiter = ';';
startRow = 1000;
endRow = 2000;
formatSpec = '%*s%*s%f%f%f%f%f%f%f%f%f%[^\n\r]';
  f1 = fopen(filename, 'r')
  dataArray = textscan(f1, formatSpec, endRow-startRow+1, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines', startRow-1, 'ReturnOnError', false);
  fclose(f1)
  Displacement = dataArray{:, 1}
  Force = (dataArray{:, 2})*1000
x = Displacement;
y = Force;
[cmin,cxix] = (min(x))
[cmax,cnix] = (max(x))
[fmin,fxix] = (min(y))
[fmax,fnix] = (max(y))
figure(j)
plot(x, y)
hold on
plot([cmin  cmax],  [fmin  fmax], '+r', 'MarkerSize', 5, 'LineWidth',1.5)
hold off
grid
text(cmin,fmin, sprintf('(%.3f, %.3f)',cmin, fmin), 'FontSize',8, 'FontName', 'Consolas', 'FontWeight', 'bold', 'VerticalAlignment','bottom', 'HorizontalAlignment','left')
text(cmax,fmax, sprintf('(%.3f, %.3f)',cmax, fmax), 'FontSize',8, 'FontName', 'Consolas', 'FontWeight', 'bold', 'VerticalAlignment','top', 'HorizontalAlignment','right')
xlabel('Displacement(mm)')
ylabel('Force(kN)')
Phi1 = @(b,I) -(b(1) .* atan(-b(2)*I+b(3)) - b(1).*I.*b(4));    % Descending
Phi2 = @(b,I)  (b(1) .* atan(-b(2)*I+b(3)) + b(1).*I.*b(4));    % Ascending
ixd = fix(length(x)/2);
vi2 = 1:ixd;
vi1 = ixd:length(x);
opts = statset('MaxIter', 5000, 'MaxFunEvals', 10000);
B1 = nlinfit(x(vi1), y(vi1), Phi1, ones(4,1),  opts )
B2 = nlinfit(x(vi2), y(vi2), Phi2, ones(4,1),  opts )
FitPhi1 = Phi1(B1,x);
FitPhi2 = Phi2(B2,x);
figure(j)
hpd1 = plot(x(vi1), y(vi1), '.b', x(vi2), y(vi2), '.b', 'LineWidth',1);
hold on
hpr1 = plot(x, FitPhi1, '-r', 'LineWidth', 1.5);
hpr2 = plot(x, FitPhi2, '-g', 'LineWidth', 1.5);
hold off
grid
legend([hpd1(1),hpr1,hpr2], 'Data', 'Descending', 'Ascending', 'Location', 'NW')
xlabel('Displacement(mm)')
ylabel('Force(kN)')
title(dinfo(j).name)


  


end

[Sadia]

I want these ascending and descending lines to fit the shape of loop and each should pass from middle of these data points. but as it can be seen from grpah its not the case. Need help plz

[Sadia]

Also in your code ascending and descending are not appropriately mentioned, upper curve should be ascending and lower should be descending, if i am not wrong? well i can fix that. I just want to be sure if I am right or not?

[Star Strider]

As always, my pleasure.

The ascending and descending limbs of the hysteresis loop are correct in my code. (I have verified this with previous data.) The ascending and descending limbs may be different (to the left or right of each other, rather than one consistently being on the left and the other on the right) for each data set. My code cannot control for this.

Your data in ‘fig 5.PNG’ appear to have several ‘force’ values for each ‘displacement’ value, so it may be best to average the ‘force’ values with respect to ‘displacement’. I am not certain how my code would handle very noisy data, except to give what it considers to be the best fit to the data it gets. Without a .mat file of your data (with an explanaztion of what the values are in it), this is the best I can do.

You most likely need to identify the ascending and descending limbs in any event, even with oisy data.

[Sadia]

I am attaching .mat file and excel file for my data set.

[Sadia]

In attachement of excel file only Force and Displacements columns are used in code, meanwhile I am trying to take average of force w.r.t Displacement, hope it may help. I highly appreciate your timely response, working on these hysteresis has consumed my time alot. I hope this time i'll get through it.

[Sadia]

I am also attaching my code where I have used a filter which some how follows the trend of loops but is only good for one data set and when I apply it to other data files it distorts them.But I prefer to follow your approach as it will work fine for all files. Following is the image of my filtered data.

[Star Strider]

Your data are definitely not easy to work with!

The data in your file compriseds one loop. There are many NaN values, so to eliminate them, I used:

dataArray = readmatrix('forceDisplacement.xlsx');
DispFrc = dataArray(:,3:4);
DispFrc = DispFrc(~any(isnan(DispFrc),2),:);
DispFrc = sortrows
Displacement = DispFrc(:,1);
Force = DispFrc(:,2);
x = Displacement;
y = Force;

It turns out that they are not neatly apportioned between the ascending and descending loops either, so simply dividing your matrix in half (that worked in applying my code to other problems) does not work with your data.

I instead used this to separate them:

ixd = gradient(x);
vi2 = ixd >  0;
vi1 = ixd <= 0;

and it appeared to work because the data you provided in your Excel file otherwise resulted in both limbs of they hysteresis loop very closely overlying each other. I could not distingush them using my existing code. I could only distinguish them using my revised ‘vi1’ and ‘vi2’ vectors based on the positive and negative gradients, and this was not optimal because the loops are supposed to approximate each other at thier extremes. Your data do not, and remain separated.

You applied my code correctly, however your data do not appear to have the characteristics my code requires to work easily with it. I am not certain what the problem is beyond that, or how to solve it.

I will do my best to work with you to solve these problems, if a solution exists.

[Sadia]

I tried to use your edited code but it gives error: Undefined function 'isnan' for input arguments of type 'cell'.

[Sadia]

Kindly see the attachment in which I have found a code for Hysteresis but I dont know how I can I process it for my code. I have also seen a discussion where I have found following comment:

"For the separation I recommend writing a small Matlab program that follows the field column and breaks up the data to sections whenever a change in direction occurs. For example: Loop through n steps and at step n check

if [field(n)-field(n-1)]*[field(n-1)-field(n-2)]<0 then ....."

May be this might give some clue to help me write some advance code. If you can help me with generating a new code?

[Star Strider]

With respect to your earlier Comment, I use the readmatrix function (R2019a and later versions). It produces a double array, not a cell array. You posted an .xlsx spreadsheet, so the easiest way to import it is with readmatrix or xlsread.

Your ‘Hysteresis1.m’ file apparently creates a hysteresis loop and then attempts to analyse it. It requires the Communications Toolbox, that I do not have. If I substitute a randn call for awgn, I find that it simulates a hysteresis loop, then calls a function called ‘Hysteresis’ that is not supplied, so I cannot run it. I was also unable to find it with an Interweb search.

I discovered that sgolayfilt can smooth your data, however you have a number of data that are not unique, and this is causing problems with the analysis. Also, you have a number of cycles in your data, not just one, and those have varying amplitudes. My code was designed to analyse one complete cycle, and for that it works correctly.

I ended up getting unique values for your data, then filtering them, then doing the analysis. I use findpeaks to isolate the last (largest) cycle, provining the most data. This is about as good as I can get your data:

dataArray = readmatrix('forceDisplacement.xlsx');
DispFrc = dataArray(:,3:4);
DispFrc = DispFrc(~any(isnan(DispFrc),2),:);
Displacement = DispFrc(:,1);
Force = DispFrc(:,2);
x = Displacement;
y = Force;

[UDispFrc,ridx] = unique(DispFrc, 'rows', 'stable');
Ux = UDispFrc(:,1);
Uy = UDispFrc(:,2);
xf = sgolayfilt(Ux, 3, 7);
yf = sgolayfilt(Uy, 3, 7);
[xpks, xlocs] = findpeaks(xf);
[ypks, ylocs] = findpeaks(yf);
xfr = xf(xlocs(end-1):xlocs(end));
yfr = yf(xlocs(end-1):xlocs(end));

Phi1 = @(b,I) -(b(1) .* atan(-b(2)*I+b(3)) - b(1).*I.*b(4));    % Descending
Phi2 = @(b,I)  (b(1) .* atan(-b(2)*I+b(3)) + b(1).*I.*b(4));    % Ascending

ixd = fix(length(xfr)/2);
vi2 = 1:ixd;
vi1 = ixd:length(xfr);

opts = statset('MaxIter', 5000, 'MaxFunEvals', 10000);
B1 = nlinfit(xfr(vi1), yfr(vi1), Phi1, ones(4,1),  opts )
B2 = nlinfit(xfr(vi2), yfr(vi2), Phi2, ones(4,1),  opts )
FitPhi1 = Phi1(B1,xfr(vi1));
FitPhi2 = Phi2(B2,xfr(vi2));

figure
hpd1 = plot(xfr(vi1), yfr(vi1), '.r', xfr(vi2), yfr(vi2), '.g', 'LineWidth',1);
hold on
hpr1 = plot(xfr(vi1), FitPhi1, '-r', 'LineWidth', 1.5);
hpr2 = plot(xfr(vi2), FitPhi2, '-g', 'LineWidth', 1.5);
hold off
grid
legend([hpd1(1),hpr1,hpr2], 'Data', 'Descending', 'Ascending', 'Location', 'NW')
xlabel('x_1')
ylabel('y_1')

That is the best I can do.

To see the problems with your data in more detail, add ‘zf’ and this plot:

zf = linspace(0, 1, numel(Ux));

figure
plot3(xf, yf, zf)
grid
xlabel('xf')
ylabel('yf')

The inconsistent data are causing serious problems with the analysis of your data in my code.

[Sadia]

Thank you so much sir for your precious time, one last request as I have multiple files to process with input argument of type "cell" not double array (R2015a), so could you please adjust my following code so that It may not give error. I have adjusted end and start rows so that i can select few cycles, you can change it. Pardon me for my lack of knoweledge in coding as that’s not within the area of my expertise. I highly appreciate your help.

projectdir = 'F:\NUST\SAMPLE 9\Measured'
dinfo = dir(fullfile(projectdir))
dinfo([dinfo.isdir]) = []    %get rid of all directories including . and ..
for j = 11:13  
  filename = fullfile(projectdir, dinfo(j).name)
  delimiter = ';';
startRow = 1000;%to get few cycles as their are 12 cycles rest are noise.
endRow = 2000;
formatSpec = '%*s%*s%f%f%f%f%f%f%f%f%f%[^\n\r]';
  f1 = fopen(filename, 'r')
  dataArray = textscan(f1, formatSpec, endRow-startRow+1, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines', startRow-1, 'ReturnOnError', false);
  fclose(f1)
  %% Allocate imported array to column variable names
 DispFrc = dataArray(:,3:4);
DispFrc = DispFrc(~any(isnan(DispFrc),2),:);
Displacement = DispFrc(:,1);
Force = DispFrc(:,2);
x = Displacement;
y = Force;
.

[Star Strider]

As always, my pleasure.

If you want to convert the cell array that textscan returns to a double array, slightly change the name of the textscan output:

dataArrayc = textscan(f1, formatSpec, endRow-startRow+1, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines', startRow-1, 'ReturnOnError', false); 	   

then:

dataArray = cell2mat(dataArrayc);

Your ‘dataArray’ variable will now be a double array.

I have no idea how to adjust your other files, unless they have essentially the same structure as the one you provided. If they do, my code will likely work, although perhaps with some minor modifications depending on the data you have.

[Sadia]

Thank you again sir, it worked.

[Star Strider]

As always, my pleasure.

[Sadia]

I have found a very easy filter for Hysteresis its as following and it works well for any shape of hysteresis loops.

sm=smooth(F,30,'loess')

[Star Strider]

Not everyone has the Curve Fitting Toolbox. (I don’t.)

[Sadia]

hmmm but R2015a version of Matlab support this command of smoothening the loops. Any way it benifited me.

[Star Strider]

MATLAB has the smooth3 function (before R2006a), and the smoothdata function was introduced in R2017a. The smooth function is only in the Curve Fitting Toolbox.

I still prefer the method I use here for hysteresis curves.

[Sadia]

Sure sir I too was using your method till I found this new one. Thank you again.

[Sadia]

see I am so happy with the outline of this hysteresis, its perfectly aligning with my data. :)

[Star Strider]

As always, my pleasure.

The nonlinear regression not only gives a smooth result, it also gives the parameters of the ascending and descending parts of the hysteresis loop.

[Sadia]

You are absolutely rights sir.

[Star Strider]

Thank you.

[Sadia]

My pleasure.