StartPoint on Curve Fitting Toolbox

I am trying to fit a linear resonance curve of a 2nd order mechanical system to get the Q-factor.

The instrument that I am using to obtain the frequency response function has an in-built curve fitting software that uses the following equation to fit the curve.

denominator = sqrt(f^2 + (Q / f0)^2 * (f^2 - f0^2)^2);

Unrecognized function or variable 'f'.

X=A * f / denominator;

(C=0) *image q2.png

A is the amplitude, 𝑄 is the quality factor and 𝑓0 is the center frequency.

My goal is to use MATLAB's Curve Fitting tool to obtain the Q-factor value; first I want to match the Q-factor value obtained from the instrument (as a validation step, so I can fit more data without requiring the instrument's curve fitting tool).

Below is the FRF obtained from the instrument--the Q factor value I obtain is 8705.9559

*image q1.png

Center Frequency (f0) = 1496.26

A = 2.4216e-3 (amplitude)

I saved this sweep and plotted and tried to recreate the this fit using the CurveFitting tool box--using the same cut-off range as above (1453 kHz to 1548 kHz)

However, I am having an issue with converging to the value for the Q-factor from the instrument--it varies as I change my StartPoint and Lower and Upper values. I understand this is the case, sensitive to inital conditions--however, I was under the assumption that it may converge to the value as I refine my StartPoint and Lower and Upper bounds---this is not the case.

For example:

With StartPoint= 0; Lower= 10; Upper= 10e3 Q=8052 (7994, 8110) with R-square of 0.9961

moving the StartPoint closer to 8000 (which is around the actual value), it drops

With StartPoint = 7000; Lower=10; Upper= 10e3 Q=8278 (8228, 8329) with R-square of 0.9973

With StartPoint= 8000; Lower= 10; Upper= 10e3 Q=8000 (7940, 8060) with R-square of 0.9961

(varying the Upper and Lower bounds also causes the Q factor value to change--and not converge--to different values that is around value obtained from the instrument)

--I was hoping to get closer and improve the goodness of fit, but the fit is very sensitive to the StartPoint---any pointers on how to get an accurate value?

Is this method, fitting this curve using the CurveFitting Toolbox the best method to get the value for Q?

I have shared my data, code and image snippets. Any pointers appreciated.

clc;clear;close all
data = load('data.mat');
f=data.d(:,1);
A=data.d(:,2);
[fitresult, gof] = createFit(f, A)
function [fitresult, gof] = createFit(f, A)
%CREATEFIT(F,A)
%  Create a fit.
%
%  Data for 'untitled fit 1' fit:
%      X Input: f
%      Y Output: A
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 12-Jun-2024 17:27:38


%% Fit: 'untitled fit 1'.
[xData, yData] = prepareCurveData( f, A );

% Set up fittype and options.
ft = fittype( '2.4216e-3*f./(sqrt(f^2+(Q/1496.26)^2*(f^2-1496.26^2)^2))', 'independent', 'f', 'dependent', 'A' );
excludedPoints = (xData < 1426.9) | (xData > 1570.4);
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = 10;
opts.MaxFunEvals = 1000;
opts.MaxIter = 1000;
opts.StartPoint = 7000;
opts.TolFun = 1e-07;
opts.Upper = 10000;
opts.Exclude = excludedPoints;

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.
figure( 'Name', 'untitled fit 1' );
h = plot( fitresult, xData, yData, excludedPoints );
legend( h, 'A vs. f', 'Excluded A vs. f', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'f', 'Interpreter', 'none' );
ylabel( 'A', 'Interpreter', 'none' );
grid on
end

TL;DR--I am trying to understand why the value of Q, obtained from my fit, does not converge to a value (and improve the fit) as I use a guess that is around the accurate value, and reduce the size of the bounds (Lower and Upper)

Answer 1

[Mathieu NOE]

hello

here a very simple code with a 2 steps approach - a relatively accurate initial guess of Aand Q, then refined with fminsearch (I don't have the CFT) - or use your prefered optimizer instead

Results obtained

A = 0.0024

Q = 8.7126e+03

Zoom on peak

code :

load('data.mat')

f = d(:,1);
df = mean(diff(f)); % frequency resolution = 0.18 Hz
y = d(:,2);

%% parameters initial guess
% peak frequency 
[y0,ind] = max(y);
f0 = f(ind);

% Q factor guess by -3dB rule (Q = f0/df @ -3dB crossing line)
[flow,fup] = find_zc(f,y,0.707*y0);
Q = f0/(fup - flow);

% initial A value (not critical)
A = y(1);

denominator = sqrt(f.^2 + (Q / f0).^2 * (f.^2 - f0^2).^2);
X= A * f./denominator;

% correct A value
A = A*mean(y./X);
X= A * f./denominator;

% correct Q value
[Xm,ind] = max(X);
Q = Q*y0/Xm;

denominator = sqrt(f.^2 + (Q / f0).^2 * (f.^2 - f0^2).^2);
X= A * f./denominator;


%% curve fit using fminsearch
% We would like to fit the function
fun = @(a,b,x) a*x./sqrt(x.^2+(b/f0).^2*(x.^2 - f0^2).^2);
obj_fun = @(params) norm(fun(params(1), params(2),f)-y);
C1_guess = [A Q];
sol = fminsearch(obj_fun, C1_guess); % 
A = sol(1)

A = 0.0024

Q = sol(2)

Q = 8.7126e+03

y_fit = fun(A, Q, f);
Rsquared = my_Rsquared_coeff(y,y_fit); % correlation coefficient


%plot
semilogy(f,y)
hold on
semilogy(f,X,f, y_fit)
title(['Fit equation to data : R² =  ' num2str(Rsquared) ], 'FontSize', 20)
xlabel('x data', 'FontSize', 20)
ylabel('y data', 'FontSize', 20)
legend('raw data','initial guess','final fit');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Cesium Modern]

Hi Mathieu!

Thanks for the great answer. I really appreciate the time you took to help me with this issue.

I retested this code with a different sample set (first mode data of the DUT from the same instrument, with the same fitting) and I got the following:

The Q factor is now 26e3, y0 is 469.78 kHz

The raw data was saved in the following file, m1data.mat, which is attached to this comment.

I was not able to replicate the same Q-factor from the instrument software's fitting (which uses the same model)

(Q value obtained is 15e3, and the final fit does not line up with the raw data)

I even try to mimic the bounds used in the software fit to see if that would improve the Q value and fit, but it did not

%bounds
upperB = 494.29e3;
lowerB = 449.09e3;
bound_idx = find(f>=lowerB & f<=upperB);
f=f(bound_idx);y=y(bound_idx);

I would like to hear your thoughts.

Also, I am trying to understand the logic behind your code, I do not follow the calculation for A--observing the model, wouldn't it be max(y)

% initial A value (not critical)
A = y(1);
denominator = sqrt(f.^2 + (Q / f0).^2 * (f.^2 - f0^2).^2);
X= A * f./denominator;
% correct A value
A = A*mean(y./X);   %?
X= A * f./denominator;

[Mathieu NOE]

hello

there many aspects in the question , and therefore I need a bit of time to best answer those

first comment : the first data set was easier to handle , as the Q factor is lower and the peak frequency exactly falls on one of the frequency bins (in the f array). So all il all , it was easier for the initial guess and also fminsearch could find the good values (there was only 2 parameters to use in the fit)

now you last data set is more difficult for the same reasons : Q factor is much higher , so frequency location of the "true" peak is going to be more delicate , and this appears in the data ,as now f0 seems to fall somewhere between two values of the f array. So now we cannot simply pick the frequency of the peak value, we have to adjust this as well, so the problem is now more difficult to solve for fminsearch as we have 3 and not 2 parameters to tweak.

to make things even more challenging , if you don't put some bounds on the parameters you can end up with similar plots but very different A and Q combinations , with a common multiplicative factor that can range from 1 to 1000, again , with very similar plots

Reducing the frequency range has also some impact on the result - so touching one or the other parameter can have a serious effect on the final result.

so , to avoid this I opted for the bounded version of fminsearch : fminsearchbnd, fminsearchcon - File Exchange - MATLAB Central (mathworks.com)

- but again , as you have the optimisation toolbox, you may find a better alternative

so the best code I can show today - for this new data - is this one :

results :

A = 0.0036

Q = 2.4991e+04

f0 = 4.6978e+05

load('m1data.mat')

f = ndata(:,1);
df = mean(diff(f)); % frequency resolution = 0.18 Hz
y = ndata(:,2);

%% parameters initial guess
% Q factor guess by -3dB rule (Q = f0/df @ -3dB crossing line)
% refined peak frequency f0 = half value of flow,fup
y0 = max(y);
[flow,fup] = find_zc(f,y,10^(-3/20)*y0);
f0 = 0.5*(flow + fup);% peak frequency 
Q = f0/(fup - flow);

% initial A value (not critical)
A = y(1);

denominator = sqrt(f.^2 + (Q / f0).^2 * (f.^2 - f0^2).^2);
X= A * f./denominator;

% correct A value
A = A*mean(y./X);
X= A * f./denominator;

% correct Q value
[Xm,ind] = max(X);
Q = Q*y0/Xm;

denominator = sqrt(f.^2 + (Q / f0).^2 * (f.^2 - f0^2).^2);
X= A * f./denominator;

%% curve fit using fminsearch
lowerB = 460e3;
upperB = 480e3;
% lowerB = f0*0.9;
% upperB = f0*1.1;
bound_idx = find(f>=lowerB & f<=upperB);

% We would like to fit the function
fun = @(a,b,c,x) a*x./sqrt(x.^2+(b/c).^2*(x.^2 - c^2).^2);
obj_fun = @(params) norm(fun(params(1), params(2), params(3),f(bound_idx))-y(bound_idx));
C1_guess = [A Q f0];
% sol = fminsearch(obj_fun, C1_guess); % unconstrained
sol = fminsearchbnd(obj_fun, C1_guess,[A*0.3 Q*0.3 f0*0.97],[A*3 Q*3 f0*1.03]); % constrained

A = sol(1)
Q = sol(2)
f0 = sol(3)

y_fit = fun(A, Q, f0, f);

Rsquared = my_Rsquared_coeff(y,y_fit); % correlation coefficient

%plot
% plot(f,y,f, y_fit)
semilogy(f,y,f, y_fit)
title(['Fit equation to data : R² =  ' num2str(Rsquared) ], 'FontSize', 20)
xlabel('x data', 'FontSize', 20)
ylabel('y data', 'FontSize', 20)
% legend('raw data','initial guess','final fit');
legend('raw data','final fit');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
% R2 correlation coefficient computation

% The total sum of squares
sum_of_squares = sum((data-mean(data)).^2);

% The sum of squares of residuals, also called the residual sum of squares:
sum_of_squares_of_residuals = sum((data-data_fit).^2);

% definition of the coefficient of correlation is
Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;


end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Mathieu NOE]

do you have also the phase and not only the magnitude ?

[Cesium Modern]

I do have the phase data as well!

Thanks for the reply, I am still reading through the new code and understanding it.

I have attached the phase data (third column in radians), could you let me know how this will be helpful in figuring out my Q-factor (which is the goal)

Again, I really appreciate your time and expertise.

[Mathieu NOE]

I have to say it's not easy ti find the optimum data

if I try to force the phase plot to match between model and your data , I get those parmeters

A = 0.0012

f0 = 4.6979e+05

Q = 6.5858e+03

notice that with such high Q factors it's easy to have some experimental problems : see the skewed shape of the modulus , the peak is not symmetrical , that happens when systems are slightly non linear or if the chirp is too fast vs the transient (decay) time (which is longer when Q increases)

if now I force Q to a higher value, the phase plots don't match anymore and I have also to correct A to keep both modulus curves as close as possible.

I wonder if your software works only on the modulus or also on the phase

wish also the frequency resolution would be increased

code

load('m1data_withphase.mat')

f = d(:,1);
df = mean(diff(f)); % frequency resolution = 0.18 Hz
modulus = d(:,2);
phase = d(:,3); % in radians
tan_phase = tan(phase); % tangent of phase data

%% parameters initial guess
% Q factor guess by -3dB rule (Q = f0/df @ -3dB crossing line)
% refined peak frequency f0 = half value of flow,fup
[y0,ind] = max(modulus);
[flow,fup] = find_zc(f,modulus,10^(-3/20)*y0);
f0 = 0.5*(flow + fup);% peak frequency 
Q = f0/(fup - flow);
beta = f/f0; % normalized frequency

% initial A value (not critical)
A = modulus(1);
%%%%%%%%%%%%%%%%%%%%
denominator = sqrt(f.^2 + (Q / f0).^2 * (f.^2 - f0^2).^2);
X= A * f./denominator;

% correct A value
A = A*mean(modulus./X);
X= A * f./denominator;

% find phase slope arounf f0 frequency
indf = ind+(-10:10);
ff = f(indf);
beta2 = ff./f0;
tan_phase_theo = -1/Q*beta2./(1-beta2.^2);
tan_phase_meas = tan_phase(indf);
tpmeas_min = min(tan_phase_meas);
tpmeas_max = max(tan_phase_meas);
tpmeas_delta = tpmeas_max-tpmeas_min;
tptheo_min = min(tan_phase_theo);
tptheo_max = max(tan_phase_theo);
tptheo_delta = tptheo_max-tptheo_min;

% correct Q value
Q = Q*tptheo_delta/tpmeas_delta;
tanphase_theo = -1/Q*beta./(1-beta.^2);

%plot
subplot 211,semilogy(f,modulus,'-*',f,X)
xlabel('freq', 'FontSize', 14)
ylabel('modulus', 'FontSize', 14)
legend('raw data','initial guess');
xlim([0.999*f0 1.001*f0])
tanphase_theo_fit = -1/Q*beta./(1-beta.^2);
subplot 212,plot(f,tan(phase),'-*',f,tanphase_theo)
ylim([-5 5])
xlabel('freq', 'FontSize', 14)
ylabel('tan phase', 'FontSize', 14)
xlim([0.999*f0 1.001*f0])

A

f0

Q

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation
zci = @(data) find(diff(sign(data))<0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end