How to use the nonlinear least square fitting to fit a transfer function？How to improve the fitting accuracy

I have a dataset Y representing the real part Y(:,1) and the imaginary part (Y:,2) of a function against the frequency fq. I want to find the coefficient of the transfer function to have the a good fit with the dataset. The transfer function is given assumed to have a fifth order, (so there are totally 10 unknowns) T(p,q,s)=p(5)*s^4+p(4)*s^3+p(3)*s^2+p(2)*s+p(1))/(s^5+q(5)*s^4+q(4)*s^3+q(3)*s^2+q(2)*s+q(1));

The code is given in the attached, I split the objective function as the imaginary and real part respectively, the objective function is (real(T)-Y(:,1))^2+(imag(T)-Y(:,2)) over all the selected points in the Y-fq curve.

The blue circle is the fitted result from the NLSQ, but they are slightly away from the original data. My question is how to improve the code to improve the fitting result.

In the attached files,"nonlsqfitting.mat" is the Y and fq data; NLSFmat is the runing code; OBJrealimag.m is the objective function, where Yout(:,1) is the real part of T and Yout(:,2) is the imaginary part of T.

[Matt J]

Your initial point,

T0 = ones(10,1); %% intial values

looks very arbitrary. You should devise a smarter initial guess.

Aside from that, I don't really understand the model OBJrealimag. It is a complicated rational function of high order, and there could very easily be a mistake in it. Is yout supposed to be the impulse response, or something else? If something else, providing that will allow a better implementation.

[Catalytic]

In your code, you have --

p5=T(10);p4=T(9);p3=T(8);p2=T(7);p1=T(6);
q5=T(7);q4=T(4);q3=T(3);q2=T(2);q1=T(1);

Surely, it is a mistake that both q5 and p2 are equal to T(7) ? If not, then T(5) goes completely unused. Why not just have 9 unknowns, instead of 10?

[yingyu]

Hi Matt J:

this objective function are the real and imaginary part from the transfer function in the frequency domain. They are from the symbolic expression. So I think the expression is correct. As to the initial guess, shall I spend a lot time on changing the initial value manually or is there any smarter way to change the initial guess?

clear all
clc
syms p5 p4 p3 p2 p1 q5 q4 q3 q2 q1 q0 w real
s=1j*w;
T=(p5*s^4+p4*s^3+p3*s^2+p2*s+p1)/(s^5+q5*s^4+q4*s^3+q3*s^2+q2*s+q1);
[n,d]=numden(expand(T));
[nr,dd]=numden(collect(expand(real(T)),w));
[ni,di]=numden(collect(expand(imag(T)),w));
nr=collect(nr,w),

nr = 

ni=collect(ni,w),

ni = 

dd=collect(dd,w),

dd = 

%% nr is the numerator of the real part
%% ni is the numerator of the imaginary part
%% dd is their denominator

[yingyu]

Hi @Catalytic:

thanks for your observation. After I changing it as q5 =T(5), the result didn't improve if using intial=ones(10,1). But after using the initial guess by @Matt J, the results converged faster.

[Matt J]

or is there any smarter way to change the initial guess?

See my answer below. But also, the problem looks very ill-conditioned, judging by the computed conditionNum. It might be worth lowering the polynomial order.

Answer 1

[Matt J]

Assuming the things @Catalytic pointed out above were mistakes, and also if Y is supposed to be the impulse response, then the following might be what you were truly after:

clear, close all, clc
load("nonlsqfitting.mat");

ft=[0.36:0.02:1]'; %% these frequency points are selected
Hreal=interp1(fq,Y(:,1),ft,"linear","extrap");
Himag=interp1(fq,Y(:,2),ft,"linear","extrap");

%% nonlinear LSQ fitting

ydata2 = [Hreal,Himag];

cofInit=Init(ydata2,2*pi*ft);  %Initial guess

conditionNum = 6.1543e+07

tol=1e-8;
options = optimoptions('lsqcurvefit','MaxIter',inf,'MaxFunEvals',inf,...
                   'OptimalityTol',tol,'StepTol',tol,'FunctionTol',tol);
[cof,resnorm,residuals,exitflag,output] = ...
    lsqcurvefit(@OBJrealimag,cofInit,2*pi*ft,ydata2,[],[],options);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

%% cof is the derived coefficient
yy=OBJrealimag(cof,2*pi*ft);
figure(4)
subplot 211
semilogx(fq,Y(:,1),'k-');
hold on
plot(ft,Hreal,'r.',ft,yy(:,1),'bo');
legend('original real part','selected data','real part after fitting',Location='best')
set(gca,'fontsize',18)
subplot 212
semilogx(fq,Y(:,2));
hold on
plot(ft,Himag,'r.',ft,yy(:,2),'bo');
legend('original imaginary part','selected data','imaginary part after fitting',Location='best')
set(gca,'fontsize',18); set(gcf,'Position',gcf().Position.*[1,1,1.5,1.5])

[cofInit,cof]

ans = 10×2

1.0e+06 * -0.0196 -0.0225 -0.0047 -0.0054 -0.0014 -0.0016 -0.0001 -0.0001 -0.0000 -0.0000 -0.1353 -0.1573 -1.4303 -1.6367 -0.4164 -0.4827 -0.0433 -0.0466 -0.0047 -0.0054

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function T=Init(ydata2,w)

    s=(1j*w(:));

    S=s.^[4,3,2,1,0];
 
    Y=ydata2*[1;1j];

    b=Y.*s.^5;
    A=[S,-Y.*S]; %[p5,....p1, q5,...q1].'
    A=[A;conj(A)];
    b=[b;conj(b)];

    conditionNum=cond(A)

    c=reshape(A\b,[],2);
    c=rot90(c,2);
    c=c(:);

    T=real(c);
    %T=T./abs(T).*abs(c);
    
end


function yout=OBJrealimag(T,w)
%% slip the real and imaginary part of objective function

    p=T(10:-1:6); p=p(:).';
    q=T(5:-1:1);  q=q(:).';

    s=(1j*w(:).');

    H=polyval(p,s)./polyval([1,q],s);

    yout=[real(H).',imag(H).'];

end

[yingyu]

Hi @Matt J:

yes that is what I am looking for. The T=Init(ydata2,w) and the "polyval" are very smart. Thanks!

Answer 2

[Star Strider]

The easiest way to do that is with the System Identification Toolbox. It was created to solve exactly this sort of problem.

Start with the frdata function, and go from there.

[yingyu]

Hi @Star Strider:

thanks for the suggestion. The code using idfrd fits well with the intial data. But how to extract the explicit expression of the identified "sysfr" for later usage?

clear, close all, clc
load("nonlsqfitting.mat");
ft=[0.36:0.02:1]'; %% these frequency points are selected
Hreal=interp1(fq,Y(:,1),ft,"linear","extrap");
Himag=interp1(fq,Y(:,2),ft,"linear","extrap");
%% use idfrd
ResponseData=Hreal+j*Himag;
sysfr = idfrd(ResponseData,ft,0);%%,'FrequencyUnits','Hz');
[resID,frID] = frdata(sysfr);
figure(10)
subplot 211
semilogx(fq,Y(:,1),'k-');
hold on
semilogx(frID,real(squeeze(resID)),'g*');
subplot 212
semilogx(fq,Y(:,2),'k-');
hold on
semilogx(frID,imag(squeeze(resID)),'g*');

[Star Strider]

You would need to map those parameters to a specific model representing the system you are identifying. It is not obvious to me what that is.

EDIT -- (01 Apr 2026 at 11:57)

The respective coefficient values would relate to the p and q values in the appropriate equations 13 and 15 ( and ) in your Comment. If you had a model for that circuit, you might even be able to estimate the component values using the derived (fitted) parameters.

[yingyu]

Hi @Star Strider: Thanks for your prompt reply. But actually the transfer function doesn't represent the impedance of a physical circuit, which may need the synthesis technique to identify the components parameters. The convolutional integral represents the retardation force from the water on the subject once it's moving in the water. I guess the reason is I didnot have the higher frequency data, which I will do later. Many thanks!

[Star Strider]

My pleasure!

Answer 3

[Torsten]

Replace

[cof,resnorm,residuals,exitflag,output] = ...
    lsqcurvefit(@OBJrealimag,T0,2*pi*ft,ydata2);

by

options = optimset('MaxFunEvals',200000,'MaxIter',20000);
[cof,resnorm,residuals,exitflag,output] = ...
    lsqcurvefit(@OBJrealimag,T0,2*pi*ft,ydata2,[],[],options);

and you'll get

clear all
clc
load("nonlsqfitting.mat");
% fq is frequency
%Y(:,1) is real part
%%Y(:,2) is imaginary part
ft=[0.36:0.02:1]'; %% these frequency points are selected
Hreal=interp1(fq,Y(:,1),ft,"linear","extrap");
Himag=interp1(fq,Y(:,2),ft,"linear","extrap");
figure(3)
subplot 211
plot(fq,Y(:,1),'k-');
hold on
plot(ft,Hreal,'ro');
legend('real part','data selected')
subplot 212
plot(fq,Y(:,2),'k-');
hold on
plot(ft,Himag,'ro');
legend('imaginary part','data selected')

%% nonlinear LSQ fitting
T0 = ones(10,1); %% intial values
ydata2 = [Hreal,Himag];
options = optimset('MaxFunEvals',200000,'MaxIter',20000);
[cof,resnorm,residuals,exitflag,output] = ...
    lsqcurvefit(@OBJrealimag,T0,2*pi*ft,ydata2,[],[],options);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

cof

cof = 10×1

1.0e+04 * 0.0281 -0.0303 -0.0272 -0.0030 0.0001 1.4104 -0.0006 -1.5035 -2.6821 -0.0693

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% cof is the derived coefficient
yy=OBJrealimag(cof,2*pi*ft);
figure(4)
subplot 211
semilogx(fq,Y(:,1),'k-');
hold on
plot(ft,Hreal,'ro',ft,yy(:,1),'bo');
legend('original real part','selected data','real part after fitting')
set(gca,'fontsize',8)
subplot 212
semilogx(fq,Y(:,2));
hold on
plot(ft,Himag,'ro',ft,yy(:,2),'bo');
legend('original imaginary part','selected data','imaginary part after fitting')
set(gca,'fontsize',8)

function yout=OBJrealimag(T,w)
%% slip the real and imaginary part of objective function
yout = zeros(length(w),2);
p5=T(10);p4=T(9);p3=T(8);p2=T(7);p1=T(6);
q5=T(7);q4=T(4);q3=T(3);q2=T(2);q1=T(1);
nr=(p5*q5 - p4)*w.^8 + (p2 - p3*q5 + p4*q4 - p5*q3)*w.^6 + (p1*q5 - p2*q4 + p3*q3 - p4*q2 + p5*q1)*w.^4 + (p2*q2 - p1*q3 - p3*q1)*w.^2 + p1*q1;
dd=w.^10 + (q5^2 - 2*q4)*w.^8 + (q4^2 + 2*q2 - 2*q3*q5)*w.^6 + (q3^2 + 2*q1*q5 - 2*q2*q4)*w.^4 + (q2^2 - 2*q1*q3)*w.^2 + q1^2;
ni=- p5*w.^9 + (p3 - p4*q5 + p5*q4)*w.^7 + (p2*q5 - p1 - p3*q4 + p4*q3 - p5*q2)*w.^5 + (p1*q4 - p2*q3 + p3*q2 - p4*q1)*w.^3 + (p2*q1 - p1*q2)*w;
yout(:,1)=nr./dd;
yout(:,2)=ni./dd;
end

[yingyu]

Hi @Torsten:

Many thanks for the suggestion. I will try the optimset. It seems T0 = ones(10,1) leads to slow convergence.

[Torsten]

I guess there is an error in your function OBJrealimag. Most probably,

p5=T(10);p4=T(9);p3=T(8);p2=T(7);p1=T(6);
q5=T(7);q4=T(4);q3=T(3);q2=T(2);q1=T(1);

should read

p5=T(10);p4=T(9);p3=T(8);p2=T(7);p1=T(6);
q5=T(5);q4=T(4);q3=T(3);q2=T(2);q1=T(1);

Answer 4

[Alex Sha]

the stable and unique solution looks like:

Sum Squared Error (SSE): 17.6992228217781
Root of Mean Square Error (RMSE): 0.15403084029212
Correlation Coef. (R): 0.999999105092174
R-Square: 0.999998210185149

Parameter	Best Estimate    
---------	-------------    
p1       	-11.8517124427641
p2       	-839.65661473504 
p3       	-1519.73716239633
p4       	-1061.21460807714
p5       	-700.396027045094
q1       	-1.82851480166084
q2       	-2.72048985167884
q3       	-5.30556239659266
q4       	-1.88516862790752
q5       	-1.53818458328205

[Matt J]

Why stable and unique? There are 3 different solution sets presented in this thread which all appear to fit very well.

[yingyu]

Hi @Matt J and @Alex Sha:

I don't know how you get the other three solutions, by changing the initial conditions? And I found the transfer function might not be the thing I am looking for. Allow me to give a bit of introduction. The hydrodynamics radiation force is in the convolution term involving the retardation function "h(t)" and the device velocity "dx" in hydrodynamics. And I need to find the state space/ transfer function for the convolution term for computing more complex dynamics. And I have the frequency magnitude and phase of the frequency H(w) of "h(t）”， so that I want to identify the transfer functions of H(w) in the 5th order term. However, after I got the parameters from NLS @Matt J, I use the impulse function to represent h(t), the response goes to infinity. So I guess, there are other solutions. The attached figure shows what the impulse response should be like.

So briefly, I am look a transfer function and its dynamics system has the impulse response as the attached figure below and the frequency real and imaginary components as shown before.

Could you @Matt J and @Alex Sha pls give some suggestions on changing the initial conditions to get a unique and stable dynamics system.

As suggested by @Torsten and @Catalytic, I have corrected the objective function. Pls the attached file.

[Matt J]

Could you @Matt J and @Alex Sha pls give some suggestions on changing the initial conditions to get a unique and stable dynamics system.

1. As I commented above, maybe degree 5 is too high.
2. You could use a broader range of frequencies for the fit
3. You could gather more ydata by measuring the responses to multiple input functions (not just an impulse input) and fit to them simultaneously.
4. Your Equation (14) above denotes the input as u(t), often the notation for a unit step. This makes me wonder if the ydata you've given is the unit step response, not the unit impulse response as the fitting model assumes.

[yingyu]

I see. I would use a wider range of frequency data, because theoretically when the frequency goes to infinity, the transfer function should goes to zero, but the current data just end at 1, did not cover larger frequency area. And then I would use the discrete impulse signal to test the new system. Thanks a lot to all!