Wie passe ich eine exponentielle Kurve an gedämpfte harmonische Schwingungsdaten in MATLAB an?

Question

Ich versuche eine Exponentialkurve an Datensätze mit gedämpften harmonischen Schwingungen anzupassen. Die Daten sind ein bisschen im Sinne kompliziert, dass die Sinusschwingungen viele Frequenzen enthalten, wie unten zu sehen:

Ich brauche die Zerfallsrate in den Daten zu finden. Die Methode, die ich verwende, kann here gefunden werden. Wie es funktioniert, ist es das Protokoll der y-Werte über dem stationären Wert annimmt und dann verwendet:

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

es zu passen.

Dies führt jedoch in der folgenden Daten passt:

Ich habe versucht, eine lineare Regressionsanpassung verwendet, die offensichtlich nicht funktioniert, weil es den Durchschnitt nahm. Ich habe auch RANSAC versucht zu denken, dass es mehr Daten in der Nähe der Spitzen gibt. Es funktionierte ein bisschen besser als die lineare Regression, aber die Methode ist fehlerhaft, da es Zeiten gibt, in denen mehr Punkte in den falschen Regionen existieren.

Kennt jemand eine gute Methode, um die Spitzen für diese Daten gerade zu passen?

Derzeit denke ich darüber nach, die 500 Datenpunkte in 10 verschiedene Regionen zu unterteilen und in jeder Region den größten Wert zu finden. Am Ende sollte ich 50 Punkte haben, die ich mit einer der oben erwähnten Exponentialanpassungsmethoden anpassen kann. Was denkst du über diese Methode?

Answer 1

Ich dachte, ich würde jedem ein Update der möglichen Lösungen geben, die funktionieren könnten. Wie bereits erwähnt, werden die Daten durch die sich ändernden Sinusfrequenzen kompliziert, so dass bestimmte Methoden aus diesem Grund nicht funktionieren. Die unten aufgeführten Methoden können abhängig von den Daten und den beteiligten Frequenzen gut sein.

Zunächst einmal, Ich gehe davon aus, dass die Daten, die die Form hat:

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

In meinem Fall ist die durchschnittliche 290 so haben wir:

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

Mit diesem wird gesagt, wir tauchen ein in die verschiedenen Methoden, die ich versucht:

findpeaks() Methode

Dies ist die Methode, die Alexander Büse vorgeschlagen hat.Es ist eine ziemlich gute Methode für die meisten Daten, aber für meine Daten, da es mehrere sinusförmige Frequenzen gibt, bekommt es die falschen Spitzen. Die roten X zeigen die Peaks.

% 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 ist gut, wenn Sie die meisten Ihrer Daten an den Spitzen haben. Sie sehen, dass in meinem wegen der vielen Frequenzen mehr Spitzen nahe der Spitze existieren. Das Problem mit meinen Daten besteht jedoch darin, dass nicht alle Datensätze so sind. Daher funktionierte es gelegentlich.

% 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 Methode

Diese Methode ist die, here verwendet. Es verwendet Lsqlin, um das System einzuschränken. Es scheint jedoch die Daten in der Mitte zu ignorieren. Abhängig von Ihrem Datensatz könnte dies wirklich gut funktionieren, wie es für die Person im ursprünglichen Beitrag getan hat.

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

Peaks in Periode findet

Dies ist die Methode, die ich in meinem Beitrag erwähnt, wo ich die Spitze in jeder Region erhalten. Diese Methode funktioniert ziemlich gut und daraus habe ich erkannt, dass meine Daten möglicherweise keine perfekte exponentielle Anpassung haben. Wir sehen, dass es nicht in der Lage ist, die großen Spitzen am Anfang anzupassen. Ich konnte das ein bisschen besser machen, indem ich nur die ersten 150 Datenpunkte benutzte und die stationären Datenpunkte ignorierte. Hier fand ich den Peak alle 25 Datenpunkte.

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

mit allen 500 Datenpunkte:

die ersten 150 Datenpunkte verwenden:

Peaks in Periode finden, mit b Constrained

Da ich es will Beginne bei der ersten Spitze, beschränkte ich den b-Wert. Ich weiß, das System ist y=290+b*e^-c*x und ich beschränke es so, dass b=y(1)-290. Dadurch muss ich nur für c wo c=(log(y-290)-logb)/x lösen. Ich kann dann den Durchschnitt oder Median von c nehmen. Diese Methode ist auch ziemlich gut, sie passt auch nicht annähernd zum Ende, aber das ist nicht so groß, da die Veränderung minimal ist.

% 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 

Hier nehme ich den Peak für alle 25 Datenpunkte und nehmen dann den Mittelwert von c

Hier nehme ich den Peak für alle 25 Datenpunkte und nehmen dann den Median c

Hier nehme ich den Peak für alle 10 Datenpunkte und nehmen dann den Mittelwert von c

Answer 2

Wenn das Hauptziel darin besteht, den Dämpfungsparameter aus der Anpassung zu extrahieren, möchten Sie vielleicht in Betracht ziehen, direkt eine gedämpfte Sinuskurve an Ihre Daten anzupassen. So etwas wie dies (mit dem Kurvenanpassungswerkzeug erstellt):

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

Vor allem, da einige Ihrer Beispieldaten wirklich nicht viele Datenpunkte in der interessanten Region hat (über dem Rauschen).

Wenn Sie jedoch direkt an Maxima der experimentellen Daten anpassen müssen, können Sie die Funktion findpeaks verwenden, um nur die Maxima auszuwählen und dann an sie anzupassen. Vielleicht möchten Sie ein wenig mit dem Parameter MinPeakProminence spielen, um ihn an Ihre Bedürfnisse anzupassen.