Biomedical Engineering Reference
In-Depth Information
is a sum of squares of two independent variables of normal distribution; hence each
frequency component of estimator
S
(
f
)
has a distribution given by:
S
χ
2
2
(
f
)
)
=
(2.35)
S
(
f
where χ
2
is a chi-square statistics of two degrees of freedom (corresponding to the
real and imaginary parts of Fourier transform) [Bendat and Piersol, 1971]. Please
note that the above expression is independent of the observation time
T
, which means
that increasing of epoch
T
does not change the distribution function describing the
error of the estimator. The increase of the time of observation
T
increases only the
number of the frequency components of the spectrum. It means that the estimator
of the spectral density obtained by means of Fourier transform is biased; also the
error of the estimator is high. Namely for a given frequency
f
1
the relative error
of
S
μ
S
f
1
. For the distribution χ
n
: σ
2
(
f
1
)
is: ε
r
=
σ
S
f
1
/
=
2
n
and
μ
=
n
,where
n
is a
number of degrees of freedom. For
n
=
2wegetε
r
=
1, which means that each single
frequency component of
S
has a relative error 100%. In consequence spectral
power calculated by means of Fourier transform is strongly fluctuating.
In order to improve statistical properties of the estimate two techniques were in-
troduced. The first one relies on averaging over neighboring frequency estimates.
Namely we calculate smoothed estimator
S
k
of the form:
(
f
)
1
l
[
S
k
=
S
k
+
S
k
+
1
+
···
+
S
k
+
l
−
1
]
(2.36)
Assuming that the frequency components
S
i
are independent, estimator
S
k
is char-
acterized by the χ
2
distribution of number degrees of freedom equal
n
=
2
l
.The
l
.
Another way of improving power spectral estimate is averaging the estimates for
successive time epochs:
relative standard error of frequency estimate will be therefore: ε
r
=
1
q
[
S
k
=
S
k
,
1
+
S
k
,
2
+
···
+
S
k
,
j
+
···
+
S
k
,
q
]
(2.37)
where
S
k
,
j
is the estimate of the frequency component
k
based on the time interval
j
.
The number of degrees of freedom in this case equals
q
; therefore the relative error
of single frequency estimate will be: ε
r
q
. This approach is known as Welch's
method and is implemented in MATLAB Signal Processing Toolbox as
pwelch
.
Both of these methods require stationarity of the signal in the epoch long enough
in order to: (i) perform averaging over frequency estimates, without excessive de-
creasing of frequency resolution or (ii) divide the epoch into segments long enough
to provide necessary frequency resolution.
An alternative method, which can be useful in case of relatively short time epochs,
is the multitaper method (MTM) described in [Thomson, 1982]. This method uses a
sequence of windows that are orthogonal to each other (discrete prolate spheroidal
=
Search WWH ::
Custom Search