sdof complex fitting with lsqcurvefit

Hello

I am trying to fit complex data of a Frequency Response measurement (i.e. frequency vs. complex response (magnitude/Phase or X/Y) to a simple sdof model, which is defined by the objective function

objfcn =  @(v,xdata)( (1i*xdata)*v(1)*v(3)/v(2) ) ./ ( -xdata.^2 + 1i*xdata*v(1)/v(2) + v(1)^2 )*exp(1i*v(4));

where v(1) is resonance frequency (in rad/sec), v(2) the Q Factor, v(3) the Amplitude at resonance and v(4) a certain (unknown) phase shift at resonance

Here is an example how this function looks like

xdata = linspace(10,20,1000)'; % frequency data
ydata = objfcn([15,100,1,50*pi/180],xdata)+0.01*rand(size(xdata)); % complex frequency response
% plot 
figure;
subplot(1,2,1); hold on; box on; grid on; legend('show'); plot(xdata,abs(ydata),'DisplayName','measurement data'); xlabel('f [rad/sec]'); ylabel('mag');
subplot(1,2,2); hold on; box on; grid on; legend('show'); plot(xdata,(angle(ydata))*180/pi,'DisplayName','measurement data');xlabel('f [rad/sec]'); ylabel('phase [deg]');

I first tried with lsqcurvefit, which gives nice results and a very robust fit.

x0 = [16,80,1,0]; %inital parameters
[vestimated,resnorm] = lsqcurvefit(objfcn,x0,xdata,ydata,[],[]); % perform the fit
% plot
subplot(1,2,1); plot(xdata,abs(objfcn(vestimated,xdata)),'DisplayName','lsqcurvefit');
subplot(1,2,2); plot(xdata,angle(objfcn(vestimated,xdata))*180/pi,'DisplayName','lsqcurvefit');

However, the coefficients are complex, which is what I want to avoid (they are real numbers in the model as well). Therefore I tried to follow the example https://ch.mathworks.com/help/optim/ug/fit-model-to-complex-data.html where it is suggested to split both the objective function and the data in complex and real parts in order to stay completely within real values.

ydata2 = [real(ydata),imag(ydata)]; % split complex frequency response
[vestimated2,resnorm2] = lsqcurvefit(@cplxreal,x0,xdata,ydata2,[],[]); % perform the fit with real values only
yout = cplxreal(vestimated2,xdata); % get the response as matrix with separated real and imag part
youtcmplx = yout(:,1) + 1i*yout(:,2); % build the complex response
% plot
subplot(1,2,1); plot(xdata,abs(youtcmplx),'DisplayName','lsqcurvefit_{real}');
subplot(1,2,2); plot(xdata,angle(youtcmplx)*180/pi,'DisplayName','lsqcurvefit_{real}');


% define the objective function 
function yout = cplxreal(v,xdata)

    yout = zeros(length(xdata),2); % allocate yout
    denom = v(2)^2*xdata.^4 + (1 - 2*v(2)^2)*v(1)^2*xdata.^2 + v(1)^4*v(2)^2; % common denominator
    a = (v(1)^2*v(2) - v(2)*xdata.^2).*v(1)*v(3).*xdata; 
    b = v(1)^2*xdata.^2*v(3);
    yout(:,1) =  -(a*sin(v(4))-b*cos(v(4)))./denom; % real part
    yout(:,2) =  (a*cos(v(4))+b*sin(v(4)))./denom; % imaginary part

end

However, the real approach is way less robust, very sensitive to inital condition and fails most of the time. I tried to play around with algorithms, tolerances, number of iterations etc., without success.

Does anybody has an Idea where the error might be?

Thank you

Tobias

P.S: I know about the TFEST function

[Alex Sha]

Hi, You may try Global Optimization Toolbox for avoiding the guess of initial start-value.

If you could post the fitting function and full data, I may have a try.

[Tobias]

Thanks for your quick reply. I just found the optimization toolbox, but I am having troubles defining the objective function such that I can select it in the toolbox.

Anyway, find attached the full data.

Best, Tobias

[Alex Sha]

Hi, Tobias, the file you uploaded seems corrupted and won't open properly, if possible, upload Excel file format please.

[Tobias]

sorry. Next try.

I was just wondering, if the small numbers (order of 1E-9) might be a problem?

Best, Tobias

[Alex Sha]

Hi, Tobias, it is easy to get unique and stable result by using 1stOpt rather than Matlab, initial start values are not required, result is perfect:

Code:

Parameter v(4);
Variable xdata, ydata[realPart],ydata[imagPart];
Function ydata=((1*i*xdata)*v1*v3/v2)/(-xdata^2+1i*xdata*v1/v2+v1^2)*exp(1*i*v4);
Data;
xdata	ydata_real	ydata_imag
268.6975941	-5.68E-09	8.97E-09
268.6976315	-5.68E-09	8.99E-09
268.697669	-5.69E-09	9.01E-09
268.6977064	-5.69E-09	9.03E-09
268.6977439	-5.69E-09	9.05E-09
268.6977813	-5.69E-09	9.07E-09
268.6978188	-5.69E-09	9.09E-09
...

Results:

Root of Mean Square Error (RMSE): 2.10638823457513E-11
Sum of Squared Residual: 8.87374278951303E-19
Correlation Coef. (R): 0.999996732864244
R-Square: 0.999993465739162

Parameter          	Best Estimate
--------------------	-------------
v1	-268.715000056698
v2	-15309.2844467811
v3	-2.35758125061554E-8
v4	4.17347605549285

[Tobias]

ok, well, this is a Matlab Forum, so is there a solution in Matlab as well?

[Alex Sha]

Then, as said previously: try Global Optimization Toolbox for avoiding the guess of initial start-value.

or try to guess different combination of initial start values until satisfactory results.