MATLAB中激光图像拟合对数正态曲线的方程求解问询

对数正态分布参数与方程说明 #

MATLAB中lognfit或fitdist返回的mu和sigma，对应的是对数变量的正态分布参数：若随机变量Y服从对数正态分布，则ln(Y) ~ N(mu, sigma²)。其概率密度函数（PDF）为：
$$
f(y) = \frac{1}{y \sigma \sqrt{2\pi}} e^{-\frac{(\ln y - \mu)^2}{2\sigma2}}, \quad y > 0
$$

如果是将像素位置x作为自变量，拟合后的强度y作为因变量（匹配你E列的拟合数据），需要将上述PDF缩放至实际强度的幅值范围，形式为：
$$
y(x) = A \cdot \frac{1}{(x - x_0) \sigma \sqrt{2\pi}} e^{-\frac{(\ln(x - x_0) - \mu)^2}{2\sigma2}}
$$
其中A是强度幅值系数，x_0是偏移量（确保x - x_0 > 0，对应分布的起始位置）。

具体实现代码 #

1. 基于像素强度的分布拟合曲线（概率密度） #

假设你已从Excel读取E列数据到变量y_fit：

% 读取数据（示例）
data = readmatrix('sample.xlsx', 'Sheet', 1, 'Range', 'E:E');
y_fit = data(~isnan(data)); % 去除空值

% 拟合对数正态分布参数
params = lognfit(y_fit);
mu = params(1);
sigma = params(2);
% 或用fitdist
pd = fitdist(y_fit, 'Lognormal');
mu = pd.mu;
sigma = pd.sigma;

% 生成拟合曲线的x轴（强度值范围）
y_range = linspace(min(y_fit), max(y_fit), 1000);
% 计算对数正态PDF
pdf_vals = lognpdf(y_range, mu, sigma);

% 绘图
figure;
histogram(y_fit, 'Normalization', 'pdf'); % 原始数据归一化直方图
hold on;
plot(y_range, pdf_vals, 'r-', 'LineWidth', 2);
xlabel('像素强度');
ylabel('概率密度');
legend('原始数据直方图', '对数正态拟合曲线');

2. 基于像素位置的强度拟合曲线 #

假设你已读取像素位置x（比如列索引，从1到N）和对应的E列拟合强度y_fit：

% 读取数据（示例）
data = readmatrix('sample.xlsx', 'Sheet', 1, 'Range', 'A:E');
x = data(~isnan(data(:,5)), 1); % 假设A列是像素位置
y_fit = data(~isnan(data(:,5)), 5); % E列是拟合后强度

% 对y_fit拟合对数正态分布，获取mu和sigma
pd = fitdist(y_fit, 'Lognormal');
mu = pd.mu;
sigma = pd.sigma;

% 生成拟合曲线的x轴（像素位置范围）
x_fit = linspace(min(x), max(x), 1000);
% 先将x转换为正数偏移（避免ln(0)或负数）
x_offset = x_fit - min(x) + 1; % 确保x_offset > 0
% 计算对数正态PDF并缩放至实际强度幅值
pdf_vals = lognpdf(x_offset, mu, sigma);
% 匹配实际强度的最大值
A = max(y_fit) / max(pdf_vals);
y_fit_curve = A * pdf_vals;

% 绘图
figure;
plot(x, y_fit, 'b.', 'DisplayName', 'SGolay拟合数据');
hold on;
plot(x_fit, y_fit_curve, 'r-', 'LineWidth', 2, 'DisplayName', '对数正态拟合曲线');
xlabel('像素位置');
ylabel('像素强度');
legend;

注意事项 #

• 若你的像素位置x本身从正数开始，可直接用x代替x_offset，无需偏移
• 强度幅值系数A的作用是将标准化的PDF曲线匹配到实际数据的强度范围，确保拟合曲线和原始拟合数据的幅值一致
• 如果需要更精准的位置-强度拟合，可使用fit函数直接拟合自定义的对数正态模型，例如：

  model = fittype('A/( (x - x0)*sigma*sqrt(2*pi) ) * exp( - (log(x - x0) - mu)^2/(2*sigma^2) )', ...
                  'independent', 'x', 'dependent', 'y');
  fit_result = fit(x, y_fit, model, 'StartPoint', [max(y_fit), min(x), mu, sigma]);
  plot(fit_result, x, y_fit);
  

内容的提问来源于stack exchange，提问作者user31729