Borrar filtros
Borrar filtros

StartPoint on Curve Fitting Toolbox

15 visualizaciones (últimos 30 días)
Cesium Modern
Cesium Modern el 12 de Jun. de 2024
Comentada: Mathieu NOE el 14 de Jun. de 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
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)

Respuesta aceptada

Mathieu NOE
Mathieu NOE el 13 de Jun. de 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 :
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
  5 comentarios
Cesium Modern
Cesium Modern el 13 de Jun. de 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
Mathieu NOE el 14 de Jun. de 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
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

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Interpolation en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by