如何擬合指數曲線來抑制MATLAB中的諧波振盪數據？

Question

我試圖擬合一個指數曲線到包含阻尼諧波振盪的數據集。該數據爲一位在這個意義上的正弦振盪包含許多複雜的頻率，如下所示：

我需要找到腐爛的數據速率。我正在使用的方法可以找到here。它是如何工作的，是否需要y值高於穩態值的對數，然後使用：

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

適合它。

然而，這會導致以下數據符合：

我試圖用線性迴歸擬合這顯然沒有奏效，因爲它採取了平均水平。我也嘗試了RANSAC，認爲峯頂附近有更多的數據。它比線性迴歸的效果好一些，但該方法存在缺陷，因爲有時在錯誤的地區存在更多的點。

有沒有人知道一個很好的方法來適合這些數據的高峯？

目前，我正在考慮將500個數據點分成10個不同的區域，並在每個區域找到最大的值。最後，我應該有50分，我可以使用上面提到的任何指數擬合方法。你對這種方法有什麼看法？

Answer 1

想到我會給每個人一個可能的解決方案的更新。如前所述，數據由於變化的正弦頻率而變得複雜，因此某些方法可能無法正常工作。根據數據和涉及的頻率，下面列出的方法可能很好。

首先，我認爲該數據有以下形式：

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

在我的情況下，平均是290，所以我們有：

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

有了這樣說，讓我們潛入不同的方法，我試圖：

findpeaks（）方法

這是AlexanderBü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=(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個數據點，然後採取ç

Answer 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參數進行調整以根據您的需求進行調整。