sdof complex fitting with lsqcurvefit

Tobias

2020 4 月 28

0 回答

2 ビュー (30 日間)

0 投票

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