# R/hilbert.r In mccreigh/EOF: EOF routines for climate data.

```## discrete time "analytic signal" via hilbert transform (via fft)

## input: real timeseries or matrix with time along columns (nspaceXntime).
##        any imaginary component is discarded.
## output: the analytic signal = input + i*HILBERT(input) for each point/timeseries,
## a complex spaceXtime result.

## The hilbert transform gives a quarter phase shift to the input timeseries.

## caveat emptor: no warranty whatsoever, use this code at your own risk.
## author: james mccreight # gmail o com

## reference:
## S. Lawrence Marple, Jr., Computing the discrete-time analytic
##     signal via FFT, IEEE Transactions on Signal Processing, Vol. 47,
##     No. 9, September 1999, pp.2600--2603.

## (this can be verified against hilbert() for univariate ts in the EMD package,
##   if that's any solace.)

hilbert <- function (input)
{

if (!is.real(input)) {
print("*** Warning: discaring imaginary part of input.")
input <- Re(input)
}

if (is.array(input)) {
nspace <- dim(as.matrix(input))[1]
ntime <- dim(as.matrix(input))[2]
} else {
ntime <- length(input)
nspace <- 1
input=t(input)
}

mult <- if (ntime%%2==0) {                        ## as per Marple's standard transform.
c(1, rep(2, ntime/2-1), 1, rep(0, ntime/2-1) )  ## even length time series
} else {
c(1, rep(2, (ntime-1)/2 ), rep(0, (ntime-1)/2 ) )  ## odd length time series
}

mult=as.array(replicate(nspace, mult), dim=c(nspace,ntime))

## mvfft computs the fft of each column rather than a 2-d fft.
## but it's annoying/i'm stupid on a single vector TS for some reason, hence the if.
if (nspace==1)
return( fft( t(mult) * fft(input), inverse=TRUE) / ntime )
else
return( t( mvfft( mult * mvfft(t(input)) , inverse=TRUE) / ntime ) )  ## inverse fft is unscaled in R.

}
```
mccreigh/EOF documentation built on May 22, 2019, 12:59 p.m.