SPECTRUM Function
Estimates the power spectrum (power spectral density) of a data sequence.
Usage
result = SPECTRUM(x, length[, overlap])
Input Parameters
x—A one-dimensional array containing the sequence from which the power spectrum will be computed.
length—The length of the subsections of x on which to operate. (Default: N_ELEMENTS(x)/5)
overlap—(optional) The number of elements the successive subsections of x will overlap. (Default: 0.5*length)
Returned Value
result—A one-dimensional array containing the power spectrum of x.
Keywords
Padding—If present and nonzero and greater than length, the subsections of x are padded with the specified number of zeros.
Print_info—If present and nonzero, returns the number and length of subsections of x.
Window_param—Window parameter for a Kaiser window.
Window_type—A scalar value indicating the type of window to use when computing the spectrum of subsections of x. See the Window_type keyword (in the FIRWIN Function on page 64) for the list of values and window types available.
Discussion
SPECTRUM computes one of several different power spectrum estimates P(f) depending on the input parameters. Specific estimates available include the periodogram, the modified periodogram, Bartlett’s method, and Welch’s method. In all cases uniform samples of the power spectrum are returned. If the parameter length = L, the frequency sample values are
fk = k/L , k = 0, 1, ..., M for real data
where M = [(L + 2)/2] for L even and M = [(L + 1)/2] for L odd for real data.
For complex data, M = L.
The periodogram is defined as:
where the frequency variable f is normalized to the Nyquist frequency of 1.0. To obtain uniform samples of the periodogram using SPECTRUM, set length equal to the length of the data x, and set Window_type to rectangular.
The modified periodogram is defined as:
where w(l) is a data window sequence. To obtain uniform samples of the modified periodogram using SPECTRUM, set length equal to the length of the data x and set Window_type to any of the windows discussed in the FIRWIN Function on page 64.
Bartlett’s method breaks the data into non-overlapping data segments represented as:
xi(n) = x(n + iL) n = 0, 1, ... , L – 1 i = 0, 1, ... , I – 1
A periodogram:
is computed for each data segment and averaged to obtain the Bartlett estimate:
To obtain uniform samples of Bartlett’s estimate using SPECTRUM, set length to be less than the data length, set overlap to 0 and set Window_type to rectangular.
The Welch method breaks the data into length L overlapping data segments represented as (n + iL).
xi(n) = x(n + iQ) n = 0, 1, ... , L – 1 i = 0, 1, ... , I – 1
where Q = (L – overlap).
A modified periodogram is computed for each data segment given by:
where:
The Welch power spectrum estimate is the average of the modified periodogram of each data segment, given by:
To obtain uniform samples of the Welch estimate using SPECTRUM, set length less than the data length, set overlap to a nonzero value and set Window_type to any of the windows discussed in the FIRWIN Function on page 64.
note | In estimating the power spectrum it is assumed that the input signal is stationary (i.e., the frequency content does not change with time). If the signal is non-stationary, the functions SPECTROGRAM and WAVELET can often provide better results than SPECTRUM. |
Example
In this example, a multiple bandpass filter is designed and the power spectral density is estimated using SPECTRUM. The results are shown in Figure 4-10: SPECTRUM Function Example, where (a) is a plot of a multiple bandpass filter and (b) is the plot of its power spectral density estimate.
!P.Multi = [0, 1, 2]
f = [0, .18, .2, .22, .38, .4, .42, .58, .6, .62, .78, .8, .82, 1]
ampl = [0, 0, 5, 10, 10, 5, 0, 0, .5, 1, 1, .5, 0, 0]
h = FIRLS(101, f, ampl, /Freqsample)
hf = FREQRESP_Z(h, Outfreq = f)
; Plot a multiple bandpass filter.
PLOT, f, 10.0d*ALOG(ABS(hf)^2), Title = 'Filter Response', $
XTitle = 'Frequency', YTitle = 'Magnitude (dB)'
n = 10000
Length = 1024
Overlap = Length/2
RANDOMOPT, Set = 31
; Generate white noise.
x = RANDOM(n, /Normal)
; Generate colored noise.
y = FILTER(h, x)
sy = SPECTRUM(y, Length, Overlap)
m = FLOAT(N_ELEMENTS(sy))
; Plot power spectral density of the multiple bandpass filter.
PLOT, FINDGEN(m)/m, 10.0d*ALOG(sy), $
Title = 'Power Spectrum Estimate', XTitle = 'Frequency',$
YTitle = 'Magnitude (dB)'
 |
See Also
For Additional Information
Oppenheim and Schafer, 1975.
