How to set an exponential curve for data of damped harmonic oscillations in MATLAB?

Question

I am trying to fit an exponential curve to datasets containing damped harmonic oscillations. The data is a little complicated in the sense that sinusoidal oscillations contain many frequencies, as shown below:

I need to find the decay rate in the data. The method I use can be found here . How it works, whether the log takes y values ​​above a steady state value, and then uses:

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

To fit that.

However, this leads to the following data:

I tried using linear regression, which obviously did not work, because it took an average value. I also tried RANSAC to think that there is more data near the peaks. It worked a little better than linear regression, but the method is wrong, as there are times when there are more points in the wrong regions.

Does anyone know of a good method to just match the peaks for this data?

Currently, I am going to divide 500 data points into 10 different regions and find the highest value in each region. In the end, I should have 50 points that I can put using any of the exponential fitting methods mentioned above. What do you think of this method?

Answer 1

Thought I'd give everyone updated information on possible solutions that might work. As mentioned earlier, data is complicated by variable sinusoidal frequencies, so some methods may not work because of this. The methods listed below may be good depending on the data and frequencies used.

First of all, I assume that the data is as follows:

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

In my case, the average is 290, so we have:

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

With that said, let's dive into the different methods I've tried:

findpeaks () method

This is the method proposed by Alexander Buz. This is a pretty good method for most data, but for my data, since there are several sinusoidal frequencies, it gets the wrong peaks. Red x indicates 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 is good if you have most of your data at peak. You see that in mine, due to several frequencies, more peaks exist near the peak. However, the problem with my data is that not all data sets are like this. Therefore, it sometimes worked.

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

This method is used here here . It uses Lsqlin to limit the system. However, they seem to ignore the data in the middle. Depending on your dataset, this may work very well, as it was for the person in the original message.

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

Find peaks in the period

This is the method that I mentioned in my post where I get a peak in each region. This method works very well, and from this I realized that my data may not have perfect exponential correspondence. We see that it cannot correspond to large peaks at the beginning. I was able to do this a little better, using only the first 150 data points and ignoring the steady state data points. Here I found a peak every 25 data points.

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

Using all 500 data points:

Using the first 150 data points:

Find peaks between b Limited

Since I want it to start from the first peak, I limited the value of b. I know that the system is y=290+b*e^-c*x , and I limit it to b=y(1)-290 . In doing so, I just need to solve for c, where c=(log(y-290)-logb)/x . Then I can take the average or average c. This method is also good, it also does not approach the value near the end, but it is not so expensive, since the change is minimal.

 % 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 

Here I take a peak for every 25 data points, and then take the average c

Here I take a peak for every 25 data points, and then take the median c

Here I take a peak for every 10 data points and then take the average c

Answer 2

If the main goal is to extract the damping parameter from the fit, you might want to consider fitting directly to the decaying sine curve of your data. Something like this (created using a curve fitting tool):

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

Moreover, some of the data in your example really does not have a large number of data points in an interesting area (above noise).

However, if you really need to fit directly to the maxima of the experimental data, you can use the findpeaks function to select only the maxima, and then approach them. You can play a little with the MinPeakProminence parameter to customize it to your needs.