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. __________________________________________________________________ See also MATERNCOV, MATERNIMP, MATERNOISE, MATERNFIT, BLURSPEC. '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