NAME¶
spectrum1d - compute auto- [and cross- ] spectra from one [or two] timeseries.
SYNOPSIS¶
spectrum1d [
x[y]file ]
-Ssegment_size] [
-C[
xycnpago] ] [
-Ddt ] [
-Nname_stem ] [
-V ] [
-W ] [
-bi[
s][
n] ] [
-bo[
s][
n] ]
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.
REQUIRED ARGUMENTS¶
- x[y]file
- ASCII (or binary, see -bi) file holding X(t) [Y(t)]
samples in the first 1 [or 2] columns. If no file is specified,
spectrum1d will read from standard input.
- -S
- 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.
OPTIONS¶
- -C
- Read the first two columns of input as samples of two
timeseries, 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 f requency response function by least squares,
such that the noise output is minimized and the coherent outpu t 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.
- -D
- dt Set the spacing between samples in the timeseries
[Default = 1].
- -N
- name_stem Supply the name stem to be used for output
files [Default = "spectrum"].
- -V
- Selects verbose mode, which will send progress reports to
stderr [Default runs "silently"].
- -W
- Write Wavelength rather than frequency in column 1 of the
output file[s] [Default = frequency, (cycles / dt)].
- -bi
- Selects binary input. Append s for single precision
[Default is double]. Append n for the number of columns in the
binary file(s). [Default is 2 input columns].
- -bo
- Selects binary output. Append s for single precision
[Default is double].
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, try
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, try
paste data.t data.g | spectrum1d
-S256
-D1.5
-Ndata
-C
SEE ALSO¶
gmt(1gmt),
grdfft(1gmt)
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.