# MATERNSPEC is the jMatern module of jLab.

``` MATERNSPEC  Fourier spectrum of the Matern random process and variations.

[F,S]=MATERNSPEC(DT,N,SIGMA,ALPHA,LAMBDA) returns the spectrum S of a
length N complex-valued Matern random process having variance SIGMA^2,
slope parameter ALPHA, and damping parameter LAMBDA.

DT is the sample interval.  Note that LAMBDA is understood to have the
same units as the inverse sample interval 1/DT.

F is an array of one-sided (positive) Fourier frequencies for a time
series of length N, F=FOURIER(N), where F is a *radian* frequency.

The lengths of the output variables F and S are N/2+1 for even N, and
(N+1)/2 for odd N.

S is the postive or negative rotary spectrum given by

S(F) = SIGMA^2 / (F^2+LAMBDA^2)^ALPHA * LAMBDA^(2*ALPHA-1)/C

where C is a normalizing constant dependent upon ALPHA.  Note that the
positive and negative spectra are identical for the Matern process.

For LAMBDA=0, the Matern spectrum reduces to the spectrum of fractional
Brownian motion.

For details on the Matern process and its spectrum, see:

Lilly, Sykulski, Early, and Olhede, (2017).  Fractional Brownian
motion, the Matern process, and stochastic modeling of turbulent
dispersion.  Nonlinear Processes in Geophysics, 24: 481--514.
__________________________________________________________________

Matrix and cell array output

[F,S]=MATERNSPEC(DT,N,SIGMA,ALPHA,LAMBDA) where N is a scalar while the
other input arguments are all either scalars or arrays of the same
length M, gives an output spectra S with LENGTH(F) rows and M columns.

[F,S]=MATERNSPEC(DT,N,SIGMA,ALPHA,LAMBDA) where N is an array of M
different lengths, returns F and S that are length M cell arrays.  Then
SIGMA, ALPHA, and LAMBDA may all either be scalars or length M arrays.

This latter format is convenient for generating sets of spectra that
do not all have the same size.

When N is an array, MATERNSPEC(...,'parallel') parallelizes the
computation of the various spectra using a PARFOR loop.  This option
requires that Matlab's Parallel Computing Toolbox be installed.

The matrix and cell array formats also work for the variations of the
Matern process described below.
__________________________________________________________________

Oscillatory Matern

[F,SPP,SNN]=MATERNSPEC(DT,N,SIGMA,ALPHA,LAMBDA,NU) with six input
arguments modifies the spectrum to have a rotation frequency NU.

This is accomplished by shifting the spectrum to be centered on F=NU
rather than F=0.  SPP and SNN are now the postive rotary and negative
rotary spectra, with the spectrum for positive frequencies +F returned
in SPP, and for negative frequencies -F in SNN.

With ALPHA=1, the oscillatory Matern becomes the complex Ornstein-
Uhlenbeck process.

Note that NU has units of radians per sample interval DT.

The oscillatory Matern is described in Lilly et al. (2017).
__________________________________________________________________

Experimental extensions

The remaining features are experimental extensions to the Matern
process.  They are not yet documented in a publication, and should
be considered as 'beta features' that are to be used with caution.
__________________________________________________________________

Extended Matern

[F,S]=MATERNSPEC(DT,N,SIGMA,ALPHA,LAMBDA,0,MU) with seven arguments
returns the spectrum of the four-parameter "extended" Matern process:

S(F) = SIGMA^2 * BESSELK(ALPHA,MU*SQRT(F^2+LAMBDA^2))
/ (MU*SQRT(F^2+LAMBDA^2))^ALPHA * C

where C is a normalizing constant dependent upon ALPHA, LAMBDA, and MU.
The additional parameter, MU, has units of time.  Here ALPHA can
take on any real value, unlike for the standard Matern case.

[F,SPP,SNN]=MATERNSPEC(DT,N,SIGMA,ALPHA,LAMBDA,NU,MU) shifts the
extended Matern spectrum to be centered at F=NU rather than F=0.

As MU becomes large with ALPHA>1/2, this becomes the standard Matern
spectrum.
__________________________________________________________________

Damped exponential

[F,S]=MATERNSPEC(DT,N,SIGMA,-1/2,LAMBDA,0,MU) with ALPHA set to -1/2
returns the spectrum of the damped exponential process, having the form

S(F) = SIGMA^2 * EXP(-MU * SQRT(F^2+LAMBDA^2)) * C

where C is again a normalizing constant, dependent upon LAMBDA and MU.

This is a special case of the extended Matern process.  As for that
process, setting NU to a nonzero value results in a shifted spectrum.
__________________________________________________________________

Composite Matern

[F,SPP,SNN]=MATERNSPEC(DT,N,SIGMA,ALPHA,LAMBDA,NU,MU,'composite')
implements the "composite" Matern spectrum having the form

SPP(F) = B * SIGMA^2 / (F^2 + MU^2)^ALPHA / [(F-NU)^2 + LAMBDA^2]
SNN(F) = B * SIGMA^2 / (F^2 + MU^2)^ALPHA / [(F+NU)^2 + LAMBDA^2]

where B is a normalizing constant discussed shortly.  This consists of
a Matern spectrum times an oscillatory Matern spectrum having ALPHA=1.

The quantity in square brackets is recognized as the transfer function
for a damped simple harmonic oscillator.  In oceanographic terms, this
composite model gives the spectrum of a damped slab model of the
surface mixed layer forced by winds having a Matern spectrum.

The interpretation of the variance SIGMA is different from the other
cases in MATERNSPEC, because an analytic form of the total variance
does not exist. Instead SIGMA^2 is an approximation to the variance
associated with the oscillatory peak at F=NU.

The additional parameter here, MU, has units of *frequency* and is the
damping parameter associated with the background process, which in this
case reprents the structure of the wind spectrum.

Here B = 2 * LAMBDA * (NU^2 + MU^2) is a normalizing constant that lets
SIGMA^2 be interpreted as an approximation to the inertial variance.
__________________________________________________________________

'maternspec --f' generates some sample figures.

Tests for MATERNSPEC can be found in MATERNCOV.

Usage:  [f,s]=maternspec(dt,N,sigma,alpha,lambda);
[f,spp,snn]=maternspec(dt,N,sigma,alpha,lambda);
[f,spp,snn]=maternspec(dt,N,sigma,alpha,lambda,nu);
[f,spp,snn]=maternspec(dt,N,sigma,alpha,lambda,nu,mu);
__________________________________________________________________
This is part of JLAB --- type 'help jlab' for more information
(C) 2013--2017 J.M. Lilly --- type 'help jlab_license' for details
```