if( !require('quantmod') ){ # I want an implementation of 'Lag'
    install.packages('quantmod')
}
library(quantmod)


my_amplitude = function(u){
    ## A function to extract the approximate amplitude of each sinusoid:
    ##

    ## Starting at the end of the series find the largest value of u(t) by walking left:
    ##
    spot = length(u)
    un = u[spot]
    unm1 = u[spot-1]
    ## While decreasing as we walk left:
    while(unm1 < un){
        spot = spot - 1
        un = u[spot]
        unm1 = u[spot-1]
    }
    ## While increasing as we walk left:
    while(unm1 > un){
        spot = spot - 1
        un = u[spot]
        unm1 = u[spot-1]
    }
    highest_value = mean(c(un, unm1))

    ## Starting at this value find the smallest value of u(t) by walking left:
    ##
    ## While decreasing as we walk left:
    while(unm1 < un){
        spot = spot - 1
        un = u[spot]
        unm1 = u[spot-1]
    }
    lowest_value = mean(c(un, unm1))

    A = 0.5 * ( highest_value - lowest_value )
    
    ##return(c(lowest_value, highest_value, A))
    return(A)
}


my_period = function(ts, u) {
    ## A function to extract the approximate period of each sinusoid:
    ##

    ## Demean the input:
    ##
    u = u - mean(u)

    ## Find the indices where the sinusoid crosses the line y=0:
    ##
    crossings = ( Lag(u, k=0) * Lag(u, k=1) < 0 )
    crossing_times = ts[crossings]
    period = 2*mean(diff(crossing_times), na.rm=TRUE)
    return(period)
}