Jak dopasować krzywą wykładniczą do tłumionych danych oscylacji harmonicznych w MATLAB?

Question

Próbuję dopasować wykładniczą krzywą do zestawów danych zawierających tłumione oscylacje harmoniczne. Dane są nieco skomplikowane, w tym sensie, że drgania sinusoidalne zawierają wiele częstotliwości, jak widać poniżej:

Muszę znaleźć szybkość rozkładu w danych. Metoda, której używam, znajduje się pod adresem here. Jak to działa, czy loguje wartości y powyżej wartości ustalonej, a następnie używa:

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

Aby to pasowało.

jednak skutkuje to wpisuje się następujące dane:

Próbowałem za pomocą regresji liniowej dopasowanie który oczywiście nie działa, ponieważ wziął średnią. Próbowałem też RANSAC, myśląc, że w pobliżu szczytów jest więcej danych. To działało nieco lepiej niż regresja liniowa, ale metoda jest wadliwa, ponieważ zdarzają się sytuacje, gdy w niewłaściwych regionach występuje więcej punktów.

Czy ktoś wie o dobrej metodzie, aby po prostu dopasować wartości szczytowe dla tych danych?

Obecnie myślę o podzieleniu 500 punktów danych na 10 różnych regionów iw każdym regionie znajduję największą wartość. Na koniec powinienem mieć 50 punktów, które umiem dopasować przy użyciu dowolnej z wyżej wymienionych metod wykładniczych. Co myślisz o tej metodzie?

Answer 1

Pomyślałem, że dam wszystkim aktualizację potencjalnych rozwiązań, które mogą działać. Jak wspomniano wcześniej, dane są komplikowane przez zmienne częstotliwości sinusoidalne, więc niektóre metody mogą z tego powodu nie działać. Wymienione poniżej metody mogą być dobre w zależności od danych i częstotliwości.

Po pierwsze, zakładam, że dane mają postać:

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

w moim przypadku, średnia wynosi 290, więc mamy:

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

Z tym mówi, niech zanurzyć różne metody, które próbowałem:

findpeaks() Metoda

Jest to metoda zalecana przez Alexandra Büse.Jest to całkiem dobra metoda dla większości danych, ale dla moich danych, ponieważ istnieje wiele częstotliwości sinusoidalnych, otrzymuje błędne wartości szczytowe. Czerwone X pokazują szczyty.

% 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 jest dobre, jeśli masz większość danych na szczyty. Widzisz to w moim, ze względu na wiele częstotliwości, więcej szczytów istnieje blisko szczytu. Jednak problem z moimi danymi jest taki, że nie wszystkie zestawy danych są takie. Dlatego od czasu do czasu zadziałało.

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

Metoda Lsqlin

Ta metoda jest jeden używany here. Używa Lsqlin do ograniczenia systemu. Jednak wydaje się, że ignoruje dane w środku. W zależności od zestawu danych może to działać naprawdę dobrze, jak w przypadku osoby w oryginalnym wpisie.

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

Znajdź Szczyty w Okresie

Jest to metoda wspominałem w moim poście gdzie mogę dostać się na szczyt w każdym regionie. Ta metoda działa całkiem dobrze iz tego powodu zdałem sobie sprawę, że moje dane mogą nie mieć doskonałego dopasowania wykładniczego. Widzimy, że nie jest on w stanie dopasować dużych szczytów na początku. Udało mi się to nieco poprawić, wykorzystując tylko pierwsze 150 punktów danych i ignorując punkty danych stanu ustalonego. Tutaj znalazłem szczyt co 25 punktów danych.

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

wykorzystując wszystkie 500 punktów danych:

Korzystanie pierwszych 150 punktów danych:

Znajdź Szczyty w Okresie Z b Ograniczone

ponieważ chcę go zacznij od pierwszego szczytu, ograniczyłem wartość b. Wiem, że system jest y=290+b*e^-c*x i ograniczam go tak, że b=y(1)-290. W ten sposób, po prostu muszę rozwiązać c, gdzie c=(log(y-290)-logb)/x. Mogę wtedy wziąć średnią lub medianę c. Ta metoda jest również dobra, nie pasuje również do wartości zbliżonej do końca, ale nie jest to tak duża oferta, ponieważ zmiana jest minimalna.

% 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 

Oto biorę szczyt na 25 punktów danych, a następnie podjąć średnio ok

Oto biorę szczyt na 25 punktów danych, a następnie podjąć medianę c

Oto biorę szczyt na każde 10 punktów danych, a następnie podjąć średnio ok

Answer 2

Jeśli głównym celem jest wyodrębnienie parametru tłumienia z dopasowania, może warto rozważyć dopasowanie bezpośrednio krzywych sinusoidalnych do danych. Coś takiego (utworzony za pomocą narzędzia dopasowania krzywej):

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

Zwłaszcza, że ​​niektóre z przykładowych danych naprawdę nie mają wiele punktów danych w regionie ciekawe (powyżej poziomu hałasu).

Jeśli jednak naprawdę potrzebujesz dopasować bezpośrednio do maksimów danych eksperymentalnych, możesz użyć funkcji findpeaks, aby wybrać tylko maksima, a następnie dopasować do nich. Możesz chcieć zagrać trochę z parametrem MinPeakProminence, aby dostosować go do swoich potrzeb.