如何在MATLAB中拟合指数曲线来抑制谐波振荡数据？

Question

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

我需要找到数据中的衰变率.我正在使用的方法可以在这里找到.它是如何工作的,它是将y值的对数高于稳态值,然后使用:

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

适合它.

但是,这会导致以下数据拟合:

我尝试使用线性回归拟合,这显然不起作用,因为它取平均值.我还尝试过RANSAC认为峰值附近有更多数据.它的工作效果比线性回归稍好一些,但方法存在缺陷,因为有时在错误的区域存在更多的点.

有谁知道一个很好的方法来适应这些数据的峰值？

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

Answer 1

我想我应该向每个人提供可能有效的潜在解决方案的最新信息。如前所述，正弦频率的变化使数据变得复杂，因此某些方法可能不起作用。根据所涉及的数据和频率，下面列出的方法可能很好。

\n\n

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

\n\n

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

\n\n

就我而言，平均值为 290，因此我们有：

\n\n

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

\n\n

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

\n\n

findpeaks() 方法

\n\n

这是Alexander B\xc3\xbcse 建议的方法。对于大多数数据来说，这是一个非常好的方法，但对于我的数据来说，由于存在多个正弦频率，因此它会得到错误的峰值。红色 x 显示峰值。

\n\n

% Find Peaks Method\n[max_num,max_ind] = findpeaks(y(ind));\nplot(max_ind,max_num,\'x\',\'Color\',\'r\'); hold on;\nx1 = max_ind;\ny1 = log(max_num-290);\ncoeffs = polyfit(x1,y1,1)\nb = exp(coeffs(2));\nc = coeffs(1);\n

\n\n

\n\n

兰萨克

\n\n

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

\n\n

% RANSAC Method\nind = (y > avg);\nx1 = x(ind);\ny1 = log(y(ind) - avg);\niterNum = 300;\nthDist = 0.5;\nthInlrRatio = .1;\n[t,r] = ransac([x1;y1\'],iterNum,thDist,thInlrRatio);\nk1 = -tan(t);\nb1 = r/cos(t);\n% plot(x1,k1*x1+b1,\'r\'); hold on;\nb = exp(b1);\nc = k1;\n

\n\n

\n\n

线性方法

\n\n

这里使用的就是这种方法。它使用 Lsqlin 来约束系统。然而，它似乎忽略了中间的数据。根据您的数据集，这可能非常有效，就像原始帖子中的人一样。

\n\n

% Lsqlin Method\navg = 290;\nind = (y > avg);\nx1 = x(ind);\ny1 = log(y(ind) - avg);\nA = [ones(numel(x1),1),x1(:)]*1.00;\ncoeffs = lsqlin(A,y1(:),-A,-y1(:),[],[],[],[],[],optimset(\'algorithm\',\'active-set\',\'display\',\'off\'));\nb = exp(coeffs(2));\nc = coeffs(1);\n

\n\n

\n\n

查找周期内的峰值

\n\n

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

\n\n

% Incremental Method 2 Unknowns\nx1 = [];\ny1 = [];\nmax_num=[];\nmax_ind=[];\nincr = 25;\nfor i=1:floor(size(y,1)/incr)\n    [max_num(end+1),max_ind(end+1)] = max(y(1+incr*(i-1):incr*i));\n    max_ind(end) = max_ind(end) + incr*(i-1);\n    if max_num(end) > avg\n        x1(end+1) = max_ind(end);\n        y1(end+1) = log(max_num(end)-290);\n    end\nend\nplot(max_ind,max_num,\'x\',\'Color\',\'r\'); hold on;\ncoeffs = polyfit(x1,y1,1)\nb = exp(coeffs(2));\nc = coeffs(1);\n

\n\n

使用全部 500 个数据点:\n

\n\n

使用前 150 个数据点:\n

\n\n

查找 b 受约束期间的峰值

\n\n

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

\n\n

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

\n\n

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

\n\n

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

\n\n

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

\n