Can't fit nonlinear filter equation by lsqcurvefit()

I want to fit a passive RLC bandstop filter equation to the following data:

f = [5500 6000 6500 7000 7500 8000 8500 9000];
y = [0.998 0.997 0.994 0.988 0.947 0.255 0.956 0.987];

Where f is the frequency and y is the value of transfer function. Here is my code:

fun = @(p, x) abs(j.*(x.*p(2)-1./(x.*p(3)))./(p(1) + j.*(x.*p(2)-1./(x.*p(3))))); % transfer function of RLC bandstop
% x, p(1), p(2), p(3) correspond to w, R, L, C respectively
x = f(1):f(end); % for plotting purpose
x0 = [1 1e-3 1e-6] % starting points
lb = [1 1e-3 1e-6];ub = [100E3 100 1000E-6];
p_fit = lsqcurvefit(fun, x0, f, y, lb, ub) % p_fit is the vector of fitted parameters
semilogx(f, y, 'ro'); hold on; 
semilogx(x, fun(p_fit, x), '--m'); grid on; % semilog plot of the transfer function
hold off;

I am not familiar with lsqcurvefit() parameter's starting points. Please, suggest any modification or any other method which could fit for any frequency range.

I have attached the image of actual plot and the transfer function.

Answer 1

[Walter Roberson]

format long g
f = [5500 6000 6500 7000 7500 8000 8500 9000];
y = [0.998 0.997 0.994 0.988 0.947 0.255 0.956 0.987];

fun = @(p, x) abs(1i.*(x.*p(2)-1./(x.*p(3)))./(p(1) + 1i.*(x.*p(2)-1./(x.*p(3))))); % transfer function of RLC bandstop
% x, p(1), p(2), p(3) correspond to w, R, L, C respectively
x = f(1):f(end); % for plotting purpose
x0 = [1 1e-3 1e-6] % starting points

x0 = 1×3

1 0.001 1e-06

lb = [1 1e-3 1e-6];ub = [100E3 100 1000E-6];
ub = [1 0.02 2e-5];
x0 = [1 0.0038733109633344873895993294337111 4e-6]

x0 = 1×3

1 0.00387331096333449 4e-06

[p_fit, fval_fit] = lsqcurvefit(fun, x0, f, y, lb, ub) % p_fit is the vector of fitted parameters

Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 2.000000e+02.

p_fit = 1×3

1 0.00342918710933376 4.51291192153049e-06

fval_fit = 
      0.000466077848839323

semilogx(f, y, 'ro', 'displayname', 'original'); hold on; 
semilogx(x, fun(p_fit, x), '--m', 'displayname', 'lsq'); grid on; % semilog plot of the transfer function

E = @(p) sum((fun(p,f)-y).^2);
[Pfc, fval_fc] = fmincon(E, x0, [], [], [], [], lb, ub)

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

Pfc = 1×3

1 0.00387331096333449 4e-06

fval_fc = 
      0.000804295477866486

[Pfs, fval_fs] = fminsearch(E, x0)

Pfs = 1×3

1.20098480142405 0.00390552292365765 3.96072393734393e-06

fval_fs = 
       0.00041975861436785

semilogx(x, fun(Pfs,x), '--g', 'displayname', 'fminsearch');
semilogx(x, fun(Pfc,x), ':b', 'displayname', 'fmincon');
hold off

legend show

[Walter Roberson]

This involved calculus and plotting, There is a solution that is nearby, with slightly higher p(2). If you calculated an error function

E = sum((fun(P,f)-y).^2)

and use 1 for P(1) and 4e-6 for P(3) and then plot for p(2) near [3.8e-3 3.9e-3] then you will see something like a w that is significantly lower than the nearby values. Take the derivative of the error function substituting in those P(1) and P(3) values, and you can vpasolve() with a starting range to get a specific value for the P(2) values that minimize. After that, you can use those as starting points to get decent results.

[Md Nure Alam Dipu]

Thank you for the answer. By the way, how did you get the assumptions of p(1), p(3) and the range of p(2)?

And have a look at this data:

f = [80 190 216 224 234 243 254 314]*1e3;
y = [999 984 897 679 205 801 938 993]*1e-3;

I can't fit this now. I have other filter equations which I want to fit. And it's not possible to manually change the starting points as the range of 'f' changes. Can you suggest any way how the program can automatically generate/guess starting points so that the curve fits the data pretty well for any range of 'f'?

[Walter Roberson]

The structure of your function is such that it becomes a constant when p(1) is 0; when you look at the slope with respect to p(1) the optimal will be at the lower bound of p(1) which is 1. That is how you get the 1 as the starting point.

Your function acts nearly constant (first 7 decimal places the same in the error) if p(3) is not pretty small compared to its permitted range.

When p(3) is pretty small your function is very steep, and the exact value of p(3) does not seem to matter much if it is sufficiently small, with variability due to p(2) causing the overwhelming majority of change. So pick a "small" p(3).

The difficulty after that is finding a good p(2). I made some fortunate guesses as to which range to look at, and found the w shape that was significantly better than the surrounding locations. Finding a useful p(2) automatically might be difficult.

[Md Nure Alam Dipu]

Thanks for the information. Could you please share the code of your following statement:

" use 1 for P(1) and 4e-6 for P(3) and then plot for p(2) near [3.8e-3 3.9e-3] then you will see something like a w that is significantly lower than the nearby values. Take the derivative of the error function substituting in those P(1) and P(3) values, and you can vpasolve() with a starting range to get a specific value for the P(2) values that minimize. "

I'm sorry that I couldn't figure that out.

[Walter Roberson]

format long g
f = [5500 6000 6500 7000 7500 8000 8500 9000];
y = [0.998 0.997 0.994 0.988 0.947 0.255 0.956 0.987];

fun = @(p, x) abs(1i.*(x.*p(2)-1./(x.*p(3)))./(p(1) + 1i.*(x.*p(2)-1./(x.*p(3))))); % transfer function of RLC bandstop
% x, p(1), p(2), p(3) correspond to w, R, L, C respectively
x = f(1):f(end); % for plotting purpose
x0 = [1 1e-3 1e-6] % starting points

x0 = 1×3

1 0.001 1e-06

lb = [1 1e-3 1e-6];ub = [100E3 100 1000E-6];
ub = [1 0.02 2e-5];
%x0 = [1 0.0038733109633344873895993294337111 4e-6]

E = @(P) sum((fun(P,f)-y).^2);

syms P [1 3]
E2 = simplify(E([1, P(2), 4e-6]))

E2 = 

dE2 = diff(E2,P(2))

dE2 = 

fplot(E2, [0.003 0.004])

fplot(E2, [0.0038 0.004])

E2f = matlabFunction(E2);
[bestP2, fval_E2] = fminsearch(E2f,50)

bestP2 = 
       0.00393867492675781

fval_E2 = 
      0.000695343018803931

x0 = [1, bestP2, 4E-6]

x0 = 1×3

1 0.00393867492675781 4e-06

[p_fit, fval_fit] = lsqcurvefit(fun, x0, f, y, lb, ub) % p_fit is the vector of fitted parameters

Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 2.000000e+02.

p_fit = 1×3

1 0.00336360777147123 4.69157914821343e-06

fval_fit = 
      9.96139569469118e-05

semilogx(f, y, 'ro', 'displayname', 'original');
hold on; 
semilogx(x, fun(p_fit, x), '--m', 'displayname', 'lsq'); grid on; % semilog plot of the transfer function

[p_fit2, fval_fit2] = lsqcurvefit(fun, x0, f, y, lb, ub)

Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 2.000000e+02.

p_fit2 = 1×3

1 0.00336360777147123 4.69157914821343e-06

fval_fit2 = 
      9.96139569469118e-05

[Pfc, fval_fc] = fmincon(E, x0, [], [], [], [], lb, ub)

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

Pfc = 1×3

1 0.00393867492675781 4e-06

fval_fc = 
      0.000695343018803941

[Pfs, fval_fs] = fminsearch(E, x0)

Pfs = 1×3

1.26924675790965 0.0039541197292072 3.9938320146631e-06

fval_fs = 
      2.13462497232316e-06

semilogx(x, fun(Pfs,x), '--g', 'displayname', 'fminsearch');
semilogx(x, fun(Pfc,x), ':b', 'displayname', 'fmincon');
hold off

legend show

The result of [bestP2, fval_E2] = fminsearch(E2f,50) shows that you do not need to start the search from anywhere close by in order to hit one of the two best regions. The one found, the second of the two dips in the w, gives an even better fit (using fminsearch) than using the first of the w dips.

[Walter Roberson]

The following was not able to find a good solution... but

format long g
f = [80 190 216 224 234 243 254 314]*1e3

f = 1×8

80000 190000 216000 224000 234000 243000 254000 314000

y = [999 984 897 679 205 801 938 993]*1e-3

y = 1×8

0.999 0.984 0.897 0.679 0.205 0.801 0.938 0.993

fun = @(p, x) abs(1i.*(x.*p(2)-1./(x.*p(3)))./(p(1) + 1i.*(x.*p(2)-1./(x.*p(3))))); % transfer function of RLC bandstop
% x, p(1), p(2), p(3) correspond to w, R, L, C respectively
x = f(1):f(end); % for plotting purpose

lb = [1 1e-3 1e-6]; %ub = [100E3 100 1000E-6];
ub = [1 0.02 2e-5];
%x0 = [1 0.0038733109633344873895993294337111 4e-6]

E = @(P) sum((fun(P,f)-y).^2);

syms P [1 3]
E2 = simplify(E([1, P(2), 4e-6]))

E2 = 

dE2 = diff(E2,P(2))

dE2 = 

%fplot(E2, [0.003 0.004])
%fplot(E2, [0.0038 0.004])

E2f = matlabFunction(E2);
[bestP2, fval_E2] = fminsearch(E2f,50)

bestP2 = 
     -1.26659870147705e-06

fval_E2 = 
         0.480361871668374

x0 = [1, bestP2, 4E-6]

x0 = 1×3

1 -1.26659870147705e-06 4e-06

Notice that search came out with a negative P(2), which is out of range.

[p_fit, fval_fit] = lsqcurvefit(fun, x0, f, y, lb, ub) % p_fit is the vector of fitted parameters

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point 
is less than the value of the optimality tolerance.

p_fit = 1×3

1 0.0105 4e-06

fval_fit = 
         0.789425744655355

Those search coordinates are in range, but the result is not great

semilogx(f, y, 'ro', 'displayname', 'original');
hold on; 
semilogx(x, fun(p_fit, x), '--m', 'displayname', 'lsq'); grid on; % semilog plot of the transfer function

[Pfc, fval_fc] = fmincon(E, x0, [], [], [], [], lb, ub)

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

Pfc = 1×3

1 0.00100169385945684 1.10465327170057e-06

fval_fc = 
         0.789396955822885

Again within range, but the results are not great

[Pfs, fval_fs] = fminsearch(E, x0)

Pfs = 1×3

0.983635708916737 -2.43204590273628e-06 6.21518622374251e-06

fval_fs = 
         0.467740743880735

and that one is out of range again

%semilogx(x, fun(Pfs,x), '--g', 'displayname', 'fminsearch');
semilogx(x, fun(Pfc,x), ':b', 'displayname', 'fmincon');
hold off

legend show

I did some grid searching, brute force. The locations I found varied a lot depending how fine a grid I used. The best I found so far was

63535.7188297698 0.710929969491426 0.001

but there was a lot of variability -- such as near 1000 0.01 1e-3 was only marginally more.

In short: there is no good fit.