find_peaks.m

function [X_peaks, Y_peaks, sums, start_positions, start_values, end_positions, end_values, detection_positions, detection_values, slope_positions, slope_values, end_bounce_positions, end_bounce_values, in_bounce_positions, in_bounce_values] = find_peaks(X, Y, thresholdY, thresholdY_second)
  X_peaks = zeros(0, 1);
  Y_peaks = zeros(0, 1);
  sums = zeros(0, 1);
  start_positions = zeros(0, 1);
  end_positions = zeros(0, 1);
  start_values = zeros(0, 1);
  end_values = zeros(0, 1);
  end_bounce_positions = zeros (0, 1);
  end_bounce_values = zeros(0, 1);
  in_bounce_positions = zeros(0, 1);
  in_bounce_values = zeros(0, 1);
  
  detection_positions = zeros(0, 1);
  detection_values = zeros(0, 1);
  slope_positions = zeros(0, 1);
  slope_values = zeros(0, 1);
  
  Y_prime = derivatives(X, Y);
  Y_second = derivatives2(X, Y);
  
  m = mean(Y(1:200));
  m1 = mean(Y_prime(1:200));
  disp(m1);
  
  L = length(Y);
  disp(L);
  
  state = 0;
  
  current_sum = 0;
  current_peak_position = 0;
  current_peak_value = 0;
  current_start_position = 0;
  current_start_value = 0;
  current_end_position = 0; 
  current_end_value = 0;
  last_bounce_position = 0;
  
 
  for i=1:L;
    switch (state)
      case 0
        if (Y(i) > thresholdY) && (Y_second(i) > thresholdY_second)
          current_sum = Y(i);
          j = i - 1;
          while (j >= 1 && X(j) > current_end_position && X(j) > last_bounce_position && abs(Y_prime(j) - m1) > 2 * abs(m1))
            current_sum += Y(j) - m;
            current_start_position = X(j);
            current_start_value = Y(j);
            j--;
          endwhile
          current_peak_position = X(i);
          current_peak_value = Y(i);
          detection_positions = [detection_positions [X(i)]];
          detection_values = [detection_values [Y(i)]];
          state = 1;
          disp("state 1");
        endif
        %break;
      case 1
        current_sum += Y(i) - m;
        if (current_peak_value < Y(i))
          current_peak_value = Y(i);
          current_peak_position = X(i);
        endif
        if (Y_second(i) < 0 || abs(Y_prime - m1) < 2 * m1)
          slope_positions = [slope_positions [X(i)]];
          slope_values = [slope_values [Y(i)]];
          state = 2;
          disp("state 2");
        endif
        %break;
      case 2 
        current_sum += Y(i) - m;
        if (current_peak_value < Y(i))
          current_peak_value = Y(i);
          current_peak_position = X(i);
        endif
        deriv1_mean = mean(Y_prime(i - 4:i));
        if (Y(i) - m < 0 || abs(deriv1_mean) <= 2 * abs(m1)) %the curve goes beyond the mean (downwards)
          current_end_position = X(i);
          current_end_value = Y(i);
          state = 3;
          disp("state 3");
        endif
        %break;
      case 3
        deriv1_mean = mean(Y_prime(i - 10: i));
        if (deriv1_mean > 2 * m1)
          in_bounce_positions = [in_bounce_positions [X(i)]];
          in_bounce_values = [in_bounce_values [Y(i)]];
          state = 4;
          disp("state 4");
        endif
      case 4
##        sum_mean = 0;
##        for k = i:-1:i+1-10;
##          sum_mean += abs(Y_prime(k) - m1);
##        endfor
##        sum_mean /= k;
        deriv1_mean = mean(Y_prime(i - 4:i));
        if (abs(current_peak_value - Y(i)) / Y(i) < 0.5 && Y(i) > thresholdY && Y_second(i) > thresholdY_second) %another peak was found
          X_peaks = [X_peaks [current_peak_position]];
          Y_peaks = [Y_peaks [current_peak_value]];
          sums = [sums [current_sum]];
          start_positions = [start_positions [current_start_position]];
          start_values = [start_values [current_start_value]];
          end_positions = [end_positions [current_end_position]];
          end_values = [end_values [current_end_value]];
          end_bounce_positions = [end_bounce_positions [X(i)]];
          end_bounce_values = [end_bounce_values [Y(i)]];
          last_bounce_position = X(i);
          current_sum = Y(i);
          j = i - 1;
          last_bounce_position = current_end_position;
          last_bounce_value = current_end_value;
          while (X(j) > current_end_position && X(j) > last_bounce_position)
            current_sum += Y(j);
            current_start_position = X(j);
            current_start_value = Y(j);
            j--;
          endwhile       
          current_peak_position = X(i);
          current_peak_value = Y(i);  
          detection_positions = [detection_positions [X(i)]];
          detection_values = [detection_values [Y(i)]];
          state = 1;
          disp("state 1");
        elseif ((abs(deriv1_mean) <= 2 * abs(m1)) || ((mean(Y_second(i-4:i)) <= 0) && (mean(Y_second(i-8:i-5)) > 0)))
          X_peaks = [X_peaks [current_peak_position]];
          Y_peaks = [Y_peaks [current_peak_value]];
          sums = [sums [current_sum]];
          start_positions = [start_positions [current_start_position]];
          start_values = [start_values [current_start_value]];
          end_positions = [end_positions [current_end_position]];
          end_values = [end_values [current_end_value]];
          end_bounce_positions = [end_bounce_positions [X(i)]];
          end_bounce_values = [end_bounce_values [Y(i)]];
          last_bounce_position = X(i);
          state = 0;
          disp("state 0");
        endif
        %break;
    endswitch
      
  endfor
endfunction