xenial (1) spectrum1d.1gmt.gz

Provided by: gmt-common_5.2.1+dfsg-3build1_all bug

NAME

       spectrum1d - Compute auto- [and cross- ] spectra from one [or two] time-series

SYNOPSIS

       spectrum1d  [  table  ]  segment_size]  [  [xycnpago]  ]  [  dt  ] [ [h|m] ] [ [name_stem ] ] [  ] [  ] [
       -b<binary> ] [ -d<nodata> ] [ -f<flags> ] [ -g<gaps> ] [ -h<headers> ] [ -i<flags> ]

       Note: No space is allowed between the option flag and the associated arguments.

DESCRIPTION

       spectrum1d reads X [and Y] values from the first [and second] columns on standard  input  [or  x[y]file].
       These  values  are  treated  as  timeseries X(t) [Y(t)] sampled at equal intervals spaced dt units apart.
       There may be any number of lines of input. spectrum1d will create file[s] containing auto- [and cross-  ]
       spectral  density estimates by Welch's method of ensemble averaging of multiple overlapped windows, using
       standard error estimates from Bendat and Piersol.

       The output files have 3 columns: f or w, p, and e. f or w is  the  frequency  or  wavelength,  p  is  the
       spectral  density  estimate,  and  e  is the one standard deviation error bar size. These files are named
       based on name_stem. If the -C option is used, up to eight files are created; otherwise only one  (xpower)
       is written. The files (which are ASCII unless -bo is set) are as follows:

       name_stem.xpower
              Power spectral density of X(t). Units of X * X * dt.

       name_stem.ypower
              Power spectral density of Y(t). Units of Y * Y * dt.

       name_stem.cpower
              Power spectral density of the coherent output. Units same as ypower.

       name_stem.npower
              Power spectral density of the noise output. Units same as ypower.

       name_stem.gain
              Gain spectrum, or modulus of the transfer function. Units of (Y / X).

       name_stem.phase
              Phase spectrum, or phase of the transfer function. Units are radians.

       name_stem.admit
              Admittance spectrum, or real part of the transfer function. Units of (Y / X).

       name_stem.coh
              (Squared)  coherency  spectrum,  or  linear  correlation  coefficient  as a function of frequency.
              Dimensionless number in [0, 1]. The Signal-to-Noise-Ratio (SNR) is coh / (1 - coh). SNR =  1  when
              coh = 0.5.

       In  addition, a single file with all of the above as individual columns will be written to stdout (unless
       disabled via -T).

REQUIRED ARGUMENTS

       -Ssegment_size]
              segment_size is a radix-2 number of samples  per  window  for  ensemble  averaging.  The  smallest
              frequency  estimated  is  1.0/(segment_size * dt), while the largest is 1.0/(2 * dt). One standard
              error in power spectral density  is  approximately  1.0  /  sqrt(n_data  /  segment_size),  so  if
              segment_size  =  256, you need 25,600 data to get a one standard error bar of 10%.  Cross-spectral
              error bars are larger and more complicated, being a function also of the coherency.

OPTIONAL ARGUMENTS

       table  One or more ASCII (or binary, see -bi) files holding X(t) [Y(t)] samples in the  first  1  [or  2]
              columns. If no files are specified, spectrum1d will read from standard input.

       -C[xycnpago]
              Read the first two columns of input as samples of two time-series, X(t) and Y(t). Consider Y(t) to
              be the output and X(t) the input in a linear system with noise.  Estimate  the  optimum  frequency
              response  function  by  least  squares,  such  that the noise output is minimized and the coherent
              output and the noise output are uncorrelated.  Optionally specify up to 8 letters from the set { x
              y  c  n p a g o } in any order to create only those output files instead of the default [all]. x =
              xpower, y = ypower, c = cpower, n = npower, p = phase, a = admit, g = gain, o = coh.

       -Ddt   dt Set the spacing between samples in the time-series [Default = 1].

       -L     Leave trend  alone.  By  default,  a  linear  trend  will  be  removed  prior  to  the  transform.
              Alternatively, append m to just remove the mean value or h to remove the mid-value.

       -N[name_stem]
              Supply  an  alternate  name  stem  to  be used for output files [Default = "spectrum"].  Giving no
              argument will disable the writing of individual output files.

       -V[level] (more ...)
              Select verbosity level [c].

       -T     Disable the writing of a single composite results file to stdout.

       -W     Write Wavelength rather than frequency in column 1 of the output  file[s]  [Default  =  frequency,
              (cycles / dt)].

       -bi[ncols][t] (more ...)
              Select native binary input. [Default is 2 input columns].

       -bo[ncols][type] (more ...)
              Select native binary output. [Default is 2 output columns].

       -d[i|o]nodata (more ...)
              Replace input columns that equal nodata with NaN and do the reverse on output.

       -f[i|o]colinfo (more ...)
              Specify data types of input and/or output columns.

       -g[a]x|y|d|X|Y|D|[col]z[+|-]gap[u] (more ...)
              Determine data gaps and line breaks.

       -h[i|o][n][+c][+d][+rremark][+rtitle] (more ...)
              Skip or produce header record(s).

       -icols[l][sscale][ooffset][,...] (more ...)
              Select input columns (0 is first column).

       -^ or just -
              Print a short message about the syntax of the command, then exits (NOTE: on Windows use just -).

       -+ or just +
              Print  an  extensive usage (help) message, including the explanation of any module-specific option
              (but not the GMT common options), then exits.

       -? or no arguments
              Print a complete usage (help) message, including the explanation of options, then exits.

       --version
              Print GMT version and exit.

       --show-datadir
              Print full path to GMT share directory and exit.

ASCII FORMAT PRECISION

       The ASCII output formats of numerical data are controlled by parameters in your gmt.conf file.  Longitude
       and  latitude  are formatted according to FORMAT_GEO_OUT, whereas other values are formatted according to
       FORMAT_FLOAT_OUT. Be aware that the format in effect can lead to loss of precision in the  output,  which
       can  lead  to  various  problems downstream. If you find the output is not written with enough precision,
       consider  switching  to  binary  output  (-bo  if  available)  or  specify  more   decimals   using   the
       FORMAT_FLOAT_OUT setting.

EXAMPLES

       Suppose data.g is gravity data in mGal, sampled every 1.5 km. To write its power spectrum, in mGal**2-km,
       to the file data.xpower, use

              gmt spectrum1d data.g -S256 -D1.5 -Ndata

       Suppose in addition to data.g you have data.t, which is topography in meters sampled at the  same  points
       as  data.g. To estimate various features of the transfer function, considering data.t as input and data.g
       as output, use

              paste data.t data.g | gmt spectrum1d -S256 -D1.5 -Ndata -C > results.txt

TUTORIAL

       The output of spectrum1d is in units of power spectral density, and so to get units of  data-squared  you
       must  divide  by delta_t, where delta_t is the sample spacing.  (There may be a factor of 2 pi somewhere,
       also. If you want to be sure of the normalization, you can  determine  a  scale  factor  from  Parseval's
       theorem:  the  sum  of  the squares of your input data should equal the sum of the squares of the outputs
       from spectrum1d, if you are simply trying to get a periodogram. [See below.])

       Suppose we simply take a data set, x(t), and compute the discrete Fourier transform (DFT) of  the  entire
       data set in one go. Call this X(f). Then suppose we form X(f) times the complex conjugate of X(f).

       P_raw(f) = X(f) * X'(f), where the ' indicates complex conjugation.

       P_raw  is called the periodogram. The sum of the samples of the periodogram equals the sum of the samples
       of the squares of x(t), by Parseval's theorem. (If you use a DFT subroutine on a  computer,  usually  the
       sum of P_raw equals the sum of x-squared, times M, where M is the number of samples in x(t).)

       Each  estimate  of  X(f)  is  now formed by a weighted linear combination of all of the x(t) values. (The
       weights are sometimes called  "twiddle  factors"  in  the  DFT  literature.)   So,  no  matter  what  the
       probability  distribution  for  the  x(t)  values  is,  the  probability distribution for the X(f) values
       approaches [complex] Gaussian, by the Central Limit Theorem. This means that the probability distribution
       for  P_raw(f)  approaches  chi-squared  with  two  degrees  of  freedom.  That  reduces to an exponential
       distribution, and the variance of the estimate of P_raw is proportional to the square of the  mean,  that
       is, the expected value of P_raw.

       In  practice  if we form P_raw, the estimates are hopelessly noisy. Thus P_raw is not useful, and we need
       to do some kind of smoothing or averaging to get a useful estimate, P_useful(f).

       There are several different ways to do this in the literature. One is to form P_raw and then  smooth  it.
       Another  is  to  form the auto-covariance function of x(t), smooth, taper and shape it, and then take the
       Fourier transform of the smoothed, tapered and shaped auto-covariance.  Another is to form  a  parametric
       model  for  the  auto-correlation  structure  in x(t), then compute the spectrum of that model. This last
       approach is what is done in what is called the "maximum entropy" or "Berg" or "Box-Jenkins" or "ARMA"  or
       "ARIMA" methods.

       Welch's  method  is  a  tried-and-true  method.  In his method, you choose a segment length, -SN, so that
       estimates will be made from segments of length N. The frequency samples (in cycles per delta_t  unit)  of
       your  P_useful will then be at k /(N * delta_t), where k is an integer, and you will get N samples (since
       the spectrum is an even function of f, only N/2 of them are really useful). If the length of your  entire
       data set, x(t), is M samples long, then the variance in your P_useful will decrease in proportion to N/M.
       Thus you need to choose N << M to get very low  noise  and  high  confidence  in  P_useful.  There  is  a
       trade-off here; see below.

       There  is  an  additional reduction in variance in that Welch's method uses a Von Hann spectral window on
       each sample of length N. This reduces side lobe leakage and has the effect of smoothing the  (N  segment)
       periodogram  as  if  the X(f) had been convolved with [1/4, 1/2, 1/4] prior to forming P_useful. But this
       slightly widens the spectral bandwidth of each estimate, because the estimate at frequency  sample  k  is
       now  a  little correlated with the estimate at frequency sample k+1. (Of course this would also happen if
       you simply formed P_raw and then smoothed it.)

       Finally, Welch's method also uses overlapped processing. Since the Von Hann window is large in the middle
       and  tapers  to near zero at the ends, only the middle of the segment of length N contributes much to its
       estimate. Therefore in taking the next segment of data, we move ahead  in  the  x(t)  sequence  only  N/2
       points.  In this way, the next segment gets large weight where the segments on either side of it will get
       little weight, and vice versa. This doubles the smoothing effect and ensures that  (if  N  <<  M)  nearly
       every point in x(t) contributes with nearly equal weight in the final answer.

       Welch's  method  of  spectral estimation has been widely used and widely studied. It is very reliable and
       its statistical properties are well understood. It is highly recommended in  such  textbooks  as  "Random
       Data: Analysis and Measurement Procedures" by Bendat and Piersol.

       In  all  problems of estimating parameters from data, there is a classic trade-off between resolution and
       variance. If you want to try to squeeze more resolution out of your data set, then you have to be willing
       to  accept more noise in the estimates. The same trade-off is evident here in Welch's method. If you want
       to have very low noise in the spectral estimates, then you have to choose N << M, and this means that you
       get  only N samples of the spectrum, and the longest period that you can resolve is only N * delta_t.  So
       you see that reducing the noise lowers the number of spectral samples  and  lowers  the  longest  period.
       Conversely,  if  you  choose  N  approaching  M,  then  you  approach  the  periodogram with its very bad
       statistical properties, but you get lots of samples and a large fundamental period.

       The other spectral estimation methods also can do a good job. Welch's method was selected because the way
       it  works,  how  one  can  code  it,  and its effects on statistical distributions, resolution, side-lobe
       leakage, bias, variance, etc. are all easily understood. Some of the other methods (e.g. Maximum Entropy)
       tend to hide where some of these trade-offs are happening inside a "black box".

SEE ALSO

       gmt, grdfft

REFERENCES

       Bendat, J. S., and A. G. Piersol, 1986, Random Data, 2nd revised ed., John Wiley & Sons.

       Welch, P. D., 1967, The use of Fast Fourier Transform for the estimation of power spectra: a method based
       on time averaging over short, modified periodograms, IEEE Transactions on Audio and Electroacoustics, Vol
       AU-15, No 2.

       2015, P. Wessel, W. H. F. Smith, R. Scharroo, J. Luis, and F. Wobbe