149 lines
5.8 KiB
Matlab
149 lines
5.8 KiB
Matlab
function [fitResults, rawCurves] = fitTwoGaussianCurvesToRadialSpectralDistribution(S_k_all, k_rho_vals, varargin)
|
|
%% fitTwoGaussianCurvesToRadialSpectralDistribution
|
|
% Fits a two-Gaussian model to multiple radial spectral curves (S(k_rho)),
|
|
% truncating each curve to user-specified k_rho and amplitude ranges.
|
|
%
|
|
% Author: Karthik
|
|
% Date: 2025-10-15
|
|
% Version: 1.0
|
|
%
|
|
% Inputs:
|
|
% S_k_all - 1xN cell array of radial spectral curves
|
|
% k_rho_vals - vector of corresponding k_rho values
|
|
%
|
|
% Optional name-value pairs:
|
|
% 'KRhoRange' - [k_min, k_max] (default: [min(k_rho_vals), max(k_rho_vals)])
|
|
% 'AmplitudeRange' - [y_min, y_max] (default: [0, Inf])
|
|
% 'ResidualThreshold', 'PositionThreshold', 'AmplitudeThreshold' as before
|
|
%
|
|
% Outputs:
|
|
% fitResults - struct array with fields:
|
|
% .pFit, .kFit, .xFit, .yFit, .kFine, .isValid, .fitMaxKRho
|
|
% rawCurves - struct array of truncated raw curves for plotting
|
|
|
|
% --- Parse optional inputs ---
|
|
p = inputParser;
|
|
addParameter(p, 'KRhoRange', [min(k_rho_vals), max(k_rho_vals)], @(x) isnumeric(x) && numel(x)==2);
|
|
addParameter(p, 'AmplitudeRange', [0, Inf], @(x) isnumeric(x) && numel(x)==2);
|
|
addParameter(p, 'ResidualThreshold', 0.3, @(x) isnumeric(x) && isscalar(x));
|
|
addParameter(p, 'PositionThreshold', 0.1, @(x) isnumeric(x) && isscalar(x));
|
|
addParameter(p, 'AmplitudeThreshold', 0.5, @(x) isnumeric(x) && isscalar(x));
|
|
parse(p, varargin{:});
|
|
opts = p.Results;
|
|
|
|
Ncurves = numel(S_k_all);
|
|
fitResults = struct('pFit',[],'kFit',[],'xFit',[],'yFit',[],...
|
|
'kFine',[],'isValid',[],'fitMaxKRho',[]);
|
|
rawCurves = struct('x', cell(1, Ncurves), 'k', cell(1, Ncurves));
|
|
|
|
for k = 1:Ncurves
|
|
% --- Truncate curve to k_rho and amplitude ranges ---
|
|
xRaw = S_k_all{k};
|
|
% --- Normalize and smooth ---
|
|
xNorm = xRaw / max(xRaw);
|
|
xSmooth = smooth(xNorm, 5, 'moving');
|
|
|
|
mask = (k_rho_vals(:) >= opts.KRhoRange(1)) & (k_rho_vals(:) <= opts.KRhoRange(2)) & ...
|
|
(xSmooth(:) >= opts.AmplitudeRange(1)) & (xSmooth(:) <= opts.AmplitudeRange(2));
|
|
|
|
x = xSmooth(mask);
|
|
kVals = k_rho_vals(mask);
|
|
|
|
rawCurves(k).x = x;
|
|
rawCurves(k).k = kVals;
|
|
|
|
if isempty(x)
|
|
warning('Curve %d is empty after truncation.', k);
|
|
fitResults(k) = fillZeroStructRadial(kVals, opts.KRhoRange(2));
|
|
continue;
|
|
elseif numel(x) < 5
|
|
warning('Curve %d has too few points after truncation, skipping.', k);
|
|
fitResults(k) = fillNaNStructRadial();
|
|
continue;
|
|
end
|
|
|
|
try
|
|
% --- Detect peaks ---
|
|
curveAmp = max(x) - min(x);
|
|
minProm = opts.AmplitudeThreshold * curveAmp;
|
|
[pk, locIdx] = findpeaks(x, 'MinPeakProminence', minProm);
|
|
kPeaks = kVals(locIdx);
|
|
|
|
%% Generic two-Gaussian
|
|
twoGauss = @(p,k) ...
|
|
p(1)*exp(-0.5*((k - p(2))/max(p(3),1e-6)).^2) + ...
|
|
p(4)*exp(-0.5*((k - p(5))/max(p(6),1e-6)).^2);
|
|
|
|
secondaryIdx = find(kPeaks > opts.PositionThreshold, 1, 'first');
|
|
if ~isempty(secondaryIdx)
|
|
kSecondary = kPeaks(secondaryIdx);
|
|
secondaryAmp = pk(secondaryIdx);
|
|
A1_guess = opts.AmplitudeRange(2);
|
|
mu1_guess = 0.0;
|
|
sigma1_guess = 0.1 * kSecondary;
|
|
A2_guess = secondaryAmp;
|
|
mu2_guess = kSecondary;
|
|
sigma2_guess = 0.1 * kSecondary;
|
|
p0 = [A1_guess, mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess];
|
|
else
|
|
p0 = [opts.AmplitudeRange(2), 0, 0.1, 0.01, opts.KRhoRange(2)/3, 0.1];
|
|
end
|
|
|
|
lb = [0, 0, 1e-6, 0, opts.PositionThreshold, 1e-6];
|
|
ub = [2, opts.KRhoRange(2)/2, opts.KRhoRange(2), 2, opts.KRhoRange(2), opts.KRhoRange(2)];
|
|
optsLSQ = optimoptions('lsqcurvefit','Display','off', ...
|
|
'MaxFunctionEvaluations',1e4,'MaxIterations',1e4);
|
|
pFit = lsqcurvefit(twoGauss, p0, kVals(:), x(:), lb, ub, optsLSQ);
|
|
|
|
if pFit(4) <= 1e-3
|
|
pFit(4) = 0;
|
|
pFit(5:6) = NaN;
|
|
end
|
|
|
|
fitMaxKRho = max(pFit(5), pFit(2));
|
|
kFine = linspace(min(kVals), max(kVals), 500);
|
|
yFit = twoGauss(pFit, kFine);
|
|
|
|
if ~isnan(pFit)
|
|
fitResults(k).pFit = pFit;
|
|
fitResults(k).kFit = kVals;
|
|
fitResults(k).xFit = x;
|
|
fitResults(k).kFine = kFine;
|
|
fitResults(k).yFit = yFit;
|
|
fitResults(k).isValid = true;
|
|
fitResults(k).fitMaxKRho = fitMaxKRho;
|
|
end
|
|
|
|
catch
|
|
warning('Curve %d fit failed completely.', k);
|
|
fitResults(k) = fillNaNStructRadial();
|
|
end
|
|
end
|
|
end
|
|
|
|
%% --- Helper: NaN struct for failed fits ---
|
|
function s = fillNaNStructRadial()
|
|
s = struct( ...
|
|
'pFit', nan(1,6), ...
|
|
'kFit', nan, ...
|
|
'xFit', nan, ...
|
|
'yFit', nan, ...
|
|
'kFine', nan, ...
|
|
'isValid', false, ...
|
|
'fitMaxKRho', nan ...
|
|
);
|
|
end
|
|
|
|
%% --- Helper: Zero struct for empty curves ---
|
|
function s = fillZeroStructRadial(kVals, MaxKRho)
|
|
s = struct( ...
|
|
'pFit', zeros(1,6), ...
|
|
'kFit', kVals, ...
|
|
'xFit', zeros(size(kVals)), ...
|
|
'yFit', zeros(size(kVals)), ...
|
|
'kFine', zeros(1,500), ...
|
|
'isValid', false, ...
|
|
'fitMaxKRho', MaxKRho ...
|
|
);
|
|
end
|