Contents Previous Next Subchapters

Estimating The Fourier Spectrum Of A Random Process
Syntax fspec(xdttaper)
See Also fft , dpss , arcov , filter

Returns an estimate for the Fourier spectrum of a random process using one or more tapers and realizations. Each column of the integer, real, double-precision, or complex column matrix x contains a realization of a random process X. The row dimension of x must be even (if its only prime factors are 2, 3, 5, and 7, an FFT is used and the calculation is quicker). The integer, real, or double-precision scalar dt specifies the time between the samples. The i-th row of x corresponds to samples at time
     t  =  (i - N / 2 - 1) dt
Each column of the integer, real, double-precision, or complex matrix taper specifies a tapering window where taper and x have the same number of rows.

If N is the row dimension of x, the return value of fspec is a double-precision column vector of length N. The k-th element of the return value is the spectral density estimate at frequency
           k - N / 2 - 1
     f  =  -------------
      k         N dt

The Fourier spectrum of X, at frequency f, is defined as the limit as T goes to infinity of the expected value of
             T/2                           2
      1  |  /                  __         |
     --- |  | X(s) exp(-2 pi \/-1 f sds |
      T  |  /                             |

The spectral estimate corresponding to the j-th column of x and the m-th column of taper is
              N                                               2
      1   | -----                            __            |
     ---- | >      taper    x    exp(-2 pi \/-1  f  t ) dt |
     N dt | -----       i,m  i,j                  k  i     |
            i = 1
                       1   -----       2
                      ---  >      taper
                       N   -----       i,m
                           i = 1
The returned value is a double-precision column vector of length N in which the k-th element is the average over all windows and realizations of the estimate above. If x is real or double-precision, and < k < N/2, the k-th element of the expression above is equal to the N-k+2-th element.

The spectral estimate corresponding to a rectangular taper is sometimes referred to as the periodogram. The following program computes the spectral estimate for 256 points of white noise using a rectangular taper. The noise is normally distributed with a mean of zero and variance of one, on a uniform time grid between time 0 and 256. It then plots the spectral estimate

N     = 256
x     = snormal(N, 1)
dt    = 1.
taper = fill(1, N, 1) 
spec  = fspec(x, dt, taper)
df    = 1. / (N * dt)
f     = (seq(N) - N / 2 - 1) * df
gtitle("Spectrum of White Noise")
gplot(f, spec, "plus")