如何在 MATLAB 中将指数曲线拟合到阻尼谐波振荡数据？

Posted 2023-03-12

【中文标题】如何在 MATLAB 中将指数曲线拟合到阻尼谐波振荡数据？

【英文标题】：How to fit an exponential curve to damped harmonic oscillation data in MATLAB?

【发布时间】：2016-08-08 06:15:51

【问题描述】：

我正在尝试将指数曲线拟合到包含阻尼谐波振荡的数据集。从正弦振荡包含许多频率的意义上说，数据有点复杂，如下所示：

我需要找出数据中的衰减率。我使用的方法可以找到here。它是如何工作的，是不是取稳态值以上的y值的对数，然后使用：

lsqlin(A,y1(:),-A,-y1(:),[],[],[],[],[],optimset('algorithm','active-set','display','off'))

为了适应它。

但是，这会导致以下数据拟合：

我尝试使用线性回归拟合，这显然行不通，因为它取平均值。我还尝试了 RANSAC，认为峰值附近有更多数据。它比线性回归效果好一点，但该方法存在缺陷，因为有时在错误区域存在更多点。

有谁知道一种适合此数据峰值的好方法？

目前，我正在考虑将 500 个数据点划分为 10 个不同的区域，并在每个区域中找到最大值。最后，我应该有 50 个点可以使用上述任何指数拟合方法进行拟合。你觉得这个方法怎么样？

【参考方案1】：

我想我会给每个人一个可能有效的潜在解决方案的更新。如前所述，数据因正弦频率的变化而变得复杂，因此某些方法可能因此而不起作用。根据所涉及的数据和频率，下面列出的方法可能很好。

首先，我假设数据具有以下形式：

y = average + b*e^-(c*x)

在我的例子中，平均值是 290，所以我们有：

y = 290 + b*e^-(c*x)

话虽如此，让我们深入了解我尝试过的不同方法：

findpeaks() 方法

这是 Alexander Büse 建议的方法。对于大多数数据来说，这是一个很好的方法，但对于我的数据，由于有多个正弦频率，它会得到错误的峰值。红色 x 表示峰。

% Find Peaks Method
[max_num,max_ind] = findpeaks(y(ind));
plot(max_ind,max_num,'x','Color','r'); hold on;
x1 = max_ind;
y1 = log(max_num-290);
coeffs = polyfit(x1,y1,1)
b = exp(coeffs(2));
c = coeffs(1);

RANSAC

如果您的大部分数据都处于峰值，那么 RANSAC 是很好的选择。你看我的，因为有多个频率，在顶部附近存在更多的峰值。但是，我的数据的问题是，并不是所有的数据集都是这样的。因此，它偶尔会起作用。

% RANSAC Method
ind = (y > avg);
x1 = x(ind);
y1 = log(y(ind) - avg);
iterNum = 300;
thDist = 0.5;
thInlrRatio = .1;
[t,r] = ransac([x1;y1'],iterNum,thDist,thInlrRatio);
k1 = -tan(t);
b1 = r/cos(t);
% plot(x1,k1*x1+b1,'r'); hold on;
b = exp(b1);
c = k1;

Lsqlin 方法

这个方法是使用here的方法。它使用 Lsqlin 来约束系统。但是，它似乎忽略了中间的数据。根据您的数据集，这可能会像对原始帖子中的人一样有效。

% Lsqlin Method
avg = 290;
ind = (y > avg);
x1 = x(ind);
y1 = log(y(ind) - avg);
A = [ones(numel(x1),1),x1(:)]*1.00;
coeffs = lsqlin(A,y1(:),-A,-y1(:),[],[],[],[],[],optimset('algorithm','active-set','display','off'));
b = exp(coeffs(2));
c = coeffs(1);

找到周期内的峰值

这是我在我的帖子中提到的方法，我在每个区域获得峰值，.这种方法效果很好，从中我意识到我的数据实际上可能没有完美的指数拟合。我们看到它无法适应一开始的大峰。通过仅使用前 150 个数据点并忽略稳态数据点，我能够使这一点变得更好。在这里，我发现每 25 个数据点的峰值。

% Incremental Method 2 Unknowns
x1 = [];
y1 = [];
max_num=[];
max_ind=[];
incr = 25;
for i=1:floor(size(y,1)/incr)
    [max_num(end+1),max_ind(end+1)] = max(y(1+incr*(i-1):incr*i));
    max_ind(end) = max_ind(end) + incr*(i-1);
    if max_num(end) > avg
        x1(end+1) = max_ind(end);
        y1(end+1) = log(max_num(end)-290);
    end
end
plot(max_ind,max_num,'x','Color','r'); hold on;
coeffs = polyfit(x1,y1,1)
b = exp(coeffs(2));
c = coeffs(1);

使用全部 500 个数据点：

使用前 150 个数据点：

在 b 约束的周期内找到峰值

因为我希望它从第一个峰值开始，所以我限制了 b 值。我知道系统是y=290+b*e^-c*x，我将其限制为b=y(1)-290。通过这样做，我只需要求解 c where c=(log(y-290)-logb)/x。然后我可以取 c 的平均值或中位数。这种方法也很好，它也不适合接近尾声的值，但这没什么大不了的，因为那里的变化很小。

% Incremental Method 1 Unknown (b is constrained y(1)-290 = b)
b = y(1) - 290;
c = [];
max_num=[];
max_ind=[];
incr = 25;
for i=1:floor(size(y,1)/incr)
    [max_num(end+1),max_ind(end+1)] = max(y(1+incr*(i-1):incr*i));
    max_ind(end) = max_ind(end) + incr*(i-1);
    if max_num(end) > avg
        c(end+1) = (log(max_num(end)-290)-log(b))/max_ind(end);
    end
end
c = mean(c); % Or median(c) works just as good

这里我取每 25 个数据点的峰值，然后取 c 的平均值

这里我取每 25 个数据点的峰值，然后取 c 的中值

这里我取每 10 个数据点的峰值，然后取 c 的平均值

【参考方案2】：

如果主要目标是从拟合中提取阻尼参数，也许您想考虑直接将阻尼正弦曲线拟合到您的数据中。像这样的东西（使用曲线拟合工具创建）：

[xData, yData] = prepareCurveData( x, y );
ft = fittype( 'a + sin(b*x - c).*exp(d*x)', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [1 0.285116122712545 0.805911873245316 0.63235924622541];
[fitresult, gof] = fit( xData, yData, ft, opts );
plot( fitresult, xData, yData );

特别是因为您的一些示例数据在感兴趣的区域（噪声之上）确实没有很多数据点。

但是，如果您确实需要直接拟合实验数据的最大值，则可以使用findpeaks 函数仅选择最大值然后拟合它们。您可能想使用MinPeakProminence 参数来调整它以满足您的需要。

【讨论】：

谢谢亚历山大。所以这个数据的棘手之处在于它不仅仅是一个正弦频率，所以使用 findpeaks 最终得到的波峰实际上并不是该区域的最大值。我用连接点的实际信号更新了我的原始帖子。我还没有尝试用曲线拟合工具箱来拟合它（我的电脑上没有），但我会在学校时尝试一下。