# StartPoint on Curve Fitting Toolbox

14 views (last 30 days)
Cesium Modern on 12 Jun 2024
Commented: Mathieu NOE on 14 Jun 2024
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
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.
%
% 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)

Mathieu NOE on 13 Jun 2024
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 :
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 on 13 Jun 2024
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 on 14 Jun 2024
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
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

### Categories

Find more on Interpolation in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!