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

采纳的回答

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