waves.Rd

% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/FitWave.R
\name{waves}
\alias{waves}
\alias{FitWave}
\alias{BuildWave}
\alias{FilterWave}
\alias{WaveEnvelope}
\title{Fourier transform functions}
\usage{
FitWave(y, k = 1)

BuildWave(
  x,
  amplitude,
  phase,
  k,
  wave = list(amplitude = amplitude, phase = phase, k = k),
  sum = TRUE
)

FilterWave(y, k, action = sign(k[k != 0][1]))

WaveEnvelope(y)
}
\arguments{
\item{y}{numeric vector to transform}

\item{k}{numeric vector of wave numbers}

\item{x}{numeric vector of locations (in radians)}

\item{amplitude}{numeric vector of amplitudes}

\item{phase}{numeric vector of phases}

\item{wave}{optional list output from \code{FitWave}}

\item{sum}{whether to perform the sum or not (see Details)}

\item{action}{integer to disambiguate action for k = 0 (see Details)}
}
\value{
\code{FitWaves} returns a a named list with components
\describe{
\item{k}{wavenumbers}
\item{amplitude}{amplitude of each wavenumber}
\item{phase}{phase of each wavenumber in radians}
\item{r2}{explained variance of each wavenumber}
}

\code{BuildWave} returns a vector of the same length of x with the reconstructed
vector if \code{sum} is \code{TRUE} or, instead, a list with components
\describe{
\item{k}{wavenumbers}
\item{x}{the vector of locations}
\item{y}{the reconstructed signal of each wavenumber}
}

\code{FilterWave} and \code{WaveEnvelope} return a vector of the same length as \code{y}
`
}
\description{
Use \code{\link[=fft]{fft()}} to fit, filter and reconstruct signals in the frequency domain, as
well as to compute the wave envelope.
}
\details{
\code{FitWave} performs a fourier transform of the input vector
and returns a list of parameters for each wave number kept.
The  amplitude (A), phase (\eqn{\phi}) and wave number (k) satisfy:
\deqn{y = \sum A cos((x - \phi)k)}
The phase is calculated so that it lies between 0 and \eqn{2\pi/k} so it
represents the location (in radians) of the first maximum of each wave number.
For the case of k = 0 (the mean), phase is arbitrarily set to 0.

\code{BuildWave} is \code{FitWave}'s inverse. It reconstructs the original data for
selected wavenumbers. If \code{sum} is \code{TRUE} (the default) it performs the above
mentioned sum and returns a single vector. If is \code{FALSE}, then it returns a list
of k vectors consisting of the reconstructed signal of each wavenumber.

\code{FilterWave} filters or removes wavenumbers specified in \code{k}. If \code{k} is positive,
then the result is the reconstructed signal of \code{y} only for wavenumbers
specified in \code{k}, if it's negative, is the signal of \code{y} minus the wavenumbers
specified in \code{k}. The argument \code{action} must be be manually set to \code{-1} or \code{+1}
if \code{k=0}.

\code{WaveEnvelope} computes the wave envelope of \code{y} following Zimin (2003). To compute
the envelope of only a restricted band, first filter it with \code{FilterWave}.
}
\examples{
\dontshow{data.table::setDTthreads(1)}

# Build a wave with specific wavenumber profile
waves <- list(k = 1:10,
              amplitude = rnorm(10)^2,
              phase = runif(10, 0, 2*pi/(1:10)))
x <- BuildWave(seq(0, 2*pi, length.out = 60)[-1], wave = waves)

# Just fancy FFT
FitWave(x, k = 1:10)

# Extract only specific wave components
plot(FilterWave(x,  1), type = "l")
plot(FilterWave(x,  2), type = "l")
plot(FilterWave(x,  1:4), type = "l")

# Remove components from the signal
plot(FilterWave(x,  -4:-1), type = "l")

# The sum of the two above is the original signal (minus floating point errors)
all.equal(x, FilterWave(x,  1:4) + FilterWave(x,  -4:-1))

# The Wave envelopes shows where the signal is the most "wavy".
plot(x, type = "l", col = "grey")
lines(WaveEnvelope(x), add = TRUE)

# Examples with real data
data(geopotential)
library(data.table)
# January mean of geopotential height
jan <- geopotential[month(date) == 1, .(gh = mean(gh)), by = .(lon, lat)]

# Stationary waves for each latitude
jan.waves <- jan[, FitWave(gh, 1:4), by = .(lat)]
library(ggplot2)
ggplot(jan.waves, aes(lat, amplitude, color = factor(k))) +
    geom_line()

# Build field of wavenumber 1
jan[, gh.1 := BuildWave(lon*pi/180, wave = FitWave(gh, 1)), by = .(lat)]
ggplot(jan, aes(lon, lat)) +
    geom_contour(aes(z = gh.1, color = after_stat(level))) +
    coord_polar()

# Build fields of wavenumber 1 and 2
waves <- jan[, BuildWave(lon*pi/180, wave = FitWave(gh, 1:2), sum = FALSE), by = .(lat)]
waves[, lon := x*180/pi]
ggplot(waves, aes(lon, lat)) +
    geom_contour(aes(z = y, color = after_stat(level))) +
    facet_wrap(~k) +
    coord_polar()

# Field with waves 0 to 2 filtered
jan[, gh.no12 := gh - BuildWave(lon*pi/180, wave = FitWave(gh, 0:2)), by = .(lat)]
ggplot(jan, aes(lon, lat)) +
    geom_contour(aes(z = gh.no12, color = after_stat(level))) +
    coord_polar()

# Much faster
jan[, gh.no12 := FilterWave(gh, -2:0), by = .(lat)]
ggplot(jan, aes(lon, lat)) +
    geom_contour(aes(z = gh.no12, color = after_stat(level))) +
    coord_polar()

# Using positive numbers returns the field
jan[, gh.only12 := FilterWave(gh, 2:1), by = .(lat)]
ggplot(jan, aes(lon, lat)) +
    geom_contour(aes(z = gh.only12, color = after_stat(level))) +
    coord_polar()

# Compute the envelope of the geopotential
jan[, envelope := WaveEnvelope(gh.no12), by = .(lat)]
ggplot(jan[lat == -60], aes(lon, gh.no12)) +
    geom_line() +
    geom_line(aes(y = envelope), color = "red")

}
\references{
Zimin, A.V., I. Szunyogh, D.J. Patil, B.R. Hunt, and E. Ott, 2003: Extracting Envelopes of Rossby Wave Packets. Mon. Wea. Rev., 131, 1011–1017, \doi{10.1175/1520-0493(2003)131<1011:EEORWP>2.0.CO;2}
}
\seealso{
Other meteorology functions: 
\code{\link{Derivate}()},
\code{\link{EOF}()},
\code{\link{GeostrophicWind}()},
\code{\link{WaveFlux}()},
\code{\link{thermodynamics}}
}
\concept{meteorology functions}