PC VI: Signal processing in spectroscopy

G. Jeschke
ETH Zürich, Laboratorium für Physikalische Chemie

[email protected]
HCI F227, Tel. 25702

December 18, 2012

Contents

1 What is spectroscopy? 9
1.1 Spectroscopy as characterization of quantum systems . . . . . 9
1.1.1 A definition of spectroscopy . . . . . . . . . . . . . . . 9
1.1.2 Description of the quantum system . . . . . . . . . . . 10
1.1.3 Description of the sample . . . . . . . . . . . . . . . . 11
1.1.4 Interaction of the sample with the environment . . . . 12
1.1.5 Interaction of the spectrometer with the sample . . . . 14

2 Linear response theory 17

2.1 Systems and their properties . . . . . . . . . . . . . . . . . . 17
2.1.1 System definition . . . . . . . . . . . . . . . . . . . . . 17
2.1.2 Linearity . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.1.3 Time invariance . . . . . . . . . . . . . . . . . . . . . 18
2.1.4 Causality . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.5 Memory . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.6 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2 Linear time-invariant systems . . . . . . . . . . . . . . . . . . 19
2.2.1 Impulse response . . . . . . . . . . . . . . . . . . . . . 20
2.2.2 Properties of LTI systems . . . . . . . . . . . . . . . . 21
2.2.3 Step response . . . . . . . . . . . . . . . . . . . . . . . 22
2.2.4 Frequency-response function . . . . . . . . . . . . . . . 22

3 Fourier transform 24
3.1 Integral transforms . . . . . . . . . . . . . . . . . . . . . . . . 24
3.1.1 Fourier series . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.2 Orthonormal systems of functions . . . . . . . . . . . 27
3.2 Definition and properties of the Fourier transform . . . . . . 28
3.2.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . 28

1
 3.2.2 Basic properties . . . . . . . . . . . . . . . . . . . . . 28
3.2.3 Convolution and deconvolution . . . . . . . . . . . . . 31
3.2.4 Correspondences . . . . . . . . . . . . . . . . . . . . . 31
3.3 Representations of complex spectra . . . . . . . . . . . . . . . 34
3.3.1 Absorption and dispersion spectrum . . . . . . . . . . 34
3.3.2 Magnitude spectrum . . . . . . . . . . . . . . . . . . . 35
3.3.3 Power spectrum . . . . . . . . . . . . . . . . . . . . . 36

4 Basics of electronics 37
4.0.4 Energy, power, amplitude, impedance . . . . . . . . . 37
4.0.5 The dB scale . . . . . . . . . . . . . . . . . . . . . . . 39
4.1 Circuit elements and rules . . . . . . . . . . . . . . . . . . . . 40
4.1.1 Circuit elements . . . . . . . . . . . . . . . . . . . . . 40
4.1.2 Rules for combining elements to a circuit . . . . . . . 43
4.2 Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2.1 RC low pass . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2.2 Ideal and real filters . . . . . . . . . . . . . . . . . . . 48
4.2.3 Butterworth filters . . . . . . . . . . . . . . . . . . . . 51
4.2.4 Chebyshev filters . . . . . . . . . . . . . . . . . . . . . 51

5 Noise 53
5.1 Noise sources and types . . . . . . . . . . . . . . . . . . . . . 53
5.1.1 Noise sources and signal-to-noise ratio . . . . . . . . . 53
5.1.2 Thermal noise . . . . . . . . . . . . . . . . . . . . . . 55
5.1.3 Shot noise . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.1.4 1/f noise and Brownian noise . . . . . . . . . . . . . . 57
5.1.5 Phase noise . . . . . . . . . . . . . . . . . . . . . . . . 58
5.2 Mathematics of noise . . . . . . . . . . . . . . . . . . . . . . . 59
5.2.1 Noise spectral density and noise figure . . . . . . . . . 59
5.2.2 Cascaded systems . . . . . . . . . . . . . . . . . . . . 62
5.2.3 Techniques for improving the signal-to-noise ratio . . . 63

6 Effect modulation and phase-sensitive detection 66

6.1 Effect Modulation . . . . . . . . . . . . . . . . . . . . . . . . 66
6.1.1 General principle of modulation . . . . . . . . . . . . . 66
6.1.2 Types of modulation and distortions . . . . . . . . . . 67
6.2 Phase-sensitive detector/ Lock-in amplifier . . . . . . . . . . . 70
6.2.1 Digital PSD . . . . . . . . . . . . . . . . . . . . . . . . 72

2
7 Time-discrete signals 73
7.1 Analog-digital conversion . . . . . . . . . . . . . . . . . . . . 73
7.1.1 Digitization . . . . . . . . . . . . . . . . . . . . . . . . 74
7.1.2 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.1.3 Aliasing . . . . . . . . . . . . . . . . . . . . . . . . . . 77
7.1.4 Oversampling . . . . . . . . . . . . . . . . . . . . . . . 78
7.2 Fourier Transformation of time-discrete signals . . . . . . . . 79
7.2.1 Computation of a continuous function in frequency
domain . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.2.2 Properties of discontinuous FT . . . . . . . . . . . . . 80
7.2.3 Discrete FT . . . . . . . . . . . . . . . . . . . . . . . . 81
7.2.4 Zero filling . . . . . . . . . . . . . . . . . . . . . . . . 83
7.2.5 Fast Fourier transformation . . . . . . . . . . . . . . . 84
7.3 Apodization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.3.1 Effects of limited observation time . . . . . . . . . . . 85
7.3.2 Apodization window . . . . . . . . . . . . . . . . . . . 87
7.3.3 Filtering for optimum signal-to-noise . . . . . . . . . . 88
7.3.4 Filtering for optimum resolution . . . . . . . . . . . . 91
7.4 Network analyzers and spectrum analyzers . . . . . . . . . . . 92

8 Stochastic signals 95
8.1 Characteristics of random variables . . . . . . . . . . . . . . . 95
8.1.1 Probability distribution . . . . . . . . . . . . . . . . . 96
8.1.2 Probability density . . . . . . . . . . . . . . . . . . . . 97
8.1.3 Mean value, variance, standard deviation, moment
analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 97
8.1.4 Gaussian distribution . . . . . . . . . . . . . . . . . . 98
8.1.5 Poisson distribution . . . . . . . . . . . . . . . . . . . 99
8.1.6 Two-dimensional autocorrelation function . . . . . . . 100
8.2 Correlation functions of stationary processes . . . . . . . . . . 101
8.2.1 Stationary processes . . . . . . . . . . . . . . . . . . . 101
8.2.2 Ergodicity theorem . . . . . . . . . . . . . . . . . . . . 101
8.2.3 Autocorrelation function . . . . . . . . . . . . . . . . . 102
8.2.4 Cross correlation functions . . . . . . . . . . . . . . . 103
8.2.5 Correlation functions of periodic signals . . . . . . . . 104
8.2.6 Felgett’s advantage . . . . . . . . . . . . . . . . . . . . 107

3
9 Alternatives to one-dimensional Fourier transformation 108
9.1 Spectral estimators and parameter estimators . . . . . . . . . 108
9.1.1 Asssumptions on the observed system . . . . . . . . . 108
9.1.2 2D spectroscopy . . . . . . . . . . . . . . . . . . . . . 109
9.1.3 Deficiencies of spectrum estimators . . . . . . . . . . . 112
9.1.4 Parameter estimators . . . . . . . . . . . . . . . . . . 113
9.2 Harmonic inversion . . . . . . . . . . . . . . . . . . . . . . . . 114
9.2.1 Linear prediction . . . . . . . . . . . . . . . . . . . . . 115
9.2.2 Singular value decomposition . . . . . . . . . . . . . . 117
9.3 Fredholm integral equations . . . . . . . . . . . . . . . . . . . 119
9.3.1 Fredholm integral equations of the first kind . . . . . . 120
9.3.2 Well-posed and ill-posed problems . . . . . . . . . . . 121
9.4 Regularization . . . . . . . . . . . . . . . . . . . . . . . . . . 122
9.4.1 Tikhonov regularization . . . . . . . . . . . . . . . . . 122

10 Digital signal processing 125

10.1 Real-time digital signal processing . . . . . . . . . . . . . . . 125
10.1.1 Digital signal processors . . . . . . . . . . . . . . . . . 126
10.1.2 Field-programmable gate arrays . . . . . . . . . . . . 126
10.2 General properties of time-discrete systems . . . . . . . . . . 127
10.2.1 Properties and description . . . . . . . . . . . . . . . . 127
10.2.2 FIR and IIR systems . . . . . . . . . . . . . . . . . . . 128
10.2.3 Causality . . . . . . . . . . . . . . . . . . . . . . . . . 128
10.2.4 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . 128
10.2.5 Frequency response . . . . . . . . . . . . . . . . . . . . 129
10.2.6 Difference equations . . . . . . . . . . . . . . . . . . . 129
10.3 z transformation . . . . . . . . . . . . . . . . . . . . . . . . . 130
10.3.1 Existence of the z transform . . . . . . . . . . . . . . . 131
10.3.2 Relation between z transformation and FT . . . . . . 131
10.3.3 Properties of the z transformation . . . . . . . . . . . 132
10.4 Transfer function . . . . . . . . . . . . . . . . . . . . . . . . . 133
10.4.1 Poles and zeros . . . . . . . . . . . . . . . . . . . . . . 134
10.4.2 Systems in series and in parallel . . . . . . . . . . . . 135
10.4.3 Frequency-response function, phase response, and group
delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
10.4.4 Relation between pole-zero plot and response functions 136
10.4.5 Minimal phase system and allpass . . . . . . . . . . . 140

4
 10.4.6 Transfer function of FIR systems . . . . . . . . . . . . 142
10.4.7 Transfer function of IIR systems . . . . . . . . . . . . 142
10.4.8 Linear phase systems . . . . . . . . . . . . . . . . . . . 142
10.5 Matlab tools for digital signal processing . . . . . . . . . . . . 145
10.5.1 Useful Matlab functions for system analysis and design 145
10.5.2 Filter design with Matlab . . . . . . . . . . . . . . . . 146

5
Introduction
Scope of this lecture course
This lecture is concerned with spectroscopic measurements and basic pro-
cessing of signals. Only those aspects are treated which are common to
different types of spectroscopy. In fact, many of the topics are also relevant
for a broad range of other measurements in physics and physical chemistry.
The lecture course discusses the procedures that are required to obtain the
most precise, most reliable and least noisy data on a system under study.
Such a general approach implies that the subject matter is rather abstract.
To relate the abstract concepts to spectroscopy, Chapter 1 is devoted to a
general discussion of spectroscopy. This Chapter clarifies what we want to
measure and how the system under study generally responds to excitation.
In Chapter 2 we go to an even higher level of abstraction by treating the
system under study as a black box. This black box receives some inputs
and responds with some outputs. We are interested in only the transfer
function that converts the inputs into outputs, not with the underlying
physics. We ask the question how we can find this function from pairs of
input (excitation) and output (signal). To simplify the treatment we assume
that the system responds linearly.
Both the general behavior of quantum ssytems during spectroscopic experi-
ments and linear response theory call for mathematical methods that relate
time and frequency domain. Such methods are integral transforms, in par-
ticular, the Fourier transform. In Chapter 3 we discuss general properties
of the Fourier transform, which are important for signal analysis.
Signals are detected and partially modified by measurement electronics.
Chapter 4 treats the most important basics of electronics, in particular, how
electronics modifies systems in time and frequency domain. The concept of
signal filtering is introduced. This sets the stage for consideration of noise in
Chapter 5. We consider how noise can be characterized, how it is averaged,
and how it can be reduced or partially supressed. Chapter 6 discusses one
particularly important method for noise reduction, effect modulation with
phase-sensitive detection.
Modern spectroscopy is based on digitization and discretization of signals.
These concepts are treated in Chapter 7, where analog-digital conversion,
sampling issues, and the discrete Fourier transformation are discussed. Fi-
nite observation time of the response may require apodization of the signal

6
by multiplication with a window function, which is introduced as well.
Spectroscopic signals may be related to stochastic processes and are super-
imposed by noise that originates from stochastic processes. To allow for
a quantitative treatment of such situations, Chapter 8 introduces tools for
characterization of stochastic signals, in particular, parameters of random
variables and correlation functions.
The relation between a spectroscopic signal and the underlying physical
information is not necessarily well described by a Fourier transformation. In
some important situations the physical data are related to the spectroscopic
data to a Fredholm integral equation, which may pose a mathematically
ill-posed problem. Such problems can be solved by Tikhonov regularization,
which is introduced in Chapter 9.
Analog filtering and phase-sensitive detection for noise reduction as dis-
cussed in Chapters 4 and 6 are increasingly replaced by equivalent processing
of digital data with digital signal processors. This is discussed in Chapter
10.

Literature
The following textbooks in German language provide good introductions of
the mathematical and theoretical concepts:

[Meyer] Martin Meyer, Signalverarbeitung: Analoge und digitale Signale,

Systeme und Filter, Vieweg + Teubner.

[Werner] Martin Werner, Digitale Signalverarbeitung mit MATLAB: Grund-

kurs mit 16 ausführlichen Versuchen, Vieweg + Teubner.

[Leon et al.] Fernando Puente Leon, Uwe Kiencke, Holger Jäkel, Signale
und Systeme, Oldenbourg.

[Kiencke & Eger] Uwe Kiencke, Ralf Eger, Messtechnik, Springer.

Among those, [Meyer] and [Werner] are most easily accessible, whereas [Leon
et al.] is most mathematically strict. The book by [Kiencke & Eger] is geared
at electrical engineers and contains some material that is treated in basic
science courses. On the other hand it also contains some information on
optimization of signal processing electronics, which is missing in the more
abstract textbooks.

7
Many, but not all important concepts are treated in English language in

[Bracewell] Ronald N. Bracewell, The Fourier Transform and its Application,

McGraw-Hill.

Non-linear system theory is not discussed in this lecture course. A good

introduction of this more advanced theory and some other advanced aspects
is given in:

[Blümich] B. Blümich, Progr. in NMR Spectrosc, 19, 331-417 (1987).

A rather comprehensive and very good introduction to digital signal pro-
cessing can be found at the homepage of Analog Devices, Inc. (sic!):

http://www.analog.com/en/processors-dsp/learning-and-development/
content/scientist engineers guide/fca.html

Acknowledgement
Helpful discussions with and corrections from Dr. Maxim Yulikov, Dr.
Yevhen Polyhach, Andrin Doll, and Rene Tschaggelar are gratefully ac-
knowledged.

8
Chapter 1

What is spectroscopy?

1.1 Spectroscopy as characterization of quantum

systems
1.1.1 A definition of spectroscopy
To talk about measurement technique in different spectroscopies on a general
level, we need an abstract definition of spectroscopy. Here it is:

• Spectroscopy is the characterization or manipulation of a

multi-level quantum system by its interaction with electro-
magnetic fields.

The multi-level quantum system is a small part of the world. The large
rest of the world is the environment of this system. A special part of this
environment is our measurement device- the spectrometer- which we can
manipulate at will and which we can use in turn to manipulate the quan-
tum system. Interaction of the measurement device and of other parts of
the environment with a single instance of the quantum system is solely by
electromagnetic fields, which can be static or oscillatory. In most cases we
work with ensembles of quantum systems and then the temperature, which
we can also manipulate, is an important parameter. Temperature generally
influences the thermodynamic equilibrium state of an ensemble of quantum
systems as well as the interaction of other parts of the environment with
this ensemble.
In this lecture course, the ensemble of quantum systems is termed sample.
The term environment excludes the spectrometer. Our world thus consists

9
 heat
Spectrometer

electro-
magnetic fields Environment
at temperature
T
electro-
Sample magnetic
(ensemble of
quantum systems) fields

Figure 1.1: Abstract view of the world of spectroscopy.

of three parts- sample, spectrometer, and environment. Interactions be-

tween spectrometer and environment are of no concern to us, with the sole
exception of temperature control (Figure 1.1). We first consider description
of the sample, as this will reveal what we can measure and how we can
measure it.

1.1.2 Description of the quantum system

A single instance of a quantum system is fully described by the Hamiltonian
b 1 which determines behavior of the system, and the wavefunction Ψ, which
H,
denotes the state of the system. The eigenvalues i and eigenfunctions ψi
of the system can be obtained by solving the time-independent Schrödinger
equation,
HΨ
b = EΨ . (1.1)

The eigenvalues i are the energies of the levels. Each level corresponds to
an eigenstate |ii = ψi .
We are interested in the behavior and state of the unperturbed system in a
mean, time-averaged state of the environment. The corresponding Hamil-
tonian is the static Hamiltonian Hb0 with eigenvalues 0,i . We can infer
energy differences between eigenstates by exciting and observing transitions
between eigenstates. The energy of transition (i, j) is

∆0,ij = ~ωij = hνij , (1.2)

1
In this lecture course, Hamiltonians are given in energy units

10
where ωij is the angular frequency of the transition, νij the frequency of the
h
transition, h Planck’s quantum of action, and ~ = 2π .
In general, the state of the system is time dependent, Ψ = Ψ(t). Only
eigenstates of H
b are time invariant. The evolution of the system is described
by the time-dependent Schrödinger equation

d i b
Ψ = − HΨ . (1.3)
dt ~

For the time-independent Hamiltonian H

b0 , Eq. (1.3) can be solved analyti-
cally. In particular, we are interested in evolution of a coherent superposition
of eigenstates |ii and |ji, as it is encountered during excitation of transition
(i, j). Such a coherence evolves according to

ψij (t) = exp (−iωij t) ψij (0) . (1.4)

Hence, the amplitude of a coherence between eigenstates is a constant of

motion, whereas the phase is proportional to ωij t.

1.1.3 Description of the sample

This description would suffice if all observed quantum systems in the sam-
ple were in the same (time-dependent) state throughout the spectroscopic
experiment. Such a pure state of an ensemble of quantum systems is an
exception. In general, the system is in a mixed state, which is described by
a density operator
N
X
ρb = pk |ψk ihψk | , (1.5)
k

where the ψk are wavefunctions of N subensembles in pure states2 and pk

is the probability for finding a system in subensemble k. The equation of
motion for this density operator ρb, disregarding interaction with the envi-
ronment, is the Liouville-von-Neumann equation

d i hb i
ρb(t) = − H(t), ρb(t) . (1.6)
dt ~

In this density operator formalism, coherence on transition (i, j) corresponds

to the complex element ρij of the matrix representation of the density op-
erator in the eigenbasis of H
b0 . The coherence amplitude is given by |ρij |
2
N can be as large as the number of quantum systems in the ensemble

11
and the phase by φij = arctan [< (ρij ) /= (ρij )]. In analogy to Eq. (1.4),
evolution of the coherence conforms to

ρij (t) = exp (−iωij t) ρij (0) . (1.7)

1.1.4 Interaction of the sample with the environment

Eqs. (1.7) applies to an ensemble of isolated quantum systems. Interaction
between the sample and the environment can be formally described in Eq.
(1.6) by considering an appropriate time-dependent Hamiltonian that arises
from time-dependent electromagnetic fields due to fluctuating electric and
magnetic dipoles and multipoles in the environment. In such a description,
H
b0 includes the time-averaged couplings within the quantum system and the
time-averaged couplings of the quantum system with external fields. The
coupling Hamiltonian H
bc quantifies fluctuation of the couplings around their
mean values and the total Hamiltonian with the spectrometer still switched
off is given by H
b=Hb0 + Hbc .
Due to the large number of degrees of freedom of the environment the fluc-
tuations are stochastic. Phenomenologically, they cause relaxation of the
ensemble towards a thermal equilibrium state ρbeq . Maximization of entropy
in the equilibrium state causes a vanishing of coherence, whereas populations
conform to a Boltzmann distribution,

exp (−i /(kB T ))

ρeq,ii = P , (1.8)
i exp (−i /(kB T ))

where kB is the Boltzmann constant. The thermal equilibrium state can be

computed by the basis-independent operator equation
 
exp −H/(k
b BT )
ρbeq = n  o , (1.9)
Tr exp −H/(k
b BT )

where Tr denotes the trace of the matrix representation.

In many situations, approach to this equilibrium state can be approximated
by first order kinetics. The polarization Pij = ρii − ρjj relaxes according to

d 1
Pij (t) = − (Pij (t) − Peq,ij ) , (1.10)
dt T1,ij

12
whereas coherence relaxes according to

d 1
ρij (t) = − ρij (t) . (1.11)
dt T2,ij

We shall see later that T2,ij determines the homogeneous linewidth. Relax-
ation of populations sets a lower limit to the rate constant 1/T2,ij , which is
given by
1 1
≥ . (1.12)
T2,ij 2T1,ij
This effect is termed lifetime broadening of the spectral line. Lifetime broad-
ening is related to Heisenberg’s uncertainty relation for time and energy
∆ (i − j ) ∆t % h.
The rate constant 1/T1,ij contains contributions from coupling of the radia-
tion field to the quantum system and from energy transfer between degrees
of freedom of the environment and the quantum system. The former contri-
bution, which includes spontaneous and stimulated emission dominates in
optical spectroscopy. In magnetic resonance spectroscopy, where frequencies
are much smaller and the strength of coupling to an incoherent radiation
field is exceedingly small, the latter contribution dominates.
An additional contribution to 1/T2,ij arises from fluctuating fields that lead
to stochastic fluctuation of ωij . This in turn leads to stochastic phase fluctu-
ations of subensemble wavefunctions, see Eq. (1.4), and thus to destructive
interference. Consequently, homogeneous linewidths may be larger than
would be expected from lifetime broadening alone.
Interaction of the ensemble with the environment may also lead to transfer
of polarization or coherence between transitions (i, j) and (k, l). Such cross
relaxation is negligible between coherences, unless the two transitions have
the same frequency. If the frequency differs, the coherence is lost due to the
stochastic nature of the transfer.
Within the approximations made in this section, coupling between the sam-
ˆ,
ple and the environment can be described by a relaxation superoperator Γ̂
with element Γij,kl describing relaxation between elements ρij and ρkl of the
density operator. This description by a stochastic Liouville equation,

d i hb i
ˆ {b
ρb(t) = − H b(t) − Γ̂
0, ρ ρ(t) − ρbeq } , (1.13)
d ~

eliminates explicit dependence on H

bc (t).

13
Since only ρb is time dependent in Eq. (1.13), this equation can be for-
mally solved. The solution shows that coherence on transition (i, j) evolves
according to    
1
ρij (t) = exp − iωij + t ρij (0) . (1.14)
T2,ij
Hence, ωij and T2,ij can be inferred from the time evolution of a signal that
is proportional to coherence on transition (i, j). This is the approach of
time-domain spectroscopy. Alternatively, the spectrum can be measured by
scanning an excitation frequency νexci in the range of expected transition
frequencies and observing absorption or emission. In such frequency-domain
spectroscopy the spectral lines are centered at νij = ωij /(2π) and their width
is related to 1/T2,ij . The information obtained by time- and frequency-
domain spectroscopy is thus equivalent.

1.1.5 Interaction of the spectrometer with the sample

When the sample is put into the spectrometer, it is in a thermal equilibrium
state described by ρbeq .3 This state is devoid of coherence and in the ensemble
average as well as in the time average there is no energy transfer between
sample and environment, here including the spectrometer, so that no signal
is measured.
The sample can be driven away from thermal equilibrium by resonant ab-
sorption of photons. Due to the coherence decay with rate constant 1/T2,ij
there is some energy uncertainty so that the photon frequency does not
need to be exactly resonant with the transition. Absorption of photons re-
quires that the transition between states |ii and |ji couples to the electric or
magnetic component of the electromagnetic field. In coherent spectroscopy,
irradiation and coupling can be described by a Hamiltonian
h i
H
b1 (t) = F1 c cos (2πνt) C
b + sin (2πνt) Sb , (1.15)

where F1 is the amplitude of the circularly polarized radiation field, c is

a coupling factor, ν the frequency of the electromagnetic radiation, and C
b
and Sb are operators that describe how the field couples to the system. The
product ω1 = F1 c has the dimension of an angular frequency. For instance, in
EPR spectroscopy F1 = B1 is the amplitude of the oscillatory magnetic field,
3
In magnetic resonance spectroscopy, the T1,ij can reach hours or days, depending on
the observed spin species and temperature. In such cases, the sample reaches thermal
equilibrium in the spectrometer only on the timescale of the longest T1,ij .

14
c = gµB /~ depends on the electron spin g factor and the Bohr magneton
µB , and C
b = Sbx as well as Sb = Sby are spin operators, assuming that the
static field is along z.
The excitation Hamiltonian H
b1 (t) is given in a laboratory-fixed frame, which
relates to the direction of field polarization. The system is not necessarily
quantized along the z axis of this frame. To see, which transitions are driven
by the electromagnetic field, H b1 (t) must be transformed into the eigen-
frame of the static Hamiltonian Hb0 . This transformation can be preceded
by transformation into a rotating or interaction frame, where H
b1 becomes
time-independent. The amplitude of the off-diagonal element Hij of the ma-
trix representation of H
b1 in the eigenframe of H
b0 is the transition amplitude
aij for transition (i, j), D E
aij = i|H
b1 |j . (1.16)

Observation of the signal is symmetric to excitation as it requires the same

type of coupling between the quantum system and the radiation field. If the
coupling to the excitation arm and detection arm of the spectrometer is the
same, the signal (line intensity) is proportional to the transition probability

wij = |aij |2 . (1.17)

The proportionality factor is an apparative constant, which together with

the noise level determines sensitivity of the spectrometer. The concept of
transition probability was introduced in a description of incoherent spec-
troscopy, where it arises from second-order time-dependent perturbation
theory (Fermi’s golden rule). It may fail in coherent spectroscopy, if exci-
tation and detection arms correspond to different polarization directions of
the field, i.e., if C
b and Sb are not the same for excitation and detection. This
problem does not, however, arise with conventional spectrometers.
During irradiation the system is characterized by the Hamiltonian H
b =
H
b0 + Hb1 . Consequently, we do not observe an unperturbed system, but
the system perturbed by H
b1 . The perturbation becomes significant when
the transition probability according to Eq. (1.17) approaches or exceeds
1/ (T1,ij T2,ij ). In this situation, relaxation cannot restore equilibrium po-
larization so that the transition is (partially) saturated. The extent of sat-
uration is maximum on resonance and decreases with increasing resonance
offset of the radiation field. Accordingly, the loss of signal intensity com-
pared to a situation without saturation is more severe on resonance than off

15
resonance, which results in line broadening. Such power broadening distorts
line shapes and may spoil resolution.
Another distortion arises in double resonance experiments when coherence
on transition (i, j) is observed while transition (k, l) is irradiated. If the
irradiation frequency is far off resonance, no transitions (i, j) are driven.
However, a Bloch-Siegert shift,

ωij
∆ωBS = ω12 2 2 , (1.18)
ωij − ωkl

arises. Accordingly, irradiation of transition (k, l) for a time texci causes a

phase shift of coherence on transition (i, j) by ∆φij = ∆ωBS texci .
Line broadening can also arise from strong coupling of coherences to the de-
tection arm of the spectrometer. This effect is most easily understood for an
NMR experiment, where coherence ρij on an allowed transition corresponds
to transverse magnetization. This transverse magnetization oscillates with
angular frequency ωij and induces a current in the detection coil. This
current in the coil in turn creates a magnetic field that also oscillates with
ωij and is thus resonant with the transition. The oscillatory field transfers
coherence ρij to polarization and thus destroys it. The resulting decay of
the coherence is termed radiation damping and is relevant only for strong
signals of narrow lines, as for instance the water signal in high-resolution
NMR.
In the following lectures we assume that power broadening and radiation
damping are avoided, so that coherences evolve according to Eq. (1.14) and
signal intensities in continuous irradiation experiments are proportional to
the transition probability wij given by Eq. (1.17).

16
Chapter 2

Linear response theory

2.1 Systems and their properties

2.1.1 System definition
In the following, a system is a device that responds to an input signal x(t)
with an output signal y(t), where t is the time (Figure 2.1). We can think
about x(t) as excitation of the sample in a spectroscopic experiment and
about y(t) as response to this excitation. However, the following theory is
much more general than that.

x(t) System y(t)

Figure 2.1: Signal processing by a system.

Mathematically, we can formulate this definition with an operator S that is

characteristic for the behavior of the system,

y(t) = S {x(t)} (2.1)

In spectroscopy, we seek information on S. We first consider time-continuous

systems, where both x(t) and y(t) are continuous signals (although they may
be zero at certain times).
All properties discussed in the following extend to systems with an input
vector x(t) and an output vector y(t).

17
2.1.2 Linearity
A simple theoretical treatment is possible if the system is linear. Linearity
is defined by the operator equation

S {c1 x1 (t) + c2 x2 (t)} = c1 S {x1 (t)} + c2 S {x2 (t)} , (2.2)

where c1 and c2 are constants. Linear systems behave according to the

superposition principle,
(N ) N
X X
S ci xi (t) = ci S {xi (t)} , (2.3)
i=1 n=1

i.e., they are additive.

Going from the sum to an integral
Z  Z
S c(τ )x(t, τ )dτ = c(τ )S {x(t, τ )} dτ (2.4)

we find as a remarkable consequence that the response to the integral of the

input signal equals the integral of the response to the input signal.
Systems in spectroscopy are not in general linear systems. However, for
low excitation power and short excitation times the systems are in a linear
regime that can be treated by linear response theory. Insights from this
treatment can often be extended beyond the linear regime by discussing
deviations, such as saturation, separately.

2.1.3 Time invariance

The next important property is time invariance, which is defined by the
operator equation

y(t) = S {x(t)} ⇒ y(t − t0 ) = S {x(t − t0 )} . (2.5)

In physical theories, time invariance is related to energy conservation. In

spectroscopy, it can be achieved by preparing and maintaining the system
under defined external conditions.

18
2.1.4 Causality
Furthermore, we characterize systems by causality. A system is causal if
and only if its response depends only on the current value and past values
of the input signal, but not on future values of the input signal. Hence, if
x1 (t) = x2 (t) for all times t ≤ t1 , it follows that y1 (t) = y2 (t) for all times
t ≤ t1 . Causality is a general property of systems that we study in science.1

2.1.5 Memory
Systems can be dynamic. The response of a dynamic system depends not
only on the current value of the input signal, but also on past values. In
other words, a dynamic system has a memory. Except for a few cases of
limiting behavior, spectroscopic systems are dynamic.

2.1.6 Stability
Systems can further be stable or instable. A stable systems responds with a
bounded output signal to a bounded input signal. In other words, if there
is a signal amplitude m < ∞ so that |x(t)| < m at all times t, there is also
a signal amplitude M < ∞ so that |y(t)| < M at all times t. An integrator
is an example for an unstable system.

2.2 Linear time-invariant systems

Linear time-invariant (LTI) systems are particulary simple to treat.

Box 1: Unit impulse function

In a Gedankenexperiment, assume that we can excite the system with an in-
finitely brief, infinitely strong unit-area impulse
δ(t) = 0 for t 6= 0
R∞
−∞
δ(t)dt = 1.
In a mathematical sense, this is an improper function, as pointed out by Dirac.
We can define it as a limiting case of a series of proper functions, for instance
2 2
δ(t) = lim 1
√ e−t /a .
a→0 a π

1
At least we assume that and science appears to work well based on this assumption

19
2.2.1 Impulse response
The response g(t) of the system S to a Dirac delta function δ(t) (unit impulse
function, see Box 2.2) is the impulse response. Any input signal can be
expressed as Z ∞
x(t) = x(t − τ )δ(τ )dτ . (2.6)
−∞

Using time invariance, Eq. (2.5), and the superposition principle, Eq. (2.3),
we find
Z ∞ 
y(t) = S {x(t)} = S x(t − τ )δ(τ )dτ
−∞
Z ∞  Z ∞
= S x(µ)δ(t − µ)dµ = x(τ )S {(δ(t − τ )} dτ
Z ∞ −∞ −∞

= x(τ )g(t − τ )dτ = x(t) ∗ g(t) , (2.7)

−∞

where the asterisk operator ∗ denotes the convolution integral. In the first
step we substituted t − τ → µ, which implies the substitution τ → t − µ.2 In
the next step we renamed µ as τ and used the superposition principle. The
results shows that the response to any input signal x(t) is the convolution of
this input signal with the impulse response of the system. In other words,
an LTI system is fully characterized by its impulse response. The meaning
of convolution can be grasped from the graphic representation in Figure 2.2.

* =

Figure 2.2: Convolution of a stick spectrum with a Gaussian line shape function.

The convolution integral is commutative,

x(t) ∗ g(t) = g(t) ∗ x(t) , (2.8)

2
R ∞substitution also proves that the convolution of g(t) and f (t), defined as g(t) ∗
This
f (t) = −∞ g(t − τ )f (τ )dτ , is commutative.

20
which follows from the substitution t0 = t − τ . It is also associative,

x(t) ∗ [g1 (t) ∗ g2 (t)] = [x(t) ∗ g1 (t)] ∗ g2 (t) . (2.9)

Because of associativity, the signal processors shown in Figure 2.3a and b

are equivalent. If a signal consecutively passes several systems, conversion
from the input to the final output signal can be described by the consecutive
convolution of the impulse responses of all systems.
The convolution integral is also distributive

x(t) ∗ [g1 (t) + g2 (t)] = x(t) ∗ g1 (t) + x(t)h2 (t) . (2.10)

Hence, if a signal is split, the two parts pass two different systems, and are
then recombined, the impulse response is the sum of the impulse responses
of the two systems.

a c h1

x(t) h1 h2 y(t) x(t) Å y(t)

h2

b d
x(t) h1 * h2 y(t) x(t) h1 + h2 y(t)

Figure 2.3: Associative and dissociative properties of the convolution integral.

The two signal processors (a) and (b) as well as the two processors (c) and (d) are
equivalent.

There is a simple relation between the derivatives of the output signal, the
input signal, and the impulse response

d d d
y(t) = x(t) ∗ g(t) = x(t) ∗ g(t) . (2.11)
dt dt dt

2.2.2 Properties of LTI systems

An LTI system is causal if and only if its impulse response vanishes at t < 0
(g(t) = 0 for t < 0). It is stable if and only if the impulse response fulfils
the condition Z ∞
|g(t)| dt < ∞ . (2.12)
−∞

21
Proofs for both propositions can be found in [Leon et al.].

2.2.3 Step response

Another important test function is the step function (Heaviside function),
(
0 for t < 0
u(t) = . (2.13)
1 for t ≥ 0

The step function, which is experimentally easier to realize, can be related

to the Dirac delta function by
Z t
u(t) = δ(τ )dτ . (2.14)
−∞

Hence, the step response is given by

Z t
a(t) = g(τ )dτ , (2.15)
−∞

which implies
d
g(t) = a(t) . (2.16)
dt
The impulse response of an LTI system is thus completely determined by the
step response. Hence, the LTI is fully characterized by the step response.
There exists a simple relation between the output signal and the derivatives
of the input signal and of the step response:

d d
y(t) = x(t) ∗ a(t) = x(t) ∗ a(t) . (2.17)
dt dt

Important characteristics of the step response are the rise time trise , the
overshoot, ringing, and the settling time. These characteristics also play a
role in pulse spectroscopies. They are illustrated in Figure 2.4.

2.2.4 Frequency-response function

Another interesting property for spectroscopic applications is the frequency-
response function.3 Any LTI system that is excited by a complex oscillation
at frequency ν,
x(t) = aei2πνt = aeiωt , (2.18)
3
Deutsch: Frequenzgang, manchmal auch Übertragungsfunktion, obwohl diese streng
genommen die Laplace-transformierte der Impulsantwort ist.

22
 overshoot

90%

10%

rise time trise

Figure 2.4: Typical step response function and main characteristics. Modified
from Wikimedia Commons.

responds with an output oscillation at the same frequency ν. The func-

tion G(ν) that characterizes the amplitude ratio between output and input
oscillation is the requency-response characteristic of system S. It can be
computed as the Fourier transform (Chapter 3) of the impulse response
g(t):
Z ∞
y(t) = x(t) ∗ g(t) = x(t − τ )g(τ )dτ
−∞
Z ∞
= aeiω(t−τ ) g(τ )dτ
−∞
Z ∞
= ae +iωt
e−iωτ g(τ )dτ . (2.19)
−∞

The frequency-response function is generally a complex-valued function,

Z ∞
G(ω) = e−iωτ g(τ )dτ = R(ω) + iI(ω) = A(ω)e−iθ(ω) , (2.20)
−∞

with a real part R(ω), an imaginary part I(ω), the amplitude spectrum
A(ω), and the phase response function θ(ω). The real part is the absorption
spectrum and the imaginary part the dispersion spectrum.

23
Chapter 3

Fourier transform

3.1 Integral transforms

In Chapter 2 we have seen that the frequency response function of an LTI
system G(ω) is the Fourier transform of the impulse response function g(t).
The equivalence of frequency-domain and time-domain representations in
spectroscopy extends beyond the linear regime, as follows from the consid-
erations in Chapter 1. It is a direct consequence of the spectrum being
determined by the solution of the time-independent Schrödinger equation
and of time evolution of the system being governed by the time-dependent
Schrödinger equation, and thus by the same Hamiltonian. In general, the re-
lation between time and frequency domain representations can be expressed
by integral transforms of the form
Z ∞
G(ν) = K(t, ν)y(t)dt . (3.1)
−∞

where K(t, ν) is the kernel of the transform. This concept can be extended
to the relation of time-domain signals to other variables instead of frequency
ν. Table 3.1 provides an overview of important integral transforms.

Table 3.1: Important integral transforms

Integral transform Second variable Kernel K(t, . . .)
Fourier transform (FT) ν e−i2πνt
FT (angular frequency) ω e−iωt
i2πkt
Fourier series k e− T
Laplace transform p ept
1
Hilbert transform τ π(t−τ )

24
Note also that integral transformations can be inverted,
Z ∞
y(t) = K −1 (t, ν)F (ν)dν , (3.2)
−∞

A scaling factor may be required to ensure that after transformation and

subsequent inverse transformation the original function is recovered. Sta-
bility of the inverse transformation with respect to small errors in the input
y(t) depends on properties of the kernel (Chapter 9). The FT is stable, since
the functions in the kernel are mutually orthogonal (Section 3.1.2).
Among the transforms given in Table 3.1, the Laplace transform is of the-
oretical interest for the characterization of relaxation time distributions.
Analytical Laplace transformation is widely used in characterization of elec-
tronic systems. In practice, direct numerical Laplace transformation from
the time domain to the relaxation rate constant domain is instable, i.e. the
problem is ill-posed problem and needs to be solved by regularization meth-
ods (see Chapter 9).
The Hilbert transform relates time-domain signals that correspond to the
absorption and dispersion spectrum, i.e., the Hilbert transform introduces
a 90◦ phase shift. Hence, the Hilbert transform when applied in angular
frequency domain allows to compute the imaginary part from the real part
of the spectrum Y (ω) and vice versa, provided that Y (ω) is the spectrum of
the response y(t) of a causal system to an input signal x(t). Consequently,
real and imaginary part of the spectrum contain the same information for
such a causal system. Note that this does not apply to time domain.
The Hilbert transform can also be used to apply a phase shift to a time-
domain signal. Assume that y(t) is purely real and is given in Matlab
variable signal, while the phase shift is given in variable phi0. The phase-
shifted signal is then computed by shifted = real(hilbert(signal) *
exp(1i*phi0)).

3.1.1 Fourier series

The Fourier series transform in Table 3.1 corresponds to a series expansion
in terms of functions with period T . Index k runs over all integers. Fourier
series can be used for analysis of periodic functions in terms of the basic
frequency and its overtones. For instance, a triangle wave with period T
can be synthesized by superposition of the sine function with period T and

25
its overtones sin(2πk/T ) with prefactors
(
1 (−1)(k−1)/2 for n odd
bk = 2 2 . (3.3)
8π k 0 for n even

If only a finite number of terms is added, the Fourier series provides an

approximation to the signal. The quality of this approximation for a triangle
wave with 1 ≤ k ≤ 5 is shown in Figure 3.1.

y(t)

Figure 3.1: Approximation of a triangle wave by a Fourier series truncated after

k = 5. Ideal triangle wave: solid black line. Approximation: dashed red line.

Such Fourier synthesis of functions with a jump discontinuity is afflicted

by the Gibbs phenomenon, which is illustrated in Figure 3.2. Near the
discontinuity an overshoot occurs. For large numbers of considered terms in
the series the amplitude of this overshoot tends to a constant value, i.e., the
overshoot does not vanish, it only decays faster.1 The Gibbs phenomenon
can lead to spectral artifacts.

a b c

Figure 3.2: Gibbs phenomenon. The overshoot at a jump discontinuity does not
vanish if a larger number of terms is used for Fourier synthesis of a rectangle wave.
(a) Terms up to k = 10. (b) Terms up to k = 50. (c) Terms up to k = 250.
Adapted from Wikimedia Commons.
1
For the infinite series, the overshoot decays infinitely fast, which ensures that in the
limit the function is perfectly reproduced.

26
In general, the Fourier series can be expressed as
∞   ∞  
a0 X 2πkt X 2πkt
y(t) = + ak cos + bk sin . (3.4)
2 T T
k=1 k=1

The real coefficients ak are all zero for odd functions, while the real coeffi-
cients bk are all zero for even functions.

3.1.2 Orthonormal systems of functions

The basis functions
1
Fk (t) = √ ei2πkt/T (3.5)
T
of the Fourier series define an orthonormal system of functions (φi (t))i∈I
with the index set I and the condition
Z t0 +T
hφi (t) , φj (t)it = φi (t)φ∗j (t)dt = δij , (3.6)
t0

where the Kronecker delta δij is defined by

(
1 for i = j
δij = . (3.7)
0 for i 6= j

This guarantees that expansion of a given function y(t) in a Fourier series is

unique. Hence, Fourier analysis (also termed harmonic analysis) is mathe-
matically stable. The system of functions is also complete, which guarantees
that any periodic function y(t) with period T can be represented by an in-
finite Fourier series. The normalization factor √1 is usually omitted.
T
Normalized Legendre polynomials,
r
2k + 1 1 dk 2
k
Pk (t) = k k
t −1 , (3.8)
2 2 k! dt

are another example for an orthonormal system of functions.

27
3.2 Definition and properties of the Fourier trans-
form
3.2.1 Definition
Only periodic signals can be analyzed with the Fourier series. To analyze
arbitray signals, we may go from a discrete set of basis functions to a contin-
uous distribution of such functions. The sum is thus replaced by an integral
and we have for the Fourier transform (FT)
Z ∞
Y (ω) = y(t)e−iωt dt , (3.9)
−∞

and the inverse transform

Z ∞
1
y(t) = Y (ω)eiωt dω . (3.10)
2π −∞

For convenience, we write

Y (ω) = F {y(t)} (3.11)

and
y(t) = F −1 {Y (ω)} (3.12)

The improper integrals converge if y(t) is absolute integrable,

Z ∞
|y(t)| dt . (3.13)
−∞

This is a sufficient, but not necessary condition for the existence of the
Fourier integral.

3.2.2 Basic properties

General mathematical properties

The FT is linear,

F {c1 y1 (t) + c2 y2 (t)} = c1 F {y1 (t)} + c2 F {y1 (t)}

= c1 Y1 (ω) + c2 Y2 (ω) . (3.14)

The conjugate complex function of the time-domain signal transforms to a

28
complex conjugated spectrum with inverted frequency axis,

F {y ∗ (t)} = Y ∗ (−ω) . (3.15)

Inversion of the time axis leads to inversion of the frequency axis

F {y(−t)} = Y (−ω) . (3.16)

Real, even functions in time domain transform to real, even functions in

frequency domain. Real, odd functions in time domain transform to imag-
inary, odd functions in frequency domain. Imaginary, even functions in
time domain transform to imaginary, even functions in frequency domain.
Imaginary, odd functions in time domain transform to real, odd functions
in frequency domain.
For non-negative signals, y(t) ≥ 0, one has |Y (ω)| ≤ Y (0).
Similarity theorem. Scaling in time domain leads to inverse scaling in fre-
quency domain,
1 ω
F {y(at)} = Y( ) . (3.17)
|a| a
Hence, a faster decaying signal leads to broader lines.
Shift theorem. A time delay introduces a phase shift that is linear in fre-
quency,
F {y(t − t0 )} = e−iωt0 Y (ω) . (3.18)

Modulation theorem. Modulation of the time-domain signal with frequency

ν0 leads to a shift by frequency ν0 ,

F eiω0 t y(t) = Y (ω − ω0 ) .

(3.19)

Derivative theorem. Differentiation of the time-domain function can be ex-

pressed as n o
F y (n) (t) = (iω)n Y (ω) . (3.20)

Differentiation of the spectrum can be expressed as

F {(−it)n y(t)} = Y (n) (ω) . (3.21)

29
Integration of the time-domain signal can be expressed as
Z t Z t1 Z tn−1 
F dt1 dt2 . . . dtn s(tn ) = (iω)−n Y (ω) . (3.22)
−∞ −∞ −∞

The area under y(t) is given by the value Y (0) and the area under Y (ω)
(integral over the whole spectrum) is given by y(0).

Parseval’s theorem

For a pair of time-domain data y1 (t) and Y2 (t) and their FTs Y1 (ω) and
Y2 (ω), the following equation holds:
Z ∞ Z ∞
1
y1 (t)y2∗ (t)dt = Y1 (ω)Y2∗ (ω)dω . (3.23)
−∞ 2π −∞

In particular, for y1 (t) = y2 (t) = y(t), we have

Z ∞ Z ∞
21
|y(t)| dt = |Y (ω)|2 dω , (3.24)
−∞ 2π −∞

which can be written in a shorter notation as a relation between the square

norm of the time-domain signal and the square norm of the spectrum

1
ky(t)k2 = ||Y (ω)| |2 . (3.25)
2π

For the relation between time-domain (t) and frequency domain (ν), the
factor 1/(2π) is dropped,

ky(t)k2 = ||Y (ν)| |2 . (3.26)

Hence, for the pair of variables (t, ν) the FT is a unitary operation.

Riemann-Lebesgue lemma

If the time-domain signal is absolute integrable, the spectrum tends to zero

when the frequency approaches −∞ or ∞,
Z ∞
|y(t)| dt < ∞ ⇒ lim Y (ω) = 0 . (3.27)
−∞ |ω|→∞

30
3.2.3 Convolution and deconvolution
Consider the FT of the convolution integral
Z ∞ 
F {y1 (t) ∗ y2 (t)} = F y1 (τ )y2 (t − τ )dτ
−∞
Z ∞ Z ∞ 
= e−iωt y1 (τ )y2 (t − τ )dτ dt
−∞ −∞
Z ∞ Z ∞ 
−iωτ −iω(t−τ )
= y1 (τ )e y2 (t − τ )e dt dτ
−∞ −∞
Z ∞
= y1 (τ )e−iωτ Y2 (ω)dτ = Y1 (ω)Y2 (ω) , (3.28)
−∞

where we have made use of the fact that a double integral does not depend
on the sequence of integration. Hence, the convolution integral in time-
domain transforms to a product in frequency domain Convolution theorem.).
Since the FT can be computed very fast by the Cooley-Tukey algorithm,
convolution can be computed conveniently by FT of the two time-domain
signals, multiplication in frequency domain, and inverse FT.
Furthermore, this property of the FT allows for solving the inverse problem
of deconvolution. Assume that all lines in a spectrum are afflicted by an
additional Gaussian broadening of the same width due to an instrumental
reason. With inverse FT z(t) of the broadening function Z(ω) and the
inverse FT y(t) of the spectrum Y (ω) we can compute a time-domain signal

y(t)
ydeconv (t) = (3.29)
z(t)

and by subsequent FT a deconvoluted spectrum Ydeconv (ω) with narrower

lines. A numerical problem arises from the fact that z(t) tends to zero at
long times, which leads to strong noise amplification at these times. This
problem can be solved by apodization of ydeconv (t) (Section 7.3).

3.2.4 Correspondences
Intuitive qualitative analysis of time-domain signals and intuitive qualitative
design of excitation waveforms with desired spectral selectivity are possible
if one knows a few correspondences between functions in time-domain and
frequency domain.

31
Constant offset

A constant offset transforms to a delta function at zero frequency,

1
F {1} = δ(ω) = δ(ν) . (3.30)
2π

Comb or sampling function

A comb of equally spaced delta impulses in time domain transforms to a

comb of delta functions in frequency domain,

∞ ∞
( )
X 1 X 2πk
F δ(t − kT ) = δ(ω − ). (3.31)
T T
k=−∞ k=−∞

This property ensures a correspondence of time-discrete data with frequency-

discrete data that is the basis of the discrete FT.

y(t) Y(w)

-3T -2T -T 0 T 2T 3T -3w0 -2w0 -w0 0 w0 2w0 3w0

t w

Figure 3.3: Fourier transform of the sampling function. Delta functions in fre-
quency domain are spaced by ω0 = 2π/T .

Rectangular pulse (box function)

A rectangular pulse transforms to a sinc function. With the definition

(
T
1 for |t| ≤ 2
rT (t) = T
. (3.32)
0 for |t| > 2

32
we find
Z ∞ Z T /2
−iωt
Y (ω) = rT (t)e dt = e−iωt dt
−∞ −T /2

1 −iωt T /2 1  −iωT /2 
= − e =− e − eiωT /2
iω −T /2 iω
 
2 ωT sin(ωT /2)
= i sin =T = T sinc (νT ) . (3.33)
iω 2 ωT /2

Note that in signal processing, the sinc function of an argument x, sinc(x),

is defined as sin(πx)/πx, such that the definite integral over all real numbers
equals 1 (normalized sinc function). The result in Eq. (3.33) is important for
understanding excitation profiles of rectangular pulses in frequency domain.
The shorter the pulse, the broader the excitation band.

Exponential decay

An exponential decay function, y(t) = e−a|t| with a > 0 transforms to a

Lorentz line centered at zero frequency.
Z ∞
Y (ω) = e−a|t| e−iωt dt
Z−∞
∞ Z ∞
−at iωt
= e e dt + e−at e−iωt dt
Z0 ∞ Z 0∞
= e(iω−a)t dt + e−(iω+a)t dt
0 0
∞ ∞
1 1
= ei(ω−a)t + − e−i(ω−a)t
(iω − a) 0 (iω + a) 0
1 1 − (iω + a) + (iω − a)
= − + =
iω − a iω + a (iω − a)(iω + a)
2a
= . (3.34)
a2 + ω 2

Gaussian decay

A Gaussian decay in time domain transforms to a Gaussian line in frequency

domain,
n o rπ
−at2 2
F e = exp−ω /(4a) (3.35)
a

33
Harmonic functions

The harmonic functions transform according to

1
F {cos(ω0 t)} = [δ(ω + ω0 ) + δ(ω − ω0 )]
2
i
F {sin(ω0 t)} = [δ(ω + ω0 ) − δ(ω − ω0 )] . (3.36)
2

Exponentially damped harmonic oscillation

Spectroscopic signals in time domain are related to evolution of coherence

on allowed transitions, as it is described by Eq. (1.14). The time-domain
function for quadrature detection of coherence order -1 can be expressed as
  
1
exp iωij − t = e−iωij t · e−|t|/T2,ij . (3.37)
T2,ij

Such a signal transforms to a Lorentz line with width 1/T2,ij centered at ω0

in angular frequency domain,
n o 2/T2,ij 2T2,ij
F e−iωij t e−|t|/T2,ij = 2 + (ω − ω )2 = 1 + T 2 (ω − ω )2 (3.38)
1/T2,ij 0 2,ij 0

Note that the factor of 2 in the denominator of the last two expressions origi-
nates from integration between −∞ and ∞ and would vanish for integration
of a causal response from 0 to ∞. In that case, however, an imaginary dis-
persion contribution would also arise.

3.3 Representations of complex spectra

3.3.1 Absorption and dispersion spectrum
With the cosine FT we obtain the absorption spectrum
Z ∞
A(ω) = y(t) cos(ωt)dt . (3.39)
−∞

With the sine FT we obtain the dispersion spectrum

Z ∞
D(ω) = y(t) sin(ωt)dt . (3.40)
−∞

The dispersion spectrum is the Hilbert transform of the absorption spectrum

34
and vice versa. The absorption spectrum is the real part and the dispersion
spectrum the imaginary part of the complex spectrum.

3.3.2 Magnitude spectrum

In the magnitude spectrum M (ω) (also absolute-value spectrum),
p
M (ω) = A2 (ω) + D2 (ω) , (3.41)

the phase information is lost. Nevertheless, the magnitude spectrum is not

a good substitute for phase correction, as is seen in Figure 3.4. The true
absorption spectrum for two overlapping Lorentz lines can be obtained if
the time-domain data are known starting at t = 0 (Figure 3.4(a)). If a
dead time t0 passes, before observation can start, the absorption spectrum
is distorted by a linear phase shift according to Eq. (3.18), as seen in (b).
For broad overlapping lines and long dead times, linear phase correction
is not successful, since the phase shift in this equation relates to the center
frequency of the lines, so that there is no proper value of the phase correction
for the range in frequency domain where the lines overlap. The problem is
not solved by taking the magnitude spectrum (c). For cases like this, data
reconstruction back into the dead time is required.

1
a 1
b 1
c
Normalized amplitude

0.8 0.8 0.8

0.6 0.6 0.6

0.4 0.4 0.4

0.2 0.2 0.2

0 0 0

-0.2 -0.2 -0.2

10 15 20 25 30 35 10 15 20 25 30 35 10 15 20 25 30 35

n (MHz) n (MHz) n (MHz)

Figure 3.4: Effects of dead time in absorption and magnitude spectra. Model
computation for two exponentially damped oscillations with frequencies ν1 = 20
and ν2 = 24 MHz, equal decay time constants of 0.3 µs, and amplitudes a1 = 1 and
a2 = 0.3. (a) Absorption spectrum obtained from time-domain data starting at
time t = 0. (b) Absorption spectrum obtained from time-domain data starting at
time t = 0.8 µs. (c) Magnitude spectrum obtained from time-domain data starting
at time t = 0.8 µs.

35
3.3.3 Power spectrum
The power spectrum is the square of the magnitude spectrum,

P (ω) = |Y (ω)|2 = A2 (ω) + D2 (ω) = M 2 (ω) , (3.42)

36
Chapter 4

Basics of electronics

4.0.4 Energy, power, amplitude, impedance

The German terminology can be inferred from F. K. Kneubühl, Repititorium
der Physik, Teubner Studienbücher, Section 5.7.
A direct current I flowing for a time t between points with a potential
difference (voltage) V leads to conversion of electrical energy

Eel = V It (4.1)

to heat in the resistance

V
R= . (4.2)
I
The power P = Eel /t is the product of current and voltage, i.e. 1 V·A = 1
W = 1 J s−1 . Hence, we have

V2
P = = RI 2 . (4.3)
R

The power is proportional to the square of the voltage.

For an alternate current or an electromagnetic field, voltage and current are
time dependent
 
Vt (t) = V0 cos (ωt − φ) = Re V0 e−iφ eiωt = Re V eiωt ,


 
It (t) = I0 cos (ωt − ψ) = Re I0 e−iψ eiωt = Re Ieiωt ,


(4.4)

where V and I are the complex voltage and current amplitudes and φ and
ψ are phases, which generally differ for voltage and current. The equivalent

37
to the resistance is the impedance

V
Z= , (4.5)
I

which is a complex quantity, since V and I are complex due to the phase
factors e−iφ and e−iψ . The imaginary part of the impedance is the reactance,
whereas the real part is the resistance.
In the following we assume that the maximal linear dimension dmax and
the minimal duration tmin of the considered electrical processes are in a
range where the speed of propagation of electrical and magnetic fields, the
velocity of light c, can be considered as infinite, i.e. tmin  dmax /c. The
corresponding theory of quasi-static currents cannot be applied in the mi-
crowave range of frequencies and generally not to antennas or delay lines.1
It is still applicable in many, but not all cases to radiofrequency circuits.
In the regime of quasi-stationary currents, impedances in series add (Z =
Z1 + Z2 ), while for impedances in parallel we have
 −1
1 1
Z= + h parallel circuit i . (4.6)
Z1 Z2

For the series circuit shown in Figure 4.1 Ohm’s law provides the voltage
divider rule Spannungsteiler-
Z2
Vout = Vin . (4.7) formel
Z1 + Z2

Z1

Vin

Z2 Vout

Figure 4.1: Voltage divider.

1
Such problems are often treated numerically by finite-element methods for solving
Maxwell’s equations, which are beyond the scope of this lecture course. For simple cases,
such as cavity resonators and waveguides, analytical solutions can be found. A good basic
introduction is given in C. P. Poole, Jr. Electron Spin Resonance, Dover Publications,
Inc., Mineola, 2nd Ed., 1983, Chapters 4 and 5

38
The instantaneous power of an alternate current in a two-terminal circuit is

P (t) = I(t)V (t) = I0 cos (ωt − ψ) V0 cos (ωt − φ)

I0 V0
= [cos (ψ − φ) + cos (2ωt − φ − ψ)] . (4.8)
2

Thus we obtain for the mean power

Z T
1 I0 V0
P̄ = P (t)dt = cos (ψ − φ) = Ieff Veff cos (ψ − φ) , (4.9)
T 0 2
√ √
where Ieff = I0 / 2 and Veff = V0 / 2 are the effective current and voltage,
respectively.

4.0.5 The dB scale

The dynamic range of quantities encountered in spectroscopic measurements
and signal processing can be very large. For instance, in pulse EPR spec-
troscopy spins are excited with powers of several hundred Watt, whereas
signals in the microwatt range are detected and amplified by orders of mag-
nitude. It is then convenient to specify amplitudes and powers on a loga-
rithmic scale.
The dB scale (decibel scale) refers to power ratios and is based on the decadic
logarithm,    
P1 P1
= 10 log10 . (4.10)
P2 dB P2
If a signal is attenuated by 10 dB, the power is reduced by a factor of 10 (one
order of magnitude), whereas 30 dB attenuation correspond to a reduction
in power by a factor of 1000. If a signal is amplified by 3 dB the power
doubles (factor ≈ 1.995, error 0.25%).
Since power is proportional to the square of voltage or current, we have
   
V1 V1
= 20 log10 ,
V2 dB V2
   
I1 I1
= 20 log10 . (4.11)
I2 dB I2

Hence, to reduce amplitude of irradiation or of a signal by a factor of two,

power needs to be attenuated by 6 dB.
The related dBm scale is used to specify absolute powers. This is achieved

39
by defining a reference point of 1 mW power at 1 dBm, so that

P
(P )dBm = 10 log10 . (4.12)
1 mW

Accordingly, 30 dBm denote a power of 1 W, while -30 dBm denote a power

of 1 µW.

Resistor Capacitor DC voltage source

R C
(IEC) +

(IEEE)

Inductor Diode AC voltage source

L
(IEC)

(IEEE)

Figure 4.2: Symbols for selected circuit elements. For resistors and inductors,
both the American (IEEE) and European (IEC) symbol are shown.

4.1 Circuit elements and rules

4.1.1 Circuit elements
Symbols of frequently encountered circuit elements are shown in Figure 4.2.
For our purposes resistors, inductors, and capacitors are the most important
elements. They form the basis of RLC networks, which can be used as res-
onance circuits and filters. Circuit diagrams represent resistors, inductors,
and capacitors by pure resistances R, inductances L, and capacitances C.
To compensate for non-ideality of physical circuit elements, it may be nec-
essary to insert more elements into the circuit diagram than are physically
present. The circuit diagram thus shows an idealized equivalent circuit of
the real device. For pure resistances, inductances, and capacitances a linear
relationship between current through the two-terminal element and voltage
between the two poles is observed.
Most linear circuits, where such a relation between current and voltage ex-
ists can be described as RLC networks, occasionally combined with linear
quadripoles ( two-port networks).

40
Resistors

A resistor is a two-terminal component with a purely Ohmic resistance R,

corresponding to a purely real impedance

ZR = R . (4.13)

Even if no explicit resistor occurs in a circuit, resistance of leads and induc-

tors may be described by an equivalent resistor. In magnetic resonance spec-
troscopies, the sample itself may correspond to a resistor (dielectric losses).
The frequency dependence of such dielectric losses can provide information
on structural relaxation processes, which can be obtained by dielectric spec-
troscopy. Resonance absorption followed by dissipation of the excitation
energy to the lattice can also be described by a (frequency-dependent) re-
sistance. From Eq. (4.2) it follows that

V (t) = RI(t) , (4.14)

so that a resistance does not introduce a phase shift between voltage and
current.

Inductors

The most simple inductor is a coil of wire. Inductance can be increased by

inserting a soft iron core or ferrite ceramics into the coil. In NMR spec-
troscopy, real coils are used for inductive coupling of the spin magnetization
to the excitation field and the detector. In EPR spectroscopy, most experi-
ments also rely on inductive coupling, with the inductance L being spatially
distributed over a resonator.
Although real coils, except for superconducting coils, have some resistance,
they are represented in circuit diagrams by a resistance-free inductance L.
The linear relation between voltage and current is

d
V (t) = L I(t) . (4.15)
dt

Since the derivative of a cosine function is the negative sine function (phase
shift -90◦ ), in an ideal inductance circuit the current lags the voltage by

41
ψ − φ = 90◦ . From Eqs. (4.4,4.5,4.15) we can derive

ZL = iωL , (4.16)

i.e., the impedance of a pure inductance is purely imaginary.

In the inductance, the energy is stored in a magnetic field. The SI unit is
the Henry, 1 H = 1 Wb A−1 , whereas the Weber (1 Wb = 1 T m2 = 1 V s
= 1 J A−1 ) is the unit of magnetic flux.

Capacitors

The most simple capacitor are two conducting plates separated by air. An
ideal capacitor has infinite resistance for direct currents. Alternate currents
R
can flow, as the charge Q = Idt can be stored in the capacitance C. The
linear relation between voltage and current is

d 1
V (t) = I(t) . (4.17)
dt C

In an ideal capacitance circuit the voltage lags the current by 90◦ , i.e.,
ψ − φ = −90◦ .
In loose terminology it could be said that inductance and capacitance are in
a inverse relationship with each other. More precisely, while an inductance
differentiates current, a capacitance integrates it
Z t
1 Q(t)
V (t) = V (0) + I(t)dt = . (4.18)
0 C C

From Eqs. (4.4,4.5,4.18) we can derive

1
ZC = (4.19)
iωC

i.e., the impedance of a capacitance is also purely imaginary.

Lumped element model

To microwave and optical components the concept of quasi-static current is

not applicable. Devices cannot constructed from elements that behave like
resistors, inductors, and capacitors. Instead, the resistance, inductance, and
capacitance are distributed in spatially overlapping regions. Such systems
may still behave like an RLC circuit and it may thus be possible to represent

42
them by an equivalent circuit diagram 2 .
This approach provides an approximation of the real system by a lumped
element model, which is equivalent to reducing the state space of the real
systems, which is continuous in time and space to a finite state space with
a finite number of parameters. For instance, a loop-gap resonator used
in EPR spectroscopy and in large-volume NMR and magnetic resonance
imaging (MRI) applications can be understood to a good approximation as
an RLC resonance circuit. Not all aspects can be understood in this model,
though. For instance, the spatial distribution of the electric and magnetic
fields in the resonator can only be modeled by renouncing a description by
quasi-static currents.

4.1.2 Rules for combining elements to a circuit

A R B C R
I1 I -I2 I1 I2 I1
+ I III -
-I3 I
I1 + I2 + I3 = 0

C V1 V2 V1 C V2

- II IV + II
II I1 I2

Figure 4.3: Two-port networks. (A) RC low pass filter with illustration of Kirch-
hoff’s current law for junction I. (B) General scheme for a two-port network. (C)
RC low pass filter.

Kirchhoff circuit laws

Kirchhoff’s current law (also Kirchhoff’s junction rule) derives from conser-
vation of charge at any junction in a circuit. The current flowing into any
junction in a circuit must equal the current flowing out of that junction. For
a junction with n branches, it can be formulated as
n
X
Ik = 0 . (4.20)
k=1
2
German: Ersatzschaltbild

43
For instance, for the LR low pass filter circuit shown in Figure 4.3A this law
must hold at junctions I and II, which have three branches each.
Kirchhoff’s voltage law (also Kirchhoff’s mesh rule) derives from conserva-
tion of energy. The directed sum of electrical potential differences around
any closed circuit must be zero. The reason is that a charge that goes full
circle in a closed circuit may not have gained or lost energy. As Kirchhoff’s
voltage law considers only the electrical field, it does not hold true in the
presence of fluctuating magnetic fields. This law thus cannot be applied
without corrections to circuits that contain inductances.

Linear two-port networks

The LR low pass filter shown in Figure 4.3A is a special case of a linear two-
port network (Figure 4.3B). The ports are terminal pairs I/II and III/IV.
They satisfy the condition that the sum of signed currents to the two port
terminals must be zero. A two-port network has an input voltage V1 , an
output voltage V2 , and input current I1 and an output current I2 . These
quantities are related by an impedance matrix
" # " #" #
V1 z11 z12 I1
= (4.21)
V2 z21 z22 I2

If the output port (III/IV) is loaded with a resistance R0 , the first port can
be considered as a linear two-terminal element. Open-circuit operation is
defined as I2 (t) ≡ 0. A two-port network can be characterized by measuring
V2 (t) as a function of V1 (t) or of I1 (t) in such open-circuit operation.

L
I1 I2
I III

V1 R V2

II IV
I1 I2

Figure 4.4: LC low pass filter as an example of a linear two-port network.

44
As an example of a linear two-port network we consider the LC low-pass
filter shown in Figure 4.4. The voltage drop of this network is due to the
inductance. With Eq. (4.15) we find

d
V2 (t) = V1 (t) − L I1 (t) . (4.22)
dt

The difference between currents I1 and I2 can be expressed with the current
flow through the resistance R. The voltage across this resistance is V2 . With
Kirchhoff’s current law, Eqs. (4.2) and (4.22) we find
 
1 L d
I2 (t) = − V1 (t) + 1 + I1 (t) . (4.23)
R R dt

For the loaded LR low pass filter we have the additional equation V2 (t) =
R0 I2 (t). Combining this with Eq. (4.22) we can solve for the relation be-
tween voltage and current at the input port
 −1
1 1 d
V1 (t) = + 0 I1 (t) + L I1 (t) . (4.24)
R R dt

For open-circuit operation we have I2 (t) ≡ 0, so that V2 must be the voltage

drop for current I1 flowing through the resistance R,

V2 (t) = RI1 (t) . (4.25)

The relation between voltage and current at the input port follows directly
from setting the left-hand side of Eq. (4.23) to zero and solving for V1 ,

d
V1 (t) = RI1 (t) + L I1 (t) . (4.26)
dt

4.2 Filters
4.2.1 RC low pass
We now consider the RC low pass filter shown in Figure 4.3C in open-
circuit operation. The current I1 is related to the voltage drop VR over the
resistance,
VR = I1 R . (4.27)

The circuit from terminal I via resistance R and capacitance C to terminal

II is closed. As no inductance exists in this circuit, we can apply Kirchhoff’s

45
voltage law and obtain
V1 = VR + VC , (4.28)

where V2 = VC is the voltage drop over the capacitance. From Eq. (4.17)
we have
dV2 I1
= , (4.29)
dt C
and thus, by substituting VR in Eq. (4.28) with the help of Eqs. (4.27) and
(4.29),
dV2
V1 = V2 + RC . (4.30)
dt
In linear response theory, we can interpret V1 (t) as the input signal x(t) and
V2 (t) as the output signal y(t). Defining the differential operator D b = d,
dt
we thus have  
x(t) = 1 + RC D
b y(t) , (4.31)

which can be solved for y(t) and gives

1
y(t) = {x(t)} = S {x(t)} . (4.32)
1 + RC D
b

The system operator S is thus the inverse of the sum of a unit operator and
the product of the differential operator with resistance and capacitance.3
We now proceed with computing the step response of the RC low pass filter.
By solving Eq. (4.31) for the derivative of y(t), we obtain a differential
equation for the response

dy 1
= (x(t) − y(t)) . (4.33)
dt RC

This equation has to be solved for a step function as the input

(
0 for t < 0
x(t) = , (4.34)
V0 for t ≥ 0

where V0 is the amplitude of the step. As the system is causal, we need to

compute the response only for t > 0, so that x(t) is constant. The solution
of the corresponding homogeneous differential equation is

y(t) = Ke−t/(RC) . (4.35)

3
Problems like this can be conveniently solved via the Laplace transform. We do not
introduce this formalism in the current lecture course. A good introduction is given in M.
Meyer, Signalverarbeitung, Vieweg+Teubner, Wiesbaden, 6. Aufl. 2011

46
In the stationary state, we must have V2 = V1 = V0 . Hence, the solution
has the form
y(t) = V0 + Ke−t/(RC) . (4.36)

The initial condition y(0) = 0 gives K = −V0 , so that we obtain

h i
y(t) = V0 1 − e−t/(RC) . (4.37)

Hence, the RC low pass filter approaches the stationary condition with a
time constant τ = RC. This is the rise time of the RC circuit.4
The formal step response function a(t) corresponds to V0 = 1, hence

a(t) = 1 − e−t/(RC) . (4.38)

We can obtain the impulse response by using Eq. (2.16),

d 1 −t/(RC)
g(t) = a(t) = e for t ≥ 0 . (4.39)
dt RC

1 0

0.8
A B
1/Ö2 3 dB
0.6
M(w)

q(w)

-45
0.4

0.2

0 -90

0 1 2 3 4 5 0 1 2 3 4 5
RCw RCw

Figure 4.5: Transfer function (frequency response) of an RC low pass filter. (A)
Amplitude function. (B) Phase response.

The frequency-response function G(ω) can in turn be computed as the FT

4
Note that this definition does not conform to the general engineering definition of the
rise time as the time passed between 10% and 90% of the stationary response.

47
of the impulse response function,
Z ∞
1
G(ω) = e−[iω+1/(RC)]t dt
RC 0
1 1 1
= =
RC iω + 1/(RC) 1 + RCiω
1 RCω
= −i , (4.40)
1 + R2 C 2 ω 2 1 + R2 C 2 ω 2

where the terms in the last line are the real and imaginary part. Hence, the
magnitude spectrum is
r
1
M (ω) = , (4.41)
1 + R2 C 2 ω 2

while the phase response function is

θ(ω) = arctan (−RCω) . (4.42)

The magnitude spectrum and phase response function of an RC low pass

filter are shown in Figure 4.5. For large frequencies the magnitude ap-
proaches zero and the phase approaches 90◦ , however, the approach is slow.
Asymptotic behavior can better be seen in the Bode plots (Figure 4.6) with
logarithmic frequency axes. The slope of the Bode gain plot (Figure 4.6A)
in the stop band is an important criterion for the performance of the filter,
as it tells how strongly unwanted frequency contributions are suppressed.
In the case at hand, the slope is -20 dB/decade of frequency change or -6
dB/octave. An octave corresponds to a doubling of the frequency.5 The di-
vision between pass band and stop band corresponds to the inflection point
in the Bode phase plot (Figure 4.6B), which corresponds to a phase shift of
45◦ . This characteristic frequency of the RC low pass filter is given by

1
ωRC = . (4.43)
RC

It is thus the inverse of the rise time in the step response.

4.2.2 Ideal and real filters

An ideal filter enables undistorted transfer of the signal in the pass band and
completely suppresses the signal in the stop band. Depending on the desired
5
This term is related to classical music.

48
 0
pass band 0

-10 A B
Magnitude (dB)

stop band
-20

q(w)
-45
-30

-40
-90
-50
-5 0 5 -5 0 5
log10(RCw) log10(RCw)

Figure 4.6: Bode plot for an RC low pass filter. (A) Magnitude function. (B)
Phase response.

pass and stop behavior, we distinguish low pass (Figure 4.7A), high pass (B),
band pass (C), and band stop (D) filters. Ideal filters cannot be realized
by analog circuits. The reason can be recognized by considering the ideal
band pass filter. Such a filter corresponds to multiplying the signal Y (ω)
in frequency domain by a rectangular function and, thus, to convolution of
the time-domain signal with a sinc function. Since the sinc function extends
to infinity towards both negative and positive times, this would imply that
the signal is delayed for an infinite time in the filter to avoid an unphysical
acausal response.

A B
M(w)

C D
M(w)

w w

Figure 4.7: Idealized magnitude spectra of different filter types. (A) Low pass
filter. (B) High pass filter. (C) Band pass filter. (D) Band stop filter

Filter design is based on a theory for normalized low pass filters. Depending
on the application, different approximations to an ideal filter are aimed for.

49
 • Butterworth filters: M (ω) is as flat as possible in the pass band.

• Chebyshev I filters: Accepting a defined ripple in the pass band, tran-

sition to the stop band is steeper than for a Butterworth filter.

• Chebyshev II filters: Accepting a defined ripple in the stop band,

transition to the stop band is steeper than for a Butterworth filter.

• Cauer filters: Defined ripple is accepted both in the stop and pass
band, transition to the stop band is even steeper than with Chebyshev
filters.

• Bessel filters: The group delay should be constant in the pass band, i.e.
different frequency components are delayed by the same time. This
preserves waveform.).

• Critically damped filters: There are no oscillations in the impulse and

step response.

For instance, the RC low pass is a critically damped filter. However, M (ω)
is not particularly flat in the pass band, the transition to the stop band is
not very steep, and the group delay in the pass band differs significantly.
Improvement of the approximation generally requires adding more elements.
This concept can be formalized via the Laplace transform H(s) of the im-
pulse response g(t). For linear, time-invariant (LTI) systems with a finite
number of lumped elements (resistances, inductances, and capacitances), the
function H(s) can be written in rational form (partial fraction expansion)
and written as

(s − sz1 )(s − sz2 ) . . . (s − szm )

H(s) = k0 (4.44)
(s − sp1 )(s − sp2 ) . . . (s − spn )

This function has m zeroes sz1 . . . szm and n poles sp1 . . . spm , which can
be real numbers or occur as complex conjugate pairs. The number of poles
is the filter order. Filters can be realized not only with passive elements
(RLC), but also with active elements, such as operational amplifiers. In the
latter case, the filter requires an external power source.

50
 1 A 0 B
-20
0.8

Magnitude (dB)
-40

0.6 -60
M(w)

-80
0.4
-100
0.2 -120

-140
0
-160
0 1 2 3 4 5 -1 0 1 2 3
w/W log10(w/W)

Figure 4.8: Butterworth filters of different order: 1 black, 2 blue, 3 magenta, 4

red, 5 green. (A) Magnitude function. (B) Bode gain plot.

4.2.3 Butterworth filters

The frequency response function of a Butterworth filter is
s
1
M (ω) = . (4.45)
ω 2n


1+ Ω

As seen in Figure 4.8 steepness of the transition from pass band to stop band
and flatness of the frequency response in the pass band increase strongly with
filter order.

4.2.4 Chebyshev filters

The frequency response function of a Chebyshev type I filter is
s
1
M (ω) = ω

, (4.46)
1+ 2 c2n Ω

where  is the ripple and cn is the Chebyshev polynomial of nth order.

51
Figure 4.9: Bode gain plot for a Chebyshev type I filter with order n = 4 and
 = 1 (3 dB ripple). Taken from Wikimedia Commons.

52
Chapter 5

Noise

5.1 Noise sources and types

5.1.1 Noise sources and signal-to-noise ratio
We consider any physical measurement on a quantum system, with spec-
troscopy being one particular class of such measurements. Noise is a random
fluctuation of the signal and can originate from

• The system itself, including quantum system and coupling to detection

arm.

• Excitation of the system.

• Signal amplification.

• Environment.

Quantum noise is related to uncertainty of measurements on quantum sys-

tems and is important mainly if the ensemble is small. Classical systems and
large ensembles of quantum systems are subject to thermal noise (5.1.2).
The system under observation is coupled in some way to the measurement
receiver, for instance, in NMR by inductive coupling to a coil. This coil is
also a noise source.
Coupling of the system to the measurement device thus generates the signal,
which is the wanted response to the excitation. The excitation waveform
usually also features unwanted fluctuations of amplitude or phase. These
fluctuations are noise superimposed on the input signal and translate into
noise superimposed on the output signal. Together with the noise of the

53
receiver and the system itself the excitation noise contributes to the noise
level of the primary output signal.
The level of the wanted signal (ideal response of the quantum system) can be
quantified by the power Psignal and the noise level can be quantified by the
power Pnoise . The quality of the primary output signal is thus characterized
by the signal-to-noise ratio (SNR)

S Psignal A2signal
= = 2 , (5.1)
N Pnoise Anoise

where Asignal and Anoise are the amplitudes of the signal and noise, respec-
tively. The signal-to-noise ratio can be given on the dB scale. In spec-
troscopy, the signal-to-noise ratio is often defined as a ration of amplitudes
of the signal and noise rather than a power ratio. Care must be taken which
definition is adhered to by a given author. The ’spectroscopic’ signal-to-noise
ratio is the square root of the ’electronic’ signal-to-noise ratio.
Usually the primary signal is too weak for direct display or digitization and
is amplified. Such amplification, and indeed any modification of the signal
by electronic components can lead to a decrease in the signal-to-noise ratio
by introducing additional fluctuations.
Finally, the receiver as well as components of the spectrometer that are not
ideally shielded pick up unwanted signals from the environment. Examples
are signals from cell phone networks, radio and TV signals, fluctutations of
the line voltage, and vibrations or fluctuating magnetic and electrical fields
from traffic.

A B
ÖáV 2ñ R

Öái 2ñ

Figure 5.1: Equivalent circuits for the description of thermal noise in resistors.
(A) Noise AC source. (B) Noise current generator.

54
5.1.2 Thermal noise
The noise in any circuit with uniform temperature T can be described by
a noise AC source (Figure 5.1(A)) in series with each resistance R. For a
small frequenccy interval dν we have

V 2 = 4kB T Rp(ν)dν , (5.2)

where kB is Boltzmann’s constant and p(ν) the Planck factor,

hν 1
p(ν) = hν . (5.3)
kB T e kB T − 1

Up to microwave frequencies at ambient temperature the Planck factor is

equal to unity, as hν  kB T . This does not apply for high microwave
frequencies at cryogenic temperatures (hν ≈ kB T for ν ≈ 100 GHz at T ≈
4.2 K) and it certainly does not apply in the infrared or optical range. The
Planck factor at T0 = 290 K is plotted as a function for the decadic logarithm
of frequency in Figure 5.2.

0.8

0.6
p(n)

0.4

0.2

0 5 10 15
log10(n/Hz)

Figure 5.2: Planck factor p(ν) at t0 = 290 K as a function of decadic logarithm of

frequency in units of Hertz.

Sometimes it is more convenient to describe noise by a current generator of

infinite impedance connected in parallel to the resistance (Figure 5.1(b)),

1
i2 = 4kB T p(ν)dν , (5.4)
R

55
or to describe the noise source in terms of available noise power (compare
Eq. (5.1),

1 V2 1 2
Pnoise = = i R
4 R 4
= kB T p(ν)dν . (5.5)

Note also that hV i = 0 for a noise source.

If Planck’s factor can be set to unity, Eqs. (5.2) and (5.5) can be integrated
over the bandwidth B of the receiver. We find

V 2 = 4kB T BR (5.6)

and
Pnoise = kB T B . (5.7)

The effective noise bandwidth B of a system with respect to the center

frequency ν0 can be computed from the frequency response function G(ν)
(or the magnitude spectrum M (ν)) by
Z ∞ Z ∞
1 1
B= |G(ν)|2 dν = M 2 (ν)dν . (5.8)
|G(ν0 )|2 0 M 2 (ν0 ) 0

The center frequency is usually defined as the frequency where M (ν) attains
its maximum. Hence, B corresponds to the width of a rectangular bandpass
filter (brick wall filter) with the same maximum throughput and the same
integrated area of the frequency response function as the actual system.
According to Eq. (5.7), noise power thus does not depend on frequency in
the range where the Planck factor is unity. Such noise is termed white noise.
As an example, consider a resistance of 50 Ω at T = 290 K in a receiver with
a bandwidth of 5 MHz. The mean square noise amplitude is V 2 = 4·10−12
V2 , corresponding to about 2 µV root mean square amplitude of the noise.
The noise power, which does not depend on the resistance, is 2 · 10−14 W or
-107 dBm.
From these considerations we can derive two basic possibilities for improv-
ing SNR. First, we could keep the receiver at low temperature. This is
an expensive way of improving SNR, but is done in radioastronomy and
liquid-state NMR (cryoprobe technology). In both cases the expense of the
additional measurement time needed with an ambient temperature receiver
is higher than the expense for a cryoreceiver. Second, bandwidth of the

56
receiver should be limited to the minimum that is required to detect the
signal of interest. Such bandwidth limitation can usually be achieved at low
expense and is commonly applied in spectroscopy.

5.1.3 Shot noise

Weak signals can be due to currents that involve a relatively small number
of electrons. Consider, for instance, a 10 ns wide echo signal in EPR spec-
troscopy with an amplitude of 1 µV at 50 Ω impedance. The underlying
current of 20 nA during a time of 10 ns corresponds to N = 1250 elec-
trons. According to Poisson statistics, the SNR, here related to the current
amplitude,1 is given by
S N √
=√ = N , (5.9)
N N
which in this case is 35.
Such noise related to the quantization of charge is termed shot noise. From
Eq. (5.9) we can derive
i2 = 2eiB , (5.10)

where e is the elementary charge. We see that shot noise is also white noise.

5.1.4 1/f noise and Brownian noise

Several important quantities in solid-state materials exhibit slow fluctua-
tions, for instance, domain structures in magnetic materials or defect con-
figurations. These fluctuations cause noise with a spectral power density
that is inversely proportional to frequency f (or ν or angular frequency ω).
In other words, the power at constant bandwidth falls with 3 dB per oc-
tave for increasing frequency. Therefore, 1/f noise is unimportant at high
frequencies.
Since the noise originating from Brownian motion has a 1/ν 2 dependence and
is called red (or brown) noise, 1/f noise is also termed pink noise. Another
synonym is flicker noise. Brownian noise is not important in electronics.
However, it may be a form of noise of the observed system, if random walk
processes influence the observed quantity.
Since the scaling behavior of 1/f noise is known and white noise (thermal
noise plus shot noise) can be computed, 1/f noise can be simply character-
ized by the corner frequency νc . This is the frequency, where for a given
1
In terms of power, the SNR is the square of the number given here.

57
system the power of 1/f noise equals the power of white noise. Above νc ,
white noise dominates, whereas below νc 1/f noise dominates. For bipolar
transistors, the corner frequency is in the kilohertz range, whereas for the
smallest metal–oxide–semiconductor field-effect transistors (MOSFETs) in
modern integrated circuits, it can reach the Gigahertz range.

-40 -20 0 20 40
n (kHz)

Figure 5.3: Phase noise. Simulation for a normal distribution of phase with stan-
drad deviation of 0.1 rad. A frequency of 10 kHz and a transverse relaxation time
of T2 = 0.2ms were assumed.

5.1.5 Phase noise

For some types of frequency sources, phase noise can be important. In
such cases, the phase of the excitation wave varies stochastically with time.
If a signal is acquired by a frequency scan (or field scan as in EPR spec-
troscopy), this leads to different phase of the signal at different frequency
(or field) offsets from the center frequency. Accordingly, different fractions
of dispersion signal are admixed to the absorption line. Such phase noise
can be recognized by its occurence only in the frequency range of the signal,
as is illustrated in Figure 5.3. Phase noise vanishes in the magnitude spec-
trum, albeit at the expense of line broadening compared to the absorption
spectrum.

58
5.2 Mathematics of noise
5.2.1 Noise spectral density and noise figure
Since white noise frequently dominates and since the power of white noise
is proportional to the effective noise bandwidth B, the noise level is often
characterized by the noise spectral density N0 , with

Pnoise
N0 = hwhite noisei . (5.11)
B

The unit of N0 is Watts/Hertz, which is equivalent to the energy unit Joule.

Specifications of amplifiers often use the dBm scale for Pnoise . For instance,
an 1 kW X-band travelling wave tube amplifier for EPR spectroscopy is
specified with 30 dBm noise power for a band between 8 and 12.4 GHz.
This corresponds to -66 dBm/Hz noise spectral density.
From Eq. (5.7) we find for the spectral density of thermal noise

Nth = kB T . (5.12)

For instance, at 290 K we find Nth = 4 · 10−21 W Hz−1 , which corresponds

to -174 dBm.

Si So
System
Ni No

Figure 5.4: If a signal passes a system, the SNR at the output (So /No ) is generally
different from the SNR at the input (Si /Ni ).

The influence of a system on the SNR (Figure 5.4) can be characterized by

the noise factor F > 1, which is the ratio betwen the signal-to-noise ratios
at the input and output of the system,

Si /Ni
F = . (5.13)
So /No

This equation can be applied without further thought only if the input noise
is thermal noise at the standard temperature T0 = 290 K, as will become
clear below. If the actual physical temperature is T0 , the noise factor can

59
be expressed by an effective noise temperature of the system,

Te = (F − 1) T0 . (5.14)

The noise factor is often specified on the dB scale (Section 4.0.5) and is then
termed noise figure,
FdB = 10 log10 F . (5.15)

Note that in large parts of the literature, the dimensionless noise factor is
also called noise figure.
Noise figures of uncooled low-power amplifiers used in detection arms of
spectrometers are typically in the range between 1.5 and 2 dB up to fre-
quencies of 6 GHz and between 2 and 5 dB up to 35 GHz. High-power
amplifiers have much higher noise figures, in particular, at high frequencies.
For instance, travelling wave tube amplifiers in EPR spectroscopy have noise
figures of about 15 dB.
Components of a system can be classified as components without gain (at-
tenuators) and components with gain G (amplifiers). Components without
gain are characterized by their attenuation2

Pi
L= . (5.16)
Po

Their noise factor is given by

(L − 1)T
F =1+ , (5.17)
T0

where T is the actual temperature of the attenuator. Hence, cooling of

attenuators reduces noise.

ÖáV 2ñ G

Figure 5.5: Equivalent circuit for a noisy amplifier.

2
Note that attenuation is usually specified on the dB scale, L = 10LdB /10 .

60
For amplifiers, the gain is specified as

Po
GdB = 10 log10 . (5.18)
Pi

For noise considerations, the amplifier can be represented by the equivalent

circuit shown in Figure 5.5, which consists of an ideal amplifier with the
same gain G = Po /Pi and noise figure F = 1 (noise temperature 0 K) and a
noise source at the effective noise temperature Te . According to Eq. (5.7),
the noise source adds a noise power Pa,0 = kB Te B. Due to the gain, the
total added noise power is Pa = GkB Te B.
Thus, we have
No = GNi + GkB Te , (5.19)

where Te is generally different from the actual ambient temperature. The

noise factor of the amplifier, as defined in Eq. (5.13), is a measure for the
performance compared to an ideal amplifier that just adds and amplifies
thermal noise at the standard temperature (Te = T0 ) and is defined given
by
Te Nadded
F =1+ =1+ , (5.20)
T0 GNth
where Nadded is the noise added by the amplifier as shown in the equivalent
circuit diagramm in Fig. 5.5 and Nth is the thermal noise at ambient tem-
perature, Ni = Nth = kB T0 . At this point we have defined the noise figure
of an amplifier with respect to thermal noise at standard temperature and
have thus restricted validity of Eq. (5.13).
Eq. (5.19) can then be written as

No = GkB T0 + GkB (F − 1)T0 = F GkB T0 = F GNi , (5.21)

which leads to
No No
= = FG , (5.22)
Ni Nth
So
As Si = G, we thus recover Eq. (5.13), as long as the input noise is thermal
noise at the standard temperature. Noise factors of amplifiers thus have to
be determined with thermal noise at the standard temperature T0 as the
input by measuring noise power at the output.

61
5.2.2 Cascaded systems
If several systems are cascaded, each of them can be characterized by its
noise factor Fj and gain Gj (Figure 5.6). On order to compute the noise
figure F of the cascade of systems with amplification, we cannot simply
apply Eq. (5.13) to the inputs and outputs of the individual systems. Such
an approach would fail, since after the first amplifier the input noise is no
longer thermal noise.
Hence, we have to compute the noise factor F by considering gains and
added noise in all steps.3 For convenience, we define the noise ratio nj after
the j th stage of the cascade as nj = Nj /Nth , where Nj is the noise power
after this stage and Nth is the thermal noise power. For the first stage, Eq.
(5.22) yields
n1 = F1 G1 . (5.23)

From Eq. (5.20), we have

Nadded = (Fj − 1)Gj Nth . (5.24)

We can now define a noise ratio nadded = Nadded /Nth . Since the two noise
ratios add, we obtain after the second stage

n2 = n1 G2 + nadded,2 = F1 G1 G2 + (F2 − 1)G2 , (5.25)

where we have considered that the input noise is amplified by a factor G2 .

We now want to express this finding by a noise factor for the two-step
cascaded amplifier with total gain G = G1 G2 and noise ratio n2 = F G =
F G1 G2 . Solving for F yields

n2 F2 − 1
F = = F1 + (5.26)
G1 G 2 G1

In adding further stages, we can use a generalization of Eq. (5.25),

j+1
Y
nj+1 = nj Gj+1 + nadded,j+1 = F (j) Gk + (Fj+1 − 1)Gj+1 (5.27)
k=1

as well as
nj+1
F = F (j+1) = Qj+1 , (5.28)
k=1 Gk
3
This derivation is based on the one in A. van der Ziel, Noise, Prentice Hall, Inc. 1954

62
where F (j) is the noise factor of the cascade after j steps. Combining Eq.
(5.27) and (5.28), we find

Fj+1 − 1
F (j+1) = F (j) + Qj . (5.29)
k=1 Gk

Complete induction then provides Friis formula

F2 − 1 F3 − 1 Fn − 1
F = F1 + + + ... + . (5.30)
G1 G1 G2 G1 G2 . . . Gn−1

Since amplifier gain factors are large numbers, this sum is usually dominated
by F1 . If several amplifiers are cascaded, the one with the lowest noise figure
should thus be the first one in the chain.

A G, F B G1, F1 G2, F2 Gn, Fn

¼
No
Nth Nth N1 N2 Nn= No

Figure 5.6: Noise figure of amplifiers and amplifier cascades. (A) The noise figure
F = No /(Ni G) of a single amplifier is determined by terminating input with a
resistance R that creates thermal noise Nth and measuring noise power No at the
output. (B) In a cascade, input noise of a given amplifier is output noise of the
previous amplifier rather than thermal noise. Noise added by the first amplifier will
dominate, as it is amplified in all further stages.

5.2.3 Techniques for improving the signal-to-noise ratio

General possibilities for noise suppression

From the above considerations it follows that the signal-to-noise ratio im-
proves at low temperatures T of the system and narrow bandwidths B. Since
noise is stochastic, it will also partially average on summation of n individ-
√
ual measurements, Nn = nN , whereas the signal simply adds, Sn = nS.
It follows that
S √ √
∝ n ∝ tm , (5.31)
N
where tm is the time spent for the accumulation of the n measurements.
With respect to 1/f noise and slow signal drifts, it is advantageous to mea-
sure at frequencies larger than 1 kHz, or, more generally, larger than the

63
corner frequency, if the quantum system permits. Whether or not the quan-
tum system permits this, may depend on relaxation times and resonance
frequencies.
To reduce spectrometer noise, components can be cooled. According to Friis’
formula, cooling the first detection amplifier is usually sufficient. Effect
modulation with phase sensitive detection can reduce bandwidth and ist
treated in the next Chapter.
Environmental noise can be reduced by shielding of components against ex-
ternal electromagnetic fields. In some cases, the external fields can also
be measured and compensated for. In NMR spectroscopy, slow drifts of the
magnetic field are compensated for by a spin lock, while in EPR spectroscopy
a field-frequency lock can offset field drifts and drifts in the resonator fre-
quency. Effect modulation also reduces the influence of environmental noise
by limiting detection bandwidth.

Impedance matching

The noise figure of an amplifier depends on source resistance Rs and fre-

quency. The first amplifier in a receiver should be selected such that the
best possible noise figure at the typical signal frequency is achieved. The
optimum source resistance Rs,opt of this amplifier at the signal frequency
may not match the true source resistance. This problem can be solved by
putting a transformer between the signal source and the amplifier input.
At low frequencies a transformer is realized by two coupled coils (induc-
tances) and is characterized by the ratio n between the number of coil wind-
ings at the amplifier side and the number of coil windings at the source
side. The impedance transforms with n2 . Hence, the winding ratio can be
computed as
Rs,opt
n2 = . (5.32)
Rs
At higher frequencies, impedance matching can be achieved with other tech-
niques. To keep things simple, components and transmission lines in high-
frequency electronics have a standardized impedance of 50 Ω. In a magnetic
resonance spectrometer, the resonance circuit for signal detection (LC res-
onance circuit in NMR, microwave resonator in EPR) may have a different
impedance, which slightly depends on sample due to dielectric losses in
the sample. Therefore, the resonator is matched to the 50 Ω impedance.
In NMR probeheads, this is achieved by variation of a capacitance. In

64
microwave resonators, coupling between transmission line and resonator is
matched by varzing the size of a coupling hole (iris coupling) or the position
of an antenna with respect to the resonator.

65
Chapter 6

Effect modulation and

phase-sensitive detection

6.1 Effect Modulation

6.1.1 General principle of modulation
Low-frequency signals are modulated for amplification and transmission to
reduce 1/f noise and, in the case of transmission, to fully use the bandwidth
of transmission channels and to allow for long-distance transmission through
the atmosphere. In spectroscopy, effect modulation can be used to effectively
reduce detection bandwidth and thus noise power.

signal modulated amplified

Modulator signal
Demodulator signal
Filter
carrier

carrier

Figure 6.1: General scheme of modulation.

The general principle of modulation is shown in Figure 6.1. The signal of

interest is superimposed on a carrier in a modulator. In effect modulation,
the modulator is the sample. The modulated signal can be amplified with
a selective amplifier at the carrier frequency, possibly transmitted via cable
or through space, and demodulated from the carrier. After demodulation,
components outside the frequency band of the signal can be suppressed by
a filter.

66
Coherent demodulation is used if a carrier with the same frequency and
phase as used in modulation is available at the receiver. If that is not the
case, incoherent demodulation can be used, which is realized in the simplest
case by a diode and a low-pass filter.

6.1.2 Types of modulation and distortions

There are many different types to modulate a signal, in particular, in trans-
mission of digital signals. In this lecture course we focus on analog modu-
lation, which is usually realized either as frequency modulation (FM) or as
amplitude modulation (AM).

Frequency modulation

In frequency modulation, the output y(t) of the modulator with carrier

frequency νc for an input x(t) is given by
 Z t 
y(t) = Vc cos 2πνc t + 2π∆ν x(τ )dτ , (6.1)
0

where Vc is the amplitude (voltage) of the carrier and ∆ν the frequency

deviation, i.e. the maximum shift with respect to the carrier frequency. The
appearance of an amplitude-modulated signal is shown in Figure 6.2A for
∆ν = 0.05νc . In many practical apllications, the relative frequency deviation
is much smaller.

A B C

Figure 6.2: Modulation types. (A) Carrier (black) and signal (blue). (B)
Frequency-modulated signal (∆ν = 0.05νc ). (C) Amplitude-modulated signal
(black line) with m = 0.8 with envelopes (dotted blue lines) corresponding to the
signal shape.

67
Amplitude modulation

Amplitude modulation is based on multiplication of the carrier oscillation

c(t) = cos(2πνc t) and the signal x(t). To see the consequences, we consider
a signal x(t) = cos(2πνs t). We find

1
c(t)·x(t) = cos(2πνc t) cos(2πνs t) = [cos(2π(νc + νs )t) + cos(2π(νc − νs )t)] .
2
(6.2)
In other words, sidebands are generated at the carrier frequency plus/minus
the signal frequency.
In practice, these sidebands are admixed to a pure carrier oscillation

m m
y(t) = Ac cos(2πνc t) + Ac cos [2π(νc + νs )t] + Ac cos [2π(νc − νs )t] ,
2 2
(6.3)
where m is the modulation index. Distortionless incoherent demodulation is
possible only for 0 < m ≤ 1.
From Figure 6.3 it is clear that one sideband carries the whole information.
With more involved modulation circuits and receivers, double sideband sup-
pressed carrier (DSBSC) and single sideband (SSB) modulation are possible.

carrier
bandwidth

0 ns nc-ns nc nc+ns n

Figure 6.3: Spectrum of a signal with maximum amplitude at frequency νs and

corresponding modulation sidebands centered at the carrier frequency νc .

Incoherent demodulation of an amplitude-modulated signal can be achieved

by rectification, a low-pass filter, and a DC block. This is illustrated in Fig- Gleichrichtung
ure 6.4(A). The amplitude-modulated signal Vi (B) is rectified by a diode to
give Vc (C). A capacitor and resistor in parallel suppress the high-frequency
components to give Vd (D). Finally, the direct-current contribution (offset)
is suppressed by a capacitor, so that the output signal Io is the original
signal x(t) supplied to the mixer in the modulator (E).

68
 A B D

Vi VC VD Vo C E

Figure 6.4: Incoherent demodulation of the amplitude-modulated signal shown in

Figure 6.3C. (A) Demodulation circuit. (B) Input voltage Vi (t) of the amplitude-
modulated signal. (C) Voltage VC after the rectifier (detector diode). (D) Smoothed
voltage VD after the parallel RC low pass filter. The dotted line denotes zero level.
(E) Output current Io after the final capacitance, which blocks the DC component.

Distortions due to non-linear characteristic curves

Most electronic components can be considered as linear systems (y(t) ∝ x(t))

only in a small range. Non-linearity can be characterized by the character-
istic curve, which is the function y = f (x). Characteristic curves are often Kennlinie
represented as the dependence of the output current on the input voltage
Io (Vi ) or as the dependence of the output power on the input power. In
detection circuits, zero input signal is usually biased to the operating point, Arbeitspunkt
which is the center of the most linear part of the characteristic curve (green
point in Figure 6.5). It is not always possible to avoid operation in the non-
linear regime. For instance, operating a high-power amplifier for excitation
well in the linear regime would impose unnecessary limitations on available
excitation power.
If components, such as amplifiers, are operated in the non-linear part of their
characteristic curve (around the red point in Figure 6.5), the signal waveform
is distorted. Such distortion generally leads to generation of overtones of
the input frequency νi . Depending on the type of distortion, even overtones
(second harmonics at frequency 2νi , fourth harmonics, and so on) or odd
overtones (third harmonics at frequency 3νi , fifth harmonics, and so on) may
dominate.
In the context of amplitude modulation and subsequent demodulation this
may introduce components at twice the original signal frequency, which can
be suppressed by a low-pass filter.

69
 Io

Vi

Figure 6.5: Signal distortion by a non-linear characteristic curve. The grey input
signal is a pure cosine wave. Due to saturation in the characetristic curve, the blue
output signal is damped near its maximum amplitude.

6.2 Phase-sensitive detector/ Lock-in amplifier

The terms phase-sensitive detector (PSD) and lock-in amplifier are synony-
mous. A PSD is a coherent demodulator, which needs to be supplied with
a reference signal (the carrier). Knowledge of the carrier frequency allows
to filter very weak signals out of a noisy background by limiting bandwidth
to a very narrow range. The PSD consists of a mixer, which creates the
difference combination of the input and the carrier frequency, and a very
low-pass filter that lets only slow variations of difference combination pass.
The effective bandwidth is determined by the bandwidth of the signal of
interest.
In spectroscopy PSDs are amply used in effect modulation experiments. A
typical example from optical spectroscopy is shown in Figure 6.6. In this
example, the irradiation of a laser is modulated by a chopper. The frequency
of this modulation is set by the chopper control. The chopper control can
thus supply the reference signal,

Vr (t) = Vr,0 cos (ωr t + φ) , (6.4)

to the PSD. The modulated laser irradiation will cause sample fluorescence
that is modulated by the same frequency, whereas the phase may differ.
This fluorescence signal is detected by a photo-detector, such as a photo-
diode or CCD camera. The signal from the sample at the output of the

70
photodetector has the form

Vs (t) = Vs,0 cos (ωr t + φs ) . (6.5)

The mixer multiplies the signal input with the reference signal

Vs,0 Vr,0
Vm (t) = cos (ωr t + φs ) cos (ωr t + φ)
V0
1 Vs,0 Vr,0
= [cos (2ωr t + φs + φ) + cos (φ − φs )] . (6.6)
2 V0

The low-pass filter eliminates the high-frequency term and any noise at
frequencies above the bandwidth of the signal from the sample. This signal
may be slowly oscillating, i.e. φs may be time-dependent at frequencies
much lower than ωr . We find for the output voltage,

1 Vs,0 Vr,0
Vo (t) = cos (φ − φs ) . (6.7)
2 V0

Since V0 is much smaller than Vr,0 the PSD acts as an amplifier. The output
π
signal is maximal for φ = φs and vanishes for φs = φ ± 2. To detect
the maximum signal, the reference signal can be properly phase shifted.
Alternatively, the reference signal can be supplied in quadrature (original
reference signal and reference signal phase shifted by π/2) and the magnitude
of the two quadrature components can be measured.

lock-in amplifier

signal signal
low-pass filter
input output

reference

Figure 6.6: Use of a PSD in an effect modulation experiment in optical spec-

troscopy. Adapted from Wikimedia Commons.

71
6.2.1 Digital PSD
Modern PSDs are based on digital signal processors. The input signal and
reference signal are first digitized and discretized by a sufficiently fast analog-
digital converter (ADC) (Chapter 7). Phase-shifting, multiplication, and
filtering can then all be performed by a digital signal processor (Section
10.1.1), which adds flexibility and allows for realizing very high-order filters.
If necessary, the output signal is converted to an analog signal by a digital
analog converter (DAC). However, in many cases the digital signal is further
processed.

72
Chapter 7

Time-discrete signals

7.1 Analog-digital conversion

The macroscopic world is continuous in time and the signals that are orig-
inally measured in spectroscopy are continuous in amplitude. These sig-
nals are ultimately processed by and stored in digital computers. Hence,
the original time-continuous analog signals must be converted to discreted
amplitude-digitized binary data. This is done by anolog-digital converters
(ADC), which are electronic devices that take samples of the analog input
signal at a given time increment ∆t and convert each of these samples to a
binary number that encodes the signal amplitude.
ADC code

Normalized input signal

Figure 7.1: Encoding of an analog signal to a 3-bit number provides 8 steps

(adapated from Wikimedia commons). Resolution of the signal level corresponds
to the least significant bit (LSB) and is 0.125 in this case.

An ADC is characterized by its vertical resolution, which is specified in bits

and by its sampling rate 1/∆t, which is specified in samples per second

73
(S/s). An n-bit ADC discretizes the signal with a resolution of 2n possible
values over the input range. For instance, a 3-bit ADC encodes the signal
level as one of eight numbers, as visualized in Figure 7.1. At the high end,
compromises must be made between resolution and sampling rate. As of
October 2011, the fastest ADC is a prototype for development applications
with a sampling rate of 56 Gs/s at 8 bit resolution (Fujitsu). At a resolution
of 12 bit, which is sufficient for most applications, devices with sampling
rates up to 3.6 GS/s (National Semiconductor) are commercially available.

7.1.1 Digitization
The vertical resolution determines the round-off error in discretization. Such
discretization corresponds to a quantization of the signal amplitude. As is
apparent from Figure 7.1 and Figure 7.2B, information on the exact shape
of the signal is lost. If the signal is constant, this loss in information cannot
be recovered by averaging. The deviation of the discretized signal from the
true signal is called quantization noise. Quantisierungs-
A different situation arises if the signal is superimposed by white noise with rauschen
an amplitude that exceeds the resolution, i.e. the change in signal level
corresponding to the least significant bit (LSB). Such a signal is shown in
Figure 7.2C. If such a signal is repeatedly sampled and the digital data are
averaged, quantization noise averages along with the white noise (Figure
7.2C). This is because the noise can push the signal above the next step or
below the previous step and the probability for such events depends on the
exact value of the signal.
Accordingly, the vertical resolution of the ADC should be selected such that
in typical measurements the noise exceeds 1 LSB. If the signal in a single
scan has a very large signal-to-noise ratio S/N> 2n , quality of the data is
limited by the ADC resolution. In such cases it may be advantageous to
switch off the first low-noise amplifier in the receiver. This will drastically
reduce S/N, which can be tolerated in such situations, if signal acquisition
can be repeated at a reasonable rate. With the low-noise amplifier switched
off and noise exceeding 1 LSB, signal quality can be improved to the required
level by averaging a sufficient number of scans.

74
 1 1

A B
Normalized signal amplitude

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

1 1

C D
Normalized signal amplitude

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
time (a.u.) time (a.u.)

Figure 7.2: Digitization noise (simulation). Dotted horizontal lines correspond

to ADC output values. (A) Original signal. (B) Signal discretized with 3 bit
resolution. (C) Original signal with noise larger than 1 LSB added. (D) Average
of 1024 discretized signals with different (psseudo-random) noise.

7.1.2 Sampling
The discretized or sampled signal s is a vector of signal values given at
the sampling points. With the sampling rate νs and the sampling period
T = ∆t = 1/νs , sampling is defined by

sk = y(kT ) k = 1 . . . N , (7.1)

where k indexes the points of the discretized signal an tmax = N T is the

total observation time.
Unlike digitization, sampling may not entail any loss of information if the
signal is bandwidth limited, which is usually the case. Specifically, we con-
sider a signal with a maximum frequency νmax . Signal loss is avoided if the

75
sampling period fulfills the sampling theorem,

π 1
T < = . (7.2)
ωmax 2νmax

The frequency νNyq = νs /2 is called Nyquist frequency of the discrete signal

processing system. In other words, the analog signal s can be exactly recon-
structed from the sk if at least two sampling points lie within one period of
the signal component with maximum frequency.
The signal is reconstructed by
∞  
X sin (πνs t − kπ)
k
y(t) = s
πνs t − kπ
νs
k=−∞
∞
sin π Tt − k



X
= sk . (7.3)
π Tt − k


k=−∞

This equation can be derived by expressing s(t) in a Fourier series, but the
derivation is lengthy.

1 1

A B
Normalized signal amplitude

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
time (a.u.) time (a.u.)

Figure 7.3: Sampling and signal reconstruction (simulation). The signal has two
frequency components ν1 = 5 and ν2 = 2.9 and is sampled with a Nyquist frequency
of 8. (A) Original signal (black), sampling times (blue lines), and sampled signal
points sk (blue circles). (B) Original signal (black dashed line), sampled signal
points Sk (blue circles) and signal reconstructed by Eq. (7.3) from the limited
number of sampled points (red line).

Sampling and signal reconstruction are illustrated in Figure 7.3. With fre-
quency components of 2.9 and 5, a Nyquist frequency of 8 is amply sufficient.
Nevertheless, signal reconstruction is not perfect in this case as the sampled
signal is known not from t = −∞ to t = ∞, but only in the range from 0 to

76
1. This leads to deviations between the reconstructed (red line) and original
signal (dashed line).

7.1.3 Aliasing
We now consider the consequence of sampling in frequency domain. The
discretized signal sd is given by
∞
X
sd = y(t)δ(t − nT ) . (7.4)
n=−∞

Hence, spectrum of the discretized signal is a convolution of the true spec-

trum S(ω) with the Fourier transform of the sampling function,

Sd (ω) = S(ω) ∗ DΩ (ω)

∞
1 X 2πk
= S(ω) ∗ δ(ω − )
T T
k=−∞
∞
1 X
= S(ω − kΩ) , (7.5)
T
k=−∞

where we have introduced the angular frequency equivalent of the sampling

frequency Ω = 2π/T and have used Eq. (3.31) to expresss the Fourier
transform DΩ of the sampling function.
The last line of Eq. (7.5) describes periodically repeated copies of the Fourier
transform S(ω) of the continuous signal. This is the reason why the original
signal needs to be bandwidth limited and the Nyquist criterion

1
T < (7.6)
2νmax

must be fulfilled. Only in that case all copies are separated and the Fourier
1 1
transform of the discretized signal in the range (−νNyq , νNyq ) = − 2T , 2T is
exactly equal to S(ω).
If the time step T = ∆t is taken to long, the copies overlap. Frequencies
with absolute values larger than the Nyquist frequency are folded back into
the accessible frequency range and lead to spurious peaks. This phenomenon
is called aliasing and should be avoided. To avoid unpleasant surprises, the
signal can be bandwidth limited by an appropriate low pass filter before
discretization.

77
 0.4
4 a
0.2 3

2
0
1

-0.2
0

0 1 2 3 4 5 0 2 4 6 8 10 12

0.4
4 b
0.2 3

2
0
1

-0.2
0

0 1 2 3 4 5 0 2 4 6 8 10 12

0.4
4 c
0.2 3

0 2

1
-0.2
0
-0.4
0 1 2 3 4 5 0 2 4 6 8 10 12
t (ms) n (MHz)

Figure 7.4: Effect of oversampling on digitization noise (simulation). Left panels

show time-domain data and right panels the corresponding absorption spectra. (A)
Signal with high amplitdue resolution (24 bit ADC) and 128 data points at time
increment ∆T = 0.391 µs corresponding to a Nyquist frequency of 12.8 MHz. (B)
Same time increment, number of data points, and Nyquist frequency as above, but
digitized with a 6 bit ADC. (C) Digitization with a 6 bit ADC, but with four times
oversampling (∆T = 0.0098 µs corresponding to a Nyquist frequency of 51.2 MHz).

7.1.4 Oversampling
If the ADC can sample faster than is required by the bandwidth limitation of
the signal, signal quality can be improved by using a time step ∆t = T that
is shorter than required by νmax . First, such oversampling guards against
aliasing of noise with frequencies larger than νmax . Such noise may exist
even if the signal has passed an analog low pass filter, since suppression
of signal beyond the cutoff frequency by the filter is not perfect. Second,

78
oversampling reduces the influence of quantization noise in the spectrum,
since digitization noise also has frequency components beyond νmax . This
effect is illustrated in Fig. 7.4. It is the reason for the use of oversampling
techniques in music players.

7.2 Fourier Transformation of time-discrete sig-

nals
7.2.1 Computation of a continuous function in frequency do-
main
In Chapter 3 we have considered the Fourier transformation of time-continuous
signals. The time-discrete signal is given as a series

sk = y(kT ) . (7.7)

We first consider the case where the signal is discrete, but known throughout
all history, presence, and future. This signal can be expressed as a time-
continuous function y∗ (t) by multiplication of the time-continuous signal
with a series of Dirac delta functions,
∞
X
y∗ (t) = y(t) · δ(t − kT ) . (7.8)
k=−∞

Using Eq. (3.9), we thus find for the FT (spectrum) of the discontinuous
(time discrete) signal
Z ∞ Z ∞ ∞
X
−iωt
Y∗ (ω) = y∗ (t)e dt = y(t)δ(t − kT )e−iωt dt
−∞ −∞ k=−∞
∞
X Z ∞ ∞
X
−iωt
= y(t)δ(t − kT )e dt = y(kT )e−iωkT
k=−∞ −∞ k=−∞
X∞
= sk e−iωkT . (7.9)
k=−∞

Note that the result is a continuous function in frequency domain.

The normalized, inverse transform from a continuous spectrum to a time-

79
discrete signal is given by
Z π
1
sk = Y∗ (Ω)eiΩkd Ω , (7.10)
2π −π

where Ω = 2π/T is again the angular sampling frequency and integration

is over the Nyquist band. This presupposes that the spectrum is identical
zero outside the Nyquist band. Conceptually, this inverse discrete-time FT
reconstructs the time-continuous signal as in Eq. (7.3) and subsequently
samples this time-continuous signal at points sk .

7.2.2 Properties of discontinuous FT

Several properties of the continuous FT have analogs in discontinuous FT.

Linearity

The DFT is linear,

n o n o n o
Fd cA s(A) (B)
n + cB sn = cA Fd s(A)
n + cB F s(B)
n

= cA SA (ω) + cB SB (ω) . (7.11)

Shift theorem

Fd {sn−n0 } = e−iωT n0 S(ω) . (7.12)

Convolution theorem
n o
(A) (B)
Fd s(A)
n ∗ s (B)
n = Sk Sk , (7.13)

i.e., convolution in time domain corresponds to the elementwise product of

the spectrum vectors in frequency domain.
Likewise,
Z π/T
n
(A) (B)
o T
Fd sn sn = S (A) (eiω̃t )S (B) (ei(ω−ω̃)t )dω̃ . (7.14)
2π −π/t

Power theorem
∞ Z π/T
X T 2
|sn |2 = S(eiωT ) dω . (7.15)
n=−∞
2π −π/T

80
7.2.3 Discrete FT
In most cases one is interested in a spectrum that is also discrete, as the
spectrum is further processed and analyzed by a digital computer. To find
out how the spectrum should be discretized, we note that in any practical
situation only a finite number of the sk in Eq. (7.9) is available. In par-
ticular, we do not know the signal at negative times k < 0 and we don’t
know the signal after the end of our observation time tmax = N T . For the
observation frequency or frequency increment, we define

1 1 νs
∆ν = = = . (7.16)
tmax NT N

In other words, the spectrum has the same number of points as the time-
domain signal.
This discrete spectrum can be obtained by the Discrete Fourier transfor-
mation (DFT), which results from defining a normalized discrete frequency
domain with νn = n∆ν and

2πνn 2πn∆ν 2πnT 2πn

Ωn = = = = . (7.17)
1/T 1/T NT N

The sum over the Ωn substitutes for the integral over ω in Eq. (7.9). This
normalization of angular frequency domain has the advantage that the com-
putation algorithm becomes independent of the actual time axis. Time
vanishes from the equation and is replaced by the running index,

N
X −1
Sn = sk e−i2πkn/N n = 0, 1, . . . N − 1 . (7.18)
k=0

DFT in this way is implemented in Matlab as a function fft.

Definition of the running index k in that way implies that the frequency
band starts at zero and is thus no longer symmetric about zero. Because of
periodicity this does not create a problem, as the discrete data can simply
be cyclically shifted. In Matlab this is done by the function fftshift.

Remarks:

• According to Eq (3.4), the first coefficient in a Fourier series must be

divided by two, as zero frequency is effectively counted twice. Eq.
(7.18) does not account for this. This adds a constant offset to the

81
 discrete spectrum, so that the integral of the spectrum is twice as large
as it should be. To avoid this, s0 needs to be divided by two before
DFT.

• The frequency axis is constructed as follows. A preliminary axis from

zero to (N − 1)/(N T ) is generated. This axis is subjected to the same
cyclic shift as the spectrum. Finally, the period in frequency domain
of 1/T is subtracted from all frequency values νk ≥ 1/(2 ∗ T ).

• The cyclic shift is performed as follows. For an odd number 2l + 1,

the number of shifted points is l + 1. For an even number 2l it is l.
This number of points is removed from the beginning of the vector and
appended at its end.

• If the input data sk are real, the spectrum is an even function. In that
case no shift is required. The spectrum and frequency axis can simply
be cut after l points for N = 2l and after l + 1 points for N = 2l + 1.

• Linearity, shift theorem, and convolution theorem of discontinuous

FT can be applied to DFT without further considerations. The power
P −1 PN −1
theorem assumes the form N 2
n=0 |sn | = N
2
k=0 |Sk | .

The Matlab function dft, programmed for this lecture course and available
on the lecture course homepage, computes the spectrum and frequency axis
from a given time-domain signal and time axis, taking these considerations
into account.
DFT can be written as a matrix equation

~ = W~s ,
S (7.19)

~ and ~s are the vectors of points Sk and sn , respectively, and W is

where S
given by
Wkn = e−i2πkn/N . (7.20)

For the inverse DFT,

1
~s = W−1 S
~ , (7.21)
N
where the inverse transformation matrix is given by

−1
Wkn = ei2πkn/N . (7.22)

82
With W = e−i2π/N , the matrix assumes the form
 
1 1 1 1 ...
 1 W W2 W3 ...
 

 
W=
 1 W2 W4 W6  .
...  (7.23)
 ... ... ... ... ... 
 

1 W N −1 W 2(N −1) ... ...

This formulation of the problem shows that computation can be set up

efficiently by recursive multiplications with W .

7.2.4 Zero filling

As it is defined, DFT does not make use of the full information content of
the signal. Let us consider a purely real signal sn with n = 0 . . . N − 1.
This signal has N real parameters. The DFT provides a complex-valued
spectrum with N points. However, the imaginary part is related to the real
part by the unique Hilbert transformation and thus does not contain any
information different from the one in the real part. Furthermore, for a purely
real signal the spectrum is an even function of frequency. This means that
half of the N points, namely the ones corresponding to negative frequencies,
do not contain any information beyond the one contained in the other half
of the points. Hence, from N real numbers we have reduced information to
N/2 real numbers. Resolution of the spectrum is worse than it should be.
An analogous argument holds for a complexed-valued signal. Now the spec-
trum is an odd function of frequency and does contain N real numbers
information. However, the input data are N pairs of numbers. Again, the
resolution of the spectrum is worse than it should be.
This problem can be solved by zero filling. For this, N zeros are appended
to the time-domain signal before DFT. The discrete spectrum then has 2N
1
values and a frequency resolution ∆ν = 2tmax . As we have added no infor-
mation, we just compute spectral amplitudes at intermediate points by Eq.
(7.18). According to Eq. (7.3) this corresponds to a sinc interpolation of
the existing points. The interpolated values correspond to a true resolution
enhancement, i.e. the better resolution is realistic.
Zero filling beyond 2N inserts further intermediate points, which are also
computed by sinc interpolation. However, this is purely cosmetic. The
spectral amplitudes at these points do not reconstruct the true spectral

83
amplitudes. For certain display and analysis purposes extended zero filling
can make sense. However, on should always be aware that the true resolution
1
si not better than ∆ν = 2tmax .
By default, the Matlab function dft applies zero filling to 2N . This behavior
can be changed by specifying in options.zerofilling either no zero filling
(options.zerofilling=0) or extended zero-filling (options.zerofilling
= M, where M specifies the number of points after zero filling and must be
larger than N ).
Reduction in number of multiplications

1200

1000

800

600

400

200

0
0 5000 10000 15000
N

Figure 7.5: Reduction factor of the number of multiplications between FFT and
DFT. Already for n = 16384 points, a factor of more than 1000 is saved in compu-
tation time.

7.2.5 Fast Fourier transformation

The generality of the FT made DFT one of the most important problems
that were solved in the early times of digital computation. Since computa-
tion speed was by several orders of magnitude slower than today, efficient
algorithms were required. A particularly efficient algorithm was found by
Cooley and Tukey1 in 1965 for the case where N is a power of 2. The
problem then becomes symmetric and one can start with a 2-point DFT,
compute the 4-point DFT as a weighted combination of two-points DFTs
and so on. The algorithm scales much more favorably with an increase in
N . Whereas a trivial algorithm according to Eq. (7.18) requires N 2 mul-
tiplications, the fast Fourier transformation (FFT) algorithm requires only
N log2 N multiplications. The computation time advantage is illustrated in
Fig. 7.5.
1
Actually they reinvented this. Gauss had the idea already in 1805.

84
This advantage is the main reason why much of digital data processing is
based on data sets whose dimensions are powers of 2. Despite the increase in
computer efficiency, the advantage of the FFT algorithm is still important,
in particular in real-time data processing and for multi-dimensional FT.

7.3 Apodization
7.3.1 Effects of limited observation time
Limited resolution

The limited acquisition time leads to limited resolution. This effect can be
understood semi-quantitatively from Eq. (3.35). If a signal exhibits Gaus-
sian decay with time-domain variance σt2 = 1/(2a), the spectrum will exhibit
an angular frequency-domain variance σω2 = 2a and thus a frequency-domain
variance σν2 = a/(2π 2 ). Hence, the product of the standard deviations is
σt σν = 1/(2π) and the uncertainty of the frequency is σν = 1/(2πσt ).
If the signal is not decayed at maximum observation time tmax , no theoretical
resolution limit is set by the observation time. For instance, a signal that is a
linear combination of a small number n of exponentially damped harmonics
is fully characterized by 2n parameters, namely the n frequencies νij and
n decay time constants T2,ij . If a noise-free signal of that type is sampled
according to the Nyquist criterion with at least N = 2n points, the full
information is available. In practice, a few more points need to be acquired
to reduce influence of noise. In that case frequencies and relaxation times
can be determined without an uncertainty that is related to tmax , if data
analysis techniques other than FT are applied (see Section 9.2). However,
such approaches break down if no theoretical model for the time-domain
signal can be assumed or if n becomes too large.

Truncation artifacts

Formally, truncation of the signal at tmax can be described by multiplication

with a box function,
 
t
ytrunc (t) = y(t) · Π , (7.24)
2tmax

85
 1 1 ms
a 50
b
0.5 40

30
0
20

-0.5 10

0
-1
0 2 4 6 8 -20 -10 0 10 20
t (ms) n (MHz)
20
30
c d
25
15
20

15 10

10
5
5

0 0

-20 -10 0 10 20 -20 -10 0 10 20

n (MHz) n (MHz)

Figure 7.6: Truncation artifacts and apodization. (a) Real part of a signal y(t) =
e(i2πν−1/T2 )t with ν = 5 MHz and T2 = 1 µs. The vertical dashed line indicates
truncation at tmax = 1 µs. The red curve is a Hamming apodization function.
(b) Spectrum obtained by DFT after zero-filling from the complete signal with
tmax = 8 µs. (c) Spectrum obtained by DFT of the signal truncated at tmax = 1 µs
after zero-filling (black) and without zero-filling (green). (d) Spectrum obtained
from the signal truncated at tmax = 1 µs by multiplication with the Hamming
window (black dashed) or a Dolph-Chebyshev window with 50 dB ripple suppresion
(red), zero-filling, and subsequent DFT.

which corresponds to
(
y(t) for t ≤ tmax
ytrunc (t) = . (7.25)
0 for t > tmax

As the FT of the box function is a sinc function, this corresponds to convo-

lution of the spectrum with the corresponding sinc function

sin (ωtmax )
Ytrunc (ω) = Y (ω) ∗ 2tmax . (7.26)
ωtmax

This convolution causes ripples in the spectrum, as can be appreciated from

86
Figure 7.6C. Although the ripples are not apparent if DFT is applied with-
out zero filling, the spectral line shape is still distorted. In Figure 7.6C
this distortion is apparent as an asymmetry and some negative intensity in
the high-frequency wing of the line in the green spectrum. Such trunca-
tion artifacts for large signals can conceal small signals and should thus be
avoided.

7.3.2 Apodization window

Dolph-Chebyshev window

Distortion by truncation artifacts can be minimized by multiplying the sig-

nal with a window function that decays smoothly towards zero in time do-
main. In such apodization the spectrum is convoluted with the FT of the
window function, which does not contain pronounced oscillations. However,
this convolution causes signal broadening beyond the digital resolution of
the spectrum obtained by DFT. Generally, better ripple suppression is re-
lated to more broadening and vice versa. For a required attenuation of the
ripples compared to the main signal, minimum broadening is obtained with
the Dolph-Chebyshev window. The Dolph-Chebyshev window function is
analytically defined in frequency domain

cos M arccos β cos πk

 

M
WDC (ωk ) = , k = 0, 1, 2, . . . , M − 1 , (7.27)
cosh [M arccosh (β)]

with  
1 α
β = cosh arccosh (10 ) , (7.28)
M
where the ripple level in dB is given by −20α. The required window func-
tion in time domain w(t) can be obained by inverse DFT. In Matlab, the
Signal Processing Toolbox provides a function chebwin that computes the
Dolph-Chebyshev window in time domain for given number of points and
ripple attenuation. The Matlab function dft.m for DFT with pre-processing,
which was written for this lecture course, allows for using Dolph-Chebyshev
apodization if the Signal Processing Toolbox is available. Ripple suppression
in dB (as a positive number) can be supplied in variable options.ripple.
As is seen in Figure 7.6D, the lineshape with Dolph-Chebyshev apodization
(red) with 50 dB ripple suppression (the default in dft.m) is symmetric and
has low ripples. Line broadening is moderate, as is apparent from compari-

87
son with Figure 7.6B.

Hamming window

Since the Dolph-Chebyshev window in time domain cannot be given in closed

form, many other apodization windows have been tried. A good all-purpose
compromise between ripple suppression and broadening is obtained with the
Hamming window
 
πt
w(t) = 0.54 + 0.46 cos . (7.29)
tmax

As is seen by comparing the dashed black line (Hamming window) with the
red line (Dolph-Chebyshev window) in Figure 7.6D, the Hamming window
is a good, simple substitute for a Dolph-Chebyshev window in many situa-
tions. In dft.m the Hamming window is implemented independently from
the Signal Processing Toolbox.

Kaiser window

The Hamming window does not allow for adjusting the compromise between
ripple suppression and line broadening. Such adjustment is possible with the
Kaiser window h p i
I0 πα 1 − t2 /t2max
w(t) = , (7.30)
I0 (πα)
where I0 is the zero-order modified Bessel function of the first kind and α de-
termines the trade-off between ripple suppression and broadening. In dft.m
the Kaiser window is implemented independently from the Signal Processing
Toolbox with a default α = 2 which provides a similar compromise between
ripple and line broadening as the Hamming window. The parameter α can
be supplied in variable options.alpha.

7.3.3 Filtering for optimum signal-to-noise

Signal truncation

Window functions can also be applied as filters. If an experimental signal

Y (ω) = Y0 (ω) + N (ω) is the sum of the wanted signal (sample response)
Y0 (ω) and noise N (ω), the wanted signal is affected in a different way by
multiplication with a window function than the noise. This is immediately

88
 1

a 50
b
0.5 40

30
0
20

-0.5
10

0
-1
0 2 4 6 8 -20 -10 0 10 20
t (ms) n (MHz)

50
c 25
1 d
40 20 0.5
0
30 15 -0.5
-1
20 10 0 2 4 6 8
t (ms)

10 5

0 0

-20 -10 0 10 20 -20 -10 0 10 20

n (MHz) n (MHz)

Figure 7.7: Signal-to-noise ratio improvement by time-domain filtering (simula-

tion). (a) Original time-domain data (blue line) of a signal y(t) = e(i2πν−1/T2 )t with
ν = 5 MHz and T2 = 1 µs superimposed by simulated noise. For noise simulation
normally distributed random numbers with standard deviation 0.1 were added. The
box window for truncation (dashed black line) and the matched filter xe (t) = e−t/T2
(red line) are also shown. (b) Spectrum obtained by DFT after zero-filling from
the original signal. The signal-to-noise ratio is 32.52 . (c) Spectrum obtained by
DFT of the signal truncated at tmax = 4 µs . The signal-to-noise ratio is 47.32 . (d)
Signal multiplied with the matched filter xe (t) tmax = 1 µs (blue line in the in set)
and spectrum obtained from this signal (black line). The signal-to-noise ratio of
the spectrum is 73.42 .

clear for a case where the wanted signal has decayed below the 
noise level
 at
t
time tdecay < tmax . Multiplication of y(t) by a box function Π 2tdecay will
set the time-domain signal to zero at times t > tdecay and will not affect the
signal at t ≤ tdecay . This operation will suppress more noise than signal and
will thus increase signal-to-noise ratio. In frequency domain, such deliberate
truncation of the signal corresponds to smoothing by convolution with a sinc
function. The effect of such multiplication with a box function is illustrated
in Figure 7.7C.

89
Matched filter

We now consider that the signal decays with an envelope function ye (t). For
damped harmonic oscillations yk (t) = exp(iωk t) exp(−t/T2,k ) this envelope
function is ye (t) = exp(−t/T2,k ). Since noise level is time independent, the
lokal signal-to-noise ratio at time t is proportional to ye (t). By multiplying
the signal once more with ye (t), each point of the time domain signal is
weighted with its own sensitivity (local signal-to-noise ratio), so that points
with larger signal-to-noise ratio dominate in the time-domain signal. The
window function
wmatched (t) = ye (t) (7.31)

is a matched filter. In general,

L(ω) = F {ye (t)} (7.32)

is the lineshape of the spectral lines. If L(ω) is known approximately or

can be inferred from the data, the signal can be detected with optimum
sensitivity by multiplying time-domain data with ye (t).
The effect of a matched filter is illustrated in Figure 7.7D. Note the broaden-
ing of the line by a factor of two and the corresponding decrease in amplitude
by a factor of two.

Box 2: Signal-to-noise ratio in discrete spectra

The signal-to-noise ratio of a peak in a numerically given discrete spectrum can
be estimated in the following way. The signal amplitude AS is approximately
the maximum amplitude of the peak. This approximation is poor for very poor
signal-to-noise ratio (SNR < 10). In such a case it is advisable to determine the
signal amplitude by fitting of the peak with an appropriate lineshape function.
The noise amplitude is defined as the root mean square amplitude of the spec-
trum in a baseline region. Assume that points k1 to k2 in the spectrum corre-
spond to baseline. Then, the noise amplitude is given by
rP
k2 2
k=k1 Yk
AN = k2 −k1 +1

Strictly speaking, signal-to-noise ratio is defined in terms of power, see Eq.

(5.1). In spectroscopy, the amplitude ratio AS /AN is often loosely termed
signal-to-noise ratio. The numbers in the caption of Figure 7.7 are thus given
as x2 , where x is the amplitude ratio.

90
7.3.4 Filtering for optimum resolution
Fourier self deconvolution

Apodization with a decaying window function, such as a matched filter

causes signal broadening. Vice versa multiplication with a rising function
in time domain causes signal narrowing. If the envelope function ye (t) is
known, the time-domain signal corresponding to infinitely narrow peaks
(Dirac delta functions) can be obtained by dividing the original signal y(t)
(or the inverse FT of the spectrum Y (ω)) by ye (t). The time-domain signal
obtained by such Fourier self deconvolution corresponds to a spectrum with
extremely poor signal-to-noise ratio, as noise is enhanced by the deconvo-
lution. Therefore, controlled line broadening is reintroduced by multiplying
with the inverse FT of the new, desired lineshape function Yshape . The
complete procedure of deconvolution line narrowing is thus described by

F −1 {Yshape (ω)} F −1 {Y (ω)}

 
Ynarrowed = F . (7.33)
ye (t)

a b

Figure 7.8: Effect of Lorentz-Gauss transformation on peak shapes in 2D FT

spectroscopy. The levels for contour plots are 2.5, 5, 10, 20, 40, and 80% of max-
imum intensity. (a) Pure absorption Lorentzian peaks. (b) After Lorentz-Gauss
transformation in both dimensions with conserved linewidth at half height.

Lorentz-Gauss transformation

In practice, it is often sufficient to change the lineshape without actually

narrowing the lines. Exponentially decaying time-domain signals correspond
to Lorentz line shapes in frequency domain, which have broad wings. This
is particularly damaging in 2D FT spectroscopy, where the crossing of the

91
wings can lead to peak-like features (Figure 7.8a) that are hard to distinguish
from true, weaker peaks. The following window function achieves a Lorentz-
Gauss transformation
2 t2 /2
wLG (t) = et/τ −σ . (7.34)

This transforms a Lorentzian line with frequency-domain linewidth 1/πτ

√
(full width at half height) to a Gaussian line with linewidth 2 ln 2/πσ.
Note also that τ = T2 for purely Lorentzian lines. In general, σ can be
chosen freely, thus trading resolution for signal-to-noise ratio. Linewidth at
√
half height is conserved for σ = 1/(τ 2 ln 2). Such a linewidth-conserving
Lorentz-Gauss transformation strongly improves the appearance of 2D spec-
tra, as can be appreciated from Figure 7.8b.
For resolution enhancement, σ = 1/(2τ ), corresponding to line narrowing
by almost a factor of two may often be a good choice. If no good estimate
for τ is known beforehand, one can fix σ at this ratio, start with τ =
tmax /3 and increase or decrease τ until the desired compromise between
resolution and linewidth is achieved. Examples for the effect of Lorentz-
Gauss transformations in 1D spectra are given in Figure 7.9.

Sinebell apodization

The phase-shifted sinebell window function is given by

 
πt
wSB (t) = sin +φ (7.35)
tmax

For φ = 0 the first point of the time-domain signal is set to zero, which im-
plies zero integral of the spectrum and thus negative ”line feet”. This effect
is reduced by introducing the phase shift φ. Larger phase shifts correspond
to less resolution enhancement, less reduction of the signal-to-noise ratio and
less negative intensity. Examples for the effect of sinebell transformations
with zero and 10◦ phase shift are given in Figure 7.9b,c.

7.4 Network analyzers and spectrum analyzers

Electronic instruments or circuits can be characterized by measuring their
frequency response. A network analyzer measures the scattering param-
eters (S-parameters) of the system, namely reflection and transmission as
a function of the frequency of a test signal (Figure 7.10). Depending on

92
 a d

b e

c f
-20 -10 0 10 20 -20 -10 0 10 20
n (Hz) n (Hz)

Figure 7.9: Resolution enhancement in 1D spectra (simulation). The spectrum

consists of a peak with relative intensity 5 at a frequency of -8 Hz and a qintuplett
with 2 Hz splitting and intensity ratio of 1:4:6:4:1 centered at 4 Hz. The uniform
linewidth is 2.1 Hz. Normally distributed random numbers with amplitude 0.001
werde added to the time-domain data to simulate noise. (a) Original spectrum
with Lorentzian peaks. (b) Sinebell apodization with zero phase shift. (c) Sinebell
apodization with 10◦ phase shift. (d) Lineshape-preserving Lorentz-Gauss trans-
formation. (e) Lorentz-Gauss transformation that reduces linewidth to 3/4 of the
original value. (f) Lorentz-Gauss transformation that reduces linewidth to 1/2 of
the original value.

the system, the transmitted wave can be attenuated or amplified and phase
shifted. Most modern network analyzers are vector analyzers, which mea-
sure the complex signal and can thus characterize not only amplitude, but
also the phase shift and in many cases the group delay. Network analyzers
provide the test signal and measure the response.
Spectrum analyzers characterize time-domain signals in terms of their spec-
trum. Analog spectrum analyzers use the super heterodyne principle. A ref-

93
erence frequency is scanned linearly by supplying a sawtooth control voltage
to a local oscillator. The signal to be characterized is mixed with the ref-
erence signal of variable frequency and then filtered by a narrow bandpass.
The sawtooth control voltage provides the horizontal axis of an oscilloscope
picture, while the amplitude of the filtered signal provides the vertical axis.
Nowadays digital spectrum analyzers are more common. In these devices,
the signal is subjected to DFT, using the FFT algorithm. Spectrum ana-
lyzers for bandwidth-limited high frequency signals are hybrid FFT-super
heterodyne analyzers, which use downconversion with a reference signal to
allow for longer sampling intervals. This makes FFT more efficient, as it
limits the number of points. Real-time spectrum analyzers process different
overlapping parts of the incoming signal in parallel and superimpose the re-
sulting spectra. Modern digital oscilloscopes usually have a simple spectrum
analyzer function without downconversion.

Test signal
System Transmission
Reflection

Figure 7.10: Scheme of measurements on a system with a network analyzer.

94
Chapter 8

Stochastic signals

8.1 Characteristics of random variables

Many physical processes are not exactly predictable. Examples are Brown-
ian motion and thermal noise, or the time between absorption and emission
of a single chromophore. Yet a signal arising from such processes contains
information, in the three mentioned cases about the diffusion constant, the
temperature, and about the lifetime of the excited state of the chromophore.
More detailed studies can reveal additional information, for instance the law
of blackbody radiation in the case of thermal noise and the presence of sev-
eral excited states with different lifetimes in the case of single chromophores.

80

60
t0 a b
2
12
40
amplitude (a.u.)

log (M)

1.5 8
20

0 4
1
-20 0

-40 0.5 -12 -8 -4 0

log(n)
-60
0
-80
0 5000 10000 15000 -0.01 0 0.01 0.02 0.03 0.04 0.05
time (a.u.) frequency (a.u.)

Figure 8.1: Stochastic signal (Brownian noise). (a) The time trace was simulated
by integrating white noise, which was in turn simulated by normally distributed
random numbers (commands randn and cumsum in Matlab). (b) Magnitude spec-
trum of the stochastic signal obtained by DFT. The inset is a double logarithmic
plot.

In this Chapter we discuss approaches for extracting information from such

95
stochastic signals. As an example, consider Brownian noise, which is the
integral of white noise (Figure 8.1a). As white noise has a constant mean
square amplitude as a function of frequency, Brownian noise has a 1/ω 2
dependence of the mean square amplitude (power spectrum). This can be
recognized by FT (Figure 8.1b) and can be seen particularly well in a dou-
ble logarithmic plot (inset). In the case at hand, the slope in the double
logarithmic plot of the magnitude spectrum is -1, as the amplitude has a
1/ω dependence if the power has a 1/ω 2 dependence. Deviations at very
low frequencies are due to the limited observation time.
If we measure the same random variable several times at time t0 after ini-
tializing (triggering) our experiment, we will obtain different values (real-
izations) x1 (t0 ), x2 (t0 ), x3 (t0 ) . . . of the random variable. If we measure (or
simulate) the Brownian noise trace several times, it will look different each
time. However, the 1/ω dependence of the amplitude will be a stable charac-
teristic. In the following, we summarize the most important characteristics
of random variables as they are used in probabilistic theory.

8.1.1 Probability distribution

In general, x can assume arbitrary values −∞ < x < ∞. However, the
values in this range are not equally probable. The probability distribution
function
F (x) = P (xm ≤ x) , (8.1)

charcaterizes the probability that any particular measured value xm is smaller

than a threshold value x. The function is nomalized and monotonically in-
creasing with

F (−∞) = 0
F (∞) = 1 . (8.2)

If the random variable is discrete, F (x) is a step function. As an example,

F (x) for an ideal dice is plotted in Figure 8.2a.

96
 1
F(x)
a b

p(x)
6

0 1 2 3 4 5 6 0 1 2 3 4 5 6
x x

Figure 8.2: Probability distribution function F (x) (a) and probability density
function p(x) (b) for an ideal dice.

8.1.2 Probability density

In physics and physical chemistry, the probability density function p(x) is
more often encountered. It is defined as

dF (x)
p(x) = . (8.3)
dx

The probability density function for an ideal dice is plotted in Figure 8.2b.
The intergal of the probability density function is unity,
Z ∞
p(x)dx = 1 . (8.4)
−∞

In many situations, we wish to know the probability P that a realization xm

of a random variable x falls in an interval (a, b). This is given by
Z b
P (a < xm < b) = p(x)dx = F (b) − F (a) . (8.5)
a

8.1.3 Mean value, variance, standard deviation, moment anal-

ysis
For any given function f (x), the k th moment is defined as
Z ∞
Mk = xk f (x)dx . (8.6)
−∞

97
The zeroth moment of a function f (x) is thus
Z ∞
M0 = f (x)dx . (8.7)
−∞

By normalization, the zeroth moment of the probability density function is

unity.
The expectation value or mean value of the random variable is
Z ∞
hxi = x̄ = M1 = xp(x)dx . (8.8)
−∞

In moment analysis, hxi is the first moment M1 ). For noise, we have hxi = 0.
The width of the distribution is characterized by the variance σ 2 , which is
the square of the standard deviation σ. This is computed as
Z ∞
2
σ = (x − hxi) p(x)dx . (8.9)
−∞

The variance is related to the second moment by

σ 2 = M2 − M12 . (8.10)

The variance is also the second central moment, with central moments µk
being defined as
Z ∞
µk = (x − hxi)k f (x)dx k = 2, 3, . . . . (8.11)
−∞

Asymmetry of the probability density can be characterized by the skewness Schiefe

γ1 , which is related to the third central moment

∞  3
x − hxi
Z
µ3
γ1 = p(x)dx = . (8.12)
−∞ σ σ3

8.1.4 Gaussian distribution

The central limit theorem of probability theory states that the mean of
a sufficiently large number of uncorrelated, identically distributed random
variables conforms to a Gaussian distribution, also called normal distribu-
tion, " #
1 (x − hxi)2
p(x) = √ exp − (8.13)
2πσ 2σ 2

98
The condition of identical distributions of the random variables can be sub-
stituted by several other, weaker assumptions. This situation is often en-
countered in physical sciences. White noise, for instance thermal noise at
low frequencies, conforms to a normal distribution, as does the error of a
measurement that is performed repeatedly under identical conditions.
Note however that in many cases where normal distributions are assumed,
the conditions of the central limit theorem are not fulfilled or at least are
not known to be fulfilled. Whenever possible, the hypothesis of a normal
distribution should be tested before it is used in data analysis. For instance,
the skewness of the normal distribution is zero. Analysis in terms of a normal
distribution will thus conceal any asymmetry of the probability distribution
with respect to the mean value.
The Gaussian function cannot be integrated in closed form. If the proba-
bility density function is a Gaussian function, the probability distribution
function is related to the error function
(x−hxi)/σ   
x − hxi
Z
1 2 1
F (x) = √ e−t dt = 1 + erf √ , (8.14)
2π −∞ 2 2σ

with Z y
2 2
erf(y) = √ e−t dt . (8.15)
π 0

In Matlab, the error function is numerically available as erf.

8.1.5 Poisson distribution

For a small number of discrete observations, the Poisson distribution is more
appropriate than the normal distribution. For instance, consider the prob-
ability p(k, r, t) of recording k counts in a time interval t if the average
counting rate is r (in counts per second). This probability is given by

e−rt (rt)k
p(k, r, t) = . (8.16)
k!

This distribution can be generalized by substituting λ = rt for the expected

number of occurences in a given interval giving

e−λ λk
p(k, λ) = (8.17)
k!

99
with the property
∞
X
p(k) = 1 . (8.18)
k=0

For large values of λ, the Poisson distribution tends to a normal distribution

of the form !
1 (k − λ)2
p(k, λ) ≈ √ exp − . (8.19)
2πλ 2λ

1
a 0.12 b
0.8

0.08
Sp(k)

0.6
p(k)

0.4
0.04
0.2

0 0

0 5 10 15 20 25 0 5 10 15 20 25
k k

Figure 8.3: Poisson distribution with λ = 10. (a) Sum of all discrete probabilities
from p(0) to p(k), corresponding to the probability distribution F (x) of a continuous
variable. (b) Discrete probability density function p(k) (black dots) and Gaussian
function corresponding to Eq. (8.19) (blue line).

8.1.6 Two-dimensional autocorrelation function

If we repeatedly measure a signal x(t) in a time interval, we can characterize
this signal by a two-dimensional autocorrelation function

Rxx (t1 , t2 ) = hx(t1 )x(t2 )i . (8.20)

This function tells us whether the signal at time t2 is correlated to the one
at time t1 (positive value), anti-correlated (negative value), or uncorrelated
(zero). The amplitude of the function quantifies the strength of the correla-
tion.

100
8.2 Correlation functions of stationary processes
8.2.1 Stationary processes
We now assume that we repeatedly measure the random variable x at a time
t0 after initializing the experiment. We shall find a probability distribution
F (x) and probability density p(x) at this given time t0 . If these functions
are independent of our choice of t0 , i.e., invariant under time translation,
the process is stationary. For instance, the process that creates thermal
noise is stationary if temperature was constant for a sufficiently long time.
Processes may be stationary at some times during our experiment, but not
at other times.
Stationary processes have time-independent mean value and variance. The
autocorrelation function depends only on the difference τ = |t1 − t2 | of the
two observation times. Therefore, the autocorrelation function is a one-
dimensional function Rxx (τ ) and it suffices to analyze it for τ ≥ 0.

8.2.2 Ergodicity theorem

Spectroscopic experiments directly provide mean values for an ensemble of
systems, if the sample contains many identical systems. This is not always
the case. In particular, in single-molecule spectroscopy ensemble mean val-
ues and variances could only be obtained by repeating the same experiment
on a large number of molecules. It may be much more convenient to observe
the same molecule over a longer time.
This raises the question whether the time average for a single particle is
the same as the average at a given time over an ensemble of particles. A
stationary stochastic process ist called ergodic if this is the case,
Z T
1
hx(t)i = lim x(t)dt = x(t) . (8.21)
T →∞ 2T −T

Precisely, the ergodic hypothesis originates from statistical thermodynamics

and states that all accessible microstates of a system are equiprobable over a
long period of time. The ergodic hypothesis is not fulfilled if parts of phase
space are kinetically inaccessible. This happens, for instance, in glasses and
solidified polymer melts. In the following, we consider ergodic processes.

101
8.2.3 Autocorrelation function
The autocorrelation function
Z T
1
Rxx (τ ) = hx(t)x(t + τ )i = lim x(t)x(t + τ )dt = x(t)x(t + τ )
T →∞ 2T −T
(8.22)
of an ergodic process is an even function,

Rxx (−τ ) = Rxx (τ ) . (8.23)

∗ (τ ).
If x(t) is complex, Rxx (τ ) is hermitian, Rxx (−τ ) = Rxx
At t = 0 the autocorrelation functions assumes the mean square value of
x(t),
Z T
2 1
Rxx (0) = hx i = lim x2 (t)dt . (8.24)
T →∞ 2T −T

In particular, for signal voltages x(t) = V (t), the value of the autocorrelation
function at t = 0 is proportional to the average power of the signal,
Z T
1 1
P̄ = lim V 2 (t)dt . (8.25)
R T →∞ T 0

This value at t = 0 is the maximum absolute value of Rxx ,

Rxx (0) ≥ |Rxx (τ )| . (8.26)

An autocorrelation function can be (partially) characterized by the correla-

tion time Z ∞
1
τc = Rxx (τ )dτ . (8.27)
Rxx (0) 0

In other words, a box function with amplitude Rxx (0) extending from τ = 0
to τ = τc has the same area as the autocorrelation function. If the auto-
correlation function is an exponential decay, τc is the time constant of this
decay.
The power spectral density (power spectrum) of a signal is the FT of its
autocorrelation function (Wiener-Khinchin theorem)

Sxx (ω) = |F {x(t)}|2 = F {Rxx (t)} . (8.28)

Likewise, the autocorrelation function can be computed from the power

102
spectra density,
Rxx (t) = F −1 {Sxx (ω)} . (8.29)

Since white noise has uniform spectral density, Eq. (8.29) predicts that the
autocorrelation function of white noise has a peak at τ = 0 and is identical
zero everywhere else. For a finite-time signal, the autocorrelation function
looks noisy for τ > 0.
For a time-discrete signal y(tk ), the autocorrelation function can be effi-
ciently computed by DFT to obtain Y (ωl ), computing the product Z(ωl ) =
Y (ωl )Y ∗ (ωl ) and subjecting this product to inverse DFT. This algorithm
is based on the convolution theorem and the fact that the autocorrelation
function is the convolution of the signal with itself.

8.2.4 Cross correlation functions

For two ergodic, and thus stationary, stochastic processes x(t) and y(t), the
cross correlation function Rxy (t) is defined by
Z T
1
Rxy (τ ) = hx(t)y(t + τ )i = lim x(t)y(t + τ )dt . (8.30)
T →∞ 2T −T

Note that this is again a convolution integral, which can be computed ef-
ficiently for time-discrete signals by DFT of both functions to frequency
domain, multiplication of one spectrum with the complex conjugate of the
other spectrum, and inverse DFT. The cross correlation function fulfils

∗
Rxy (τ ) = Ryx (−τ ) . (8.31)

Cross correlation functions are useful for instance for characterizing the cor-
relation between input and output of a device or circuit. As a particu-
larly simple example, we consider the voltage divider shown in Figure 4.1
with x(t) = Vin and y(t) = Vout with purely real impedance Z1 = R1
and Z2 = R2 . We assume an autocorrelation function of the input signal
Rxx (τ ) = σ 2 e−k|τ | . Since

R2
y(t) = x(t) , (8.32)
R1 + R2

we have
R22
y(t)y(t + τ ) = x(t)x(t + τ ) , (8.33)
(R1 + R2 )2

103
so that the autocorrelation function of y(t) is

R22 R22
Ryy (τ ) = Rxx (τ ) = σ 2 e−k|τ | . (8.34)
(R1 + R2 )2 (R1 + R2 )2

Analogously we find for the cross correlation function

R2 R2
Rxy (τ ) = Rxx (τ ) = σ 2 e−k|τ | . (8.35)
R1 + R2 R1 + R2

In this computation we have neglected that the resistors add noise to the
signal. Thus, the result is only valid for the case where the input signal is
much larger than the added thermal noise.

8.2.5 Correlation functions of periodic signals

We first consider the function

x(t) = A cos(ω0 t) . (8.36)

By trigonometric relations, we have

1 1
x(t)x(t + τ ) = A2 cos(ω0 τ ) + A2 cos(2ω0 t + ω0 τ ) . (8.37)
2 2

Because of periodicity, it suffices to integrate over one period T = 2π/ω0 .

The integral of the second term on the right-hand side of Eq. (8.37) vanishes
and because of the normalization we find

A2
Rxx (τ ) = cos(ω0 τ ) . (8.38)
2

For
y(t) = B sin(ω0 t) , (8.39)

we have

1 1
y(t)y(t + τ ) = B 2 cos(ω0 τ ) − B 2 sin(2ω0 t + ω0 τ ) . (8.40)
2 2

and find
B2
Ryy (τ ) = cos(ω0 τ ) . (8.41)
2
We find average powers Rxx (0) = A2 /2 and Ryy (0) = B 2 /2 and in both
cases the period of the autocorrelation function is the same as the one of

104
the input function x(t).
This is a general property, as can be shown by expressing an arbitrary
function x(t) with period T by a Fourier series. With Eq. (3.4) we find
∞ ∞ 2
a20 X a2k
  X  
2πkτ bk 2πkτ
Rxx (τ ) = + cos + cos
4 2 T 2 T
k=1 k=1
∞  
X 2πkτ
+ ak bk sin . (8.42)
T
k=1

The term on the last line of Eq. (8.42) originates from cross correlation of
sine and cosine terms of the series, which can be computed by integrating
ak bk sin(ω0 t) cos(ω0 (t + τ )) over one period.

500
a 400
b
5
300

200
0
100

-5 0

-100
-10
-200
0 2 4 6 8 0 10 20 30 40 50 60

300 7000

250
c 6000
1 d
200 5000 0

150 4000
-1
100 3000
0 2 4 6 8
2000 time (a.u.)
50

1000
0

0 2 4 6 8 0 10 20 30 40 50 60
time (a.u.) frequency (a.u.)

Figure 8.4: Detection of a noisy signal x(t) via autocorrelation. (a) Cosine wave
with frequency ν = 3 a.u. and unity amplitude superimposed by white noise with
S/N = (1/3)2 (simulation with 1024 data points). (b) Spectrum of the noisy signal
obtained by DFT. (c) Autocorrelation function Rx x(τ ) of the noisy signal. (d)
Spectrum of the autocorrelation function. The approximate frequency ν0 = 3.05
a.u. can be determined from the clearly visible peak. The inset shows the cross
correlation function Rzx (τ ) obtained with a subsidiary signal z(t) = cos(2πν0 t).
Note that phase and amplitude of Rzx (τ ) have noise-related deviations from the
corresponding parameters of x(t).

105
This property of the autocorrelation function can be used to detect a small
oscillatory signal superimposed by strong noise. For example, the oscilla-
tory signal with 1024 data points and a period 0.33 a.u. at signal-to-noise
ratio of (1/3)2 is not at all reccognizable by looking at Figure 8.4a. In the
spectrum obtained by DFT (Figure 8.4b), a peak at frequency 3 a.u. can be
recognized, but at this signal-to-noise ratio it is hard to be certain that this
is due to a true signal. Surprisingly, the oscillation is rather well seen in the
autocorrelation function of the signal (Figure 8.4c) and from the DFT of the
autocorrelation function (Figure 8.4d) the frequency can be extracted with
rather good precision. Cross correlation of the original signal with a sub-
sidiary function with this frequency provides an estimate for the oscillatory
function, which has approximately 20% amplitude error and 45◦ phase error
in this particular case. Larger errors can occur for this noise level (the errors
are themselves stochastic) and the frequency estimate can also have a larger
error, leading to beating in the cross correlation function Rzx . However, the
presence of the signal is always reliably detected.1 .
This can be understood by considering the total signal y(t) as a sum of the
wanted signal x(t) and noise n(t),

y(t) = x(t) + n(t) . (8.43)

The autocorrelation function can be expanded as

Ryy (τ ) = Rxx (τ ) + Rxn (τ ) + Rnx (τ ) + Rnn (τ ) . (8.44)

Since x(t) and n(t) are uncorrelated, the terms Rxn (τ ) and Rnx (τ ) vanish.
The term Rnn (τ ) consists of a large peak at τ = 0 (in the center of the plot
in Figure 8.4c), and much smaller noise at other times which stems from
finite-time observation. An estimate ν0 for the frequency of the oscillation
can be directly read off, determined by fitting, or taken from the spectrum
of the autocorrelation function.2 This estimate may be slightly erroneous
due to the noise in the DFT spectrum.
We can now construct a subsidiary signal z(t) = cos(2πν0 t). The cross
correlation function
1
Ryz (τ ) = x(τ ) (8.45)
T
1
This is no longer true for 1024 data points at S/N = (1/4)2
2
It is as well recognizable in the power spectrum.

106
is an estimate of the original signal x(t). Due to the error in ν0 and finite
observation time, amplitude and phase of Ryz (τ ) are slightly erroneous.

8.2.6 Felgett’s advantage

We have seen in the previous section that an oscillatory signal is more easily
recognized in frequency domain than in time domain. This is related to
an improvement of signal-to-noise ratio due to the integration in FT, which
averages noise. An oscillatory signal will contribute to the Fourier integral
proportionally to the length of the time-domain signal, whereas noise is
only proportional to the square root of this length. It can be shown that
this fact leads to an advantage in sensitivity when measuring and Fourier
transforming time-domain data rather than performing a scan of frequency,
wavelength, or field, provided that the whole spectrum can be excited at
once. This multiplex advantage was first observed by Felgett and is thus
also termed Felgett’s advantage. In NMR spectroscopy, Felgett’s advantage
is usually large, so that continuous-wave frequency scans are no longer used.
The same advantage arises in FT-IR spectroscopy, where the spectrum can
be computed by FT from an interferogram. If the spectrum cannot be
completely excited, as is usually the case in EPR spectroscopy, continuous-
wave sweep techniques may be more sensitive than FT spectroscopy.

107
Chapter 9

Alternatives to
one-dimensional Fourier
transformation

9.1 Spectral estimators and parameter estimators

9.1.1 Asssumptions on the observed system
To this point we have only assumed that the system to be characterized is
linear and time invariant. In this situation, the system operator is com-
pletely determined if we know the impulse response, or the step response,
or the spectrum. In general, systems of quantum states do not respond
linearly to excitation. For systems involving more than two levels, non-
linear response can provide additional information on the system. Such
experiments are not, however, simple to generalize and are therefore dis-
cussed in the context of a particular kind of spectroscopy. Examples are
two-dimensional (2D) and 3D NMR spectroscopy, which provides informa-
tion on correlations between nuclear spin transitions, 2D IR spectroscopy,
which correlates molecular vibrations, and second harmonic generation in
optical spectroscopy, which provides information on second-order electric
susceptibility. Two-dimensional experiments are best analyzed in terms of
two-dimensional spectra, which can be obtained by 2D FT.

108
9.1.2 2D spectroscopy
Data from some non-linear spectroscopic experiments can be processed by
two-dimensional FT, a concept that can be extended to multi-dimensional
FT. We consider a spectroscopic experiment where coherence on a transition
(m, n) evolves for time t1 with angular frequency ωmn and is then transferred
to a transition (p, q) where it evolves for time t2 with angular frequency ωpq
before being detected. Two types of coherence transfer allow for recognizing
the correlation between the two transitions from the dependence of the signal
on times t1 and t2 .

Phase modulation

In the first transfer type, the phase of the coherence on transition (p, q)
immediately after the transfer is given by

φ1 (t1 , 0) = ωmn t1 + φ(0) (0)

mn + φpq , (9.1)

(0)
where φmn is the initial phase of the coherence on transition (m, n) at time
(0)
t1 = 0 and φpq is a constant phase shift introduced by the transfer. After
evolution, at detection time t2 , the phase of the coherence is given by

φ2 (t1 , t2 ) = ωmn t1 + φ(0) (0)

mn + ωpq t2 + φpq . (9.2)

If detection causes a phase shift, we can include this by redefining φ0pq . The
signal eiφ2 can be written as
n  o
hp (t1 , t2 ) = exp i ωmn t1 + φ(0)
mn + ω pq t2 + φ (0)
pq
n  o n  o
= exp i ωmn t1 + φ(0)
mn exp i ω pq t2 + φ (0)
pq

= c1 c2 − s1 s2 + ic1 s2 + is1 c2
= c1 (c2 + is2 ) + is1 (c2 + is2 ) . (9.3)

with
   
c1 = cos ωmn t1 + φ(0)
mn , s 1 = sin ωmn 1t + φ (0)
mn
   
c2 = cos ωpq t1 + φ(0)
pq , s 2 = sin ω t
pq 1 + φ (0)
pq . (9.4)

109
If required, coherence decay according to Eq. (1.14) can be accounted for
by multiplication with the decay factor
   
t1 t2
d(t1 , t2 ) = exp − exp − . (9.5)
T2,mn T2,pq

Amplitude modulation

In the second transfer type, the phase of the coherence on transition (p, q)
(0)
after the transfer is independent of t1 (φpq ), whereas the amplitude is given
by c1 . In this case, the signal without decay is given by

ha (t1 , t2 ) = c1 c2 + ic1 s2 = c1 (c2 + is2 ) , (9.6)

and coherence decay can again be accounted for by the factor defined in Eq.
(9.5).

2D Fourier transformation

The signals hp (t1 , t2 ) and ha (t1 , t2 ) are processed by applying FT subse-

quently along the two time dimensions,

H(ω1 , ω2 ) = F {h(t1 , t2 )} = F (1) F (2) {h(t1 , t2 )}

Z ∞ Z ∞
= exp(−iω1 t1 )dt1 exp(−iω2 t2 )h(t1 , t2 )dt2 .(9.7)
−∞ −∞

The result is most easily understood for amplitude modulation, ha (t1 , t2 ).

In the first FT performed along dimension t2 , c1 is a constant. Hence, we
obtain the spectrum along frequency domain ω2 , as we would have obtained
it from a one-dimensional experiment. This spectrum is a Dirac peak at
angular frequency ωpq or, if the decay factor is included, a Lorentzian peak
at ωpq with width 2/T2,pq . The same spectrum is obtained at different times
t1 , multiplied by the amplitude factor c1 (t1 ). We also obtain an imaginary
part, which contains a dispersion line modulated along t1 with a cosine. This
imaginary part can be dropped, since c1 is an even function of ωmn . Hence,
we do not have information on the sign of ωmn in the t1 dimension and can
be content with a cosine FT instead of a complex FT.
During the second FT along dimension t1 , ω2 is a constant. The FT of the
cosine function are Dirac peaks at frequencies ωmn and −ωmn . Because of
this symmetry, the half plane with ω1 < 0 can be discarded. The spectrum

110
then contains a single peak at (|ωmn |, ωpq ). Information on the sign of ωmn
is lost.

a b

c d

Figure 9.1: Phase-twisted lineshape in 2D FT spectra with phase modulation. (a)

2D absorption peak A1 A2 . (b) 2D dispersion peak D1 D2 . (c) Phase-twisted peak
A1 A2 − D1 D2 . (c) Interference effects for two overlapping phase-twisted peaks.
Adapted from A. Schweiger, G. Jeschke, Principles of pulse electron paramagnetic
resonance, Oxford University Press, Oxford, 2001.

For phase modulation, the same peak is obtained, but according to the last
line of Eq. (9.3), we have to consider an additional term is1 (c2 + is2 ).
With respect to the t2 dimension, this term behaves the same way as the
other one. Hence, before FT along the t1 dimension we have an imaginary
part that is a sine modulated spectrum with a single peak at frequency
ωpq . A problem arises, since the first FT also generates an imaginary part,
which contains a dispersion peak at frequency ωpq . If Eq. (9.7) is strictly
adhered to, the product of this dispersion spectrum with is2 is a real input
for the FT along the t1 dimension. This leads to superposition of a negative
dispersion component D1 D2 with the wanted positive absorption component
A1 A2 , which is termed a phase-twisted lineshape (Figure 9.1). Dropping
the imaginary part after the first FT does not help in this case, since the
dispersion lineshape along t1 remains admixed. The problem can only be
circumvented if a mirror image h0p (t1 , t2 ) can be generated in which ωmn is
replaced by −ωmn . Sometimes this is possible by appropriate phase changes
of the excitation pulses. If so, the sum of phase-modulated signal and mirror
image correspond to an amplitude-modulated signal.

111
9.1.3 Deficiencies of spectrum estimators
Both one- and two-dimensional FT do not presuppose that the signal is a
linear combination of damped harmonics in time domain and thus a lin-
ear superposition of Lorentzian peaks in frequency domain. The numerical
versions DFT and 2D DFT provide estimates for the frequency response
(spectrum) of the system, whatever shape this response happens to have.
Such spectral estimators work equally well for spectra with a small, known
number of resolved transitions and spectra with a large unknown or even
infinite number of overlapping transitions.
The price paid for this universality is a fixed relation between relative fre-
quency resolution and the number of points acquired in time domain. The
sampling interval T is determined by the required bandwidth (Nyquist cri-
terion, Eq. (7.6)) and the frequency resolution by the number of data points
within the bandwidth, which is only twice as large as the number of data
points in time domain.1 This becomes a problem in multi-dimensional spec-
troscopy, where the size of the data array and the acquisition times for
indirect dimensions scale with the product of the number of time-domain
data points in all dimensions.2
A further disadvantage of FT arises from the assumption that the data are
known starting at time t = 0. Time-domain spectrometers tend to have a
certain deadtime that passes after high-power excitationn, before the small
response of the system can be detected. According to the shift theorem,
such a time shift leads to a phase shift. The phase shift is proportional
to frequency ωk throughout the whole lineshape of a peak centered at ωk .
Thus, it cannot be perfectly eliminated by a phase correction that is linear
with respect to the frequency variable ω. If broad peaks overlap, destructive
interference may occur. In that case, a small absorption peak superimposed
on the wing of the lineshape of a large absorption peak may appear as a hole
in the lineshape even in a magnitude spectrum (Figure 9.2a). The problem
can be understood as an admixture of dispersive cross terms Dm Dn to the
lineshape during computation of the magnitude spectrum, similar to the
admixture of a dispersive lineshape in 2D FT of phase-modulated signals.
In the one-dimensional case this admixture can be systematically averaged
1
Further zero filling does not improve resolution, but merely interpolates intermediate
data points in frequency domain.
2
This section and the following sections have profited from reading V.A. Mandelshtam,
Progr. Nucl. Magn. Reson. 38, 159-196 (2001) and from personal discussions with
Vladimir Mandelshtam.

112
by summing power spectra with different time delays and computing the
magnitude spectrum as the square root of the sum (Figure 9.2a). In practice,
the first power spectrum P0 (ω) is computed from the original signal with the
first data point at deadtime t0 = td . For each of the N − 1 subsequent power
spectra, with N being the number of original data points, a new time-domain
signal yk (t) is constructed by removing the first point of the current signal
and appending a zero. This corresponds to an initial time t0 = td + kT for
the k th data set. The cross term-averaged magnitude spectrum is obtained
as v
uN −1
uX
Mca =t Pk . (9.8)
k=0

a b

0 5 10 15 20 25 0 5 10 15 20 25
n (MHz) n (MHz)

Figure 9.2: Destructive lineshape interference due to deadtime (simulation). A

magnitude spectrum obtained by DFT is shown. Two absorption peaks at fre-
quencies of 7 and 14 MHz with transverse relaxation times of 0.5 and 0.15 µs and
amplitudes of 0.2 and 10 were simulated. Data were acquired with a sampling in-
terval of 18 ns and the first point was assumed to be missing (deadtime of 18 ns).
(a) Magnitude spectrum obtained by DFT. The small peak at 7 MHz appears as a
hole in the lineshape. (b) Magnitude spectrum obtained by cross-term averaging.
The small peak at 7 MHz is absorptive.

9.1.4 Parameter estimators

If we have additional information on our system, a spectrum estimator may
not be the best way to process the signal, as it ignores this information. In
high-resolution spectroscopy, the signal is of the form

N   
X 1
y(t) = dk exp iωk − t , (9.9)
T2,k
k=1

113
where the number of lines N may be known or unknown. The total infor-
mation content in this signal are 3N real parameters, namely the amplitude
dk , angular frequencies ωk , and widths 2/T2,k of the N lines. If the sig-
nal is noiseless, 3N data points will theoretically suffice to determine the
parameters with arbitrary resolution, independently of the maximum obser-
vation time. This requires a signal processing algorithm that directly fits
the parameters rather than first computing a spectrum and determining the
lineshape parameters from this spectrum. Such an algorithm is a parameter
estimator.
Depending on the problem at hand, parameter estimators can suffer from
numerical instability in the presence of noise. In general, such instability is
encountered when the fit problem is non-linear. At first sight, the problem
defined by Eq. (9.9) appears to be non-linear. However, that problem can
be recast in linear form.

9.2 Harmonic inversion

It is easy to compute the time-domain signal y(t) in Eq. (9.9) from the
known parameters dk , ωk , and T2,k . In contrast, it is hard to solve the
inverse problem of computing the parameters from the known time-domain
signal. This is a general feature of signal processing: it is usually easy
to compute the response of a system with known parameters and hard to
estimate the parameters from the system response, in particular, in the
presence of noise. For the system defined by Eq. (9.9) parameter estimation
is called the harmonic inversion problem (HIP), since the individual terms of
the sum are harmonic functions. Harmonic inversion is more elaborate than
FT, but mathematically it is not a particularly hard problem.3 In practice,
when significant noise is present, the technique works rather well for data
extrapolation, but not so well for direct computation of the parameters.
In the following, we thus focus on data extrapolation, which can solve the
deadtime and resolution problems and assume that the extrapolated signal
is further processed by FT.
3
The following sections are based on P. Koehl, Progr. Nucl. Magn. Reson. 34, 257–299
(1999)

114
9.2.1 Linear prediction
We first consider a time-domain signal that arises from only a single term in
Eq. (9.9), sampled at a constant time interval T . This signal can be written
as
(k)
yl = Ae−lT /T2,k eiωk lT = AZkl , (9.10)

with
Zk = e(iωk −1/T2,k )T . (9.11)

Hence, for l > 0 we have

(k) (k)
yl = Zk yl−1 . (9.12)

The next point in the signal can be linearly predicted if we know Zk (note
(k)
that the yl are complex). For a noiseless signal we can compute Zk af-
(k) (k) (k) (k)
ter measuring two points, for instance y0 and y1 , since Zk = y1 /y0 .
Hence, if we know beforehand that the signal is a single damped harmonic
function and is noiseless, we can determine frequency and width with ar-
bitrary precision from two complex data points. The complex amplitude
(magnitude and initial phase φ0 ) is also encoded in the data points, but
does not influence Zk .
This property of damped harmonic functions can help to solve the deadtime
problem. Assume that we do not know data points y0 through y3 , but we do
(k) (k) (k) (k) (k)
know y4 and y5 . As Zk = y5 /y4 and, from Eq. (9.12), yl−1 = yl /Zk ,
we can iteratively predict all the missing points.
In practice, data are not noiseless, hence Zk is determined with a certain er-
ror. This leads to errors in the parameter estimates and to an extrapolation
error in linear prediction of missing parts of the signal that increases with
increasing time difference from the last known point. If we have more than
two points, we can reduce the error due to noise, since Zk can be computed
from any pair of points and in the mean of these values, noise will partially
average.
The concept carries over to the sum in Eq. (9.9), as the problem cast in this
form is linear. If we know N noiseless data points, we can predict earlier
and later data points, using

N
X
yl = bm yl−m , (9.13)
m=1

where the complex coefficients bm take the role previously taken by Zk and,

115
together with the known data points, encode the full information on ampli-
tudes, frequencies, widths, and phases of all N lines.4 In fact, the concept
of Eq. (9.13) can be applied to other problems than HIP, as long as they
are autoregressive.
As the signal is not noiseless, we shall need more data points to obtain a
reasonable quality of the linear prediction and reasonable estimates of the
line parameters. Hence, we define

M
X
yl = bm yl−m + el (M ≥ N ) , (9.14)
m=1

where the el account for the noise, which is not autoregressive. Note that the
el are not the noise amplitudes at point l themselves, but only the deviation
of the actual noise amplitude from a linear prediction from noise of the
preceding M points.
This equation can be restated in a matrix form appropriate for linear least-
square solution,
Tb = a − e , (9.15)

where b is the vector of length M of the linear prediction coefficients with,

a is a vector of length L − M of the measured data points from point M up
to point L − 1 and e is a vector of length L − M of noise coefficients from
point M up to N − 1. T is an (L − M ) × M Toeplitz matrix
 
yM −1 yM −2 . . . y0
 
 yM yM −1 . . . y1 
T=
 ...
 . (9.16)
 ... ... ... 

yL−2 yL−3 . . . yL−1−M

the index ranges arise from the condition that M should be as large as
possible and imply M = L/2, so that the length of vectors a and e is
also M . Since the rows of T are not necessarily linearly independent, the
equation cannot generally be solved by direct matrix inversion. In other
words, the system of equations may not have a unique solution. We can
ensure a unique solution by requiring that the sum of squares of the noise
4
The N complex coeffeicients bm encode 2N real informations, the same applies to the
N complex data points.

116
coefficients e attains a minimum,

bLSQ = minb ||Tb − a||2

. (9.17)

9.2.2 Singular value decomposition

If we do not know the number of lines N beforehand, we do not know the
number of linearly independent rows in Eq. (9.15) and thus, we do not know
the rank of matrix T. According to Eq. (9.13) only N coefficients bm encode
spectral information, whereas the remaining ones correspond to noise. With
an algorithm that can reveal the rank of T, we will obtain an estimate of the
number of lines with significant amplitude (as compared to the noise level)
and we can distinguish between signal and noise based on the property of
signal to be lineraly predictable.
The mathematical standard approach for revealing rank and null space of a
matrix is singular-value decomposition (SVD), which is a matrix factoriza-
tion
T = USV∗ , (9.18)

where U (left singular vectors) and V (right singular vectors) are orthogonal
matrices, ∗ denotes the conjugate transpose, and S is a diagonal matrix of
singular values. SVD is related to the eigen problems of matrices TT∗
and T∗ T which provide the left and right singular vectors, respectively,
as eigenvectors. The non-zero eigenvalues in both eigenproblems are the
squares of the singular values.
With SVD, the solution of Eq. (9.17) is obtained as

bLSQ = VS+ U∗ a , (9.19)

where S+ is the pseudo-inverse of S,

(
1
for Skk 6= 0
S+
kk =
Skk
. (9.20)
0 for Skk = 0

For a noiseless signal y, the rank of T is N , so that there exist only N

non-zero singular values and, correspondingly, N non-zero coefficients bm .
Hence, we can redefine M = N for this case. In principle the frequencies and
relaxation rates of the N lines can be obtained from the bm by the extended
Prony method. The Zk defined by Eq. (9.11) are the complex roots of a

117
 3
Skk
0
4 8
k
12 16
a
c

b
d
0 1 2 3 4 5 -20 -10 0 10 20
time (a.u.) frequency (a.u.)

Figure 9.3: Linear prediction with singular value decomposition (simulation). The
inset is a plot of the singular values. (a) Original noisy time-domain data trace with
256 points. The real part is blue, the imaginary part red. (b) Points 17–256 were
predicted from coefficients bm that were computed from points 0–15 only. The sin-
gular value matrix S was truncated at N = 7. (c) Absorption spectrum computed
by DFT from the original data trace. (d) Absorption spectrum computed by DFT
from a data trace where points 16–255 were replaced by the linear prediction.

characteristic polynomial P (z) with

M
X
P (z) = z M
− bm z M −m , (9.21)
m=1

where b0 = 1. In computing a peak list, one can reject all complex roots
that correspond to exponentially rising rather than decaying signals, i.e. all
those one that have abs(Zk ) > 1.
The complex amplitudes ck eiφ0,k are obtained by least squares computation
from
M
X
yl = ck eiφ0,k zkl + wl , (9.22)
k=1

where the wl are the true noise contributions and ||wl ||2 is to be minimized.
However, this approach is sensitive to the ratio between sampling interval
and transverse relaxation times and, generally, producing a peak list in that
way from noisy data is not a very stable procedure.
If the signal contains noise, singular values beyond N approach a plateau
value larger than zero. The number of singular values that are significantly
larger than this plateau value is then an estimate for the number of lines
that have magnitudes significantly larger than the noise level. If the number
of lines N can be estimated in that way, signal-to-noise ratio of the predicted

118
points can be improved by setting all diagonal values Skk of S with k > N −1
to zero5 before computing the pseudo-inverse S+ and the bm from Eq. (9.19).
Such truncation of S does not reduce the number of non-zero prediction
coefficients bm .
Truncation does, however, reduce the noise in the predicted points. This
is illustrated in Figure 9.3. Linear prediction coefficients bm (m = 0 . . . 15)
were computed from a slightly noisy time-domain signal (a) corresponding
to the spectrum shown in (c). Matrix S was truncated by setting diagonal
elements Skk with k > 6 to zero. In this case, the original spectrum clearly
reveals that there are only seven lines. Points 16–255 of the original data
set were then replaced by the linear prediction according to Eq. (9.13). The
spectrum obtained by DFT from this data set (d) looks less noisy, although
the amplitudes, frequencies, and widths of the lines may not necessarily be
better defined.

9.3 Fredholm integral equations

In solid-state magnetic resonance spectroscopy, frequency, amplitude and
possibly relaxation rate and phase may depend on sample orientation. In
this case, the number of lines N in Eq. (9.9) is infinite. Such cases may also
be encountered in optical spectroscopy due to inhomogeneous line broaden-
ing (lattice strain). The LP-SVD technique further become impractical, if N
is larger than the number of data points that can be measured. This number
is limited by the maximum observation time tmax , where the signal has de-
cayed below noise level, and the sampling interval T , at which the spectrum
is fully contained in the Nyquist band. Roughly speaking, LP-SVD becomes
impractical if the homogeneously broadened lines of individual transitions
are not resolved.
In such situations, spectra may still be described by a small number of
parameters. As an example, we consider the dipole-dipole interaction be-
tween two spins at distance r in a macroscopically isotropic sample. The
normalized time-domain signal corresponding to the separated dipole-dipole
5
In the equations above, the first index is always 0.

119
interaction has the form
Z π/2
y(t)
3 cos2 θ − 1 ωdd t sin θdθ


= cos
y(0) 0
Z 1
3x2 − 1 ωdd t dx ,


= cos (9.23)
0

with
µ0 ~ 1 CAB
ωdd = γA γB 3 = 3 . (9.24)
4π r r
Since the permittivity of vacuum µ0 , Planck’s quantum of action ~, and the
gyromagnetic rations γA and γB are known constants, the only unknown is
the distance r. The spectrum Y (ω) that results by FT from Eq (9.23) is a
Pake pattern.
A distribution P (r) of distances r cannot easily be extracted from Y (ω). As
the time-domain signal depends only on r, it must be possible to transform
y(t)/y(0) to P (r). This corresponds to inversion of the integral equation
Z ∞ Z 1  
y(t) 2

CAB
= P (r) cos 3x − 1 t dxdr . (9.25)
y(0) 0 0 r3

9.3.1 Fredholm integral equations of the first kind

Transformation from frequency to distance domain is an example for a Fred-
holm integral equation of the first kind,
Z b
g(t) = K(t, s)p(s)ds (9.26)
a

where we can make the identifications g(t) = y(t)/y(0), p(s) = P (r) and
where K(t, s) is the kernel. Inferring p(s) from g(t) is termed an inverse
problem.
In practical situations, g(t) is given as a vector g of discrete values. In
analogy to FT, the length of this vector and the sampling interval will have
some relation to the resolution and range of values in s domain. Hence, p(s)
is also well represented by a vector p, so that K(t, s) can be replaced by
a matrix. The discrete form of the Fredholm integral equation of the first
time thus reads
g = Kp . (9.27)

Apparently, we only need to invert the kernel K to solve for p. By proper

definition of the sampling intervals for p we can ensure that the kernel is a

120
square matrix, which can be inverted.

9.3.2 Well-posed and ill-posed problems

It turns out that solution of discrete Fredholm equations of the first kind by
simple inversion of the kernel matrix is unstable with respect to noise. Very
small changes in g lead to very large changes in p. Such behavior is typical
for an ill-posed problem. Stable solutions by kernel inversion are obtained
for well posed problems, with the conditions

• Existence. For every input vector g ∈ G, there exists an output vector

p = F (g), where p ∈ P . Here and in the following, G and P represent
the time and s domain, respectively.

• Uniqueness. For any pair of input vectors g, z ∈ G, it follows that

F (g) = F (z) if and only if g = z.

• Continuity (stability). The mapping is continuous, that is, for any

 > 0 there exists ξ = ξ() such that the condition dG (g, z) < ξ
implies that dP (F (g), F (z)) < , where d(·, ·) represents the distance
metric between two arguments in their respective spaces.6

The relation between time and frequency domain in Fourier transformation

is an example for a well-posed problem. The relation between time and dis-
tance domain for the dipole-dipole interaction is a case where the conditions
of existence and uniqueness are fulfilled. However, it turns out that one
cannot construct a kernel matrix for the discrete problem that fulfils the
continuity condition.
The degree of ”ill-posedness” of a discrete form of the Fredholm equation can
be quantified by the condition number of the kernel matrix, which measures Konditionszahl
the sensitivity of solution of a system of linear equations to errors in the
input data. The condition number is the ratio between the largest and
smallest singular value of K and can thus be computed by SVD (Matlab
has a dedicated function cond for this purpose). Condition numbers near
unity indicate a well-posed discrete problem.
6
qPIn the simplest case, the distance metric for two vectors g and z is given by dG (g, z) =
2
k (gk − zk )

121
9.4 Regularization
Regularization techniques stabilize the solution of ill-posed problems by in-
troducing additional information, in particular, additional constraints on
the solution. A typical restraints is smoothness of the solution, which can
be defined as a small square norm of the derivative or second derivative. If
the solution is known to be a probability density function p(s), the addi-
tional constraint p(s) ≥ 0 can be imposed, since probability densities cannot
be negative. A non-negativity constraint often stabilizes the solution very
strongly, since part of the problems with continuity arise from the fact that
negative contributions in s domain can partially offset positive contributions
at nearby points.

9.4.1 Tikhonov regularization

The conventional, but unstable least-squares solution of Eq. (9.27) is

pLSQ = minp ||Kp − g||2 = minp (E(p)) ,

(9.28)

where E(p) = ||Kp − g||2 is the function to be minimized. Smoothness can

be imposed by redefining E as

E(p, α) = ||Kp − g||2 + α||Dp||

b 2, (9.29)

where D
b is a linear differential operator. Depending on the problem at
hand, D
b can be the first or second derivative. For the example of time-
domain signals from dipole-dipole interactions, the second derivative D
b =
d2 /dr2 performs best. The derivative or second derivative is a measure for
roughness, hence it minimization maximizes smoothness.
The regularization parameter α in Eq. (9.29) determines the relative weight-
ing between the square deviation of the simulated signal from the measured
signal (first term) and roughness (second term). In general, the optimum
choice of α depends on the expected error in g and the expected smooth-
ness of the data, for example, the width of the features that are expected
in the distance distribution P (r). Whereas the noise level in g can often be
estimated, the resolution of p(s) is usually not known beforehand.

122
 -12
1
a a = 0.01 b
0.9 -16

log h
0.8
-20 aopt = 100
0.7
a = 1000

-24
0.6
0 0.5 1 1.5 2 -3 -2.8 -2.6 -2.4 -2.2
t (µs) log r

c d
P(r)

1.5 2 2.5 3 3.5 4 4.5 1.5 2 2.5 3 3.5 4 4.5

r (nm) r (nm)

Figure 9.4: Tikhonov regularization (simulation). Dipole-dipole interaction be-

tween two electron spins at a mean distance of 3 nm was assumed. The distance
has a standrad deviation of 0.3 nm. The additional constraint P (r) > 0 was im-
posed. (a) Original noisy time-domain data trace (black) and simulation for the
optimum regularization parameter aopt = 100 (green). (b) L curve. Regularization
parameters α = 0.01 (red, overfitting), αopt = 100 in the corner of the L curve
(green), and α = 1000 (blue, oversmoothing) are indicated. (c) Original distance
distribution (black) superimposed by the optimum result of regularization (green).
(d) Original distance distribution (black) superimposed by the results obtained
with regularization parameters corresponding to overfitting (red) and oversmooth-
ing (blue).

L curve method

A general method for estimating the optimum regularization parameter αopt

can be derived from the following plausible considerations. As long as α is
too small, solution will have a large roughness

b α ||2
η(α) = ||Dp (9.30)

and a small square deviation between simulation and experimental data

ρ(α) = ||Kpα − g||2 . (9.31)

123
Accordingly, on increasing α we expect a strong decrease in η, but only a
moderate increase in ρ. Vice versa, if α is too large, the roughness η will be
very small, but the square deviation between simulation and experimental
data ρ will be large. A further increase in α will then lead to only a moderate
decrease in η, but to a strong increase in ρ. A parametric plot of log η versus
log ρ will thus have an L shape, with the corner of the L corresponding
to an optimum regularization parameter αopt that balances errors due to
oversmoothing a errors due to overfitting the data by reproducing noise by
features in p(s).

124
Chapter 10

Digital signal processing

Spectrometers always perform some analog signal processing, since the de-
tected signal is transmitted and usually amplified before digitization and the
transmission lines and amplifiers have a certain bandwidth. Hence, some ef-
fective bandpass filter is invariably applied to the signal. In addition, it
may be necessary to apply a lowpass filter before analog-digital conversion
to avoid aliasing of high-frequency noise. Such aliasing would otherwise oc-
cur if the Nyquist frequency corresponding to the ADC sampling rate is
lower than the upper cutoff frequency of the effective bandpass filter.
Additional signal processing is more flexible and less expensive when it is
performed as digital signal processing. For instance, very high filter orders
can be realized digitally, while they would require a lot of electronic com-
ponents in analog realization. Digital signal processing is limited only by
computing times and this limitation is relevant only if the processing needs
to be performed in real time.

10.1 Real-time digital signal processing

In some applications, signals need to be processed at the rate at which they
are produced. Obvious examples are hearing aids or telecommunication
systems. In a spectrometer such a need can arise if the amount of data pro-
duced has to be reduced before storage. In such applications it is acceptable
and unavoidable to introduce a delay of a few sampling intervals between
input and output, although the maximum length of this delay may also be
dictated by the application. In most of this chapter we shall assume data
post-processing, where processing rate is of minor importance. However,

125
the same algorithms can be implemented on the digital signal processors
and field-programmable gate arrays introduced in the next two sections.

10.1.1 Digital signal processors

Digital signal processors (DSP) are the most widely used hardware for real-
time signal processing. The main principle that distinguishes DSP from
general central processing units (CPUs) is parallel execution of several op-
erations in a single processor cycle in a way that is optimized for signal
processing. For instance, by separating the memory for commands and data
it is possible to load a command and a datum in the same cycle (Harvard
architecture). Cyclical data buffers allow for fast addressing of the current
signal point without software overhead and for parallel loading and access
to data. Instructions for multiplication and addition in a single cycle en-
hance the speed of those operations that are typical for the implementation
of digital filters.
Depending on the application, DSP may work with fixed-point or floating-
point arithmetics. Fixed-point arithmetics processors are somewhat cheaper
and more energy efficient and are thus used in consumer products that are
sold in large quantities at rather low prices. Their disadvantage is more
complicated programming and the introduction of significant quantization
noise (round-off errors) in each step of signal processing. For spectroscopy,
floating-point signal processors are better suited.
Typical DSPs from the super-Harvard single-chip computer (SHARC series)
of Analog devices work at rates of 40 MIPS (Mega instructions per second)
corresponding to a clock rate of 25 ns and, due to parallelization, can perform
up to 120 MFLOPS (Mega floating-point operations per second) at peak
rates. Typical tasks performed by such DSPs are filtering and FFT. If
data are generated at sampling intervals much shorter than 25 ns, real-time
processing with DSPs is currently not possible.

10.1.2 Field-programmable gate arrays

At high sampling rates or large numbers of operations per sample field-
programmable gate arrays (FPGA) may have better performance in digital
signal processing than DSPs. An FPGA is a very large scale integrated
circuit that contains a large number of repeated small logic blocks. Pro-
gramming of such a device is more flexible, but also more complicated than

126
for typical microprocessors or digital signal processors. This flexibility does,
however, allow for more extensive parallelization than with DSPs and thus
for higher processing rates.
An example for an FGPA is the Virtex-5 processor of Xilinx, for which a
DSP development kit exists. This development kit interfaces to Matlab,
which can be used for filter design.

10.2 General properties of time-discrete systems

10.2.1 Properties and description
In the previous chapters we assumed that the signal is time-continuous,
although we did already discuss consequences of discretization and digitiza-
tion. In the following we take it for granted that a discrete, digital signal
with sampling interval T contains all the information that we need. We
can then replace the system that converts a continuous input signal to a
continuous output signal by a time-discrete system that converts an input
series xn = x(n) to an output series yn = y(n). The time-discrete system is
described by the operator equation

yn = S {xn } , (10.1)

where yn and xn denote the complete series. In the follwoing we shall use
the index k, as in xk to denote a single data point in the series.
We assume that the system is linear and time-invariant, and, unless noted
otherwise, causal. Most systems are also dynamic, i.e., their response de-
pends not only on the current input signal xk but also on previous input
values xl with l < k. We are most interested in stable systems. For a stable,
time-discrete system any bounded allowed time-discrete input xn results in
a bounded response yn .
Many of the concepts that we introduced for time-continuous LTI systems
can be easily extended to time-discrete LTI systems. For instance, we can
define a time discrete Dirac impulse as
(
1 for k = 0
δk = . (10.2)
0 for k 6= 0

The value yk of a series is thus the convolution of the series yn with a Dirac

127
impulse at time k,
∞
X
yk = yl δl−k = yn ∗ δn |n=k (10.3)
l=−∞

Convolution of two series is defined as

∞
X ∞
X
x n ∗ yn = xn−k yk = xk yn−k . (10.4)
k=−∞ k=−∞

Accordingly, the time-discrete impulse response is

gn = S{δn } . (10.5)

The impulse response completely characterizes the time-discrete system.

The response to any input signal xn can be computed as

yn = xn ∗ gn (10.6)

10.2.2 FIR and IIR systems

Time-discrete systems can be classified by the length of their impulse re-
sponse. If the impulse response is finite, i.e., exactly zero after a certain
number of points, the system is called an FIR system (finite impulse re-
sponse). If the impulse response is infinite, the system is an IIR system
(infinite impulse response).

10.2.3 Causality
A time-discrete system is causal if and only if the impulse response vanishes
for negative indices,
gk = 0 for k < 0 (10.7)

10.2.4 Stability
A time-discrete system is stable if and only if the impulse response is abso-
lutely summable,
∞
X
|gk | < ∞ . (10.8)
k=−∞

128
10.2.5 Frequency response
If a time-discrete LTI system is excited by a harmonic oscillation

xk = Aeiω0 kT , (10.9)

it responds with an output signal of the same angular frequency ω0 . The

proportionality factor between input and output signal is the Fourier trans-
form G∗ (ω) of the time-discrete impulse response gn at frequency ω0 . The
Fourier transform G∗ (ω) is called the frequency response of the time-discrete
system S.

10.2.6 Difference equations

In Chapter 4 we had seen that the response of (analog) electronic circuits can
be described by differential equations, i.e. equations that contain differential
operators. Likewise, a time-discrete system can be described by a linear
difference equation
n2
X m2
X
aν yk−ν = bµ xk−µ , (10.10)
ν=n1 µ=m1

where n1 ≤ n2 and m1 ≤ m2 . The coefficients aν and bµ completely char-

acterize the system. In a causal system, all indices n1 , n2 , m1 and m2 are
larger than or equal to zero.

Figure 10.1: Structure of a digital IIR filter consisting of multiplier/accumulators

and delay elements (from Wikimedia Commons).

We now consider the case n1 = m1 = 0 and n2 = m2 = M , which is of

particular interest for digital IIR filters. All coefficients can be normalized
by a0 , as this leads only to multiplication of the signal by a constant. Such
a difference equation can be easily realized by the structure shown in Figure

129
10.1 for M = 2, where the boxes z −1 represent delay by the sampling interval
T , the triangles multiplication with the corresponding coefficient, and the
plus characters in circles addition (accumulation). Because of feeding output
signal back via the multipliers with ak coefficients, the filter is recursive.
This leads to the infinite impulse response.
A finite impulse response is thus obtained for n1 = n2 = 0 and m1 = 0, m2 =
M . Again, we can normalize all coefficients by a0 . The structure of such an
FIR filter is shown in Figure 10.2, where a slightly different representation
with Σ characters for accumulation is used.

Figure 10.2: Structure of a digital FIR filter consisting of multiplier/accumulators

and delay elements (from Wikimedia Commons).

10.3 z transformation
When discussing filters in Section 4.2.2, Laplace transformation was shortly
mentioned as a tool in filter theory. The transfer function of a filter was
written based on the poles and zeros in Eq. (4.44). For digital filters, we
consider this theory in more detail. The equivalent of Laplace transforma-
tion for time-discrete signals is the z transformation,
∞
X
Y (z) = Z {yn } = yk z −k . (10.11)
k=−∞

An acausal series yn can be split into a causal part y+ (k) with y+ (k) = 0
for k < 0 and an anticausal part y− (k) with y− (k) = 0 for k > 0. The

130
corresponding z transforms are
∞
X
Y+ (z) = yk z −k and
k=0
X0
Y− (z) = yk z −k (10.12)
k=−∞

with the two-sided z transform given by

Y (z) = Y+ (z) + Y− (z) − y(0) . (10.13)

Note that z is a complex argument and that the z transform is a continuous

function in the complex plane, even though the input data yn are discrete.

10.3.1 Existence of the z transform

The z transform exists for a given z and input series yn only if the sum in
Eq. (10.11) converges. The causal part Y+ exists outside a circle with center
0 and radius r+ = 1/R (|z| > r+ ) with

r+ = lim sup|yn |1/n , (10.14)

n→∞

whereas the anticausal part Y− exists inside a circle with center 0 and radius
r− = r (|z| < r− ) with

r− = (lim sup|y−n |1/n )−1 , (10.15)

n→∞

Hence, Y (z) exists in a circular ring between radii r− and r+ . We are mostly
concerned with strictly causal series, so that the z transform exists outside
a circle with radius r+ .

10.3.2 Relation between z transformation and FT

If the convergence region of the z transform contains the unit circle, the
time-discrete Fourier transform also exists. The Fourier transform G(ω)
can be computed from the z transform by setting z = eiωT . Note the
subtle difference between the time-discrete Fourier transform and the result
of DFT- the time-discrete Fourier transform is a continuous function. Note
also that the unit circle is a mapping of the Nyquist band to z domain, i.e.

131
the argument ω can be considered to run from −Ω/2 to Ω/2 with Ω = 2π/T ,
so that ωT runs from −π to π.

10.3.3 Properties of the z transformation

The z transformation is linear. The z transform of a time-inverted series is
 
1
Yb (z) = Y . (10.16)
z

Shift theorem

For a general series, a time shift by l points leads to

Yb (z) = z −l Y (z) (shif t theorem) . (10.17)

The time-shift expression has to be used with caution for causal series,
which would become acausal for l < 0. This is corrected by eliminating the
contribution from the acausal part

l−1
X
Z {yn+l } = z l Y+ (z) − yk z l−k for l ≥ 0 . (10.18)
k=0

Convolution theorem

The most important property relates to convolution, which was defined for
series in Eq. (10.4). We have

Z {xn ∗ yb } = X(z) · Y (z) , (10.19)

i.e. as for the FT, convolution is replaced by a multiplication. The conver-

gence region is given by

max {rx+ , ry+ } < |z| < min {rx− , ry− } . (10.20)

The convolution theorem is important for cascaded filters.

132
Correlation functions

Further important properties relate to cross correlation and autocorrelation.

Cross correlation of two time series is defined as
∞
X
Rxy (k) = xl+k yl
l=−∞
X∞
= xν y−(k−ν) = xn ∗ y−n . (10.21)
ν=−∞

The z transform of correlation functions can thus be computed with the

rules for time inversion and convolution. We find for cross correlation
 
1
Sxy (z) = Z {Rxy (n)} = X(z)Y (10.22)
z

and for autocorrelation

 
1
Sxx (x) = X(z)X . (10.23)
z

The z transform of a cross correlation function converges for

   
1 1
max rx+ , < |z| < min rx− , . (10.24)
ry− ry+

10.4 Transfer function

An LTI system is completely characterized by the impulse response and the
response to any input signal can be computed by Eq. (10.6). This relation
simplifies by z transformation

Y (z) = X(z) · G(z) , (10.25)

where the z transform G(z) of the impulse response is the transfer function Systemfunktion
of the system. If the convergence region is also given, the transfer function
is a complete description of the system.
The transfer function can be easily obtained by z transformation of the
difference equation, Eq. (10.10). Using the shift theorem we obtain
n2
X m2
X
aν z −ν Y (z) = z −µ bµ X(z) , (10.26)
ν=n1 µ=m1

133
which can be inserted into Eq. (10.25) to give
Pm2 −µ
Y (z) µ=m1 bµ z
G(z) = = Pn2 −ν
. (10.27)
X(z) ν=n1 aν z

To simplify this equation, we choose the same index ranges for both sums,
namely M1 = min {n1 , m1 } and M2 = max {n2 , m2 }, with missing coeffi-
cients being set to zero. Hence,
PM2 −µ
µ=M1 bµ z
G(z) = PM2 . (10.28)
ν=M1 aν z −ν

Furthermore, we can shift all indices by −M1 (the system is time invariant)
and expand the fraction by z M with M = M2 − M1 to have non-negative
exponents of z. Thus, we obtain
PM M −µ
µ=0 bµ z P (z)
G(z) = PM = . (10.29)
ν=0 aν z M −ν Q(z)

Hence, the transfer function of a time-invariant system described by a dif-

ference equation is a rational function with numerator polynomial P (z) and
denominator polynomial Q(z).

10.4.1 Poles and zeros

The numerator and denominator polynoms can be written in product form
using their M roots z0µ and z∞ν ,
QM
µ=1 (z − z0µ )
G(z) = G0 QM . (10.30)
ν=1 (z − z∞ν )

Note that, depending on the degree of the polynomials before reindexing,

there may be less than M zeroes z0µ or less than M poles z∞ν .
The transfer function in sum form,

M
X Aν
G(z) = A0 + , (10.31)
z − z∞ν
ν=1

can be obtained by partial fraction decomposition. From this, the impulse

response is computed by an inverse z transformation, which is in general a
non-trivial procedure. For the case at hand, correspondence tables such as

134
Table C.3 in [Leon et al.]1 give

M
X
n−1
gn = A0 δn + Aν z∞ν u(n − 1) , (10.32)
ν=1

where u(n) is a time-discrete step function analogous to the one defined in

Eq. (2.13), which is given by
(
0 for n < 0
u(n) = . (10.33)
1 for n ≥ 0

With Eq. (10.32) we can consider the stability criterion given in Eq. (10.8).
We find that a system is stable, if and only if the condition
∞
X
|z∞ν |n−1 < ∞ (10.34)
n=1

is fulfilled for all poles with indices ν. This leads to the stability criterion

|z∞ν | < 1 , ν = 1, . . . , M . (10.35)

In other words, all poles of the transfer function G(z) must be positioned
inside the unit circle for the system to be stable.

10.4.2 Systems in series and in parallel

If two systems S1 with transfer function G1 (z) and S2 with transfer function
G2 (z) are cascaded, i.e. passed by the signal in series, the transfer function
of the combined system is the product of the individual transfer functions

G(z) = G1 (z) · G2 (z) cascaded systems . (10.36)

If the input signal is split to two systems and the outputs of these two
systems are additively combined (parallel systems), the transfer function is
the sum of the individual transfer functions

G(z) = G1 (z) + G2 (z) parallel systems . (10.37)

1
Fernando Puente Leon, Uwe Kiencke, Holger Jäkel, Signale und Systeme, Oldenbourg.

135
10.4.3 Frequency-response function, phase response, and group
delay
The continuous frequency response characteristic G∗ (ω) can be computed
by the z transform at z = iωT , as indicated in Section 10.3.2. This can
be seen by computing G∗ (ω) as the FT of the impulse response (see Eq.
(2.19)),
∞ ∞
X X
−n
G∗ (ω) = gn e−iωnT = gn eiωT
n=−∞ n=−∞
iωT
= G(z = e ). (10.38)

This relation is applicable to any stable system, as in such a case the poles
are inside the unit circle, while the z = iωT describe a line on the unit circle.
Note also that this implies the limited frequency band of the discrete signal,
which we have already encountered in Section 7.1.3. For some discussions, it
is thus convenient to define a normalized angular frequency Ω = ωT , which
falls into an interval of width 2π.
The amplitude response M (ω), corresponding to the magnitude spectrum
of the time-cintinuous system, is given by
q
|G∗ (z)|2 = G(z)G∗ (z)
p p
M (ω) = = G(z)G(z −1 ) .
z=eiωT z=eiωT
(10.39)
The phase response is defined as

φ(ω) = ∠G∗ (ω) = ∠G(z = eiωT ) . (10.40)

Finally, we are interested in the group delay

1 d d
τg = − φ(ν) = − φ(ω) . (10.41)
2π dν dω

If the group delay is not constant, different frequency components are de-
layed by the system for different times.

10.4.4 Relation between pole-zero plot and response func-

tions
According to Eq. (10.29) the transfer function is determined by the poles
and zeros up to a multiplicative constant G0 . Hence, the frequency response

136
 imaginary W=p/2

z01 P1
w0
z¥ 1 Q1

W0 W=0
Q2 1 real

P2
a2
z¥ 2
b2 e
z02

Figure 10.3: Pole-zero plot of a notch filter with two zeros on the unit circle at
Ω = ±π/2 and two zeros (1 − )e±iπ/2 . The value of the frequency response at
angular frequency ω0 , corresponding to normalized angular frequency Ω0 , can be
computed from the lengths P1 , P2 , Q1 , and Q2 of the zeros Z01 and z02 and poles
z∞1 and z∞2 by also taking into account angles α2 and β2 as well as the analogously
defined angles α1 and β1 .

and phase response can also be described in terms of the poles and zeros,
QM iωT − z


iωT µ=1 e 0µ
G∗ (ω) = G(z = e )= G 0 QM . (10.42)
iωT − z
ν=1 (e ∞ν )

In a pole-zero plot in the complex plane, the factors in the numerator and
denominator correspond to vectors between a zero or pole, respectively, and
the point on the unit circle that corresponds to current frequency ω0 . These
vectors are most conveniently described in polar coordinates, as illustrated
in Fig. 10.3. In such plots, a zero is marked by a circle, while a pole is
marked by a cross. The polar coordinates are given by

eiω0 T − z0µ = Pµ · eiβµ , Pµ ≥ 0 ,

eiω0 T − z∞ν = Qν · eiαν , Qν ≥ 0 ,
G0 = C · eiγ , C ≥ 0 . (10.43)

137
With these definitions we find
M
Y M
Y PM PM
βµ −i
iω0 T
Q−1 iγ+i αν


G z=e =C Pµ ν e
µ=1 ν=1 . (10.44)
µ=1 ν=1

The advantage of the polar coordinate representation becomes apparent

when computing frequency and phase response,

M
Y M
Y
M (ω) = C Pµ Q−1
ν , (10.45)
µ=1 ν=1

M
X M
X
φ(ω) = γ + βµ − αν . (10.46)
µ=1 ν=1

The group delay τg is obtained as the negative derivative of Eq. (10.46) with
respect to ω0 .
Note than in a Bode diagram (see Section 4.2.1) of log M (ω) vs. ω the
components from the individual poles and zeros in Eq. (10.45) become
additive. A zero z0k on the unit circle leads to a factor Pk = 0 and thus
to complete suppression of the frequency corresponding to this point on the
unit circle. A pole on the unit circle leads to an infinite signal.

Example: Notch filter

2 A notch filter is a special type of band-stop filter that suppresses a narrow Kerbfilter
frequency band around an undesired frequency ωA . From discussion of poles
and zeros it is clear that complete suppression can be achieved by a zero on
the unit circle at ΩA = ωA T . To generate the narrow stop band without
making the filter unstable, a pole should be placed inside the unit circle
close to the zero. The principle of this compensation of the zero by the pole
can be appreciated from Fig. 10.3. Although frequency ω0 is not far from
the zero at ΩA = π/2, the length of vectors P1 and Q1 are almost equal,
corresponding to no suppression at all. It is also clear from this geometrical
picture that frequencies with P1 ≈ Q1 are the closer to ΩA the closer the
pole is put to the zero. Limitations are put by numerical stability of the
filter and by the long rise and decay times that result for a very narrowband
filter.
2
This example is adapted from [Leon et al.]. A detailed discussion, including analytical
treatment can be found there.

138
Based on these considerations, a symmetric notch filter can be designed
with two zeros at e±ΩA and two poles at (1 − )e±ΩA , as shown in Fig. 10.3
for ΩA = π/2. The notch frequency follows with the sampling interval T as
νA = ΩA /(2πT ). We now consider the transfer function of such a notch filter
with ΩA = π/2. In this case the zeros are z01 = i and z02 = −i, whereas the
poles are z∞1 = (1 − )i and z∞2 = −(1 − )i. Thus, the transfer function
is given by
(z − i) (z + i)
G(z) = . (10.47)
(z − (1 − )i) (z + (1 − )i)
In the complex plane, z01 has the coordinates (<{z01 }, ={z01 } = (0, 1)
whereas the point ω0 = Ω0 /T has the coordinates cos(Ω0 ), sin(Ω0 ). Hence,
q
P1 = (0 − cos(Ω0 ))2 + (1 − sin(Ω0 ))2
q
= cos2 (Ω0 ) + 1 − 2 sin(Ω0 ) + sin2 (Ω0 )
p p
= 2 − 2 sin(Ω0 ) = 2 − 2 cos(Ω0 − π/2) . (10.48)

Analogous computations give

p
P2 = 2 − 2 cos(Ω0 + π/2)
p
Q1 = 1 − 2(1 − ) cos(Ω0 − π/2) + (1 − )2
p
Q2 = 1 − 2(1 − ) cos(Ω0 + π/2) + (1 − )2 . (10.49)

Thus, the amplitude characteristic according to Eq. (10.45) becomes

P1 · P2
A(Ω) =
Q1 · Q2
p
2 − 2 cos(Ω0 − π/2)
= p
1 − 2(1 − ) cos(Ω0 − π/2) + (1 − )2
p
2 − 2 cos(Ω0 + π/2)
·p . (10.50)
1 − 2(1 − ) cos(Ω0 + π/2) + (1 − )2

This amplitude characteristic is shown in Fig. 10.4 as a black line for  =

0.2667, corresponding to the pole-zero plot shown in Fig. 10.3. Often,
analytical expressions are not required and it is sufficient to compute Eq.
(10.45) numerically. The result of such a numerical computation (Matlab
script notch filter.m) for ΩA = π/4 and  = 0.05 is shown as a blue
line in Fig. 10.4. Note that signal components at frequencies far from the
notch frequency are slightly amplified and that this amplification is larger for

139
broader stop bands. The two-pole notch filter as defined here has an even
amplitude characteristic and an odd phase characteristic. The strongest
phase shifts occur near the notch frequency.

1.4

1.2

0.8
A(W)

0.6

0.4

0.2

0
-p -p/2 0 p/2 p
W

Figure 10.4: Amplitude characteristics of two-pole notch filters with normalized

angular notch frequencies ΩA = π/2 (black,  = 0.2667, corresponding to Fig. 10.3)
and ΩA = π/4 (blue,  = 0.05).

10.4.5 Minimal phase system and allpass

For a given amplitude characteristic, the phase characteristic becomes min-
imal if all zeros are within the unit circle,

|z0µ | ≤ 1 (µ = 1 . . . M ) minimal phase system . (10.51)

For any time-discrete system S we can find a corresponding minimal phase

system SM with the same amplitude characteristic by partioning S into SM
and an allpass, which is defined as a time-discrete stable system with rational
transfer function and A(ω) = 1 for all angular frequencies ω. An equivalent
definition can be based on the pole-zero plot: An allpass is a stable system
with zeros having mirror symmetry to poles.
For partitioning we factorize the numerator polynom P (z) into two poly-

140
noms P1 (z) and P2 (z), with

M1
Y
P1 (z) = (z − z0µ ) , |z0µ | ≤ 1 ,
µ=1
M
Y
P2 (z) = (z − z0µ ) , |z0µ | > 1 . (10.52)
µ=M1 +1

The polynom

M M
" !#
Y Y 1
P20 (z) = ∗ ∗



1 − zz0µ = −z0µ · z− ∗ (10.53)
z0µ
µ=M1 +1 µ=M1 +1

∗ that feature mirror symmetry with respect to the roots

has roots 1/z0µ
z0µ of P2 . Hence P2 /P20 is an allpass. A system with the original transfer
function G(z) can be obtained by cascading this allpass with a system that
has additional zeros P20 (z),

P (z) P1 (z) · P2 (z) P1 (z)P20 (z) P2 (z)

G(z) = = = · 0 . (10.54)
Q(z) Q(z) Q(z) P2 (z)

Here,
P1 (z)P20 (z)
GM = (10.55)
Q(z)
is the transfer function of a minimal phase system with the same amplitude
characteristic as the one of the original system. This follows from the unity
amplitude characteristic of the allpass and from fact that the mirrored roots
of P20 must all reside inside the unit circle, as the roots of P2 are all residing
outside the unit circle.
This principle can be used, for instance, to compensate for the amplitude
characteristic and construct a pure phase shifter. For that, we can compute
GM as shown above and invert it. The time-discrete system with transfer
function 1/GM is stable, since all its poles are inside the unit circle (they are
the zeros of the minimal phase system GM . By cascading the original system
with a system with transfer function 1/GM , the amplitude characteristic is
compensated and a phase shifter with the phase characteristic of the allpass
is obtained.

141
10.4.6 Transfer function of FIR systems
For an FIR sytem the current output value yn is determined by M + 1
earlier input values xn . . . xn−M , where we neglect the delay introduced by
processing. This is expressed by the difference equation

M
X
a0 yn = bµ xn−µ . (10.56)
µ=0

The transfer function is obtained by z transformation

M
X
0
a0 z Y (z) = bµ z −µ X(z) , (10.57)
µ=0

which gives PM −µ
Y (z) µ=0 bµ z
G(z) = = . (10.58)
X(z) a0
We divide all coefficients bµ by a0 to obtain

M
X
G(z) = bµ z −µ . (10.59)
µ=0

Hence, the FIR system is fully specified by M + 1 coefficients bµ .

10.4.7 Transfer function of IIR systems

From Eq. (10.29) we find by dividing all coefficients by a0 ,
PM M −µ
µ=0 bµ z
G(z) = . (10.60)
zM + M M −ν
P
ν=1 aν z

Hence, the IIR system is fully specified by M + 1 coefficients bµ and M

coefficients aν .

10.4.8 Linear phase systems

For a filter, we would prefer zero phase shift, i.e. φ(ω) = 0 for all ω within
the Nyquist band. Such zero phase-shift filters cannot be realized by causal
systems. In many applications, the optimum solution is then a constant
group delay, where all frequency components are delayed by the same time.
According to Eq. (10.41) this corresponds to a constant derivative of the

142
phase response with respect to frequency and thus to a linear dependence
of the phase shift on frequency,

φ(ω) = cω , c ∈ R . (10.61)

The transfer function of such a linear phase system can be expressed as

G(z) = |G(z)| · z −k , k ∈ R . (10.62)

This can be easily seen by considering the frequency response (z = eiωT ),

G∗ (ω) = |G∗ (ω)| · eiωkT = A(ω) · eiωkT , (10.63)

which gives
φ(ω) = ∠G∗ (ω) = −ωT k , (10.64)

corresponding to c = −kT . Note that −c = kT is the group delay.

We now express the normalized group delay k as the sum of an integer
number nν (Matlab floor(k)) and a number α smaller than one (Matlab
k-floor(k)),
k = nν + α , nν ∈ Z , 0 ≤ α < 1 . (10.65)

Hence,
φ(ω) = −ωT nn u − ωT α . (10.66)

It can be shown ([Leon et al.]) that a real-valued system can have linear
phase in sections only if
∞
X
gn sin [ωT (n − nν − α) + φ0 ] = 0 , (10.67)
n=−∞

where we have introduced an additive constant that allows for phase jumps
by 2π.

Design of a linear phase FIR filter

We assume that the impulse response of our filter is purely real (gn ∈ <)
and vanishes outside the interval [0, M T ]. Furthermore, we take M as even
and choose the normalized group delay k so that

M
nν = . (10.68)
2

143
From Eq. (10.65) it then follows that α = 0. With these assumptions we
can design linear-phase FIR filters with symmetric or antisymmetric impulse
response. We choose the case of symmetric impulse response,
 
M
gn = gM −n , n = 0, 1, . . . , . (10.69)
2

The condition for linear phase response, Eq. (10.67), thus becomes3

M/2     
!
X M M
0 = gn sin ωT n − + gM −n sin ωT −n
2 2
n=0
M/2   
X M
(gn − gM −n ) sin ωT n − . (10.70)
2
n=0

This condition is fulfilled indeed, since (gn − gM −n ) follows from Eq. (10.69).
We can also demonstrate the linear phase response by computing the fre-
quency response

M M
M
−iωT M
gn e−iωT (n− 2 )
X X
−iωT n
G∗ (ω) = gn e =e 2

n=0 n=0
 
M/2−1
M M M
gn e−iωT (n− 2 ) + gM −n e−iωT ( 2 −n) 
X
= e−iωT 2 g M
2
n=0
 
M/2−1 h i
M M M
e−iωT (n− 2 ) + e−iωT (n− 2 ) 
X
= e−iωT 2 g M gn
2
n=0
 
M/2−1  
M X M 
= e−iωT 2 g M + 2 gn cos ωT n −
2 2
n=0
M
= e−iωT 2 Rg (ω) , (10.71)

where Rg (ω) is the amplitude response, as it is a pure cosine.4

For this symmetric impulse response filter, the amplitude response is an
even function of frequency. An odd frequency response can be obtained
analogously with an antisymmetric impulse response, gn = −gM −n .
3
We double count the term for n = M/2. This does not matter, however, as the sine
function is zero for that term.
4
In fact, Rg can become negative. This can be cured by introducing phase jumps of
±π.

144
10.5 Matlab tools for digital signal processing
10.5.1 Useful Matlab functions for system analysis and de-
sign
The following (selected) Matlab functions are convenient for numerical im-
plementation of the concepts treated in this Chapter

• roots finds zeros of a polynom; given coefficients bµ and aν in Eq.

(10.59) or Eq. (10.60), poles and zeros of the system can be obtained

• poly finds polynom coefficients from zeros; given poles and zeros of the
system, the transfer function according to Eq. (10.59) or Eq. (10.60)
can be found

• filter filters a time series x and provides output y, when given coef-
ficients aν and bµ

• residuez Z-transform residue partial fraction expansion, directly finds

poles from given coefficients aν and bµ or coefficients from poles and
residues

• prony given the impulse response and the orders of the numerator and
denominator polynom, finds coefficients aν and bµ

• freqz returns the digital filter complex frequency response, when given
coefficients aν and bµ

• invfreqz returns real coefficients aν and bµ for a filter that has least
square deviation from a requested frequency response function, given
also the orders of numerator and denominator polynom (numbers of
poles and zeros); a discrete frequency axis must be specified

• zplane creates a pole-zero plot when given filter coefficients aν and bµ

• grpdelay returns a vector of group delays, when given filter coefficients

and a discrete normalized frequency axis

• impz returns the impulse response, when given filter coefficients

Of these functions, roots, poly, and filter are contained in the basic
Matlab program. The remaining functions require the Signal Processing
Toolbox.

145
10.5.2 Filter design with Matlab
The Signal processing toolbox in Matlab also contains tools with a graphical
user interface for filter analysis and design. The analysis and visualization
tool fvtool is convenient if one has a set of filter coefficients aν and bµ . Co-
efficients for Butterworth and Chebyshev type I and II filters mentioned in
Section 4.2.2 can be computed by functions butter, cheby1, and cheby2, re-
spectively. Cauer filters are also called elliptical filters and can be computed
by ellipt. Additional functions, such as buttord help in specification, by
returning the lowest order of a filter that fulfills certain specifications for the
stop and passband and for the edge frequencies of stopband and passband.
Digital filters are not restricted to the analog types discussed in Section 4.2.2.
It is possible to specify maximum deviations δp and δs for the passband and
stopband, such that

1 − δp ≤ Y ei Ω 1 + δp , 0 ≤ Ω ≤ Ωp


(10.72)

and
Y ei Ω


≤ δs , ΩS ≤ Ω ≤ π . (10.73)

Note that Ω is a normalized angular frequency in the range between −π and

π with ω = Ω/T . The normalized angular frequencies Ωp and ΩS are edge
frequencies of the passband and stop band.
Filters with such specifications can be designed with fdatool. The δp and
δs are specified on a dB scale. In fdatool the response type can be specified
as lowpass, highpass, bandpass, bandstop (a generalized notch filter) and
several specialized responses, such as a differentiator or a system for Hilbert
transformation. For each filter, IIR and FIR types and, within those types,
different design methods are provided. From the specifications, fdatool
computes the filter coefficients and displays the various response functions,
the pole-zero plot, and even round-off noise power, which corresponds to
quantization noise introduced by the digital filter.
With the concepts learned in this lecture course, fdatool is fairly intuitive
to use. As output, the tool can also provide a Matlab .m file that sets the
specifications and computes the filter coefficients.

146