help with deleting zeros from an array
조회 수: 3 (최근 30일)
이전 댓글 표시
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
댓글 수: 0
채택된 답변
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
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
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 Center 및 File Exchange에서 Interpolation에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!