Find peak positions & plot vline in ggplot

30 views Asked by At

I am plotting some centered signal distribution and I observe peaks that I would like to automatically detect as well as the range.

Here is a snippet of my data with the plot result.

x <- c(-505, -492.278481012658, -479.556962025316, -466.835443037975, 
         -454.113924050633, -441.392405063291, -428.670886075949, -415.949367088608, 
         -403.227848101266, -390.506329113924, -377.784810126582, -365.063291139241, 
         -352.341772151899, -339.620253164557, -326.898734177215, -314.177215189873, 
         -301.455696202532, -288.73417721519, -276.012658227848, -263.291139240506, 
         -250.569620253165, -237.848101265823, -225.126582278481, -212.405063291139, 
         -199.683544303797, -186.962025316456, -174.240506329114, -161.518987341772, 
         -148.79746835443, -136.075949367089, -123.354430379747, -110.632911392405, 
         -97.9113924050633, -85.1898734177215, -72.4683544303797, -59.746835443038, 
         -47.0253164556962, -34.3037974683544, -21.5822784810126, -8.86075949367086, 
         3.86075949367091, 16.5822784810126, 29.3037974683544, 42.0253164556963, 
         54.746835443038, 67.4683544303798, 80.1898734177215, 92.9113924050633, 
         105.632911392405, 118.354430379747, 131.075949367089, 143.79746835443, 
         156.518987341772, 169.240506329114, 181.962025316456, 194.683544303797, 
         207.405063291139, 220.126582278481, 232.848101265823, 245.569620253165, 
         258.291139240506, 271.012658227848, 283.73417721519, 296.455696202532, 
         309.177215189873, 321.898734177215, 334.620253164557, 347.341772151899, 
         360.063291139241, 372.784810126582, 385.506329113924, 398.227848101266, 
         410.949367088608, 423.670886075949, 436.392405063291, 449.113924050633, 
         461.835443037975, 474.556962025316, 487.278481012658, 500)

y <- c(1.08962117485998, 1.04114060875037, 0.996776905965624, 0.960646929830632, 
          0.93686754367026, 0.929555610809385, 0.94282799457288, 0.980801558285618, 
          1.04759316527247, 1.14729099315313, 1.27932537816759, 1.43335917393974, 
          1.5978253344836, 1.76115681381319, 1.91178656594254, 2.03814754488565, 
          2.12867270465656, 2.17179499926927, 2.15634540312282, 2.08322610260535, 
          1.96729487385438, 1.82418687657192, 1.66953727045998, 1.51898121522055, 
          1.38815387055563, 1.29269039616724, 1.24822595175736, 1.2690426113437, 
          1.3536928278639, 1.49058091162332, 1.66794203721677, 1.87401137923908, 
          2.09702411228505, 2.32521541094951, 2.54682044982726, 2.75007440351314, 
          2.92421034174006, 3.06404525918371, 3.16633736840561, 3.22784683098126, 
          3.24533380848613, 3.21555846249569, 3.13528095458543, 3.00126144633082, 
          2.81026468113987, 2.56361890762165, 2.27577932458706, 1.96354702910295, 
          1.6437231182362, 1.33310868905365, 1.04850483862219, 0.806712664008658, 
          0.624533262279933, 0.518563841282694, 0.493168255653025, 0.533748662550571, 
          0.624076105373563, 0.747921627520233, 0.889056272388814, 1.03125108337753, 
          1.15827710388463, 1.25390537730833, 1.30264258298335, 1.30220153657616, 
          1.26171800857497, 1.19070441506743, 1.09867317214123, 0.995136695884031, 
          0.889607402383517, 0.791597707727357, 0.710620028003226, 0.655184811066718, 
          0.625845182720423, 0.619359642721832, 0.632463321306988, 0.661891348711933, 
          0.704378855172708, 0.756660970925356, 0.81547282620592, 0.877549551250441
)

test <- as.data.frame(cbind(x, y))
ggplot(test, aes(x = x,y =  y)) +
  geom_line()+
  theme_bw()

I am currently testing the pracma package which has a findPeaks() function. It looks like the function manage to find the right peak summits however the position of the peaks are odd, I don't see yet how to transform them so I can plot vline on the above plot. See my result :

library(pracma)
peaks <- as.data.frame(pracma::findpeaks(y, npeaks = 3))

ggplot(test,  aes(x = x,y =  y)) +
  geom_line()+
  theme_bw()+
  sapply(peaks$V2, function(xint) geom_vline(aes(xintercept = xint)))

I would understand if the position of peaks was recorded from 0:N positions in the vector but the distance between peaks don't make sense to me.

What I am missing ?

1

There are 1 answers

0
Carl On BEST ANSWER

In peaks V1 gives the height and V2 the row number index. So you could left_join those to test and plot the geom_vline (and/or hline and points).

library(tidyverse)
library(pracma)

test <- as.data.frame(cbind(x, y))

peaks <- as.data.frame(findpeaks(y, npeaks = 3))

test_with_peaks <- test |> 
  mutate(row = row_number()) |> 
  left_join(peaks, join_by(row == V2)) |> 
  mutate(
    peak_x = if_else(!is.na(V1), x, NA),
    peak_y = if_else(!is.na(V1), V1, NA)
    )

ggplot(test_with_peaks,  aes(x = x,y =  y)) +
  geom_line()+
  theme_bw()+
  geom_vline(aes(xintercept = peak_x), alpha = 0.2) +
  geom_hline(aes(yintercept = peak_y), alpha = 0.2) +
  geom_point(aes(peak_x, peak_y), size = 5, colour = "red")

Created on 2024-03-29 with reprex v2.1.0