MATLAB_Utilities/fit_power_curve.m at master · Campbell-Muscle-Lab/MATLAB_Utilities · GitHub

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
28
29
30
31
32
33
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
function [x0,a,b,r_squared,x_fit,y_fit,max_power,rel_x_at_max_power] = ...
    fit_power_curve(x,y,varargin);
% Fits a curve of the form y = x*b*(((x0+a)/(x+a))-1)

% Defaults
params.x0_guess = [];
params.a_guess = [];
params.b_guess = [];
params.x0_min = [];
params.x0_max = [];
params.x_fit = [];
params.figure_display = 0;


% Update
params=parse_pv_pairs(params,varargin);

% Some defaults
if (isempty(params.x_fit))
    params.x_fit = linspace(min(x),max(x),100);
end

% Deduce some starting values
if (isempty(params.x0_guess))
    params.x0_guess=max(x);
end

if (isempty(params.a_guess))
    params.a_guess = 0.2*max(x);
end

if (isempty(params.b_guess))
    params.b_guess = 1e-3;
end

if (isempty(params.x0_min))
    lower_bounds=[0 min(x)+eps 0];
else
    lower_bounds=[params.x0_min min(x)+eps 0];
end
if (isempty(params.x0_max))
    upper_bounds=Inf*ones(3,1);
else
    upper_bounds=[params.x0_max Inf Inf];
end

% Set p
p = [params.x0_guess params.a_guess params.b_guess];

[p,~,status]=fminsearchbnd(@power_error, ...
    p, ...
    lower_bounds, upper_bounds, ...
    optimset('MaxFunEvals',5000), ...
    x,y,params.figure_display);

% Calculate y_fit
x0=p(1);
a=p(2);
b=p(3);
for i=1:numel(x)
    y_fit(i) = power_value(x(i),x0,a,b);
end
r_squared=calculate_r_squared(y,y_fit);

% Calculate fit curve
x_fit = params.x_fit;
for i=1:length(x_fit)
    y_fit(i)=power_value(x_fit(i),x0,a,b);
end

% Create a function handle and interpolate to get the max power
fh = @(x)negative_power_value(x,x0,a,b);
[force_at_max_power,neg_power]=fminbnd(fh,min(x_fit),max(x_fit));
max_power = -neg_power;
rel_x_at_max_power = force_at_max_power/x0;

end


function error_value = power_error(p,x,y,figure_display)

    for i=1:length(x)
        y_fit(i) = power_value(x(i),p(1),p(2),p(3));
    end
    [r1,c1]=size(y);
    [r2,c2]=size(y_fit);
    if (r1>r2)
        y=y';
    end
    error_value = sum((y-y_fit).^2);

    if (figure_display)
        figure(figure_display);
        clf;
        plot(x,y,'bo');
        hold on;
        x_plot=linspace(min(x),max(x),100);
        for i=1:length(x_plot)
            y_plot(i)=power_value(x_plot(i),p(1),p(2),p(3));
        end
        plot(x_plot,y_plot,'r-');
        drawnow;
    end
end

function y=power_value(x,x0,a,b)
    y = x*b*(((x0+a)/(x+a))-1);
end

function y = negative_power_value(x,x0,a,b);
    y = -power_value(x,x0,a,b);
end