---
title: "Vertical Jump Data Processing Example"
author: "Colin Bond"
date: "3/22/2022"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

New R users need to install the packages they require. Packages only need to be installed one time, unless a new version of the package is released and you wish to update to the new version. 

signal: Used to filter the ground reaction force data. https://cran.r-project.org/web/packages/signal/index.html
ggplot2: Used to plot results. https://cran.r-project.org/web/packages/ggplot2/index.html
All other functions used in this demonstration are available in Base R, which is what you get with installing R or RStudio.

```{r}
install.packages("signal") #Only run if you have not installed.
install.packages("ggplot2") #Only run if you have not installed.
```

For every R session, you need to "load" the packages that you want to use, assuming they are already installed. You only have to load packages one time per session. 
```{r}
library(signal) #Load package.
library(ggplot2) #Load package.
```

The next step is to import data. The below line creates a data frame called "data" by reading a csv file found at the given path with the first row in the csv file as a header row containing variable names. Note that for windows machines you need to use double "\\" between components in the path. The first 5 rows of data are hten printed to visualize the structure of the data Fz is vertical ground reaction force and the variable of interest in this analysis. 

```{r}
data <- read.csv("C:\\Users\\703354911\\Desktop\\NLACSM Presentation\\F1.csv", header=TRUE) #Loads data from directory. Make sure to change directory to your own.
head(data) #Prints first 6 rows of data.
```

The below lines creates new variables in the data frame called "sample" and "time" to make visualizing and plotting data easier.

```{r}
data$sample <- seq.int(nrow(data)) #Creates a sequential column from 1 to end called "sample".
Hz = 1/1000 #The sampling rate from the force plate used.
data$time <- data$sample*Hz #Creates a sequential column from 0.001 to end called "time".
```

Plot Fz over time to visualize data. Note that you have to tell R to "show" you the plot by printing it on the screen.

```{r}
plot1 <- ggplot(data, aes(x=time, y=Fz)) + geom_line() #Plots Fz over time.
print(plot1) #Prints plot to sceen.
```

Filter the Fz signal to remove some of the noise using the "signal" package. A 4th order 10Hz lowpass Butterworth filter is used to create new variable called "Fz_filt" in the data frame. 

```{r}
bf <- butter(4, 0.1, type="low") #Sets up 4th order 10Hz lowpass Butterworth filter
data$Fz_filt <- filter(bf, data$Fz) #Applies filter.
```

Plot the filtered signal to verify it looks appropriate compared to the unfiltered signal.

```{r}
plot2 <- ggplot(data, aes(x=time, y=Fz)) + geom_line(color="red") + geom_line(inherit.aes=FALSE,data=data,aes(x=time, y=Fz_filt), color="blue") #Plots Fz over time in red and Fz_filt over time in blue.
print(plot2) #Prints plot to screen.
```

Determine the subjects body weight in newtons (BW) using the quite stance period at the beginning of the trial. Then subtract BW from Fz_fil to remove impulse caused by body weight.

```{r}
print(BW <- mean(data[100:500,"Fz_filt"])) #BW is the mean of Fz_filt from 0.1 to 0.5 seconds.
data$Fz_filt_norm <- data$Fz_filt - BW #Subtracts BW from Fz_filt and creates new variable called Fz_filt_norm in data frame.
```

Plot Fz_filt_norm over time to visually inspect it.

```{r}
plot3 <- ggplot(data, aes(x=time, y=Fz_filt_norm)) + geom_line() #Plots Fz_filt_norm over time. 
print(plot3) #Prints plot to screen.
```

Plot the same data, but apply different colors depending on if the "Fz_filt_norm" signal is greater or less than 0, indicating positive and negative impulse, respectively. 

```{r}
plot4 <- ggplot(data, aes(x=time, y=Fz_filt_norm,ymin = 0,ymax = Fz_filt_norm)) + geom_line() + geom_linerange(data = data, aes(colour = ifelse(Fz_filt_norm <0, "Negative Impulse", "Positive Impulse")),stat = "identity", position = "identity",size=0.5) + scale_color_discrete(name = "Impulse") #Plots Fz_filt_norm over time and applies different colors depending on if it is greater or less than 0.
print(plot4) #Print plot to screen.
```

Mark the rows in the data frame where certain events occur. 

```{r}
begin_negative_impulse = min(which(data[100:nrow(data),"Fz_filt_norm"] < -20))+100 #The beginning of the counter movement. Looks beginning at the 100th sample (0.1 s) and finds the first sample where Fz_filt_norm goes below -20 N. 20 or -20 N is commonly used threshold to distinguish an event, but this could be adjusted depending on the inherit noise in a given force plate.
begin_negative_impulse #Print sample to screen.
begin_positive_impulse = min(which(data[begin_negative_impulse:nrow(data),"Fz_filt_norm"] > 20))+begin_negative_impulse #The end of the counter movement. Looks beginning at the sample identified as begin_negative_impulse and finds the first sample where Fz_filt_norm goes above 20 N.
begin_positive_impulse #Print sample to screen.
end_positive_impulse = min(which(data[begin_positive_impulse:nrow(data),"Fz_filt_norm"] < 0))+begin_positive_impulse #The end of the positive impulse. Looks beginning at the sample identified as begin_positive_impulse and finds the first sample where Fz_filt_norm goes below 0 N.
end_positive_impulse #Print sample to screen.
```

Impulse (Ns) = sum of [Fz_filt_norm*t] from the beginning of positive impulse to the end of positive impulse where t equals our sampling rate (1000).

```{r}
Impulse = sum(data[begin_positive_impulse:end_positive_impulse,"Fz_filt_norm"]*Hz) #Calculates the impulse, or integral of force, from the beginning to end of the positive impulse.
Impulse #Prints impulse to screen.
```

Takeoff_Velocity (m/s) = Impulse /(BW/9.81)

```{r}
Takeoff_Velocity = Impulse/(BW/9.81) #Calculates the take off velocity.
Takeoff_Velocity #Prints take off velocity to screen.
```

Jump Height (m) = 0.5 * (Takeoff_Velocity^2/9.81)

```{r}
Jump_Height = 0.5 * (Takeoff_Velocity^2/9.81) #Calculated jump height in m.
Jump_Height #Prints jump height in m to screen.
Jump_Height_in = Jump_Height*39.37 #Converts jump height in m to jump height in inches.
Jump_Height_in #Prints jump height in inches to screen.
```

