Estimating The Fourier Spectrum Of A Random Process

Subchapters

Estimating The Fourier Spectrum Of A Random Process

Syntax	`fspec(`x`,` dt`,` taper`)`
See Also	fft , dpss , arcov , filter

Description
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

      i

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 s) ds |

      T  |  /                             |

            -T/2

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

     ----------------------------------------------------------

                             N

                       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 2 < k < N/2, the k-th element of the expression above is equal to the N-k+2-th element.

Example
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






clear

#

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

gyaxis("log")

gtitle("Spectrum of White Noise")

gplot(f, spec, "plus")