NTIA Technical Report TR-09-458 

Statistical Considerations for Noise and 

Interference Measurements 

Roger Dalke 
 
 

NTIA Technical Report TR-09-458 

Statistical Considerations for Noise and 

Interference Measurements 

Roger Dalke 

U.S. DEPARTMENT OF COMMERCE 

Carlos M. Gutierrez, Secretary 

Meredith A. Baker, Acting Assistant Secretary 

for Communications and Information 

November 2008 

CONTENTS

Page

1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2. PRELIMINARY REMARKS ABOUT STOCHASTIC PROCESSES . . . . . 3

2.1 Stationarity and Ergodicity . . . . . . . . . . . . . . . . . . . . . . . 5

3. ESTIMATION OF THE MEAN . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.1 Stationary Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.2 Cyclostationary Processes . . . . . . . . . . . . . . . . . . . . . . . . 13

4. ESTIMATION OF THE AMPLITUDE PROBABILITY DENSITY . . . . . 21

4.1 Probability Coordinates for the APD . . . . . . . . . . . . . . . . . . 21

4.2 Accuracy of the APD Estimate . . . . . . . . . . . . . . . . . . . . . 22

5. SECOND-ORDER STATISTICS . . . . . . . . . . . . . . . . . . . . . . . . 25

5.1 Estimation of the Covariance Function . . . . . . . . . . . . . . . . . 26

5.2 Estimation of the Power Spectrum . . . . . . . . . . . . . . . . . . . 27

6. CONCLUDING REMARKS . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

7. REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

iii

STATISTICAL CONSIDERATIONS FOR NOISE AND INTERFERENCE

MEASUREMENTS

Roger Dalke∗

The study of noise and interference in the radio environment is essential to

the development of efficacious communications systems. Many of the characteristics of radio noise and interference of interest to radio system designers

can be expressed in terms of first- and second-order statistics. These statistics

are necessarily calculated using measured data. In this report, we discuss uncertainties that arise when measured data sets are used to calculate statistics

of radio noise and interference.

Key words: amplitude probability distribution; average power; measurement uncertainties;

power spectrum; radio interference; radio noise

1. INTRODUCTION

The performance of wireless communications systems is adversely affected by a variety of

environmental processes which from the point of view of a victim receiver are time-varying

stochastic processes. These processes can be either additive or multiplicative. Examples of

additive processes are natural and man-made noise and interference. Multiplicative processes

involve radio channel effects that result in signal fading such as multipath. In what follows,

we will simply use the term noise to describe these random processes.

Since communications systems need to resist noise, it is important that the salient characteristics of these processes be well understood. For many applications of interest to communications engineers, these characteristics can be expressed as properties of their first- and

second-order statistics. Such statistics are necessarily determined via measurements. Measurements, however, only provide finite data sets which introduce uncertainties in calculated

statistics.

The intent of this report is to discuss some of the statistical uncertainties that arise when

measured data are used to characterize radio noise and interference. In particular, we discuss

uncertainties for estimates of mean power, statistical distribution functions (viz. the amplitude probability distribution), covariance functions and power spectral density functions. It

should be noted that while the concepts described here are not new, they are presented in

a format that is intended to provide a useful reference for communications engineers performing noise and interference measurements. For a more detailed discussion of the topics

∗The author is with the Institute for Telecommunication Sciences, National Telecommunications and

Information Administration, U. S. Department of Commerce, Boulder, Colorado 80305.
described in this report, the interested reader is directed to the list of publications in the

reference section.

Section 2 of this report gives a brief overview of relevant concepts from the mathematical

theory of stochastic processes. This section also gives mathematical expressions for functions

and definitions used throughout the report. Estimators and related uncertainties for commonly measured first order statistics are described in Sections 3 and 4. Section 5 describes

the statistical properties of commonly used estimators for the covariance and power spectral

density functions.

2
2. PRELIMINARY REMARKS ABOUT STOCHASTIC PROCESSES

The study of noise in the radio environment is essential to the development of efficacious communications systems. From the standpoint of a victim receiver, noise processes are generally

random functions in time and space. The development of appropriate noise measurement

techniques and the proper interpretation of results is greatly enhanced by an understanding

of the mathematical theory of stochastic processes. The synergism of theory and measurements allows engineers to develop a keen understanding of the probabilistic structure of

important noise processes and provides a basis for analyzing the performance of radio systems. In this section we provide a brief discussion on pertinent aspects of the mathematical

description of stochastic processes.

Formally a stochastic process is defined in terms of a probability space (Ω, S,P) and a

parameter set T representing, for example, time variation (see [1] for a complete description

of this and what follows). Ω is the space of elementary events ω (e.g., possible outcomes of

an experiment), S is a σ-field of sets in Ω, and P is a probability measure. A stochastic

process is a finite function ξ(t, ω) which, for every fixed t ∈ T, is a measurable function of

ω ∈ Ω.

With regard to the use of the independent variable ω, we will adopt practices commonly

found in the literature. Specifically, the independent variable ω will be suppressed unless it

is required for calculations or it provides clarification. At times it will be useful to denote

the realization of a random variable or process based on a particular elementary event.

For example, ξ(t, ω1) denotes the realization (or sample function) resulting from the event

{ω = ω1}. Throughout this report ξ(t) is used to denote a random process. The specific

type or class of random process represented by the function ξ(t) will be defined as required

in each section or subsection.

For a fixed t = t1 there is a random variable ξ(t1) = ξ(t1, ω) with a definite probability

distribution denoted by

F(x;t1) = P{ξ(t1) ≤ x}

that is, F(x;t1) is the probability of the event that the random experiment yields a value

of ξ(t1) which does not exceed the real number x. For an arbitrary finite set of t-values,

t1, . . . , tn, the random variables ξ(t1), . . . , ξ(tn) have the joint distribution function

F(x1, . . . , xn;t1, . . . , tn) = P{ξ(t1) ≤ x1, . . . , ξ(tn) ≤ xn}.

All joint probability distributions for all n = 1, 2, . . . and all possible values of ti constitutes

the family of finite-dimensional distributions needed to completely characterize the process

ξ(t, ω). It is clear that completely characterizing many of the random noise processes of

interest in communications engineering would be quite difficult if not impossible. Evidently

a random process is described in terms of a denumerable number of random variables. In this

context, the statistical properties of one random variable (e.g., F(x1;t1)) are usually referred

to as first-order statistics and the joint statistics of 2 random variables (e.g., F(x1, x2;t1, t2))

are called second-order statistics.

3
Fortunately, it is often the case that first- and second-order statistics provide sufficient

information for predicting how noise affects the performance of communication systems. In

addition to the mean,

μ(t) = E {ξ(t)}

where E denotes expected value, the usual first-order statistic of interest is the amplitude

probability distribution (APD) which is the compliment of the distribution function, i.e.,

1 − F(x;t), for the process amplitudes.

Second-order statistics of interest are the correlation function

r(t1, t2) = E {ξ(t1)ξ

∗

(t2)} ,

the covariance function

cov(ξ(t1), ξ(t2)) = γ(t1, t2) = E {ξ(t1)ξ

∗

(t2)} − E {ξ(t1)} E {ξ

∗

(t2)} , (1)

the variance

var(ξ(t)) = cov(ξ(t), ξ(t))

and the Fourier transform of the covariance function

Γ(f1, f2) = F{γ(t1, t2)}.

The superscript ∗ used above denotes complex conjugate.

If the process of interest is discrete, the time parameter is usually denoted by an integer

subscript (e.g., ξk(ω)). The extension of the previous equations to discrete processes is

obvious, for example the discrete covariance is

γij = E

n

ξiξ

∗

j

o

− E {ξi} E nξ

∗

j

o

. (2)

At this point it is useful to introduce some standard terminology used in the text. Let M

denote the average of N random samples of the random variable ξ where E {ξ} = μ (i.e., M

is an estimator for the mean). If E {M} = μ the estimator is called unbiased. If E {M} 6=

μ, but we have E {M} → μ as N → ∞ then the estimator is said to be asymptotically

unbiased. It is also important to describe convergence properties of the estimator. Meansquare convergence means that var(M) = E {(M − μ)

2} → 0 as N → ∞. If in addition we

have P{|M − μ| < } = 1 for any  > 0 then we will say that M → μ with probability

one or in probability. An estimate that converges in probability is also called a consistent

estimate. Note that convergence in probability is a stronger condition than mean-square

convergence (see for example the random sampling variance calculation in Section 3.2.4).

We will sometimes use the standard symbol O when describing the order of magnitude of a

function when the independent variable tends to some limit (usually infinity). For example

var(M) = O(N −1) means that Nvar(M) remains bounded as N → ∞.

4
2.1 Stationarity and Ergodicity

The family of processes that are invariant under a translation in time (or space) are called stationary processes. This property greatly simplifies measurement procedures. Different classes

of stationarity are obtained by imposing more or less stringent conditions of invariance. A

process is strictly stationary when the whole family of finite-dimensional distributions is

invariant under a translation in the parameter of interest (e.g., time). Communications engineers are often concerned with the broader class of processes that have time invariant firstand second-order statistics. This class is usually referred to as wide-sense stationary. In the

sequel, these will simply be called stationary processes. Thus, for stationary processes we

have the following relations

μ(t) = μ F(x;t) = F(x) r(t1, t2) = r(t1 − t2) γ(t1, t2) = γ(t1 − t2) Γ(f) = F{γ(t1 − t2)}.

Since the variance of a stationary process plays an important role in characterizing uncertainties we will at times simply call it σ

2

, i.e., σ

2 = γ(0). Note that in the case of the covariance

and correlation function, we have used the same symbol to represent both a function of two

variables and a single variable which is the difference of the two independent variables. We

will frequently use such simplified notation when the relationship between the two functions

is obvious as with stationary processes.

While, from the standpoint of performing experiments, the attribute of stationarity is highly

desirable, many of the noise processes encountered in radio communications are nonstationary. For example, radio noise originating from human activity can vary significantly within

the hour as well as hour-to-hour. However, it is usually possible to adequately characterize

these processes using measurements over time scales during which the process is approximately stationary. For longer time scales, the first- and higher-order statistics are usually

treated as random variables.

In principle, the statistical properties of a noise process can be calculated from an ensemble

of sample functions obtained from a sufficiently large number of random events. Implementing the experiments required to construct such an ensemble is usually not practical and

may not even be possible. An obvious and more practical procedure is to try and obtain

the desired statistics from a sample function ξ(t) = ξ(t, ω1) corresponding to a single random experiment. Such a procedure is valid if each sample function is representative, when

viewed over sufficiently long times, of the statistical properties of the entire ensemble—and

conversely. Random processes with this property are known as ergodic processes. In cases

where it is only possible to measure a single sample function, and the process is known to be

stochastic, the usual practice is to assume the process is ergodic and construct the desired

statistics from the measured sample function.

Formally, for a strictly stationary discrete time process ξk(ω) with a bounded first moment

(E {|ξk|} < ∞), a point in the basic probability space ω = (ω0, ω1, ω2, . . .) can be thought of

as a possible time history of an event, i.e., ξk(ω) = ωk. If T is a time shift transformation,

(T ω)k = ωk+1, and if for every invariant set A ∈ Ω (i.e., T

−1A = A) P{A} = 0 or 1, the

5
ergodic theorem [2] states that

lim

N→∞

1

N

N

X−1

k=0

ξ0(T

kω) = E {ξ0}

with probability one.

More generally, this includes time averages of bounded (in the first moment) integrable

functions f of the process, i.e.,

lim

N→∞

1

N

N

X−1

k=0

f(T

kω) = E {f} .

For example, if the first-order distribution function F(x) is estimated by a function Fx which

gives the fraction of time that the observed sample function is less than a particular level x,

then ergodicity implies that Fx → F(x) as N → ∞.

2.1.1 A Strictly Stationary and Ergodic Process

Consider the process ξ0, ξ1, ξ2 . . . where the {ξk} are independent and identically distributed

(iid) random variables with E {|ξk|} < ∞ and form the time average

M =

1

N

N

X−1

k=0

ξk.

By the strong law of large numbers |M − E {ξk} | → 0 with probability one.

2.1.2 A Strictly Stationary Process that is Not Ergodic

Let ω be uniformly distributed over the unit interval (ω ∈ [0, 1)) and define the random

variable X as a real function that can be represented as the Fourier series

X(ω) =

+

X∞

n=−∞

cne

i2πnω

where i =

√

−1. For any λ ∈ [0, 1) and transformation T ω = ω + λ (mod 1) obtain the

random process

ξk(ω) = X(T

kω) k = 0, 1, 2, . . . .

Since T is a measure preserving transformation (P{T A} = P{A} for any set A ∈ Ω), the

process is stationary and E {ξk} = c0. The time average is

M =

X

nλ∈Z

cne

i2πωn +

1

N

X

nλ /∈Z

cne

i2πωn 

e

i2πNnλ − 1

e

i2πnλ − 1

!

,

6
where Z is the set of integers. In the limit N → ∞ we have

M →

X

nλ∈Z

cne

i2πωn

.

If λ is a rational number, the process is periodic and there are non-empty invariant sets that

do not have probabilities of 0 or 1. Hence, the transformation is not ergodic. Clearly the

time average M does not in general converge to c0.

If, on the other hand, λ is irrational, the process is not periodic and it can be shown that

the invariant sets have probability 1 or 0. The transformation is ergodic and M → c0 with

probability one [2].

7
3. ESTIMATION OF THE MEAN

Perhaps the first order statistic that is most widely used to characterize radio noise is the

mean power. In this section, we will examine the statistical properties of the commonly used

power estimator. Expressions for the mean and variance of the estimator will be obtained.

We will see that the variance depends on correlation properties of the data and hence data

sampling strategies. By way of an example, these results will be applied to narrowband

Gaussian radio noise.

When the distribution of the estimator is known or can be reasonably approximated, a

more precise statistical description of uncertainties can be obtained. An example of such a

procedure, the method of confidence intervals, is also given in this section.

Finally, we examine the mean and variance of the estimator when the underlying process is

cyclostationary. This class of nonstationary processes often occurs in the analysis of communications systems. This is followed by an application of the results to a pulse modulated

process.

3.1 Stationary Processes

In this subsection, we describe uncertainties associated with estimating the mean of a stationary process. Let ξ(t, ω) be such a process with mean μ = E {ξ(t)} . The usual procedure

is to calculate the time average of a measured sample function ξ(t, ω1) based on a particular

experiment ω1. If we obtain N discrete samples equally spaced ∆t time units apart and

denote ξk(ω1) = ξ(k∆t, ω1), the estimate of the mean is

M(ω1) = 1

N

N

X−1

k=0

ξk(ω1). (3)

Note that the estimator M is a function on the probability space consisting of elementary

events ω ∈ Ω, i.e., it is a random variable. Hence, the characteristics of the estimate are

properly described in terms of its statistics.

The mean of the estimate is

E {M} =

1

N

N

X−1

k=0

E {ξk} = μ , (4)

hence M is an unbiased estimate. The variance is

var(M) = 1

N2

N

X−1

j=0

N

X−1

k=0

E {ξj ξ

∗

k} − |μ|

2 =

1

N2

N

X−1

j=0

N

X−1

k=0

γjk , (5)

where γjk is given by Equation 2.

8
For a real stationary process γjk = γn, where n = k − j, so that Equation 5 simplifies to

var(M) = σ

2

N

+

2

N2

X

N

n=1

(N − n)γn

where σ

2 = γ0.

If the process is stationary and

1

N

N

X−1

k=0

γk → 0

the estimator is mean-square convergent. This is sometimes referred to as mean-square

ergodicity. If in addition γ(t) = O(|t|

−β

) for some β > 0 as |t| → ∞, then M → μ

with probability one (i.e., the estimate is consistent) [1]. It should also be noted that if

γn = 0 for n 6= 0, the variance of the estimator is simply

var(M) = σ

2

N

. (6)

3.1.1 Sample Correlation Effects

As we have seen, when the samples are uncorrelated, the variance of the estimator decreases

in proportion to the number of data samples. In this subsection we will compare this result

with what would be obtained using a fixed length record of densely sampled data. In fact,

we will assume that the number of samples is so large that the summation in Equation 3

can be replaced with an integral.

To obtain uncorrelated samples we need to know or otherwise estimate the minimum sample

spacing τ0 that yields (at least approximately) uncorrelated samples. This sample spacing is

often referred to as correlation time and we assume that γ(τ ) ≈ 0 when τ ≥ τ0. A practical

definition of correlation time is

τ0

2

=

1

γ(0) Re Z∞

0

γ(τ )d τ. (7)

As might be anticipated, there is a relationship between correlation time and the bandwidth

of the power spectrum. If the equivalent bandwidth of the power spectrum Γ(f) is defined

as

Beq =

1

Γ(0)

Z∞

−∞

Γ(f)df

then

2τ0 =

1

Beq

.

9
In the case of a triangular function, 2τ0 is exactly equal to the width of the normalized

covariance function. More generally, the quantity 2τ0 can be considered as a measure of the

width of the covariance function.

If the time increment is greater than the correlation time, we can simply calculate the

variance using Equation 6. Note that evaluating the variance requires a knowledge of σ

2

which is often not directly available. When N is large, a common practice is to replace σ

2

with the unbiased estimate of the variance calculated directly from the measured data

S

2 =

1

N − 1

N

X−1

k=0

(ξk − M)

2

(8)

which converges to σ

2

in probability.

If a real stationary and ergodic noise process is densely sampled over a time period T, the

variance of the estimate can be calculated via integration

var(M) = 1

T

2

Z

T

0

Z

T

0

γ(s − t) dt ds =

1

T

Z

T

−T

 

1 −

|τ |

T

!

γ(τ )d τ. (9)

If τ0 is finite and T >> τ0 we have

var(M) ≈

2

T

Z∞

0

γ(τ )d τ =

γ(0)

T/τ0

. (10)

But T/τ0 is approximately equal to the maximum number of uncorrelated samples contained

in the sample time period. Hence, the variance for dense sampling is not much different than

what would be obtained for sparse sampling using a sample spacing of τ0.

3.1.2 An Example – Narrowband Noise

In this subsection we will use the results obtained so far to calculate the variance of mean

power estimators assuming densely and sparsely sampled narrowband Gaussian noise.

Let n(t) be a stationary narrowband Gaussian noise process with E {n(t)} = 0, center

frequency fc, and baseband representation ˆn(t), i.e.,

n(t) = Re ˆn(t)e

−i2πfct

.

Let ξ(t) = |nˆ(t)|

2 and assume that we wish to estimate the mean power μ = E {ξ(t)} using

Equation 3. We know from Equation 4 that the estimate is unbiased. The variance of the

estimate depends on the covariance of ξ(t), so we will first obtain an expression for the

covariance in terms of the baseband noise process. From Equation 1 we immediately obtain

γ(s − t) = E {nˆ(t)ˆn

∗

(t)ˆn(s)ˆn

∗

(s)} − μ

2

10
This result can be simplified using the following identity for stationary zero mean Gaussian

processes [3]

E {nˆ(t)ˆn

∗

(s)ˆn

∗

(t)ˆn(s)} = E {nˆ(t)ˆn

∗

(s)} E {nˆ

∗

(t)ˆn(s)} + E {nˆ(t)ˆn(s)} E {nˆ

∗

(t)ˆn

∗

(s)} + μ

2

to obtain

γ(τ ) = |E {nˆ(t)ˆn

∗

(s)} | 2 + |E {nˆ(t)ˆn(s)} |2

where τ = s − t.

Often the noise equivalent bandwidth is known (at least approximately) so it is useful to

formulate the variance of the estimator in terms of this parameter. We begin by noting that

narrowband Gaussian noise can be represented as bandpass filtered white noise z(t) with

spectral density kT where k is Boltzmann’s constant (1.3807 × 10−23JK−1) and T is the

temperature in Kelvin. If hˆ is the baseband representation of the bandpass filter, we have

nˆ(t) = Z∞

−∞

z(s)e

i2πfcshˆ(t − s)ds.

Since E {z(t)z(s)} = kT δ(s − t) we find

E {nˆ(t)ˆn

∗

(s)} = kT

Z∞

−∞

hˆ(τ + s)hˆ∗(s)ds

and

E {nˆ(t)ˆn(s)} = kT e

i4πfct

Z∞

−∞

e

−i4πfcshˆ(τ + s)hˆ(s)ds.

Using the notation Hˆ (f) = F{hˆ(t)}, the integral given above can readily be evaluated as

Z∞

−∞

Hˆ (2fc − g)Hˆ (g)e

−i2πgτdg ≡ 0

and hence E {nˆ(t)ˆn(s)} ≡ 0 so that

γ(τ ) =













kT

Z∞

−∞

hˆ(τ + s)hˆ∗(s)ds













2

.

Usually, the filter is specified so that the maximum gain is 1 (i.e., max|Hˆ (f)|

2 = 1. The

mean and variance of ξ(t) can then be expressed in terms of the filter’s equivalent bandwidth

B(H)

eq =

R ∞

−∞ |Hˆ (f)|

2df as follows

E {ξ(t)} = μ = (kT B

(H)

eq )

11
and

γ(0) = (kT B

(H)

eq )

2

.

Using Equation 7 and Parseval’s theorem the correlation time can be written as

τ0 =

1

(B

(H)

eq )

2

Z∞

−∞

|Hˆ (f)|

4

df

from which we find that τ0 ≤ 1/B(H)

eq . In the case of dense sampling over a time period T,

Equation 10 gives

var(M) ≤

(kT B(H)

eq )

2

T B(H)

eq

.

In the case of sparse uncorrelated samples we have from Equation 6

var(M) ≈ μ

2

/N = (kT B

(H)

eq )

2

/N.

3.1.3 Confidence Intervals

If the sampling distribution of the estimate M is known, it is possible to calculate the endpoints of an interval that with some specified probability contains the value of the population

mean. Such a procedure is known as the method of confidence intervals.

Given a small number ε > 0 and confidence coefficient p = 1−ε, let c1(M) and c2(M) be the

endpoints of an interval so that P{c1(M) ≤ μ ≤ c2(M)} = p. This is equivalent to the known

probability P{η1(μ) ≤ M ≤ η2(μ)} if we choose c1(M) = η

−1

2

(M) and c2(M) = η

−1

1

(M).

For example, if a large number of independent and identically distributed random samples

are used, M is approximately normal with mean μ and standard deviation σ/√

N. For large

N it follows that

P{M − λε

σ

√

N

< μ < M + λε

σ

√

N

} = P{μ − λε

σ

√

N

< M < μ + λε

σ

√

N

} = p

where λε is the 100ε % value of the standard normal deviate Z = (M − μ)/(σ/√

N), i.e.,

ε = P{|Z| > λε} =

2

√

2π

Z∞

λε

e

−t

2/2

dt.

Now, for a given estimate M and confidence coefficient p we can determine the endpoints

M ± λεσ/√

N of the interval that with probability 100p % contains μ.

As before, if the value of σ is not known and N is large, a common practice is to use the

unbiased estimate obtained from Equation 8.

12
3.2 Cyclostationary Processes

Many nonstationary processes of interest in communications exhibit periodically time-variant

statistics. Examples include pulse-amplitude, -width, and -position modulation, amplitude

modulation, phase modulation, and periodic noise pulses. A process ξ(t) is said to be

cyclostationary if the mean and correlation are periodic. If the period is T we have

μ(t) = μ(t + T)

r(t, s) = r(t + T, s + T),

(11)

and if R T

0

|μ(t)|

2dt < ∞, the mean can be represented as a Fourier series

μ(t) = X∞

n=−∞

μne

i2πnt/T

.

The Fourier coefficient μn is sometimes referred to as a cyclic mean and can be defined as

μn = lim

Z→∞

1

Z

Z/2

Z

−Z/2

μ(t)e

−i2πnt/T dt. (12)

The covariance function γ(t, s) = r(t, s) − μ(t)μ(s) is also periodic. It is sometimes useful

to define the covariance in terms of clock time t and lag τ = s − t as follows

γ ̃(t, τ ) = γ(t − τ/2, t + τ/2)

and if R T

0

|γ ̃(t, τ )|

2dt < ∞ as a Fourier series

γ ̃(t, τ ) = X∞

n=−∞

γ ̃n(τ )e

i2πnt/T

.

with Fourier coefficients

γ ̃n(τ ) = limZ→∞

1

Z

Z/2

Z

−Z/2

γ ̃(t, τ )e

−i2πnt/T dt. (13)

As before, it is useful to define a correlation time. Since the process is not stationary, the

average correlation time can be defined in terms of clock time averages as follows

τ0

2

=

h

R ∞

0

γ ̃(t, τ )dτ i

hγ ̃(t, 0)i

=

R ∞

0

γ ̃0(τ )dτ

γ ̃0(0) .

13
3.2.1 Dense Sampling

Equation 12 indicates that an estimate of the cyclic mean μn can be obtained from a sample

function as follows

Mn(Z, ω1) = 1

Z

Z/2

Z

−Z/2

ξ(t, ω1)e

−i2πnt/T dt.

The expected value of the estimate is asymptotically unbiased

E {Mn(Z)} = μn +

X

k6=n

μk sinc(π(k − n)Z/T) → μn (14)

with variance

var(Mn(Z)) = 1

Z2

Z/2

ZZ

−Z/2

γ(t, s)e

−i2πn(t−s)/T dt ds. (15)

If limZ→∞ var(Mn(Z)) = 0, then the process is mean-square cycloergodic at the frequency

n/T.

3.2.2 Sparse Uncorrelated Samples

Suppose that ξ(t) is cycloergodic at all frequencies of interest and that γ(t, s) ≈ 0 when

|s − t| ≥ τ0. We can then use approximately uncorrelated samples to obtain an estimate of

the cyclic mean as follows

Mn =

1

N

N

X−1

k=0

ξ(tk)e

−i2πntk/T

where |tk − tk±1| ≥ τ0. If we let

tk = kT(J + ε) (16)

where T J ≥ τ0, for some positive integer J, and 0 < ε < 1, the estimate can be written as

Mn =

1

N

N

X−1

k=0

ξ(tk)e

−i2πnkε

with expected value

E {Mn} =

1

N

N

X−1

k=0

X∞

j=−∞

μj e

i2π(j−n)kε =

X

jε∈Z

μj+n +

X

jε /∈Z

μj+n

1

N

e

i2πNjε − 1

e

i2πjε − 1

.

In the limit we find that the estimate is in general biased

lim

N→∞

E {Mn} =

X

jε∈Z

μj+n.

14
However, if ε is irrational or if μ(t) is bandlimited so that μj+n = 0, |j| ≥ b1/εc, the estimate

is asymptotically unbiased.

For practical applications, we assume that the increment is a small rational number which

can be written as ε = 1/K, where K is a positive integer. For this case, the limiting expected

value of the estimate can also be written as a discrete Fourier transform

lim

N→∞

E {Mn} =

1

K

K

X−1

k=0

μ(kT/K)e

i2πnk/K.

We see that, for sparse uncorrelated samples, the expected value of Mn is equivalent to

estimating μn using the discrete Fourier transform as an approximation to Equation 12. It

should be noted that we must also have K >> n when estimating higher order cyclic means.

The variance of the estimate is

var(Mn) = 1

N2

N

X−1

k=0

γ(tk, tk). (17)

Since the covariance is periodic, we can use Equation 16 to obtain

1

N

N

X−1

k=0

γ(kT ε, kT ε) = 1

N

N

X−1

k=0

γ ̃(kT ε, 0).

As before, if we let K = 1/ε and N is large, we need only average the covariance function

over a period and therefore

1

N

N

X−1

k=0

γ(kT ε, kT ε) →

1

K

K

X−1

k=0

γ ̃(kT/K, 0).

Note that 1

K

PK−1

k=0 γ ̃(kT/K, 0) is the discrete Fourier transform approximation of  ̃γ0(0) (see

Equation 13) so that when ε is sufficiently small the variance can be approximated as

var(Mn) ≈

γ ̃0(0)

N

. (18)

This result is similar to that obtained for uncorrelated samples of a stationary process

(Equation 6) with the zero-lag stationary process covariance function replaced by the zerolag time average of the cyclostationary process covariance function.

As with stationary processes, calculating the uncertainty in the sparse sampling estimate

requires additional knowledge of the process, in this case  ̃γ0(0), which is often not directly

available. When N is large,  ̃γ0(0) can be replaced by an estimate of the form

S ̃

0(0) = 1

N − 1

N

X−1

k=0

|ξ(tk) − M0|

2

.

15
3.2.3 Random Sampling

In this subsection we will show that a cyclostationary process can be converted into a stationary process by random sampling. The average of the stationary process is an unbiased

estimate of μ0 and the results of Section 3.1 can be used to calculate uncertainties for this

estimate.

By random sampling, we mean that the observation times are random and uniformly distributed over a period T. This gives a new process defined as follows

X(t) = ξ(t + θ), θ uniform over [0, T).

The new process is stationary since the mean is constant

E {X(t)} = E {E {ξ(t + θ)| θ}} =

1

T

Z

T

0

μ(θ)dθ = μ0 (19)

and the correlation function depends only on the lag τ = s − t

E {X(t)X

∗

(s)} = E {E {ξ(t + θ)ξ

∗

(s + θ)| θ}}

=

1

T

Z

T

0

r(t + θ, s + θ)dθ

=

1

T

Z

T

0

r(θ, s − t + θ)dθ = r(τ ).

Note that the mean and correlation functions are time averages of the cyclostationary statistics.

The covariance function can be written in terms of the correlation function and the mean as

follows

γ(τ ) = r(τ ) − μ

2

0

.

3.2.4 An Example – Pulse Modulation

In this subsection, we will calculate the variance of the mean power estimator for a baseband

pulse modulated process using the different sampling schemes described above. The process

used in this example has the form

nˆ(t) = X∞

n=−∞

Zn p (t − nT)

where p(t) is real and is zero outside of the interval [0, T) and Zn are iid complex random

variables. The instantaneous power can be written as

ξ(t) = X∞

n=−∞

Wn p

2

(t − nT)

16
where Wn = |Zn|

2

.

To calculate the variance of the estimator we first need to obtain expressions for the process

mean and correlation functions. Using the notation μw = E {Wn} and σ

2

w = E {W2n} − μ

2

w,

we have

E {WnWm} = σ

2

wδnm + μ

2

w. (20)

The process mean is simply

μ(t) = E {ξ(t)} = μw

X∞

n=−∞

p

2

(t − nT).

and the correlation function is

r(t, s) = E {ξ(t)ξ(s)} =

X∞

n=−∞

X∞

m=−∞

E {WnWm} p

2

(t − nT) p

2

(s − mT).

Using Equation 20 gives

r(t, s) = μ(t)μ(s) + σ

2

w

X∞

n=−∞

p

2

(t − nT)p

2

(s − nT)

and hence

γ(t, s) = σ

2

w

X∞

n=−∞

p

2

(t − nT)p

2

(s − nT).

Clearly, this process is cyclostationary and has the periodicity as described in Equation 11.

We now assume that the estimator is obtained from a densely sampled record of length Z

and calculate the variance of the estimator. To minimize the bias (Equation 14) we require

Z >> T. It is useful to approximate Equation 15 in terms of the minimum number of

intervals of length T in the record. Using the floor function notation define

2N =



Z

T



and hence

var(Mn(N)) ≈

σ

2

w

(2NT)

2

X∞

m=−∞















Z

NT

−NT

p

2

(t − mT)e

−i2πnt/T dt















2

.

We can sum the series (approximately) and write the expression for the variance in terms of

the Fourier transform of the pulse squared

var(Mn(N)) ≈

σ

2

w

2NT2













Z∞

−∞

p

2

(t)e

−i2πnt/T dt













2

.

17
Applying the commonly used rectangular pulse

p(t) = prect (t) = (

1 0 ≤ t < T0, T0 < T

0 else (21)

gives

var(Mn(N)) ≈

σ

2

wT

2

0

2NT2

sinc2(πnT0/T). (22)

Next we will assume that the estimator is obtained from N sparse uncorrelated samples.

Using Equation 17 gives

var(Mn) = σ

2

w

N2

N

X−1

k=0

X∞

n=−∞

p

4

(kT ε − nT).

For N >> K and ε = 1/K << 1, we can use Equation 18 and Equation 13 to obtain

var(Mn) ≈

γ ̃0(0)

N

=

σ

2

w

NT

Z

T

0

p

4

(t)dt.

For the rectangular pulse of Equation 21

var(Mn(N)) ≈

σ

2

wT0

NT

. (23)

So far we have compared the variance of the mean estimator when dense and sparse uncorrelated sampling schemes are used. Often we only need to estimate the time average mean

power μ0. Such an estimate is readily obtained by using random sampling. In what follows,

we will calculate the variance of the estimator M0 obtained from random samples of the

pulse modulated processes.

We begin with the modified pulse modulated process

X(t) = X∞

n=−∞

Zn p

2

(t − nT + θ)

where θ is uniformly distributed over the period T. Using the results of Section 3.2.3 we can

calculate the process mean and correlation, viz.,

E {X(t)} =

μw

T

Z∞

−∞

p

2

(θ)dθ = μ0

and the correlation as a function of lag τ = s − t

r(τ ) = r1(τ ) + r2(τ )

18
where

r1(τ ) = σ

2

w

X∞

n=−∞

E

n

p

2

(t + θ − nT)p

2

(s + θ − nT)

o

r2(τ ) = μ

2

w

X∞

n=−∞

X∞

m=−∞

E

n

p

2

(t + θ − nT)p

2

(s + θ − mT)

o

.

These equations reduce to

r1(τ ) = σ

2

w

T

Z∞

−∞

p

2

(θ)p

2

(τ + θ)dθ

r2(τ ) = μ

2

w

σ

2

w

X∞

n=−∞

r1(τ − nT).

Note that r2(τ ) is a scaled periodic extension of r1(τ ). Using prect(t) defined by Equation 21

gives

r1(τ ) = (

σ

2

w(T0 − |τ |)/T |τ | ≤ T0

0 else

i.e., the correlation function consists of a single triangle function of height σ

2

wT0/T and base

2T0 centered at τ = 0 and the periodic extension of a triangle function of height μ

2

wT0/T

and base 2T0. The covariance function is

γ(τ ) = r(τ ) − μ

2

w



T0

T

2

.

Clearly, we have

lim

Z→∞

1

Z

Z/2

Z

−Z/2

γ(τ )dτ = 0

so the process is mean-square ergodic. However, since the covariance function itself does not

tend to a limit, the time average does not tend to the mean with probability one.

If we choose the special case T0 = T then

r1(τ ) = (

σ

2

w(1 − |τ |/T) |τ | ≤ T

0 else

r2(τ ) = μ

2

w

The process correlation time is τ0 = T and the time average tends to the mean with probability one. If we obtain N uncorrelated samples at time intervals of τ0 (or greater)the variance

of the estimate is

var(M) = σ

2

w

N

.

19
In the case of dense sampling over the same time window (i.e., TN = Nτ0), we have (Equation 9)

var(M) = 2σ

2

w

TN

Z

TN

0



1 −

τ

TN



T − τ

T

dτ =

2σ

2

w

N

Z

1

0



1 −

τ

N



(1 − τ ) dτ

which reduces to

var(M) = σ

2

w

N



1 −

1

3N



.

Hence even for small N, the variance is essentially the same as that for sparse uncorrelated

samples using a sample spacing of τ0.

It is useful to compare this result with the estimator variance given in Equation 22 and

Equation 23 for M0 and T0 = T. Clearly, the random sampling variance and the sparse

sampling variance are the same (assuming N is large). To compare this with the variance

of the dense sampling estimator we note that 2N in Equation 22 is essentially the same as

TN /T given above. If the sample spacing is equal to the correlation time τ0 = T we find that

the variance is approximately the same for all the sampling schemes. It should be noted,

however, that the random sampling estimator is unconditionally unbiased (Equation 19)

while the other estimators are only asymptotically unbiased.

20
4. ESTIMATION OF THE AMPLITUDE PROBABILITY DENSITY

Quantile statistics play an important role in characterizing radio noise processes and in

predicting the performance of radio receivers. Typically, radio engineers are interested in

the percent of time that noise amplitudes exceed a specified threshold. The threshold may

represent, for example, the point at which the symbol error rate exceeds the capability

of error correction algorithms. Radio engineers often refer to such statistics as amplitude

probability distributions. The APD is the compliment of the first-order distribution function

of the complex baseband noise envelope expressed in decibels. In this section, we define the

APD estimator and derive expressions for calculating quantile confidence intervals.

An estimate of the APD can be obtained from an experiment which yields discrete data

samples from a complex baseband noise process ˆn(t, ω). Here, we assume that the samples

used to construct the APD are approximately independent. This is usually accomplished

in practice by designing the experiment so that the time between data samples exceeds the

process correlation time. We have then a random sample consisting of a sequence of N

observed values of the noise envelope

ξk(ω1) = |nˆ(tk, ω1)| k = 1, 2, · · · , N |tk+1 − tk| ≥ τ0.

If the process is stationary and ergodic, the ξk can be considered as iid random variables

with distribution F(x) and P{ξk > x} = 1 − F(x).

If we denote by ν the number of sample values that are greater than some specified value

x, we obtain a step-function estimate of the APD with an incremental step height of 1/N.

The probability that x is exceeded is simply the frequency of the event {ξk > x} in the

sequence of N observations. In practice, the APD estimate is easily calculated by ordering

the samples ξ[1] ≥ ξ[2] ≥ · · · ≥ ξ[N] where ξ[ν] denotes the ν

th order-statistic and assigning

the probabilities

P{ξk > ξ[ν]} =

ν

N

.

Obviously the total number of samples needed is predicated on the desired accuracy of the

probability estimate. This will, however, vary with the quantile, being most accurate for

the median and least accurate for the tails. For example, it would seem quite unlikely that

P{ξk > ξ[N]} = 1 which indicates that the endpoints of the APD estimate are dubious as

is shown in the quantitative analysis which follows.

4.1 Probability Coordinates for the APD

Since the amplitude of a complex Gaussian process (e.g., radio receiver noise) is Rayleigh

distributed, the usual engineering practice is to graph APDs on probability paper that

linearizes Rayleigh distributed amplitudes. This allows the engineer to determine, at a

glance, if the measured process is Rayleigh distributed and if not, how and why it deviates

from a Rayleigh distribution. To determine the scales for Rayleigh graph paper we begin

21
with the distribution of the amplitude (squared) of a stationary complex Gaussian process

P{|nˆ(t)|

2 ≤ v2

} = e

−v

2/v2

0

where v

2

0 = E {|nˆ(t)|

2}. Engineers prefer quantities in decibels, e.g., W(t) = 20 log |nˆ(t)| in

which case the distribution becomes

P = P{W(t) ≤ w} = e

−10w/10/v2

0 .

If we take for the abscissa scale x = −10 log ln 1/P and the ordinate scale y = w we have

the desired linear relationship x = −y + w0, where w0 = 20 log v0.

4.2 Accuracy of the APD Estimate

As in the case of the mean, it is important to make quantitative remarks about the accuracy

of the APD estimate. In this section we will discuss both convergence of the estimate and a

method for calculating confidence intervals (ref [4]). It is important to note that the following

concepts require that the data ξk are iid observations.

To begin, we define the quantile xp as

P{ξk > xp} = 1 − F(xp) = p

and the related random variables

ηk =

(

1 if ξk > xp

0 else and ν =

X

N

k=1

ηk.

The random variable ν is just the number of sample values greater than xp in N observations.

The probability that exactly m observations are greater than xp and N − m observations

are less than or equal to xp is p

mqN−m where q = 1 − p. The number of such sample sets is

 

N

m

!

=

N!

m!(N − m)!

hence, the probability that a sample set contains exactly m observations greater than xp,

irrespective of the order, is

P{ν = m; N} =

 

N

m

!

p

mqN−m.

By De Moivre’s Theorem [5] ν/N is asymptotically normal with E {ν/N} = p and var(ν/N) =

p q/N. We see immediately that when the ξk are iid, as N → ∞, the estimate converges in

probability, i.e., P{|ν/N − p| < } = 1 for any  > 0.

To calculate confidence intervals for a population quantile xp we can use the following relations

P{ν ≥ m; N} =

PN

k=m



N

k



p

k

q

N−k = b(m)

P{ν < m; N} = 1 − b(m).

22
Now let m− and m+ be the indices for two ordered samples so that ξ[m−] > ξ[m+]. Then it

follows that

P{ξ[m−] ≥ xp ≥ ξ[m+]; m+ > m−} = P{ν ≥ m−} − P{ν ≥ m+} = b(m−) − b(m+)

hence, we can say that with (b(m−) − b(m+)) · 100% confidence the population quantile xp

is in the interval [ξ[m+], ξ[m−] ]

For example, given N ≥ 2 samples, the probability that the median falls between the smallest

and the largest sample is given by

P{ξ[1] ≥ θ1/2 ≥ ξ[N]} = 1 −

1

2

N−1

.

4.2.1 Approximation for Large N

If Npq >> 1 and |m − Np|

3/(Npq)2 << 1 the following approximation for the binomial

distribution can be used [6]

P{ν ≥ m} =

X

N

k=m

 

N

k

!

p

k

q

N−k → 1 − Φ

 

m − Np

√

Npq !

where

Φ(x) = 1

√

2π

Zx

−∞

e

−t

2/2

dt.

We can now specify a desired confidence c · 100% and calculate the indices of the order

statistics that bound the confidence interval. Obviously we require that

b(m−) − b(m+) = c.

In addition we want the mean value of ν (Np) for the quantile xp to be near the center of

the interval. This is accomplished by setting b(m+) + c/2 = 1/2 or

b(m+) = (1 − c)/2.

Applying these two conditions gives the indices of the order statistics that bound the confidence interval

m+ = dNp +

√

Npq Φ

−1



1+c

2



e

m− = bNp −

√

Npq Φ

−1



1+c

2



c .

The median m0 is given by P{ξ[m0] < x1/2} = 1/2, hence m0 = [N/2].

As an example, assume that we wish to find the 90% confidence interval for the median

given 2000 iid samples. In this case, Φ−1



1+c

2



= 1.64 and hence,

m+ = d(1000 + 36.7)e = 1037

m− = b(1000 − 36.7)c = 963

m0 = 1000

23
therefore

P{X[963] ≥ x1/2 ≥ X[1037]} ≥ 0.9 .

24
5. SECOND-ORDER STATISTICS

Second-order statistics provide information about the statistical dependence of a noise process at two different times. This information can be used, for example, to predict the

duration of undesirable correlated events and for the development of algorithms to mitigate

these effects. Fourier transforms of second-order statistics provide information about the

characteristics of the frequency content of the process.

An important and easily measured second-order statistic is the covariance. Specifically, the

covariance function gives a measure of the correlation of the process at two different times. If

the process is stationary, the covariance only depends on the time difference. Furthermore, if

all of the joint distribution functions of the process are normal, the covariance and the mean

are all that is required to completely characterize the process. If the process is stationary,

the Fourier transform of the covariance function is

Γ(f) = Z∞

−∞

γ(τ )e

−i2πfτdτ ,

which evidently shows how the mean power of the process is distributed over frequency, e.g.,

γ(0) = Z∞

−∞

Γ(f)df.

Hence, the function Γ(f) is usually referred to as the power spectrum.

In this section we will describe covariance and power spectrum estimators that are commonly

applied to noise measurement data. Mathematical expressions for the expected value and

covariance of these estimators will be given. Finally, we will obtain asymptotic statistical

distributions for discrete estimates of the power spectrum based on the discrete Fourier

Transform.

Before proceeding, it is useful to define some standard notation used in the remainder of this

report. First, the notation N(μ, σ2) used in conjunction with a real random variable ξ means

that ξ is normally distributed with mean μ and variance σ

2

. If ξ is complex normal we will

use the notation Nc(μ, σ2), where in general the mean is complex. Note that for identically

distributed independent complex random variables, Re ξ and Im ξ are N(Re E {ξ}, var(ξ)/2)

and N(Im E {ξ}, var(ξ)/2). We will also adopt the standard notation χ

2

n

to represent a chisquare random variable with n degrees of freedom, i.e., if ξk are a sequence of real N(0, 1)

independent random variables,

χ

2

n =

Xn

k=1

ξ

2

k

. (24)

Finally, the symbol Γn will be used to denote the discrete Fourier transform of the covariance

function defined as

Γn =

N

X−1

k=0

γke

−i2πnk/N .

25
5.1 Estimation of the Covariance Function

Suppose that ξ(t, ω) is a zero mean complex baseband stationary and ergodic noise process.

A commonly used method for estimating γ(τ ) = cov(ξ(t), ξ(t + τ )) is to obtain a sample

function ξ(t, ω1) of time duration T and calculate the time average

c(τ, ω1) = ( 1

T

R T −τ

0

ξ

∗

(t, ω1)ξ(t + τ, ω1)dt if 0 ≤ τ ≤ T

0 if τ > T (25)

and c(−τ ) = c

∗

(τ ). This estimator is a random process, hence, we need to determine its

mean and covariance.

Since the underlying process is stationary, the mean of the estimator is easily calculated

E {c(τ )} =

1

T

T

Z−τ

0

E {ξ

∗

(t)ξ(t + τ )} dt =

(

γ(τ )



1 −

τ

T



if 0 ≤ τ ≤ T

0 if τ > T

and E {c(−τ )} = E {c

∗

(τ )}. This result shows that the estimator is biased, however it is

asymptotically unbiased (i.e., E {τ} → γ(τ ) as T → ∞).

Calculating the covariance of the estimator is somewhat involved and is formally obtained

from the following expression

cov(c(τ1), c(τ2)) = 1

T

2

T

Z−τ1

0

T

Z−τ2

0

cov (ξ

∗

(t)ξ(t + τ1), ξ∗(s)ξ(s + τ2)) dt ds

with 0 ≤ τ1 ≤ T and 0 ≤ τ2 ≤ T.

To simplify this expression, the usual procedure is to invoke the fourth-order cumulant

κ = cum(ξ1, ξ2, ξ∗

3

, ξ∗

4

), and apply the following formula

cov(ξ1ξ2, ξ3ξ4) = cov(ξ1, ξ3)cov(ξ2, ξ4) + cov(ξ1, ξ4)cov(ξ2, ξ3) + κ.

It is a common practice to assume that the contribution due to κ is small and can be ignored

[7]. This assumption is apparently based on the fact that κ = 0 for normal processes and

that the contribution due to κ can be neglected when ξ(t) is a non-Gaussian linear process.

If the contribution of the fourth order cumulant is neglected, we obtain

cov(c(τ1), c(τ2)) = 1

T

2

T

Z−τ1

0

T

Z−τ2

0

h

E {ξ

∗

(t)ξ(s)} E {ξ(t + τ1)ξ

∗

(s + τ2)}

+ E {ξ

∗

(t)ξ

∗

(s + τ2)} E {ξ(t + τ1)ξ(s)}

i

dt ds. (26)

If T is large, the following approximate result can be used [7]

cov(c(τ1), c(τ2) ≈

1

T

Z∞

−∞

h

γ

∗

(t)γ(t + τ2 − τ1) + γ

0∗(t + τ2)γ0

(t − τ1)

i

dt

26
where γ

0

(τ ) = E {ξ(t)ξ(t + τ )}. For a complex baseband stationary process, γ

0

(τ ) = 0 and

hence

cov(c(τ1), c(τ2) ≈

1

T

Z∞

−∞

γ

∗

(t)γ(t + τ2 − τ1)dt.

This result shows that for ergodic processes, the variance of the estimator decreases as the

record length increases.

In practice, the sample function observations yield a discrete time sequence ξk. The corresponding discrete estimator of the covariance function is given by

cm =

1

N

N−Xm−1

k=0

ξ

∗

k

ξk+m m = 0, 1, · · · , N − 1.

The covariance of the estimator is approximately

cov(cm, cn) ≈

1

N

X∞

k=−∞

{γ

∗

kγk+n−m}

where now the process covariance function is discrete and is defined as a sequence

γm = E {ξ

∗

k

ξk+m} .

5.2 Estimation of the Power Spectrum

The estimate of the power spectrum is obtained by taking the Fourier transform of a sample

function. Assuming ξ(t) is a zero mean sample function we have

C(f) = 1

T













Z

T

0

ξ(t)e

−i2πf tdt













2

=

1

T

Z

T

0

Z

T

0

ξ(t)ξ

∗

(s)e

−i2πf(t−s)

dt ds.

A change of variables u = t − s and v = s yields

C(f) = Z

T

0





1

T

T

Z−u

0

ξ

∗

(v)ξ(v + u)dv



 e

−i2πfudu +

Z

0

−T





1

T

Z

T

−u

ξ

∗

(v)ξ(v + u)dv



 e

−i2πfudu.

Using the definition of c(u) from Equation 25, the first integral above is simply

Z

T

0

c(u)e

−i2πfudu.

After a change of variables, the second integral can be written as

Z

0

−T







1

T

T −|u|

Z

0

ξ

∗

(v + |u|)ξ(v)dv





 e

−i2πfudu =

Z

0

−T

c

∗

(|u|)e

−i2πfudu.

27
Recalling that c

∗

(|u|) = c(−|u|) gives

C(f) = Z

T

−T

c(u)e

−i2πfudu

which is the Fourier transform of the covariance function estimator. Thus the covariance

function estimator and the sample spectrum are Fourier transform pairs

C(f) = Z

T

−T

c(u)e

−i2πfudu c(u) = Z∞

−∞

C(f)e

i2πfudf.

The expected value of the spectral estimator of an ergodic stationary process is readily

calculated as

E {C(f)} =

Z

T

−T

E {c(u)} e

−i2πfudu =

Z

T

−T

γ(u)

 

1 −

|u|

T

!

e

−i2πfudu.

Hence, the spectral estimator is asymptotically unbiased.

Since the measured data is usually digitized we need to examine the properties of the spectral

estimator based on the discrete Fourier Transform. For equally spaced samples ξk and spacing

increment ∆t, the discrete version of the spectral estimator is

C(fn) = Cn =

∆t

N











N

X−1

k=0

ξke

−i2πnk/N











2

(27)

where fn = n/(N∆t). In the sequel, we will refer to the Cn as spectral components.

5.2.1 Statistical Characteristics of the Discrete Spectral Estimate

The expected value of the discrete estimator is readily calculated as follows

E {Cn} = ∆t

N

X−1

m=−N

 

N − |m|

N

!

γm e

−i2πnm/N

which is asymptotically unbiased.

We have seen that if the process is ergodic, the covariance function estimator converges to

the covariance function as the data record length increases. This is not true, however, for

the spectral estimator. This state of affairs can readily be demonstrated when the ξk are

real iid normally distributed samples and the spectral estimator is calculated via the discrete

Fourier transform.

28
We begin by assuming that each sample is N(0, σ2

) and consider the expected values

E {Cn} =

∆t

N

N

X−1

k=0

N

X−1

m=0

E {ξkξm} e

i2πn(k−m)/N = σ2∆t

and

E {Cn1Cn2 } =

(∆t)

2

N2

X

j,k,`,m

E {ξj ξkξ`ξm} e

i2π(n1(j−k)+n2(`−m))/N .

Since the process is iid normal, the expectation of the four samples is zero when the indices

are all different. Futhermore E {ξ

4

k} = σ

4 and E {ξj ξj ξkξk} = E

n

ξ

2

j

o

E {ξ

2

k} = σ

4 when k 6= j.

This can be written in terms of Kronecker deltas as follows

E {ξj ξkξ`ξm} = σ

4

(δjkδ`m + δj`δkm + δjmδk`)

and hence

E {Cn1Cn2 } =

σ

4

(∆t)

2

N2



N

2 +

N

X−1

j=0

N

X−1

k=0

e

i2π(j−k)(n1+n2)/N +

N

X−1

`=0

N

X−1

m=0

e

i2π(`−m)(n1−n2)/N



or

E {Cn1Cn2 } = σ

4

(∆t)

2



1 + 

sin π(n1 + n2)

N sin π(n1 + n2)/N !2

+

 

sin π(n1 − n2)

N sin π(n1 − n2)/N !2



 .

The covariance is obtained from cov(Cn1, Cn2) = E {Cn1Cn2 } − E {Cn1 } E {Cn1 } yielding

cov(Cn1, Cn2) = σ

4

(∆t)

2





 

sin π(n1 + n2)

N sin π(n1 + n2)/N !2

+

 

sin π(n1 − n2)

N sin π(n1 − n2)/N !2





and

cov(Cn1, Cn2) = (

σ

4

(∆t)

2 n1 = n2

0 n1 6= n2

.

This result demonstrates that the variance of the estimator does not decrease with the

number of samples. Furthermore, adjacent samples are uncorrelated and, for a fixed record

length, as N increases, the time between uncorrelated samples decreases. As a consequence,

the rapidity of the fluctuations of the estimator increase and the variance does not decrease.

General results for the variance have been obtained for non-Gaussian non-white stationary

processes with specific mixing assumptions (i.e., ξk and ξ` become statistically independent

as |k − `| → ∞) [8]. For example,

cov(Cn1, Cn2) = (∆tΓn1)

2

δn1n2 + O(N

−1

).

29
5.2.2 Averaged Spectral Estimates

A variety of methods have been developed to mitigate the unwieldy character of the spectral

estimator described above. In particular, a simple average of independent spectral components is an effective way to reduce the variance. This is sometimes called the Welch Method

[9]. The asymptotic statistical characteristics of such an estimator are described below.

We begin with the Fourier transform of a windowed segment of data

Xn =

N

X−1

k=0

w(k/N)ξke

−i2πkn/N .

Here, it is assumed that the window, w(t), is of bounded variations and w(t) = 0 for t < 0

and t > 1. Also, the ξk are complex samples of a stationary zero-mean noise process with

mixing properties as described in [8].

According to a central limit theorem for time series [8] the complex random variables

Xn are asymptotically independent Nc(0, Nβn) variates where βn = Γn

R

w

2

(t)dt. The

|Xn|

2 are asymptotically independent Nβnχ2

2

/2 random variables (see Equation 24) with

E {|Xn|

2} = Nβn and var(|Xn|2

) = (Nβn)

2

. Hence, we will assume that for large N, the calculated spectral components (Equation 27) can be treated as independent ∆tβnχ

2

2

/2 random

variables with E {Cn} = ∆tβn and var(Cn) = (∆tβn)

2

.

Again we see that merely increasing the length of the time series does not reduce the variance.

We can, however, make the variance small by averaging a number of the independent spectral

components. An obvious way to accomplish this is to use spectral components obtained from

each of M windowed segments of a measured time series, i.e.,

C

(`)

n =

∆t

N











N

X−1

k=0

w(k/N)ξk+(`−1)N e

−i2πnk/N











2

` = 1, . . . , M.

The asymptotic statistics of this average are easily calculated. First we note that

Cˆ(M)

n =

1

M

X

M

`=1

C

(`)

n =

∆tβn

2M

X

2M

`=1

u

2

`

where the u` are N(0, 1). The sum of the u

2

`

is a χ

2

2M variate with E {χ

2

2M} = 2M, var(χ

2

2M) =

4M and probability distribution

P{0 ≤ χ

2

2M ≤ x} = P(x| 2M) = 1

2M(M − 1)!

Zx

0

t

M−1

e

−t/2

dt.

It follows immediately that the distribution of Cˆ(M)

n

is given by

P{0 ≤ Cˆ(M)

n ≤ x} = P

 

2M

∆tβn

x











2M

!

.

30
For 2M > 30, we can use the normal approximation

P(y| n) ≈ Φ





(y/n)

1/3 − (1 − 2/9n)

q

2/9n





where

Φ(x) = 1

√

2π

Zx

−∞

e

−t

2/2

dt.

Note that the asymptotic estimate of each spectral component Cˆ(M)

n

is just an estimate of

the mean as described in Section 3. Furthermore, the samples are independent and we can

use the results presented in Section 3.3.1. The expected value of the spectral component

estimate is

E

n

Cˆ(M)

n

o

= ∆tβn = ∆tΓn

Z

w

2

(t)dt.

This estimate is biased. The contribution of the window can be removed by normalization, in

which case the expected value of the estimate is simply equal to the discrete Fourier transform

of the covariance function γ(u) = cov(ξ(t), ξ(t + u)) using N samples and a time increment

∆t. Because the time series must satisfy a mixing condition [8], the expected value of each

component can be made to approach the Fourier transform of the continuous covariance

function by increasing the number of samples N and decreasing the time increment.

The variance of the estimate is

var(Cˆ(M)

n

) = (∆tβn)

2

M

.

In practical applications, we can use Equation 8 to obtain an approximation for the variance,

i.e.,

var(Cˆ(M)

n

) ≈

1

M − 1

X

M

`=1

(C

(`)

n − Cˆ(M)n

)

2

.

Since N is large, it is a straightforward matter to calculate confidence intervals using the

normal distribution approximation (as described in Section 3.1.4).

31
6. CONCLUDING REMARKS

We have described some of the statistical uncertainties that arise when measured data are

used to characterize radio noise and interference. Mathematical expressions for uncertainties

in estimates of the mean, statistical distribution functions (viz. the amplitude probability

distribution), covariance functions and power spectral density functions are presented. In

addition to calculating uncertainties, the material presented in this report can be used to

develop measurement strategies that minimize data sampling errors.

32
7. REFERENCES

[1] H. Cramer and M. R. Leadbetter, Stationary and Related Stochastic Processes, New

York, NY: John Wiley and Sons, Inc., 1967.

[2] L. Breiman, Probability, Philadelphia, PA: Addison-Wesley Publishing Co., Inc., 1968,

p. 113.

[3] E. Parzen, Stochastic Processes, San Francisco, CA: Holden-Day, Inc., 1962, p. 93.

[4] J. Walsh, Handbook of Nonparametric Statistics, Princeton, NJ: D. Van Nostrand Company, Inc., 1962, p. 193.

[5] H. Cramer, Mathematical Methods of Statistics, Princeton, NJ: Princeton University

Press, 1954, p. 198.

[6] W. Feller, An Introduction to Probability and Its Applications, Vol 1, New York, NY:

John Wiley & Sons Inc., 1968, p. 184.

[7] G. M. Jenkins and D. G. Watts, Spectral Analysis and Its Applications, San Francisco,

CA: Holden-Day, 1968.

[8] D. R. Brillinger, Time Series Data Analysis and Theory, San Francisco, CA: Holden-Day

Inc., 1981.

[9] A. V. Oppenheim and R. W. Schafer, Digital Signal Processing, Englewood Cliffs, NJ:

Prentice-Hall, 1975.

33
 
FORM NTIA-29 U.S. DEPARTMENT OF COMMERCE 

(4-80) NAT’L. TELECOMMUNICATIONS AND INFORMATION ADMINISTRATION 

BIBLIOGRAPHIC DATA SHEET 

1. PUBLICATION NO. 

TR-09-458 

2. Government Accession No. 3. Recipient’s Accession No. 

5. Publication Date 

Nov. 2008 

4. TITLE AND SUBTITLE 

Statistical Considerations for Noise and Interference Measurements 6. Performing Organization 

NTIA/ITS.T 

7. AUTHOR(S) 

Roger Dalke 

9. Project/Task/Work Unit No. 

3105011-300 

8. PERFORMING ORGANIZATION NAME AND ADDRESS 

Institute for Telecommunication Sciences 

National Telecommunications & Information Administration 

U.S. Department of Commerce 

325 Broadway 

Boulder, CO 80305

10. Contract/Grant No. 

11. Sponsoring Organization Name and Address 

National Telecommunications & Information Administration 

Herbert C. Hoover Building 

14th & Constitution Ave., NW 

Washington, DC 20230

12. Type of Report and Period 

Covered 

14. SUPPLEMENTARY NOTES 

15. ABSTRACT 

The study of noise and interference in the radio environment is essential to the development of efficacious 

communications systems. Many of the characteristics of radio noise and interference of interest to radio 

system designers can be expressed in terms of first- and second-order statistics. These statistics are 

necessarily calculated using measured data. In this report, we discuss uncertainties that arise when 

measured data sets are used to calculate statistics of radio noise and interference. 

16. Key Words 

amplitude probability distribution; average power; measurement uncertainties; power spectrum; radio 

interference; radio noise 

18. Security Class. (This report) 

Unclassified 

20. Number of pages 

 33 

17. AVAILABILITY STATEMENT 

 F UNLIMITED. 

19. Security Class. (This page) 

Unclassified 

21. Price: 

NTIA FORMAL PUBLICATION SERIES

NTIA MONOGRAPH (MG)

A scholarly, professionally oriented publication dealing with state-of-the-art research or 

an authoritative treatment of a broad area. Expected to have long-lasting value. 

NTIA SPECIAL PUBLICATION (SP)

Conference proceedings, bibliographies, selected speeches, course and instructional 

materials, directories, and major studies mandated by Congress. 

NTIA REPORT (TR)

Important contributions to existing knowledge of less breadth than a monograph, such as 

results of completed projects and major activities. Subsets of this series include: 

NTIA RESTRICTED REPORT (RR)

Contributions that are limited in distribution because of national security 

classification or Departmental constraints. 

 NTIA CONTRACTOR REPORT (CR)

Information generated under an NTIA contract or grant, written by the contractor, 

and considered an important contribution to existing knowledge. 

 JOINT NTIA/OTHER-AGENCY REPORT (JR)

This report receives both local NTIA and other agency review. Both agencies’ 

logos and report series numbering appear on the cover. 

NTIA SOFTWARE & DATA PRODUCTS (SD)

Software such as programs, test data, and sound/video files. This series can be used to 

transfer technology to U.S. industry. 

NTIA HANDBOOK (HB)

Information pertaining to technical procedures, reference and data guides, and formal 

user's manuals that are expected to be pertinent for a long time. 

NTIA TECHNICAL MEMORANDUM (TM)

Technical information typically of less breadth than an NTIA Report. The series includes 

data, preliminary project results, and information for a specific, limited audience. 

For information about NTIA publications, contact the NTIA/ITS Technical Publications Office at 

325 Broadway, Boulder, CO, 80305 Tel. (303) 497-3572 or e-mail info@its.bldrdoc.gov. 

This report is for sale by the National Technical Information Service, 5285 Port Royal Road, 

Springfield, VA 22161,Tel. (800) 553-6847.