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

Question

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

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

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

适合它。

然而，这导致以下数据适合：

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

有人知道适合此数据峰值的好方法吗？

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

Answer 1

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

[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 函数仅 select 最大值然后拟合它们。您可能想尝试一下 MinPeakProminence 参数以根据您的需要进行调整。

Answer 2

我想给大家更新可能有效的潜在解决方案。如前所述，数据因正弦频率的变化而变得复杂，因此某些方法可能无法正常工作。根据所涉及的数据和频率，下面列出的方法可能很好。

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

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 来约束系统。但是，它似乎忽略了中间的数据。根据您的数据集，这可能非常有效，就像它对原始 post.

中的人所做的那样

% 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);

查找周期中的峰值

这是我在 post 中提到的方法，我在其中获得了每个区域的峰值。这种方法效果很好，由此我意识到我的数据实际上可能没有完美的指数拟合。我们看到它一开始无法适应大峰。通过仅使用前 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=(log(y-290)-logb)/x 处的 c。然后我可以取 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的平均值

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

How to fit an exponential curve to damped harmonic oscillation data in MATLAB?

matlab

statistics

signal-processing

regression

curve-fitting