Program src/ts/wf/foutra
Short description
Applies a spectral analysis of time series by computing the Fourier transform and providing either one of the following output quantities:
 Power spectral density
 Average signal power in a finite relative bandwidth
 Fourier amplitude spectrum
 Amplitudes of harmonic signals
Location
Repository location: src/ts/wf
Contents
Program description
Purpose
foutra applies spectral analysis to time series data. Its primary purpose is the computation of power spectral density (power). Additionally it offers two other modes of scaling. Values of the Fourier amplitude spectrum can be written to file (amplitude) as well as values of a spectral representation, where peak values in the spectrum represent amplitudes of harmonic signals present in the time series (harmonic). If none of these modes is selected, a tapered version of the input time series is written to the output file.
Scaling of Fourier transforms
The Fourier transform presented in libfourier is scaled appropriately to represent coefficients of the Fourier integral transform (see libfourier documentation).
Scaling factor for Hanning taper when applied in computation of powerspectral density
Tapering stationary noise with a Hanning taper reduces total signal energy
an thus the signal power obtained through an FFT. In his tutorial Agnew
recommends to scale each sample s_k
of the series by w_k/W
,
where w_k
is the k
th sample of the taper and
W^2 = \frac{1}{N} \sum\limits_0^{N1} w_k^2
is a measure for the loss in total signal power due to application of the taper. This applies only for staionary noise.
Walter Zürn implements a correction for loss of signal power in his Matlab code as follows.
If data
(not used here)
is the time series vector and y
is the Hanning taper of appropriate
length, the program computes
w = sqrt(sum(y.*y)/ndata);
and
y=y/w;
where ndata
is the length of data
and y
.
Factor
w_k = \sin^2(k \pi/(N1)).
of a Hanning taper will be applied to sample k
of a time series of N
samples in total.
Thus
W^2 = \frac{1}{N} \sum\limits_0^{N1} \sin^4(k\pi/(N1))
gives the fraction of signal energy which remains after tapering.
From Gradshteyn and Ryzhik (eq. 1.321 3)
\sin^4(k\pi/(N1))
= \frac{1}{8} \left( \cos(4 k\pi/(N1))
 4 \cos(2 k\pi/(N1))
+ 3 \right).
Within the sum over k=0\ldots N1
the contribution of both costerms will
vanish, since both are averaged over one and two periods, respectively. Thus
W^2 = \frac{1}{N} \sum\limits_0^{N1} \frac{3}{8} = \frac{3}{8}.
Because foutra is not scaling the taper but scaling the power spectrum, it must apply the factor 8/3 to the result of power spectrum computation in order to compensate the loss in signal energy.
This factor 8/3=2.66667 was tested against the value for W^2
, when
explicitely derived from a Hanning taper time series by the above formula.
Makefile.foutra
in the collection of test
cases provides a section with test code
for foutra.
This offers instantaneous testing of PSD scaling in foutra.
Scaling for harmonic signals
For harmonic signals the FFT normalized to the duration of the available time window has an amplitude of the value of the maximum of the Fourier transform of the window function times half the amplitude of the time domain signal. To obtain a spectral representation with peaks of the amplitude of the time domain signal, the FFT must be scaled accordingly.
The Fourier transform
\tilde{f}(\omega)=\int\limits_{\infty}^{+\infty}
f(t) \,e^{i\omega t} \textrm{d}{t}
corresponds to the time domain signal
f(t)=\int\limits_{\infty}^{+\infty}
\tilde{f}(\omega) \,e^{i\omega t} \frac{\textrm{d}{\omega}}{2\pi}.
Application of a time domain window function w(t)
to the function
f(t)
results in the tapered function
g(t)=f(t) w(t)
and its Fourier transform
\tilde{g}(\omega)=\int\limits_{\infty}^{+\infty}
\tilde{f}(\omega')\,\tilde{w}(\omega\omega')\,\textrm{d}\omega.
Application to harmonic signals
Let
f(t)=A\,\cos(\omega_0 t+\phi)
be the harmonic signal under investigation. Then
f(t)=A\left\{\cos(\omega_0 t)\cos(\phi)
\sin(\omega_0 t)\sin(\phi)\right\}
and
\tilde{f}(\omega) =
\frac{A}{2}\left\{
\cos(\phi)
\left[
\delta(\omega\omega_0)
+
\delta(\omega+\omega_0)
\right]
+i\sin(\phi)
\left[
\delta(\omega\omega_0)

\delta(\omega+\omega_0)
\right]
\right\},
where \delta(\omega)
is Dirac's delta function with
\delta(\omega)=\left\{
\begin{array}{ll}
\infty & \textrm{if $\omega=0$ and}\\
0 & \textrm{otherwise}
\end{array}
\right.
and
\int\limits_{\infty}^{+\infty}
\delta(\omega)\,\textrm{d}\omega=1
such that
\tilde{f}(\omega)=\int\limits_{\infty}^{+\infty}
\tilde{f}(\omega')\,\delta(\omega\omega')\,\textrm{d}\omega.
This way
\tilde{g}(\omega)=
\frac{A}{2}\left\{
e^{i\phi}\tilde{w}(\omega\omega_0)
+
e^{i\phi}\tilde{w}(\omega+\omega_0)
\right\}
is the Fourier transform of the tapered function.
Ignoring interference with sidelobes from the negative frequency
\omega_0
and sidelobes of potential other harmonics at nearby
frequencies, we can approximate
\tilde{g}(\omega_0)\approx\frac{A}{2}e^{i\phi}\tilde{w}(0)
and
\left\tilde{g}(\omega_0)\right\approx
\frac{A}{2}\left\tilde{w}(0)\right.
Hence
A=\frac{2\,\left\tilde{f}(\omega_0)\right}{\left\tilde{w}(0)\right}
is the amplitude of the corresponding time domain harmonic signal with
frequency \omega_0
.
The maximum w_{\textrm{max}}=\tilde{w}(0)
of the transform of the Boxcar taper is T
, where T
is the
duration of the window. Similarly the maximum of the Fourier transform of the
Hanning taper is 0.5*T
, where the duration of the window is T
.
Boxcar taper function
The boxcar taper is defined as
w(t) =
\left\{
\begin{array}{ll}
1 & \textrm{if $t\leq T/2$ and}\\
0 & \textrm{otherwise}
\end{array}
\right.
with
\tilde{w}(\omega)
= T \frac{\sin(\omega T/2)}{\omega T/2}
and
w_{\textrm{max}}=\tilde{w}(0)=T.
Hanning taper function
The Hanning taper is defined as
w(t) =
\left\{
\begin{array}{ll}
\cos^2(\pi t / T))
= \frac{1}{2}\left[
1 + \cos(2\pi t / T)
\right]
& \textrm{if $t\leq T/2$ and}\\
0 & \textrm{otherwise}
\end{array}
\right.
with
\tilde{w}(\omega)
= T/2 \frac{\sin(\omega T/2)}{\omega T/2}
+ T/4 \frac{\sin(\omega T/2+\pi)}{\omega T/2+\pi}
+ T/4 \frac{\sin(\omega T/2\pi)}{\omega T/2\pi}
(see Blackman, R.B. and Tukey, J.W. 1958. The measurement of power spectra. Dover Publications, Inc., New York. Section B.5) and
w_{\textrm{max}}=\tilde{w}(0)=\frac{T}{2}.
Usage
FOUTRA V1.10 Fourier transforms
usage: foutra [v] [o] [type type] [Type type] [D] [ASCII[=base]]
[amplitude] [power] [logascii[=n]]
[boxcar] [avg[=n]] [rbw[=n]] [avgascii]
[demean[=n]] [dtrend[=n]] [scalerbw[=n]]
[divisor[=n]] [rms] [harmonic] [pad n]
[derivative n]
[nsegments n]
outfile infile [t:T] [infile [t:T] ...]
or: foutra helph
or: foutra xhelp
outfile output filename
infile input filename
t:T select traces T, where T may be any range
specification like '34' or '5,6,712,20'
help prints this help text
xhelp print information concerning supported data formats
v be verbose
D debug mode
o overwrite output
type type select input file type
Type type select output file type
amplitude compute Fourier amplitude spectrum
power compute power spectral density (or power spectrum)
harmonic scale output appropriate fro harmonic signals
useful for twotonetests of linearity (see below)
rms report rms values of input data
boxcar apply boxcar taper (i.e. no taper; default is Hanning)
demean[=n] remove average from input time series
(determined from n samples)
detrend[=n] remove trend from input time series
(determined from n samples)
divisor[=n] FFT becomes very inefficient if the factorization
of the number of samples includes large prime numbers.
This option removes the least number of samples to
the total number of samples a multiple of "n"
pad n pad time series with zeroes; n gives the integer factor
for the number of samples; the raw amplitude spectrum
has to be understood as the spectrum of the whole
series including the padded zeroes; PSD and harmonic
signals are scaled to represent the taper time window
only, such that padding is a means of smoothing only
nsegments n subdivide input time series in "n" overlapping
sequences (overlap is 50%); spectral analysis is
done for each sequence individually; the result
is the mean of all sequences; this applies to PSD
and harmonic signal analysis only; it is particularly
useful for twotonetest where spectral smoothing
of background noise is anticpated, while maintaining
the full resolution for harmonic peaks
avg[=n] smooth power spectrum by averaging over n samples
values still specify power spectral density
rbw[=n] smooth power spectrum by averaging over n decades
values still specify power spectral density
derivative[=n] effectively take nth derivative of time series
computation is done in the Fourier domain and has no effect on
time series
scalerbw[=n] scale to mean value in n decades
output value then are not power spectral density but average
signal power in the specified bandwidth
ASCII[=base] write result to twocolumn ASCII files with basename 'base'
logascii[=n] write ASCII data on logarithmic frequency axis with
one value per 'n' decades
avgascii only average values for output to ASCII file
this option speeds up calculation together with
scalerbw which increases computation time
with the square of frequency
foutra applies spectral analysis to time series data. Its primary purpose is
the computation of power spectral density (power). Additionally it offers two
other modes of scaling. Values of the Fourier amplitude spectrum can be
written to file (amplitude) as well as values of a spectral representation,
where peak values in the spectrum represent amplitudes of harmonic signals
present in the time series (harmonic). If non of these modes is selected, a
tapered version of the input time series is written to the output file.
Output is written to a file in any of the time series data formats (Type)
supported for output. Samples are then values along increasing frequency,
sampling interval is given in Hertz. The first sample in the file is at 0 Hz.
The last is the sample at Nyquist frequency.
As an option (ASCII) output to plain ASCII is supported, where the first
column in each file specifies frequency in Hertz and the second column
specifies the corresponding spectral values. This format supports output with
logarithmic sampling along the frequency scale (logascii).
Power spectral density

Option: power
If input units are K, then the output units of power spectra will
be K*K/Hz. The units for amplitude spectra then are K/Hz. If scaling
to the mean in a relative bandwidth is used (only applies for power
spectra; switch scalerbw) the output units are K*K.
The integral over the power spectral density calculated by foutra
over the total bandwidth (over all frequencies from 0 Hz to Nyquist
frequency) provides the total power of the signal (i.e. the
variance of the input signal, i.e. the square of the rms value).
This is called the onesided power spectral density.
Option scalerbw selects the computation of average signal power in a finite
relative bandwidth (rather than power spectral density).
Amplitudes of harmonic signals

Option: harmonic
The Fourier transformation does not exist for harmonic signals. The option
"harmonic" supports the analysis of a time limited portion of a harmonic
signal. If this option is set, the output is scaled such that the peak values
in the spectrum equal the amplitude of the corresponding harmonic time domain
signal in units of K.
Details of the scaling used for this option are given in doxygen formatted
comments at the end of the source code file.
The input signal can be extended by padding with zeroes (pad). This mainly is
used to obtain a smoother representation of the corresponding spectral display
where the amplitude of narrow peaks might not fall on a spectral node when
analysing harmonic signals. The spectral representation consequently is scaled
to the length of the applied taper function (and not to the full length of the
padded time series) for harmonic signal analysis and power spectral densities,
since the signal under investigation is understood as being infinite in time.
For amplitude spectra of transient signals (the default) no scaling to a time
window takes place, since padding with zeroes is understood as not altering
the signal.
data formats supported for reading:

libdatrwxx (version 20180406)
data formats supported by ianystream:
sff: Stuttgart File Format
hpmo: HPMO data format defined by W. Grossmann (BFO)
pdas: PDAS100 (i.e. DaDisp)
mseed: MiniSEED (SeisComP, EDL, etc.)
bonjer: K2 ASCII data format (defined by K. Bonjer?)
sac: SAC binary format
gse: raw GSE format
tsoft: TSOFT format
tfascii: ASCII format of T. Forbrigers any2ascii
su: SeismicUn*x format
seife: seife format (E. Wielandt)
thiesdl1: Thies DL1 pluviometer data at BFO
ascii: simple single column ASCII data
bin: binary data
data formats supported for writing:

libdatrwxx (version 20180406)
data formats supported by oanystream:
sff: Stuttgart File Format
gse: raw GSE format
su: SeismicUn*x format
seife: seife format (E. Wielandt)
ascii: simple single column ASCII data
bin: binary data