# Using loop inside loop curve fit and extract x intercept_Data Attached

21 views (last 30 days)
Darpan Verma on 11 Mar 2019
Commented: Darpan Verma on 11 Mar 2019
Hi I would appreciate any helps on code for building a logic for this problem. I am trying to write a code which will only scan data between the two lines as shown in figure below 4.3 to 5.1. and do the curve fit (only for the left side portion of curve) (linear portion) of all these graphs and give me its x intercepts. Basically summarizing properly:
1. Take curve data values between 4.3>x>5.1
2. Do the curve fitting (only for the left side linear portion of the curve: the curve fitting is a straight line) for all the graphs.
3. Give me the values of x intercept for all the fitted lines.
I HAVE ATTACHED THE DATA AS WELL FOR YOUR REFERENCE I have my code like this:
while x < 5.1 | x > 4.3 %Take curve between 4.3>x>5.1
%Curve fitting for all the data
for i = 1:1:42
[fitresult, gof] = createFit(x, y) % Call the fitting function
coeffs = coeffvalues(fitresult);
m = coeffs(1);
b = coeffs(2);
Xintercept = -b/m;
end
end
FITRESULT FUNCTION CREATED:
function [fitresult, gof] = createFit(x, y)
%CREATEFIT(X,Y)
% Create a fit.
%
% Data for 'untitled fit 1' fit:
% X Input : x
% Y Output: y
% Output:
% fitresult : a fit object representing the fit.
% gof : structure with goodness-of fit info.
%
% Auto-generated by MATLAB on 06-Mar-2019 21:38:52
%% Fit: 'untitled fit 1'.
[xData, yData] = prepareCurveData( x, y );
% Set up fittype and options.
ft = fittype( 'poly1' );
excludedPoints = excludedata( xData, yData, 'Indices', [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72] );
opts = fitoptions( 'Method', 'LinearLeastSquares' );
opts.Exclude = excludedPoints;
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
% Plot fit with data.
figure( 'Name', 'untitled fit 1' );
h = plot( fitresult, xData, yData, excludedPoints );
legend( h, 'y vs. x', 'Excluded y vs. x', 'untitled fit 1', 'Location', 'NorthEast' );
% Label axes
xlabel x
ylabel y
grid on
I am not able to device an idea for doing this.
Any help would be appreciated.
Darpan
Darpan Verma on 11 Mar 2019
Actually Walter this condition which I wrote in my code for a while loop might not be right.
Actually I what I am doing is doing curve fitting for the linear portion of these curves. Since I can't do it one by one as there would be whole lot of data, I am working on a code which will do that for me and give the x intercept value for all the fitted curves.

Raghunandan V on 11 Mar 2019
Edited: Raghunandan V on 11 Mar 2019
I have done it using a for loop. You can also try without for loop too. Please do optimise the code with functions for your particular application.
clear all;
clc
% find the indices between 4.5 qand 5
Index1 = (Num>4.5);
Index = (Num<5.1);
Index = Index1==1 & Index==1;
Index = Index(1,:);
X = Num(1, Index);
Y = Num(3, Index);
Correlation = 0;
%initialize slope =0
m=0; %slope
for i = 3: length(X)
%Calculate slope between two adjacent points
m(i) = (Y(i) - Y(i-1))/(X(i) - X(i-1));
% calculate correlation only if slope is positive beacuse you have only
% asked for left side of the graph
if (m(i) > 0)
%initialize variables for calculate=ing correlation between last 3
%points
C = [Y(i-2),Y(i-1), Y(i)];
B = [X(i-2),X(i-1), X(i)];
Product = C .* B;
% calculate correlation
Numerator = (length(C)* sum(Product) - sum(C) * sum(B));
Denominator = sqrt(((length(C) * sum (C.^2)) - ((sum(C))^2)) * ((length(B) * sum (B.^2)) - ((sum(B))^2)));
Correlation(i) = Numerator / Denominator;
end
end
[MaxCorr, Ind] = max(Correlation)
display('The best correlation points are')
display(Y(Ind: Ind+2) )
This code only provides you the points where correlation is highest between 3 points. If required the number of points can be made dynamic too. Please continue the code to draw the line and find X intercept. To draw the line follow this link . To know about correlation follow this link.
Please let me know if you do not understand the code.
Darpan Verma on 11 Mar 2019
Thanks a lot Raghunandan. I ll try this and see if its working.

R2018b

### Community Treasure Hunt

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

Start Hunting!