español - ¿Cómo adaptar una curva exponencial a los datos de oscilación armónica amortiguada en MATLAB?

polyfit matlab español

Question

Estoy tratando de ajustar una curva exponencial a los conjuntos de datos que contienen oscilaciones armónicas amortiguadas. Los datos son un poco complicados en el sentido de que las oscilaciones sinusoidales contienen muchas frecuencias como se ve a continuación:

Necesito encontrar la tasa de descomposición en los datos. El método que estoy usando se puede encontrar aquí . Cómo funciona, es que toma el registro de los valores de y por encima del valor de estado estable y luego utiliza:

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

Para encajarlo

Sin embargo, esto da como resultado los siguientes ajustes de datos:

Traté de usar un ajuste de regresión lineal que, obviamente, no funcionó porque tomó el promedio. También probé RANSAC pensando que hay más datos cerca de los picos. Funcionó un poco mejor que la regresión lineal, pero el método es defectuoso ya que hay momentos en que existen más puntos en las regiones incorrectas.

¿Alguien sabe de un buen método para ajustar los picos de estos datos?

Actualmente, estoy pensando en dividir los 500 puntos de datos en 10 regiones diferentes y en cada región encontrar el valor más grande. Al final, debería tener 50 puntos que pueda encajar utilizando cualquiera de los métodos de adaptación exponencial mencionados anteriormente. ¿Qué piensas de este método?

Answer 1

Si el objetivo principal es extraer el parámetro de amortiguación del ajuste, tal vez desee considerar ajustar directamente una curva sinusoidal amortiguada a sus datos. Algo así (creado con la herramienta de ajuste de 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 );

Especialmente porque algunos de sus datos de ejemplo realmente no tienen muchos puntos de datos en la región interesante (por encima del ruido).

Sin embargo, si realmente necesita ajustarse directamente al máximo de los datos experimentales, puede usar la función findpeaks para seleccionar solo los máximos y luego ajustarlos. Es posible que desee jugar un poco con el parámetro MinPeakProminence para ajustarlo a sus necesidades.

Answer 2

Pensé que le daría a todos una actualización de posibles soluciones que pueden funcionar. Como se mencionó anteriormente, los datos son complicados por las frecuencias sinusoidales variables, por lo que ciertos métodos pueden no funcionar debido a esto. Los métodos enumerados a continuación pueden ser buenos según los datos y las frecuencias involucradas.

Primero, supongo que los datos tienen la forma:

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

En mi caso, el promedio es 290, entonces tenemos:

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

Con eso dicho, profundicemos en los diferentes métodos que probé:

Método findpeaks ()

Este es el método que sugirió Alexander Büse. Es un método bastante bueno para la mayoría de los datos, pero para mis datos, dado que hay múltiples frecuencias sinusoidales, obtiene los picos incorrectos. Las x rojas muestran los picos.

% 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 es bueno si tiene la mayoría de sus datos en los picos. Ya ves que en el mío, debido a las múltiples frecuencias, existen más picos cerca de la parte superior. Sin embargo, el problema con mis datos es que no todos los conjuntos de datos son así. Por lo tanto, ocasionalmente funcionó.

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

Método Lsqlin

Este método es el que se usa aquí . Utiliza Lsqlin para restringir el sistema. Sin embargo, parece ignorar los datos en el medio. Dependiendo de su conjunto de datos, esto podría funcionar muy bien como lo hizo para la persona en la publicación original.

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

Encuentra picos en el período

Este es el método que mencioné en mi publicación donde obtengo el pico en cada región,. Este método funciona bastante bien y de esto me di cuenta de que mis datos pueden no tener un ajuste exponencial perfecto. Vemos que no es capaz de adaptarse a los grandes picos al principio. Pude mejorar esto un poco utilizando solo los primeros 150 puntos de datos e ignorando los puntos de datos de estado estacionario. Aquí encontré el pico cada 25 puntos de datos.

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

Usando todos los 500 puntos de datos:

Usando los primeros 150 puntos de datos:

Encontrar picos en el período con b restringido

Como quiero que comience en el primer pico, restringí el valor b. Sé que el sistema es y=290+b*e^-c*x y lo constriño tal que b=y(1)-290 . Al hacerlo, solo necesito resolver c donde c=(log(y-290)-logb)/x . Entonces puedo tomar el promedio o la mediana de c. Este método también es bastante bueno, no se ajusta al valor cerca del final, pero eso no es tan importante, ya que el cambio es mínimo.

% 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

Aquí tomo el pico por cada 25 puntos de datos y luego tomo el promedio de c

Aquí tomo el pico por cada 25 puntos de datos y luego tomo la mediana de c

Aquí tomo el pico por cada 10 puntos de datos y luego tomo el promedio de c