Come adattare una curva esponenziale ai dati di oscillazione armonica smorzata in MATLAB?

Question

Sto provando ad adattare una curva esponenziale a serie di dati contenenti oscillazioni armoniche smorzate. Il dato è un po 'complicato, nel senso che le oscillazioni sinusoidali contengono molte frequenze, come si vede qui sotto:

ho bisogno di trovare il tasso di decadimento nei dati. Il metodo che sto usando può essere trovato here. Come funziona, è che ci vuole il registro dei valori y al di sopra del valore di regime e quindi utilizza:

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

per adattarlo.

Tuttavia, questo si traduce nei seguenti dati si inserisce:

Ho provato ad utilizzare una forma di regressione lineare che ovviamente non ha funzionato perché ci sono voluti la media. Ho anche provato RANSAC a pensare che ci siano più dati vicino ai picchi. Ha funzionato un po 'meglio della regressione lineare, ma il metodo è imperfetto perché ci sono momenti in cui più punti esistono nelle regioni sbagliate.

Qualcuno sa di un buon metodo per montare i picchi per questi dati?

Attualmente, sto pensando di dividere i 500 punti dati in 10 diverse regioni e in ogni regione trovare il valore più grande. Alla fine, dovrei avere 50 punti che posso adattare utilizzando uno dei metodi di adattamento esponenziale sopra menzionati. Cosa ne pensi di questo metodo?

Answer 1

Pensavo di dare a tutti un aggiornamento delle possibili soluzioni che potrebbero funzionare. Come accennato in precedenza, i dati sono complicati dalle frequenze sinusoidali variabili, quindi alcuni metodi potrebbero non funzionare a causa di ciò. I metodi elencati di seguito possono essere buoni a seconda dei dati e delle frequenze coinvolte.

Prima di tutto, suppongo che i dati ha la forma:

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

Nel mio caso, la media è 290 quindi abbiamo:

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

Detto questo, cerchiamo di tuffarsi i diversi metodi che ho provato:

findpeaks() Metodo

Questo è il metodo suggerito da Alexander Büse.È un metodo abbastanza buono per la maggior parte dei dati, ma per i miei dati, dato che ci sono più frequenze sinusoidali, si ottengono i picchi errati. Le x rosse mostrano le vette.

% 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 è buono se si dispone di maggior parte dei dati in corrispondenza dei picchi. Lo vedi nel mio, a causa delle frequenze multiple, più picchi esistono vicino alla cima. Tuttavia, il problema con i miei dati è che non tutti i set di dati sono così. Quindi, occasionalmente ha funzionato.

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

Metodo Lsqlin

Questo metodo è quello utilizzato here. Usa Lsqlin per vincolare il sistema. Tuttavia, sembra ignorare i dati nel mezzo. A seconda del set di dati, questo potrebbe funzionare molto bene come ha fatto per la persona nel post originale.

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

Trova Picchi nella Periodo

questo è il metodo che ho citato nel mio post dove ho il picco in ogni regione,. Questo metodo funziona abbastanza bene e da questo ho capito che i miei dati potrebbero non avere una perfetta corrispondenza esponenziale. Vediamo che non è in grado di adattarsi ai grandi picchi all'inizio. Sono riuscito a migliorarlo un po 'usando solo i primi 150 punti dati e ignorando i punti dati a stato stazionario. Qui ho trovato il picco ogni 25 punti dati.

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

utilizzando tutti 500 punti dati:

utilizzando i primi 150 punti dati:

Trova Picchi in epoca con b vincolata

Dal momento che voglio che iniziare al primo picco, ho limitato il valore b. So che il sistema è y=290+b*e^-c*x e lo vincolo in modo tale che b=y(1)-290. In tal modo, ho solo bisogno di risolvere per c dove c=(log(y-290)-logb)/x. Posso quindi prendere la media o mediana di c. Questo metodo è abbastanza buono, non si adatta anche al valore vicino alla fine, ma non è un grosso problema dato che il cambiamento è minimo.

% 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 

Qui prendo il picco per ogni 25 punti di dati e poi prendere la media dei c

Qui prendo il picco per ogni 25 punti di dati e poi prendere il mediano di c

Qui prendo il picco per ogni 10 punti di dati e poi prendere la media dei c

Answer 2

Se l'obiettivo principale è estrarre il parametro di smorzamento dall'adattamento, forse si desidera considerare l'adattamento diretto di una curva sinusoidale ai dati. Qualcosa di simile (creato con lo strumento di montaggio di curva):

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

Tanto più che alcuni dei vostri dati di esempio in realtà non hanno molti punti di dati nella regione interessante (al di sopra del rumore).

Se tuttavia è necessario adattarsi direttamente ai valori massimi dei dati sperimentali, è possibile utilizzare la funzione findpeaks per selezionare solo i massimi e quindi adattarli. Potresti voler giocare un po 'con il parametro MinPeakProminence per adattarlo alle tue esigenze.