You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Removing noise from Hysteresis loops
10 views (last 30 days)
Show older comments
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
Accepted Answer
Star Strider
on 6 Jul 2019
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.
25 Comments
Sadia
on 6 Jul 2019
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
on 6 Jul 2019
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
on 6 Jul 2019
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
on 6 Jul 2019
Edited: Sadia
on 6 Jul 2019
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
on 6 Jul 2019
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
on 6 Jul 2019
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
on 7 Jul 2019
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
on 7 Jul 2019
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
on 7 Jul 2019
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
on 7 Jul 2019
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.
Star Strider
on 8 Jul 2019
As always, my pleasure.
Star Strider
on 19 Sep 2019
Not everyone has the Curve Fitting Toolbox. (I don’t.)
Star Strider
on 21 Sep 2019
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.
Star Strider
on 21 Sep 2019
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.
Star Strider
on 23 Sep 2019
Thank you.
More Answers (0)
See Also
Categories
Find more on Visual Exploration in Help Center and File Exchange
Tags
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)