help with deleting zeros from an array
5 次查看(过去 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!