fittting data using nonlinear optimization fminsearch

I have data as below that i want to fit to nonlinear equation. in theory it says that this data can be fit by taking square envelope and will obey equation

my aim is to get optimum Qc based of the data. I wrote a code such below, but the result seems unreasonable (compare to other studies).

clear all; clc; close all;

% Use current working directory
input_dir = pwd; % Directory containing ACF files
output_dir = pwd; % Directory to save outputs

% Get the list of ACF files
acf_files = dir(fullfile(input_dir, 'cBPCen1_s2D1D_ACF_HGSD_HHZ_.HHZ-.HHZ_2023-12_2023-12-01T00_00.txt'));

% Model function
freq=1;
model = @(Qc, tACF) (1 ./ tACF.^2) .*exp(-2 * pi * freq * tACF ./ Qc);

% Nonlinear optimization using fminsearch
objective_function = @(Qc, t, y_obs) sum((y_obs - model(Qc, t)).^2);

% Loop through all ACF files
for i = 1:length(acf_files)
    % Load ACF file
    acf_file = fullfile(input_dir, acf_files(i).name);
    data = load(acf_file);
    
    % ACF signal
    y_obs_ACF = data((length(data)/2):(end-2000))';
    t_ACF = 1:length(y_obs_ACF);
    
    % Compute the envelope
    envelope = abs(hilbert(y_obs_ACF)).^2; %square enveloped of ACF
    t_envel = t_ACF;
    y_obs_envel = envelope;
    
    % Initial guess for Qc
    Qc_initial = 10;

    % Fit ACF using fminsearch
    Qc_fitted_ACF = fminsearch(@(Qc) objective_function(Qc, t_ACF, y_obs_ACF), Qc_initial);
    
    % Calculate RMS error for ACF
    rms_error_ACF = sqrt(mean((y_obs_ACF - model(Qc_fitted_ACF, t_ACF)).^2));
    
    % Fit envelope using fminsearch
    Qc_fitted_envel = fminsearch(@(Qc) objective_function(Qc, t_envel, y_obs_envel), Qc_initial);
    
    % Calculate RMS error for envelope
    rms_error_envel = sqrt(mean((y_obs_envel - model(Qc_fitted_envel, t_envel)).^2));

    % Save the envelope to a text file
    envelope_file = fullfile(output_dir, ['envel_', acf_files(i).name]);
    save(envelope_file, 'envelope', '-ascii');
    
    % Plot both the ACF and the envelope fits
    figure;
    subplot(2, 1, 1);
    plot(t_ACF, y_obs_ACF, 'bo', 'DisplayName', 'ACF Data');
    hold on;
    plot(t_ACF, model(Qc_fitted_ACF, t_ACF), 'r-', 'LineWidth', 1.5, 'DisplayName', 'Fitted ACF');
    xlabel('Time (t)');
    ylabel('Amplitude');
    legend('Location', 'northeast');
    title(sprintf('ACF Fit - Qc = %.4f, RMS Error = %.4f', Qc_fitted_ACF, rms_error_ACF));
    grid on;

    subplot(2, 1, 2);
    plot(t_envel, y_obs_envel, 'bo', 'DisplayName', 'Envelope Data');
    hold on;
    plot(t_envel, model(Qc_fitted_envel, t_envel), 'r-', 'LineWidth', 1.5, 'DisplayName', 'Fitted Envelope');
    xlabel('Time (t)');
    ylabel('Amplitude');
    legend('Location', 'northeast');
    title(sprintf('Envelope Fit - Qc = %.4f, RMS Error = %.4f', Qc_fitted_envel, rms_error_envel));
    grid on;
end

this is my result. I don't know weather my i made mistake in my code or not.

I wonder why i get unreasonable Q is, and that i try to plot how the exuation should look like if i vary the Qc up to reasonable value (again based on reserach paper it has range between 100 to 200 for freq =1, and it looks like below. As I increase my Qc, it does even slightly approach my envelope.

so my question is

• does it because the eq can not really fit my data that make i get big value for Qc or do I make mistake in my code?
• what other suggestion that help to get optimum Qc based on my data is welcome.

[Sam Chak]

Hi @nirwana

Was your MATLAB code generated by a chatbot? If not, please share the link so we can look into the problem.

By the way, do you want to continue with this thread or the newly posted one?

https://www.mathworks.com/matlabcentral/answers/2172326-help-me-to-do-nonlinear-fitting-use-lsqcurvefit?s_tid=srchtitle

Answer 1

[Matt J]

I see two things that might be a problem. First, your model equation might be reasonable for the envelope of the ACF curves you are showing, but you are giving the entirety of the ACF data to the fitting routine, rather than just the data for its envelope. Second, you are choosing seemingly some random number Qc=10 for the initial guess. The more methodical way to develop an initial guess is to solve log(E(t))=.... for Qc.

[nirwana]

related with "First, your model equation might be reasonable for the envelope of the ACF curves you are showing, but you are giving the entirety of the ACF data to the fitting routine, rather than just the data for its envelope." I already triying to cut the ACF data up to 500 point (even i tried also 100 point start from second half data), but still Qc seems unreaasonable.

Answer 2

[Mathieu NOE]

hello

seems to me the data looks like a traditionnal exp decaying pulse - I don't undertsand how you can have the 1/t² amplitude term in the equation (and it would go to infinity at t = 0 ?? .)

A simple code without any optimisation can give you already some ggod fit (could be further refined if needed with your prefered optimizer)

NB : we don't know the sampling rate so time increment is 1 by default

y = readmatrix('cBPCen1_s2D1D_ACF_HGSD_HHZ_.HHZ-.HHZ_2023-12_2023-12-01T00_00.txt');

%% LINEAR DAMPING identification
[val,ind] = max(y); 
y = y(ind:ind+800); % second half of signal after the peak
t = 1:numel(y);

[Yp,Xp,Wp,Pp] = findpeaks(y,'NPeaks',6,'MinPeakHeight',val/50,'MinPeakDistance',10); % all positive and negative peaks
Xpp = t(Xp);
Ypp = y(Xp);
n_peaks=numel(Xpp);  


%---Calculate Logarithmic Decrement,  undamped natural frequency, damped natural frequency, damping ratio
Log_Dec = zeros(length(n_peaks));  
for nn = 1:n_peaks-1 %
    Log_Dec(nn)= log(Ypp(nn)/Ypp(nn+1));  % computed with n = 1 periods apart
end  
Mean_dec = mean(Log_Dec); %Calculate Average Logarithmic Decrement  
wn = 2*pi/mean(diff(Xpp)); %Calculate Average Pulsation

%Damping  
dzeta = 1/sqrt(1+((2*pi/(Mean_dec))^2)); % normalized damping coefficient
Q = 1/(dzeta);

yf = val*exp(-wn*t/Q); % enveloppe
yfs = yf.*cos(wn*t) ; % enveloppe * sinusoid

Rsquared = my_Rsquared_coeff(y(:),yfs(:));

figure(1),plot(t,y,t,yf,t,yfs,Xpp,Ypp,'dr')
title(['Data fit , R² = ' num2str(Rsquared)]);

xlabel('time [s]'); ylabel('amplitude')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Mathieu NOE]

Added a second layer with fminsearch , but the IC obtained with my first code was already giving good estimates

Q (ic) = 8.4160

Q (final) = 7.7196

here fminsearch tends to improve the fit on the large amplitudes (first 3 periods) , but the error is getting bigger afterwards

y = readmatrix('cBPCen1_s2D1D_ACF_HGSD_HHZ_.HHZ-.HHZ_2023-12_2023-12-01T00_00.txt');

%% LINEAR DAMPING identification
[val,ind] = max(y); 
y = y(ind:ind+800); % second half of signal after the peak
t = (1:numel(y))';

[Yp,Xp,Wp,Pp] = findpeaks(y,'NPeaks',6,'MinPeakHeight',val/50,'MinPeakDistance',10); % all positive and negative peaks
Xpp = t(Xp);
Ypp = y(Xp);
n_peaks=numel(Xpp);  


%%---Calculate Logarithmic Decrement,  undamped natural frequency, damped natural frequency, damping ratio
Log_Dec = zeros(length(n_peaks));  
for nn = 1:n_peaks-1 %
    Log_Dec(nn)= log(Ypp(nn)/Ypp(nn+1));  % computed with n = 1 periods apart
end  
Mean_dec = mean(Log_Dec); %Calculate Average Logarithmic Decrement  
wn = 2*pi/mean(diff(Xpp)); %Calculate Average Pulsation

%Damping  
dzeta = 1/sqrt(1+((2*pi/(Mean_dec))^2)); % normalized damping coefficient
Q = 1/(dzeta);

yf = val*exp(-wn*t/Q); % enveloppe
yfs = yf.*cos(wn*t) ; % enveloppe * sinusoid

Rsquared = my_Rsquared_coeff(y(:),yfs(:));

figure(1),plot(t,y,t,yf,t,yfs,Xpp,Ypp,'dr')
title(['Data fit , R² = ' num2str(Rsquared)]);

xlabel('time [s]'); ylabel('amplitude')
legend('data','enveloppe fit','signal fit','peaks');

%% step 2 : refine with fminsearch 
f = @(a,b,c,x) a.*exp(-b*x/c).*cos(b*x);
obj_fun = @(params) norm(f(params(1), params(2), params(3),t)-y);
IC = [val,wn,Q];
sol = fminsearch(obj_fun, IC);
val = sol(1);
wn = sol(2);
Q = sol(3);

yfit = f(val,wn,Q,t);
Rsquared = my_Rsquared_coeff(y(:),yfit(:));

figure(2),plot(t,y,t,yfit)
title(['Data fit , R² = ' num2str(Rsquared)]);
xlabel('time [s]'); ylabel('amplitude')
legend('data','signal 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%