help with deleting zeros from an array

조회 수: 6 (최근 30일)
Teshan Rezel
Teshan Rezel 2022년 2월 7일
편집: Mathieu NOE 2022년 2월 8일
Hi folks,
I am trying to use polyfit and polyval to get the line of best fit across two sets of data, but am getting a very wierd graph as follows:
I think this is because the regression is fitting to the zero values on the top left of the graph. Therefore, I am trying to remove the zeros, but with no real success! I have below a snippet of my code for review. Any ideas on how to fix this please?
The graph is myLog (y axis) versus RT (x-axis).
Thanks!
for i = 1:numSheets
test = data{1, i};
DWC{i} = test.("Derivative Weight Change");
Time = test.("Time");
Temperature{i} = test.("Temperature");
Weight = test.("Weight");
[peakDWC,locationDWC,~, ~] = findpeaks(DWC{i}, 'SortStr','descend','NPeaks',1);
peakPositions = ischange(DWC{i}, "linear");
Max_MassLoss = peakDWC;
INDEX_MaxMassLoss = locationDWC;
AVG_MassLoss = mean(DWC{i}, "omitnan");
INDEX_Ignition(i) = 3000+find(DWC{i}(3000:end)>0.1, 1);
Time_MaxMassLoss = Time(INDEX_MaxMassLoss);
Time_Ignition = Time(INDEX_Ignition(i));
Temperature_Ignition(i) = 273+Temperature{i}(INDEX_Ignition(i));
INDEX_burnoutTemp(i) = find(peakPositions, 1, 'last');
burnoutTemp(i) = 273+Temperature{i}(INDEX_burnoutTemp(i));
widthDWC(i) = INDEX_burnoutTemp(i) - INDEX_Ignition(i);
massDiff = Weight(1) - Weight(end);
for j = 1 : widthDWC(i)
alpha(j) = (Weight(1)-Weight(INDEX_Ignition(i)+j))/massDiff;
kelvin(j) = 273+Temperature{i}(INDEX_Ignition(i)+j);
logBottom(j) = kelvin(j).^2;
logTop(j) = log(1-alpha(j));
logTop(isinf(logTop)) = 0;
myLog(j, i) = log(-logTop(j)./logBottom(j));
RT(j, i) = 1/(gasConstant*kelvin(j));
end
myLog(myLog==0) = [];
RT(RT==0) = [];
mySize = size(RT(:, i), 1);
coefficients = polyfit(RT(:, i), myLog(:, i), 1);
xFit{i} = linspace(min(RT(:, i)), max(RT(:, i)), mySize);
yFit{i} = polyval(coefficients , xFit{i});
[fitresult, gof] = fit(RT(:, i), myLog(:, i), 'poly1', 'Normalize', 'on');
intercept(i) = fitresult.p2;
slope(i) = fitresult.p1;
R_Squared(i) = gof.rsquare;
Energy(i) = -slope(i)*100;
K0_first = exp(intercept(i));
K0(i) = Energy(i)*K0_first/gasConstant;
end

채택된 답변

Mathieu NOE
Mathieu NOE 2022년 2월 7일
hello
could not really figure out what in your plot is weird but at least if you really want to remove zero or almost zero values , this is how I would do it : do a conditionnal statement with < threshold , and not == 0 because this is not very robust due to floating point rounding error
ind = (myLog<eps); % eps is a very small positive value
myLog(ind) = [];
RT(ind) = [];
  댓글 수: 6
Mathieu NOE
Mathieu NOE 2022년 2월 8일
oops, sorry
my last suggestion work only if the data is a vector not a 2D array.
By the way I noticed that among the 56 data sets, they have all different amount of zeros.
I still don't know exactly what is the goal of the fit ? does it have to go through the origine (no constant term ) or do you want with a non zero constant term ?
why the need to remove the zero data in the fit ?
what are the data you want to store ? the fitted data ?
Mathieu NOE
Mathieu NOE 2022년 2월 8일
편집: Mathieu NOE 2022년 2월 8일
i tried another approach
I guessed that what was the issue is maybe the "weird" looking straigth portion of the curve that goes from the origin (0,0) to the first experimental point
I fitted then the experimental data (without the zeros) to a 1st order polynomial with polyfit0 , so the constant term of the polynomial is forced to zero here (again, use regular polyfit otherwise) .
as the xfit vector is the same for all data sets (from zero to max value in all array RT) , its fairly easy to concatenate the fitted y data (if you wish)
clc
clearvars
load myLog.mat
load RT.mat
maxx = max(RT,[],'all');
xfit = (linspace(0,maxx,100))';
for ci = 1:56
x = RT(:,ci);
y = myLog(:,ci);
% remove data below 1e-6 threshold (RT)
ind = (x<1e-6);
x(ind) = [];
y(ind) = [];
% fit 1st degree polynomial fit with zero forcing (the curve must go through
% origin (0,0)
[slope,~] = polyfit0(x,y,1);
yfit(:,ci) = slope*xfit;
figure(ci);
plot(x,y,xfit,yfit(:,ci));
end

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Interpolation에 대해 자세히 알아보기

제품


릴리스

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by