Как установить экспоненциальную кривую для данных затухающих гармонических колебаний в MATLAB?

Question

Я пытаюсь подогнать экспоненциальную кривую к наборам данных, содержащим затухающие гармонические колебания. Данные немного сложнее в том смысле, что синусоидальные колебания содержат много частот, как показано ниже:

Мне нужно найти скорость распада в данных. Метод, который я использую, можно найти here. Как это работает, он берет журнал значений у выше стационарного значения, а затем использует:

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() метод

Это метод, предложенный Александром Бюсом.Это очень хороший метод для большинства данных, но для моих данных, поскольку есть несколько синусоидальных частот, он получает неправильные пики. Красные 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. Я знаю, что система y=290+b*e^-c*x, и я ограничиваю ее таким, что b=y(1)-290. Таким образом, мне просто нужно решить для c, где 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 точек данных, а затем взять среднее с

Здесь я беру пик на каждые 25 точек данных, а затем взять медиану с

Здесь я беру пик на каждые 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, чтобы настроить его на свои нужды.