MSPEC Multitaper power and cross spectra. _______________________________________________________________________ _______________________________________________________________________ MSPEC implements spectral and cross-spectral analysis using the multi- taper method for real or complex-valued data, and along any dimension. MSPEC is to be run after calling SLEPTAP to compute the multitapers. Confidence intervals can be computed by calling MCONF. _______________________________________________________________________ Real-valued time series [F,S]=MSPEC(X,PSI) returns the power spectrum of the real-valued time series X at positive frequencies using data tapers PSI. The spectrum at negative frequencies, which is not returned, is identical. X may be an array with an arbitrary number of dimensions and with time along it first dimensions. Let M be its number of rows, M=SIZE(X,1). PSI is a then matrix of K data tapers having M rows and K columns. F is an array of frequencies with FlOOR(M/2)+1 rows, while the spectral matrix S has FlOOR(M/2)+1 rather than M rows, and the same size as X along all of its other dimensions. The spectrum can also be computed with time oriented along a different dimension, as described below. By default, MSPEC removes the mean from each time series before computing the spectra. This is suppressed by MSPEC(...,'nodemean'). ______________________________________________________________________ Cross-spectra of real-valued data [F,SXX,SYY,SXY]=MSPEC(X,Y,PSI) computes the cross-spectrum of two real-valued time series or sets of time series. Here SXX and SYY are the one-sided spectra of X and Y, while SXY is their cross spectrum. See TWOSPECPLOT for plotting SXX and SYY simultaneously. ______________________________________________________________________ Rotary spectra of complex-valued data [F,SPP,SNN,SPN]=MSPEC(Z,PSI) where Z is complex-valued computes the so- called "rotary spectra". Here SPP and SNN are the positively-rotating and negatively rotating spectra, and SPN is the rotary cross spectrum. the one-sided spectra of X and Y, while SPN is their cross spectrum. Note that the rotary spectra are defined such that SXX+SYY=SPP+SNN. The rotary spectra SPP and SNN are normalized such that the sum of SPP over all frequencies plus that of SNN approximates the variance of Z. See TWOSPECPLOT for plotting SPP and SNN simultaneously. ______________________________________________________________________ Periodogram MPSEC can be used to form the naive spectral estimator, known as the periodogram. Although this is not generally a good way to estimate the spectrum, it can be useful as a comparision. MSPEC(X,[]) or MSPEC(X,Y,[]) with PSI empty uses the default, or boxcar taper, normalized to unit energy. This returns the periodogram. ______________________________________________________________________ Sample rate [F,S]=MSPEC(DT,...) specifies the sample interval to be used in the calculation of the frequency array F. DT defaults to unity. Spectral values depend linearly upon the sample rate in order that the integral of the spectra over frequency approximate the variance. ______________________________________________________________________ Spectra along arbitary dimension MSPEC(...,DIM) computes the spectrum with time oriented along dimension DIM, with the default behavior corresponding to DIM=1. Let M be the length of the input X, Z, or X and Y along dimension DIM. PSI is again a matrix of K data tapers having M rows and K columns. F will be again an array of frequencies with FlOOR(M/2)+1 rows, while the output spectral matrices will the same size as the input arrays, but with FlOOR(M/2)+1 rather than M elements along dimension DIM. ______________________________________________________________________ Normalizations By default, MSPEC uses *radian* frequency as in cos(f t). Optionally MSPEC(,...,'cyclic') will use *cyclic* frequency, as in cos(2 pi f t). MSPEC is normalized to approximately recover the time series variance. For the MSPEC periodogram, this recovery is exact, although the expressions are complicated somewhat by the use of one-sided spectra. For simplicity, the normalizations will be explained for the case of a single time series with M oriented along rows, that is, with DIM=1. Real-valued data [F,S]=MSPEC(DT,X,[]) where X is a real-valued time series of length M recovers the variance of X, STD(X,1).^2, as follows: 2*(1/2/pi)*(F(2)-F(1))*SUM(S(2:end)) -- M odd 2*(1/2/pi)*(F(2)-F(1))*(SUM(S(2:end-1))+S(end)/2) -- M even where the initial factor of two accounts for the fact that the spectrum at negative frequencies is the same as that at positive frequencies. Note that the zero frequency is omitted in the summation, and for even time series length, the power at the Nyquist S(end) must be divided by two to avoid double-counting by the one-sided spectrum. The "1" in the argument of STD forces STD to use an N rather than N-1 normalization. Complex-valued data [F,SPP,SNN]=MSPEC(DT,Z,[]) where Z is a complex-valued time series of length M recovers the variance of Z, STD(Z,1).^2, as follows: (1/2/pi)*(F(2)-F(1))*(SUM(SPP(2:end))+SUM(SNN(2:end))) -- M odd (1/2/pi)*(F(2)-F(1))*(SUM(SPP(2:end))+SUM(SNN(2:end-1))) -- M even Again the modification for even M prevents the power at the Nyquist from being double-counted. This modification is necessary because the negative rotary spectrum duplicates the Nyquist when M is even. ______________________________________________________________________ Cross-spectra of complex-valued data To compute the cross-spectra of two complex-valued time series or sets of time series Z1 and Z2, run MSPEC repeatedly. [F,SP1P1,SP2P2,SP1P2]=MSPEC(Z1,Z2,PSI); [F,SN1N1,SN2N2,SN1N2]=MSPEC(CONJ(Z1),CONJ(Z2),PSI); The first call returns the spectra and cross-spectra of Z1 and Z2 at positive frequencies, while the second returns their spectra and the *conjugate* of the cross-spectrum at negative frequencies. Finally [F,SP1P1,SN2N2,SP1N2]=MSPEC(Z1,CONJ(Z2),PSI); [F,SN1N1,SP2P2,SN1P2]=MSPEC(CONJ(Z1),Z2,PSI); returns the so-called outer or complementary cross-spectra. ______________________________________________________________________ Adaptive spectra MSPEC(...,LAMBDA,'adaptive'), where LAMBDA contains the eigenvalues of the tapers as computed by SLEPTAP, alternately uses the "adaptive" multitaper method of Thomson (1982). This implementation follows that of Park et al. (1987a), JGR. For cross-spectra or for rotary spectra, the weights appearing in the adaptive spectra are derived for the total spectrum of each signal compoment, that is for SXX+SYY or SPP+SNN as appropriate. Then the separate spectra and co-spectra are computed using identical weights. ______________________________________________________________________ Cell array input / output MSPEC generates cell array output given cell array input. Let's say one has P different time series, X1, X2,..., XP. Put these into a cell array X{1}=X1, X{2}=X2, ..., X{P}=XP, and then use "[psi,lambda]=sleptap(cellength(x))" to make a cell array of tapers. [F,S]=MSPEC(X,PSI) then returns cell arrays F and S corresponding to the Fourier frequencies and spectra of the P arrays. The other argument forms given above also work. In particular, specifiying the sample time through MPSEC(DT,...) works, with DT either a scalar or an array of the same length as the cell array X. The spectra can then be plotted with CELLPLOT(F,S), or TWOSPECPLOT for a pair of output spectra. ______________________________________________________________________ Parallelization MSPEC(..., 'parallel') when the input fields X, X and Y, or Z are cell arrays, parallelizes the spectral estimation by looping over the cells with a PARFOR loop. This requires Matlab's Parallel Computing Toolbox. ______________________________________________________________________ Example The example at the top of this help file shows clockwise (left) and counterclockwise (right) rotary spectra from moored current meter measurements of the ocean currents in the Labrador Sea. The periodogram is in gray, and blue and red are multitaper spectra with P=4 and P=32, respectively. The local Coriolis frequency is marked with a dashed line. Tidal and inertial peaks are apparent. The main point of this figure is to show that increasing P increases the degree of frequency-domain smoothing. ______________________________________________________________________ 'mspec --t' runs some tests. 'mspec --f' generates the above sample figure from Bravo mooring data. See also: SLEPTAP, MCONF, HERMFUN, MSVD, TWOSPECPLOT. Usage [f,s]=mspec(x,psi); [f,s]=mspec(dt,x,psi); [f,s]=mspec(dt,x,psi,dim); [f,spp,snn,spn]=mspec(z,psi); [f,sxx,syy,sxy]=mspec(x,y,psi); [f,sxx,syy,sxy]=mspec(x,y,psi,dim); _________________________________________________________________ This is part of JLAB --- type 'help jlab' for more information (C) 2000--2020 J.M. Lilly --- type 'help jlab_license' for details