Power Spectrum: 6.1 Outline

Download as pdf or txt
Download as pdf or txt
You are on page 1of 7

Chapter 6

Power Spectrum
The power spectrumanswers the question Howmuch of the signal is at a frequency
?. We have seen that periodic signals give peaks at a fundamental and its
harmonics; quasiperiodic signals give peaks at linear combinations of two or more
irrationally related frequencies (often giving the appearance of a main sequence and
sidebands); and chaotic dynamics give broad band components to the spectrum.
Indeed this later may be used as a criterion for identifying the dynamics as chaotic.
Examples are shown in demonstration 1. These are all statements about the ideal
power spectrum, if innitely long sequences of continuous data are available to
process. In practice there are always limitations of restricted data length and
sampling frequency, and it is important to investigate how these limitations affect
the appearance of the power spectrum.
6.1 Outline
If we did not have to worry about limitations in the datai.e. we have a continuous
time series y(t ) innite in lengththe power spectrumof the signal would be given
simply by the Fourier transform: P
ideal
() | y()|
2
with
y() =
_

y(t )e
it
dt. (6.1)
On the other hand we usually have a signal y(t ) measured over a nite interval
0 t T and with a nite sample rate so that we have N values of y is measured
at intervals t = int eger (so that T = N). Then to estimate the power
1
CHAPTER 6. POWER SPECTRUM 2
spectrum of the signal we calculate the Fourier series
y
k
=
N1

j=0
y
j
exp
_
2ijk
N
_
=
N1

0
y(t
j
) exp
_
i
k
t
j
_
(6.2)
where in the latter expression the discrete frequencies and times
k
= 2k/T
and t
j
= j are introduced. (For a discrete time system of course the dynamics
is given in terms of the index j.) For concreteness we take N to be even in the
following. The power spectrum estimate is then (see ref. [1] for details)
P()
_

_
N
2
| y
0
|
2
for = 0
N
2
_
| y
k
|
2
+| y
Nk
|
2
_
for = 2k/T, k = 1, 2, . . . (
N
2
1)
N
2

y
N/2

2
for = N/T = /
(6.3)
where we have used | y
k
| = | y
k
| for a real signal and y
k
= y
Nk
from (6.2). We
will use P only for positive frequenciesthe reason for the combination appearing
for frequencies away from the two end points. The normalization is chosen so that
the sumof P() over the
N
2
+1 frequencies
k
is the mean square amplitude of the
signal y
j
. Other techniques for estimating the power spectrum also exist [1]. The
question immediately arises of howthe estimate P() relates to the ideal spectrum
P
ideal
().
To understand the full properties of this estimate of the spectrum it is useful to
identify (6.2) as the Fourier transform of the modied time series y(t ) dened by
a rather complicated expression that separates out the different imperfections in
the measurement:
y(t ) = [(y(t ) H(t, T )) S(t, T )] S(t, ) . (6.4)
Here denotes a convolution, H(t, ) is the top hat function, with value 1 for
0 < t < and zero elsewhere, and S(t, ) is the periodic spike function zero
everywhere except at t = n with n any integer. Note that the convolution with
S(t, T ) in (6.4) simply repeats the function y(t )H(t, T ) periodically with period
T which is what the Fourier series construction (6.2) effectively doesand the
multiplication by S(t, ) is the sampling at the discrete intervals. The ideal
power spectrum would be given simply by the Fourier transform y() of y(t ).
Now we need three facts (see section 6.2 for further details):
CHAPTER 6. POWER SPECTRUM 3
(i) The Fourier transformof a product is the convolution of the Fourier transforms:
if Y(t ) = y
1
(t ) y
2
(t ) then

Y() = y
1
() y
2
() .
(ii) The Fourier transformof a convolution is the product of the Fourier transforms:
if Y(t ) = y
1
(t ) y
2
(t ) then

Y() = y
1
() y
2
() .
(iii) The Fourier transform of a periodic set of spikes with separation is also a
periodic set of spikes and with separation in frequency
2

; i.e. the Fourier


transform of S(t, ) with respect to t is just proportional to S(,
2

). This
can be seen by noting the analogy to the diffraction pattern of a diffraction
grating. Alternatively note that since S(t, ) is periodic in t with period ,
its Fourier transform is nonzero only at the frequency
2

and its harmonics,


and also since S(t, ) is very non-sinusoidal (each peak has amplitude at all
frequencies) we would expect a large amplitude in each of the harmonics.
So now we have
y()
__
y()

H(,
2
T
)
_
S(,
2
T
)
_
S(,
2

) (6.5)
where


H(,
2
T
)

is just the familiar sinch function


_
sin
1
2
T
_
/
1
2
, which is
maximum at = 0, and has a width
1
T
, but with oscillating tails falling
off only as
1
. Figure 6.1 shows the effect of

H(,
2
T
) on a pure frequency
s
that is not an exact multiple of
2
T
, so does not t with the measurement period
T . The line shows the functional form of the convolution in the [] in (6.5) for
y() = (
s
) i.e.


H(
s
)

, and the points show the values at the discrete


frequencies n
2
T
. Note that the oscillations do not show up because the same
scale
2
T
determines the oscillations and the discrete frequencies. (The values used
in this plot are
s
= .126, T = 8 so that
2
T
=

8
, and = 0.25 so that the
Nyquist frequency / is 4.)
So we are told in (6.5) to take the ideal Fourier transform y(), convolve
(broaden) it with the resolution function

H, sample the result at = int eger
2
T
and then superimpose the result repeatedly displaced by = int eger
2

. The
last step means that we can restrict the range of to the Nyquist range

<

since other ranges are just duplications of this range, and that amplitudes
CHAPTER 6. POWER SPECTRUM 4
Figure 6.1: Power spectrumof exp (i
s
t ) withfrequency
s
= 0.126 constructed
from a time series of length T = 8 with steps = 0.25. (A complex signal is
used for clarity. A real signal cos
s
t would be the superposition of this curve with
the same curve centered around
s
.)
in the ideal y() outside this range will be folded back or aliased into the range
(e.g. an amplitude at =
3
2
will appear at =

2
due to the copy shifted to
=
2

).
To get a better estimate of the spectrumit is customary to replace the discontin-
uous function H by a window function W, such as a tent, parabola or other shape
(i.e. multiply the data by such a function). Since the high frequency tails of the
Fourier transform are determined by the order of the derivative in which a discon-
tinuity rst appearsthe power p of the
p
tail is the order of this derivative
plus 1this pushes the discontinuities to higher order derivatives and so makes
the Fourier transform

W fall off more quickly with . Windowing is investigated
in demonstration 2.
It is also important to make sure that the Nyquist frequency / is large
enough (i.e. the sampling rate high enough), so that there is not signicant power
beyond this frequency, to minimize the effects of the aliasing. A nal point is
CHAPTER 6. POWER SPECTRUM 5
that for nonperiodic signals the estimate y
k
is a very noisy estimate of the power
spectrum i.e. different choices of which time interval T to measure will lead to a
power spectrum rather different in the details. (You can see this from the gure:
the largest value in the apparent power spectrum depends how close
s
is to some
integer multiple of 2/T ). In fact the variance of y
k
is equal to the mean! To reduce
the scatter take a number of power spectra (from different sets of measurements
over intervals of length T ) and average the power spectra.
The dependence of the appearance of the power spectrumon various parameters
is illustrated in demonstration 3.
6.2 Details
This section is for reference only, and can be skipped on a rst reading.
The convolution is dened as
y
1
(t ) y
2
(t ) =
_

y
1
(t

)y
2
(t t

) dt

, (6.6)
and we dene the Fourier transform as
y() =
_

y(t )e
it
dt
y(t ) =
1
2
_

y()e
it
d
. (6.7)
The periodic spike function S(t, ) is dened formally as
S(t, ) =

n=
(t n)
with the Dirac Delta function. The Fourier transform of S(t, ) with respect to t
is just
2

S(,
2

) since we have

S() = lim
M
_

n=M
(t n) e
it
dt = lim
M
M

n=M
e
in
. (6.8)
giving

S() = lim
M
sin(M +
1
2
)
sin
1
2

(6.9)
CHAPTER 6. POWER SPECTRUM 6
and then using the representation of a periodic sequence of delta functions
lim
M
sin(M +
1
2
)x
sin
1
2
x
= 2

n=
(x 2n) .
You can see this latter result by noting the value is very large, 2M+1, at x = 2n
where the denominator goes to zero, falling to zero over the narrow distance /M
and the integral is
1
2
_

sin(M +
1
2
)x
sin
1
2
x
dx = 1 .
(Plot the function for M = 10 if you do not believe this!) Also

H(,
2
T
) is:

H(,
2
T
) =
_
T
0
e
it
dt = e
iT/2
sin
1
2
T
1
2

.
January 16, 2000
Bibliography
[1] Numerical Recipes by W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T.
Vetterling, (Cambridge University Press, 2nd Edition) p. 542
7

You might also like