Time Varying Narrowband Communications Channels: Analysis and Implementation

Niklas Grip and Götz E.
Pfander
Time Varying Narrowband Communications

Channels: Analysis and Implementation
Technical Report No. 12

October 2007
School of Engineering and Science

Time Varying Narrowband
Communications Channels:
Analysis and Implementation1
Niklas Grip Götz E. Pfander
Department of Mathematics School of Engineering and Science

Luleå University of Technology Jacobs University Bremen
SE-971 87 Luleå Campus Ring 12
Sweden 28759 Bremen
Germany
E-Mail: [email protected], [email protected],

http: // math. jacobs-university. de/ pfander
Summary
We derive and describe a Matlab implementation of an efficient numerical algo-
rithm for the analysis of certain classes of Hilbert–Schmidt operators that naturally
occur in models of wireless radio and sonar communications channels.
A common short-time model of these channels writes the channel output as a
weighted superposition of time- and frequency shifted copies of the transmitted
signal, where the weight function is often called the spreading function of the
channel operator.
It is often believed that a good channel model must allow for spreading func-
tions containing Dirac delta distributions. However, we show that many narrow-
band finite lifelength channels such as wireless radio communications can be well
modelled by smooth and compactly supported spreading functions.
We derive a fast algorithm for computing the matrix representation of such
operators with respect to well time-frequency localized Gabor bases, such as pulse
shaped OFDM bases. Hereby we use a minimum of approximations, simplifica-
tions, and assumptions on the channel. The primary intended target application
1
The work on this project was supported by the German Research Foundation (DFG) Project
PF 450/1-1 as part of the DFG priority program TakeOFDM. During the final preparations of
the manuscript, the first mentioned author was supported by the Swedish Research Council,
project registration number 2004-3862.
is comparisons of how different parameter settings and pulse shapes affect the
diagonalization properties of an OFDM system acting on a given channel.
Finally, we provide and demonstrate the use of a Matlab implementation
that is fast enough for the number of carrier frequencies typically in use today.
A shortened and improved version of the main results in this report is published
in the journal paper [GP07], which excludes some discussions, proofs and the
software documentation contained here.
3
Contents
1 Introduction 2
2 Preliminaries 4
2.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Frequency localization of compactly supported functions . . . . . . 6
2.3 Gabor analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4 Hilbert–Schmidt operators . . . . . . . . . . . . . . . . . . . . . . . 9
3 The channel model 11

3.1 Single path frequency response . . . . . . . . . . . . . . . . . . . . . 12
3.2 Narrowband signals . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Finite lifelength channels . . . . . . . . . . . . . . . . . . . . . . . . 17
4 Discretization of the channel model 19

4.1 Gabor bases for near-diagonalization . . . . . . . . . . . . . . . . . 20
4.2 Discretization tools . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.3 Computing the channel matrix . . . . . . . . . . . . . . . . . . . . . 24
5 The algorithm and its implementation 31

5.1 The algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.2 Refinements and complexity . . . . . . . . . . . . . . . . . . . . . . 32
5.3 Parameters and window functions . . . . . . . . . . . . . . . . . . . 34
6 Applications 35
6.1 Mobile phone communications . . . . . . . . . . . . . . . . . . . . . 36
6.2 Satellite communications . . . . . . . . . . . . . . . . . . . . . . . . 38
6.3 Underwater sonar communications . . . . . . . . . . . . . . . . . . 39
6.4 Further spreading function examples . . . . . . . . . . . . . . . . . 40
7 Conclusions 42
A Fourier transform decay vs differentiability 44
B Matlab implementation: overview and function reference 46

Ω00
B.1 BPFS: Compute SΩ 00 (ν, t) . . . . . . . . . . . . . . . . . . . .
c
. . . . 47
B.2 ChannelSetupAndAnalysis: The main program . . . . . . . . . . . 50
B.3 CoeffOpMat: Compute the Coefficient Operator Matrix . . . . . . . 53
B.4 CutToEssSupp: Truncate a function to its essential support . . . . . 65
B.5 figsize: Set the orientation and size of the current figure . . . . . 67
Ωcd
,Ω
B.6 FTBPFS: Compute SH (·, t)(t0 ) . . . . . . . . . . . . . . . . . . . 69
B.7 GaborDual: Compute the dual of a Gabor Riesz basis . . . . . . . . 71
B.8 H: Compute samples (Hg)(kTγ ) using (4.12) . . . . . . . . . . . . . 76
B.9 l2Innerprod: Computes the l2 inner product (4.8a) . . . . . . . . . 78
B.10 multiplot: Command for saving the current figure in .eps and .emf
format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
B.11 PlotSpreadingFunction: Plots bandpass filtered spreading function 82
B.12 resample: Compute new samples from Nyquist frequency samples . 85
B.13 SetParameters: Set system parameters for different applications . . 87
B.14 SetPlot: Changes of line thickness, font size etc in plots . . . . . . 102
B.15 SincOmega: Computes the function sincω (x) . . . . . . . . . . . . 104
B.16 TFshift: Time-frequency shifts (integer time-shifts, fast) . . . . . . 105
B.17 TFshiftAndResample: Time-frequency shifts . . . . . . . . . . . . . 107
C Functions for two of the plots 109

C.1 FourierTranformDecay: Matlab-file for producing Figure 1, p. 8 . . 109
C.2 NrOfNonzeroNyquistFreqSamples: Produce Figure 6, p. 33 . . . . 114
5
Contents
List of Figures
1 Fourier transform decay examples. . . . . . . . . . . . . . . . . . . . 8
2 Multipath propagation . . . . . . . . . . . . . . . . . . . . . . . . . 12
3 System function supports and decays . . . . . . . . . . . . . . . . . 18
4 Convolution of compactly supported functions . . . . . . . . . . . . 20
5 Fourier transform supports of the synthesis Gabor system . . . . . . 26
6 Number of nonzero Nyquist frequency samples for B-splines . . . . 33
7 Approximative Gaussian pulse shape and truncation errors . . . . . 36
8 OFDM channel example . . . . . . . . . . . . . . . . . . . . . . . . 37
9 Satellite channel example . . . . . . . . . . . . . . . . . . . . . . . . 39
10 Underwater channel example . . . . . . . . . . . . . . . . . . . . . . 40
11 Spreading functions with more correlated coefficients . . . . . . . . 41
12 Off-diagonal decay comparisons . . . . . . . . . . . . . . . . . . . . 42
13 Program structure for the channel matrix computations . . . . . . . 46
14 Gabor Gram matrix indexing . . . . . . . . . . . . . . . . . . . . . 71
1
1 Introduction
Channel–dependent customization is expected to provide considerable performance
improvements in time-varying systems such as future generations of wireless com-
munications systems. Consequently, the idea of shaping the transmission pulses
in order to minimize the InterCarrier and InterSymbol Interference (ISI and ICI)
in Orthogonal Frequency Division Multiplexing (OFDM) communications is an
active research area in the applied harmonic analysis and signal processing com-
munities (see [BS01, MSG+ 05] and references therein). Even though some insights
can be gained from careful mathematical modelling and analysis, there remains a
need for fast algorithms and implementations aimed at the numerical evaluation of
performance improvements through pulse shaping. In this report, we discuss two
closely connected topics that we regard of vital importance to fulfill this demand.
1. We review the most important physical properties of wireless channels and

show how these lead naturally to a model of the short-time behavior of
a channel as an operator H that maps an input signal s to a weighted
superposition of time and frequency shifts of s, that is,
Z
Hs(·) = SH (ν, t) ei2πν(t−t0 ) s( · − t) d(ν, t) , K compact. (1.1)
K×[A,∞)
This model is well-known and the coefficient function SH is usually called the
spreading function of H. The model given by (1.1) is mostly used either under
the assumption that SH is square-integrable or that SH is a tempered distribu-
tion. The latter, weaker assumption suggests that the input signal s should be
a Schwartz class function and requires the use of distribution theory in the anal-
ysis of H, both of which we shall try to avoid. Therefore, we derive (1.1) using
some refinements of the standard multipath propagation model of the channel.
Our analysis implies that the short-time behaviour in many communications ap-
plications can be completely described by a smooth SH with rapid decay ensuring
“essentially compact” support2 . This model has the big advantage that it allows
for both Fourier analysis and numerical evaluation of the performance of OFDM
procedures without the need of deviating into distribution theory.
2. We employ the channel model described above to derive an efficient algo-

rithm for the numerical evaluation of ISI and ICI in pulse shaped OFDM
systems. Although we, for simplicity, justify our channel model only for sig-
nals s ∈ L2 (R), we derive our algorithm for a straightforward extension to
functions s ∈ L2 (Rd ). We expect this extension to be useful, for example,
for combining measurements from several antennas to a large set of samples
2
With essentially compact we mean that the function decays fast enough to assure that in
any practical application, the function values outside some “reasonably small” compact set are
very small compared to the overall noise level and therefore negligible (see also Section 6.1).
2
s(xn , tk ) of s in an (irregular) sampling grid of points xn in space and tk in
time, as is sometimes done in radio astronomy. Similarly for wireless com-
munications, both on the transmitter and receiver side, space diversity is one
of the most popular forms of diversity [Rap02]. The main idea there is to use
antenna arrays with antennas preferably placed several tens of wavelengths
apart for improving the chance to always have good transmission between at
least two of the antennas (see, e.g., [Rap02, BAB+ 01, JYGR74] for a more
complete coverage).
We shall now motivate and describe the principles of our discretization in some
detail:
For multicarrier modulation systems in general, the aim is the joint diagonal-
ization of a class of possible channel operators in a given environment. That is, we
try to find a transmission basis (gi ) and a receiver basis (filters) (e γj ) with the prop-
erty that all coefficient mappings that correspond to channels in the environment
have matrix representations Gi,j = hHgi , γ ej i that are as close to diagonal as possi-
ble, that is, |Gi,j | decays fast with |i − j|. In general, an easily computable inverse
of this coefficient mapping wouldP allow us to regain the transmitted coefficients
(ci ) in the input signal s = ci gi , and, therefore, the information embedded in
these coefficients, from the inner products hHs, γ ei i which are calculated on the
receiver side.
In wireline communications, the problem described above has a well accepted
solution, namely OFDM (also called Discrete MultiTone or DMT) with cyclic pre-
fix. Here, the transmission basis (gi ) and the receiver basis (γi ) are so-called Gabor
bases, that is, each basis consists of time and frequency shifts of a single proto-
type function which is often referred to as window function. Diagonalization of the
channel operator using Gabor bases with rectangular prototype function is then
possible since wired channels are assumed to be time-invariant. This allows us to
model such channel operators as convolution operators with complex exponentials
ei2πωt as “eigenfunctions”. This cyclic prefix procedure applies if the channel has
finite lifelength and is explained in more detail in Section 4.1 and with further
references in [Gri02]. The superiority of Gabor bases in comparison to Wavelet
and Wilson bases for wireline communications is examined in detail in [KPZ02].
Wireless channels are inherently time-varying. The generality of time-varying
channel operators and, in particular, the fact that they do not commute in general,
implies that joint diagonalization of classes of such channels cannot be achieved
as in the general case, so approximate diagonalization becomes our goal. In many
cases, for example in mobile telephony, the channel varies only “slowly” with
time. Hence, we use the results for time-invariant channels as a starting point and
consider in this report only the use of Gabor bases as transmitter and receiver
bases.
For such slowly time-varying systems, Matz, Schafhuber, Gröchenig, Hartmann
and Hlawatsch conclude that excellent joint time-frequency concentration of the
windows g and γ is the most important requirement for low ISI and ICI [MSG+ 05].
3
There, it is shown how to compute a γ (or an orthogonalization of the basis (gi ))
that diagonalizes the coefficient mapping in the idealized borderline case when the
channel is the identity operator (Gi,j = δi,j ). They show that both γ and the
corresponding orthogonalized basis inherit certain polynomial or subexponential
time-frequency decay properties from g. They also derive exact and approximate
expressions for the ISI and ICI and present an efficient FFT-based modulator and
demodulator implementation.
For multicarrier systems with excellent joint time-frequency localization of g
and γ, we derive, starting from our channel model, a procedure for the numerical
computation of the matrix entries Gi,j = hHgi , γ ej i under a minimum of assump-
tions, simplifications, or approximations. These properties make our approach
different from and complementary to a number of papers that use discrete Gabor
bases (sometimes under the name BEM or Basis Expansion Model) for time-
varying channels and statistical applications, such as [BLM04, BLM05, LM04,
LZGk03, MGk02, MGk03b, MGk03a, MLGk05, SA99, TL04].
The paper is organized as follows: The Notation and some mathematical pre-
liminaries are described in Section 2. We derive a channel model in Section 3, and
use it to derive formulas for the matrix elements in Section 4. In Section 5 we
describe a Matlab implementation of these formulas and give suggestions on how
to do the necessary parameter and window/pulse shape choices. In Section 6 we
provide typical system-dependent parameters for and demonstrate our software
on some example mobile phone communications, satellite communications and
underwater sonar communications applications. Finally, our conclusions follow in
Section 7.
2 Preliminaries
For completeness and easy availability we collect our notation in Section 2.1 and
give an overview of the mathematical tools that we shall use. In Section 2.2 we shall
discuss the availability of functions that are compactly supported and “essentially
bandlimited”, in particular, we explain how compactly supported functions can
be designed to have subexponential decay. Section 2.3 covers the Gabor system
expansions which are used to obtain diagonal dominant coefficient mappings of
channel operators. Finally, in Section 2.4 we discuss the Hilbert–Schmidt operator
theory and the integral representation of channel operators in terms of system
functions such as the spreading function and the time-varying impulse response.
2.1 Notation
We assume the reader to know some basic tools and notation from functional
analysis and measure theory, which otherwise can be found in [Fol99, Rud87].
The conjugate of a complex number z is denoted z. We use boldface font for
def def def
elements in Rd , write Rd+ = (0, ∞)d = R+ ×R+ ×· · ·×R+ and Zd+ = Zd ∩Rd+ . The
4
R
Fourier transform of a function f is formally given by fb(ξ) = Rd
f (t)e−i2πhξ,ti dt
def
for ξ ∈ Rd and l2 = l2 (Zd ×Zd ) is the Hilbert space of sequences (cq,r ) for which
the l2 -norm is given by
 1/2
def
X
k(cq,r )k2 =  |cq,r |2  < ∞.
q,r∈Zd
Throughout the paper we use Roman and Greek letters for variables that have a
physical interpretation as time or spacial variable and frequency, respectively. For
A, B, C, x, y, t, ν, ω ∈ Rd and r ∈ R we use the following shorthand notation:
def def ¡ ¢T
[A, B] = [A1 , B1 ] × · · · × [Ad , Bd ], 1 = 1 1 ··· 1 ,
def def ¡ ¢T
hx, yi = x1 y1 + x2 y2 + · · · + xd yd , xy = x1 y1 · · · xd yd
x def ¡ x1 x2 ¢T def ¡ ¢T
= r r · · · xrd , xr = xr1 xr2 · · · xrd ,
r
x def ¡ x1 x2 ¢T def
= y1 y2 · · · xydd , |x| = |x1 x2 · · · xd | ,
y
def def
Tt g = g(· − t), Mν g = g(·)ei2πhν,·i ,
· ¸ d
Y
def B B def sin(πωj xj )
IC,B = C − , C + and sincω (x) = .
2 2 j=1
πxj
Here, sincω is extended continuously to Rd and we shall frequently use that

Z
e−i2πhξ,xi dξ = e−i2πhC,xi sincB (x). (2.1)
IC,B
For ² > 0 we define the ²-essential support of a bounded function f : Rd → C to

be the closure of the set {x : ² ≤ |f (x)| / supx |f (x)| }. For an almost everywhere
defined function f , supp f denotes the intersection of the supports of all repre-
sentatives of f (and similarly for ²-essential support). For any set I, χI is the
characteristic function χI (x) = 1 if x ∈ I and χI (x) = 0 otherwise. The sets of n
times, respectively infinitely many, times continuously differentiable functions are
denoted C (n) and C (∞) , respectively.
We denote by Lp = Lp (Rd ) the Banach space of complex-valued measurable
functions f with norm
µZ ¶1/p
def p
kf kp = |f (x)| dx < ∞.
Rd
def R
L2 (Rd ) is a Hilbert space with inner product hf, gi = Rd f (x)g(x) dx. We say
that two sequences (fn ) and (gn ) of functions are biorthogonal if hfm , gn i = 0
whenever m 6= n and hfn , gn i = 1 for all n. The Wiener amalgam space
5
W (A, l1 ) = S0 (Rd ) (also named the Feichtinger algebra) consists of the set of
all continuous f : Rd → C for which
X
k(f (·)ψ(· − n))b k1 < ∞
n∈Zd
P
for some compactly supported ψ such that ψb ∈ L1 (Rd ) and n∈Zd ψ(x − n) = 1.
We write S00 for the space of linear bounded functionals on S0 . S0 is also a so-called
modulation space, described at more depth and with notation S0 = M 1,1 = M 1
and S00 = M ∞,∞ = M ∞ in [Grö00, FZ98].
A real-valued, measurable and locally bounded function w on Rd is said to be
a weight function if for all x, y ∈ Rd ,
w(x) ≥ 1 and w(x + y) ≤ w(x)w(y). (2.2)
For weight functions w we define L1w = L1w (Rd ) to be the family of functions
f ∈ L1 (Rd ) such that
def
kf k1,w = kf wk1 < ∞.
2.2 Frequency localization of compactly supported func-

tions
The Gabor window g in the introduction needs to be compactly supported in a
time interval short enough to satisfy typical maximum delay restrictions, such as
25 ms for voice communications. Moreover, its Fourier transform gb has to decay
fast enough to allow for reasonably high transmission power (which determines the
signal-to-noise ratio) without exceeding standard regulations on the allowed power
leakage into other frequency bands. In other words, g should have good, joint
time-frequency localization, which also is of great importance for achieving low
ISI and ICI [MSG+ 05]. For this reason, we seek to know to what extent compact
support of a function can be combined with good decay of its Fourier transform.
This classical question was first answered by Beurling [Beu38, Theorem V B] and
generalized from functions on R to functions on locally compact abelian groups,
such as Rd , by Domar [Dom56, Theorem 2.11]. Domar’s results are explained in
much more detail in [RS00, Ch. 6.3 + appendices]. We will describe a partial
result that describes the speed of decay of a Fourier transform by specifying for
how fast growing weight functions w it belongs to L1w (R). For describing the
asymptotic decay, we only need to consider continuous w with nondecreasing w(ξ)
and w(−ξ) for positive ξ. Suppose for one such w that there is some compactly
supported f ∈ L2 (R) such that fb ∈ L1w (R), then for any open set U ∈ R we can
choose a dilation and translation f (a · −x) of f with support in U . Since also
\
f (a · −x) ∈ L1w (R), the following theorem then provides the so-called logarithmic
integral condition (2.3) on the decay of w:
6
Theorem 1. Let w be a continuous weight function such that w(ξ) and w(−ξ)
are nondecreasing for positive ξ. Suppose that for every open set U ⊆ R there is
a non-zero function f ∈ L2 (R) such that supp f ⊆ U and fb ∈ L1w (R). Then
Z
log(w(ξ))
dξ < ∞. (2.3)
R 1 + ξ2
Proof. This theorem was proved in [RS00, Theorem A.1.13] for functions defined
on locally compact abelian groups, without our assumption about continuous and
nondecreasing w(±ξ) and with (2.3) replaced by the conclusion
X log(w(nν))
<∞ for all ν ∈ R. (2.4)
n∈Z
1 + n2
However, by our additional assumptions on w, the sum criterion for integrals
implies that (2.3) holds if and only if (2.4) holds for ν = 1. This concludes the
proof.
The logarithmic integral condition (2.3) limits the decay of both the am-
plitude and “the area under the tail” of fb. For example, the Fourier trans-
form of a compactly supported function f cannot be either O(e−α|ξ| ) nor fb(ξ) =
P α|n|
n∈Z φ(e (ξ − n)) for any φ ∈ C (∞) with support supp φ ⊆ [0, 1], because in
both cases, fb ∈ L1w (R) for w(ξ) = ea|ξ| and a < α but w does not satisfy (2.3). This
fact rules out the existence of compactly supported functions f with exponentially
decaying fb. However, Dziubański and Hernández [DH98] have shown how to use
a construction by Hörmander [Hör03, Theorem 1.3.5] to construct a compactly
supported function f whose Fourier transform is subexponentially decaying. That
is, they construct f such that for every 0 < ε < 1 there exists Cε > 0 such that
¯ ¯
¯b ¯ 1−ε
¯f (ξ)¯ ≤ Cε e−|ξ| , ∀ξ ∈ R.
From their example and standard techniques such as convolution with a charac-
teristic function, it is then easy to design for any compact set K a compactly
supported function f such that f (x) = 1 for x ∈ K, and fb is subexponentially
decaying.
Note however, that subexponential decay is not everything. For example,
− 1
the function f (x) = e 1−x2 χ[−1,1] (x) is a compactly supported C (∞) function,
so that fb(ξ) = O(1 + |ξ|)−n for all n ∈ N, whereas the function g(x) = (1 +
cos(πx))4 χ[−1,1] (x) has Fourier transform decay gb(ξ) = O((1 + |ξ|)−α ) for α = 9
but not for α > 9 (see Corollary 2, page 44). However, Figure 1 shows that
gb decays much faster down to amplitude thresholds such as the power leakage
restrictions described above (see also Figure 7, page 36). Thus it can be an im-
portant design issue to choose functions and forms of decay that are optimal for
a given application.
However, for simplicity and a clear presentation in this report, we shall con-
sistently claim subexponential decay although also other forms of decay are rapid
enough for all of our results to hold.
7
|g (ξ)|
|f (ξ)|
∧
∧
−10 −10
10 10
0 20 40 0 20 40
ξ ξ
(a) (b)
Figure 1: The Fourier transform decay after normalizing the following posi-
def − 1
tive functions to have integral 1: (a) f (x) = e 1−x2 χ[−1,1] (x). (b) g(x) =
(1 + cos(πx))4 χ[−1,1] (x).
2.3 Gabor analysis

Here, we give a brief review of some basic Gabor frame theory that is needed to
understand the relevance of the coefficient mappings that we introduce in (2.7)
below. For a more complete and general coverage of this subject, see, for exam-
ple, [Chr02, Gri02, Grö00].
A Gabor (or Weyl-Heisenberg) system with window g and lattice constants a
and b is the sequence (gq,r )q,r∈Zd of translated and modulated functions
def
gq,r = Tra Mqb g = ei2πhqb,x−rai g(x − ra).
The corresponding synthesis or reconstruction operator

def
X
Rg : l2 → L2 , Rg c = cq,r gq,r
q,r∈Zd
is defined with convergence in the L2 -norm if and only if its adjoint, the so-called
analysis operator
Rg∗ : L2 → l2 , Rg∗ f = (hf, gq,r i)q,r∈Zd ,

P
is bounded, i.e., if and only if q,r∈Zd |hf, gq,r i|2 ≤ B kf k22 for some B ∈ R+ and
all f ∈ L2 (Rd ) [Gri02, p. 14].
We call (gq,r )q,r∈Zd a Gabor frame for L2 (Rd ) if there are frame bounds A, B ∈
R+ such that for all f ∈ L2 (Rd ),
° °2
A kf k22 ≤ °Rg∗ f °2 ≤ B kf k22 . (2.5)
8
def
It follows from (2.5) that the frame operator Sg = Rg Rg∗ is invertible. We call a
def
frame with elements geq,r = Tra Mqb ge a dual Gabor frame if for every f ∈ L2 ,
X X
f= hf, geq,r i gq,r = hf, gq,r i geq,r (2.6)
q,r∈Zd q,r∈Zd
with L2 -norm convergence of both series. There may exist (infinitely) many dif-
ferent dual windows ge for g. However, we shall always consider the canonical dual
window, which is the minimum L2 –norm dual window [Jan98, p. 51]. The dual
frame has frame bounds A−1 , B −1 and the coefficients in (2.6) are not unique in l2 ,
but they are the unique minimum l2 -norm coefficients. It follows also from (2.5)
def
and (2.6) that Rge∗ picks coefficients from Cge = Rge∗ (L2 (Rd )) ⊆ l2 and that Rg is a
bounded invertible mapping of Cge onto L2 with bounded inverse Rg−1 = Rge∗ . By
this isomorphism and the usual definition of operator norms we can use two Ga-
bor frames (gq,r ) and (γq,r ) (possibly with different lattice constants) to obtain an
isomorphism of the family of linear bounded operators H : L2 (Rd ) → L2 (Rd ) with
the coefficient mappings G = Rγ∗ HRg , as illustrated in the following commutative
diagram.
H / L2 (Rd )
L2 (Rd )
O
Rg Rγ∗ (2.7)
G=Rγ∗ HRg ²
Cge / Cγ
We will provide an explicit expression for G in Section 4.

The frame (gq,r )q,r∈Zd is called a Riesz basis if Cge = l2 . Then, the coefficients
in (2.6) are truly unique and, as a consequence, (gq,r ) and (e gq,r ) are biorthogonal.
2.4 Hilbert–Schmidt operators

The mathematical framework for the use of Hilbert–Schmidt operators acting on
functions defined on locally compact abelian groups has been developed in great
generality in harmonic and functional analysis [FK98]. For the basic theory, see,
for example, [Con00, RS80] or [Fol95, Appendix 2], where the following classic
definition is used.
Definition 1 (Hilbert–Schmidt operator). Let H be a separable Hilbert

P space. A
bounded operator H : H → H is called a Hilbert–Schmidt operator if n∈Z kHen k2H <
∞ for some orthogonal basis (en )n∈Z of H. The Hilbert–Schmidt norm of H is
given by
Ã ! 12
X
kHkH.S. = kHen k2H .
n∈Z
9
This norm is independent of the choice of the orthonormal basis (en )n∈Z [Con00,
RS80] and can be calculated in a number of ways, one of which is given below.
Proposition 1 ([RS80, Theorem VI.23]). Let {M, µ} be a measure space and

H = L2 (M, dµ). Then a linear bounded operator H : H → H is Hilbert–Schmidt
if and only if there is a function
κ ∈ L2 (M × M, µ ⊗ µ)
such that Z
(Hf )(x) = κ(x, y)f (y) dµ(y) a.e. x . (2.8)
Moreover, Z
kHk2H.S. = |κ(x, y)|2 dµ⊗µ(x, y).
Within this report, we shall only consider Hilbert–Schmidt operators on the

Hilbert space H = L2 (Rd ). There are many similar versions of Proposition 1,
some of which can be obtained by applying partial unitary Fourier transforms to
the function κ and replacing (2.8) with corresponding mappings relating f or fb to
d as done in (2.11) below. For example, application of the leftmost
either Hf or Hf
unitary Fourier transform in (2.11a) gives the following corollary:
Corollary 1. A linear bounded operator H : L2 (Rd ) → L2 (Rd ) is Hilbert–Schmidt

if and only if there is a function SH ∈ L2 (Rd × Rd ) such that for all s ∈ L2 (Rd ),
Z
(Hs)(t0 ) = SH (ν, t)s(t0 − t)ei2πhν,x−ti d(t, ν). (2.9)
Rd ×Rd
The integral in (2.9) is defined in a weak sense. In fact, for s, g ∈ L2 (Rd ),

we have g(·)s(· − t) ∈ L1 , so that the short-time Fourier transform Vg g of g with
window s is well-defined as the function
Z
def
Vs g(ν, t) = g(t0 )s(t0 − t)e−i2πhν,t0 −ti dt0 (2.10)
Rd
on L2 (Rd × Rd ). Hence, Hs is defined to be the unique L2 (Rd )-function with
hHs, giL2 (Rd ) = hSH , Vs giL2 (Rd ×Rd ) .
Many similarly obtained system functions are known under a rich plethora
of different names in the literature, ranging back to a first systematic study by
Zadeh and Bello [Zad50, Bel63, Bel64] (see also [Ric03] for an overview). The
integral representations of importance in this text describe H in terms of the
spreading function SH , the kernel κH , the time-varying impulse response h, the
10
Kohn-Nirenberg symbol σH and the bifrequency function BH . These system func-
tions are related via the following partial Fourier transforms:
F t→ξ
κH (t0 , t0 − t) = h(t0 , t) Â / σH (t0 , ξ)
_ _
F t0 →ν F t0 →ν (2.11a)
² ²
Â F t→ξ / BH (ν, ξ)
e−i2πhν,ti SH (ν, t)
For κH being smooth and compactly supported, we apply the Fubini–Tonelli the-
orem, (2.11a) and Plancherel’s theorem to (2.8) to get
Z
(Hs)(t0 ) = κH (t0 , t)s(t) dµ(t)
Z Z
= SH (ν, t)ei2πhν,x−ti dνs(t0 − t) dt (2.11b)
Rd Rd
Z
= h(t0 , t)s(t0 − t) dt (2.11c)
Rd
Z
= s(ξ)ei2πht0 ,ξi dξ
σH (t0 , ξ)b (2.11d)
Rd
Z Z
= BH (ν − ξ, ξ)b s(ξ) dξei2πht0 ,νi dν. (2.11e)
Rd Rd
Note that the validity above extends to general Hilbert–Schmidt operators via a
density argument. Certainly, the convergence of the integrals is considered in the
L2 sense as generally done in L2 -Fourier analysis. In this case, the equalities above
hold for almost every t0 .
It follows naturally from (2.11) to view h(t0 , t) as the impulse response at t0
to an impulse at t0 − t and to view σH (t0 , ξ) as the frequency response at t0 to a
complex exponential with frequency ξ.
A Hilbert–Schmidt operator H is usually called underspread if its spreading
function is contained in a rectangle with area less than one and overspread oth-
erwise. Underspread operators have the important property that they are iden-
tifiable [KP06, PW06], which means that the operator H can be computed from
its response to a selected single input function. The most well-known example of
identifiability is the fact that linear time-invariant channels are completely char-
acterized by their action on a Dirac delta distribution, that is, by their impulse
response.
3 The channel model

An important and uniting property for radio and sonar communications in air and
water, respectively, is multipath propagation, which means that due to reflections
on different structures in the environment, the transmitted signal reaches the
11
receiver via a possibly infinite number of different wave propagation paths, as
illustrated in Figure 2 (see, for example, [Rap02]).
In Sections 3.1–3.2, we examine the multipath propagation model at some
depth under the standard assumptions that the electric field component at the
receiver is the superposition of the contributions from all signal paths leading there,
and that the action of the channel on a transmitted signal is the superposition of
the action on all complex exponentials in a Fourier expansion of the signal. For this
we use a standard and straightforward linear extension (H(u + iv) = Hu + iHv)
of H to complex valued functions. Initially, we do also allow the channel to be of
infinite lifelength, which is necessary for our class of modelled channels to include,
for example, the identity operator, which has a Dirac delta distribution spreading
function.
For communications applications, however, only finite lifelength channels are
important. We show in Section 3.3 that this subclass of channels can be completely
described by very smooth spreading functions with “essentially compact” support.
3.1 Single path frequency response

Most wireless communication channels change their characteristics slowly com-
pared to the rate at which transmission symbols are sent. Significant changes
either require a long time-period to evolve, or they are caused by abrupt changes
in the environment, for example, when a mobile telephone user drives into a tun-
nel. The standard countermeasure is to regularly make new estimates of the
channel. In OFDM based methods this is usually done by sending pilot sym-
bols, pilot tones or scattered pilots [GHS+ 01]. For a more general treatment,
see [GHS+ 01, KP06, LPW05, MMH+ 02, PW06].
Thus, from now on we shall only consider the short-time behaviour of the
channel during time intervals I that are short enough to assume a fixed collection
Figure 2: The transmitted signal reaches the receiver along a continuum of dif-
ferent signal paths. Each path P has time-varying length lP (t).
12
of signal paths with the length lP (t) of path P being a linear function of the time.
That is, we assume the length and prolongation-speed of each path to be such
that for some T0 ∈ I and all t ∈ I,
lP (t) = LP + VP · (t − T0 ) with |VP · (t − T0 )| ¿ LP . (3.1)
Physical constraints on the speed of antenna and reflecting object movements give
some upper bound Vmax for |VP |. We will assume Vmax to be smaller than the wave
propagation speed Vw , so that
Vw > Vmax ≥ |VP | for all paths P.
Hence, if a simple harmonic ei2πξt is sent along the path P without any attenuation
or perturbations, then the received signal would be
³ ´ ³³ ´ ´ ³ ´³ ´
l (t) V L −V T V L −VP T0
i2πξ t− P i2πξ 1− VP t− P V P 0 i2πξ 1− VP t− VP −V
e V w =e w w =e w w P = TtP MνP ξ ei2πξ(·)
(3.2a)
where the time and frequency shifts tP and νP ξ satisfy
· ¸
LP − VP T0 VP Vmax Vmax
tP = and νP = − ∈ − , ⊂ (−1, 1). (3.2b)
Vw − VP Vw Vw Vw
This mapping from (VP , LP ) to (νP , tP ) is invertible with inverse
VP = −νP Vw and LP = Vw (tP (1 + νP ) − νP T0 ) . (3.2c)
By (3.2b), (3.1) and (3.2c), there is a compact set K ⊂ (−1, 1) and some A ∈ R
such that (νP , tP ) ∈ K × [A, ∞) for all paths P.
Now fix some arbitrary (ν, t) ∈ K × [A, ∞) and some path P with frequency
response parameters (νP , tP ) = (ν, t) in (3.2a). The channel operator action on a
complex exponential sξ (t0 ) = ei2πξt0 consists of the following components:
1. A multiplication by a transmitter amplitude gain GT (P).

If we identify the path with the angular direction in which it leaves the
transmitter, then we can integrate (or sum) over all
R P and 2note that for
energy conservation reasons the total power gain |GT (P)| dP must be
finite.
2. The time-frequency shift by (νP ξ, tP ) that is given in (3.2a).
3. Attenuation with a factor3 Aξ (P) ∈ R that for free space transmission has
size O(L−2
P ) for large LP [Rap02, Reu74].
However, the decay is usually much faster and exponential decay O(e−aξ LP )
can be argued for if we assume some fixed minimum attenuation every time
3
We allow Aξ (P) to be negative to include potential sign-changes caused by reflections.
13
a signal is reflected [Str06]. Even for radio signals propagating through
the atmosphere without reflections (line-of-sight propagation, see Figure 2),
frequency selective absorption causes exponential decay with faster decay for
higher frequencies [Reu74, Section 2.1.7]. From this and (3.2c) we get that
for some aξ , C > 0,
|Aξ (P)| ≤ Ce−aξ LP ≤ Cξ e−αtP χ[A,∞) (tP ) (3.3)
with Cξ = sup|ν|<1 Ceaξ Vw νT0 < Ceaξ Vw |T0 | and α = inf ξ inf |ν|<1 aξ Vw (1 − ν) >
0.
4. Multiplication by a receiver amplitude
R gain GR (P), which for any kind of
practical use must also satisfy that |GR (P)|2 dP < ∞.
Altogether, the above steps add up to the following single path frequency response:
def Transmitter
sξ (·) = ei2πξ(·) −−−−−−−−→ GT (P)sξ
TF-shift (3.2a)
−−−−−−−−−−→ GT (P)TtP MνP ξ sξ
(3.4)
Attenuation
−−−−−−−−→ GT (P)Aξ (P)TtP MνP ξ sξ
Receiver
−−−−−→ GT (P)Aξ (P)GR (P)TtP MνP ξ sξ
Now set
def
Bξ (P) = Aξ (P)eαξ tP
for all paths P, so that by (3.3), |Bξ (P)| ≤ Cξ for all P. Further, we set Pν,t =
{P : (νP , tP ) = (ν, t) }.
As usual for electromagnetic waves, we expect the electric field component
measured at the receiver to be the superposition of the electric field components
received from the different paths P (and similarly for sonar waves), or written as
a formal integration4
Z ÃZ !
(Hsξ )(t0 ) = GT (P)Bξ (P)GR (P) dP e−αξ t (Tt Mνξ sξ )(t0 ) d(ν, t).
K×[A,∞) Pν,t
(3.5)
If we denote the inner integral
Z
def
rξ (ν, t) = GT (P)Bξ (P)GR (P) dP, (3.6a)
Pν,t
then it follows from the Hölder inequality, the bound |Bξ (P)| ≤ Cξ and items 1
and 4 above that
Z Z
|rξ (ν, t)| d(ν, t) ≤ |GT (P)Bξ (P)GR (P)| dP
K×[A,∞) ∪(ν,t)∈K×[A,∞) P(ν,t) (3.6b)
≤Cξ · kGT k2 · kGR k2 < ∞.
4
R
Here again, · · · dP is shorthand notation for the integration over the different angles in a
polar coordinate system.
14
Both for the actual transmitted real-valued signals and for our linear extension
of H to complex-valued signals, the gain and attenuation factors are all real-
valued. We shall however, without any extra effort, allow for complex-valued r
in our mathematical model. Moreover, for inclusion of some important idealized
borderline cases such as r being a Dirac delta distribution, and for avoiding some
computational distribution theory technicalities, we choose to model the integrals
def R
ρξ (U × V ) = U ×V rξ (ν, t) d(ν, t) over sets U × V ⊆ K × [A, ∞) to be a complex
Borel measure ρξ with finite total variation, that is
Z
ρξ (U × V ) = dρξ (ν, t) with |ρξ | (K × [A, ∞)) < ∞. (3.7a)
U ×V
Thus, in this mathematical model, (3.5) takes the form

Z
(Hsξ )(t0 ) = e−αξ t (Tt Mνξ sξ )(t0 ) dρξ (ν, t). (3.7b)
K×[A,∞)
Note that this model includes, for example, Dirac measures and thus also the
identity operator.
3.2 Narrowband signals

We shall call the transmitted signal s narrowband if sb is well-localized enough to
justify the approximations
νξ ≈ νξ0 , e−αξ t ≈ e−αξ0 t and ρξ ≈ ρξ0 (3.8)
in the computations leading to (3.9b) below. We will primarily assume this nar-
rowband assumption to hold for the same ξ0 in the entire transmission frequency
band. In the remark on page 29, we will show that this assumption holds true for
some radio communications examples and discuss a refined model with different ξ0
for different basis functions that is necessary in underwater sonar communications.
Suppose now that the physical channel has the property that its action on a
signal s is the superposition of its action on each complex exponential in a Fourier
expansion of s, that is,
Z Z
i2πξ(·)
Hs(·) = H sb(ξ)e dξ = sb(ξ)Hei2πξ(·) dξ. (3.9a)
R R
Then, at least for bandlimited and thus continuous narrowband L1 -signals s, we
15
can apply (3.7b), (3.8) and the Fubini–Tonelli theorem to obtain
Z Z
Hs(t0 ) = sb(ξ) e−αξ t ei2πνξ(t0 −t) ei2πξ(t0 −t) dρξ (ν, t) dξ
R K×[A,∞)
Z Z
≈ sb(ξ) e−αξ0 t ei2πνξ0 (t0 −t) ei2πξ(t0 −t) dρξ0 (ν, t) dξ
R K×[A,∞)
Z Z
−αξ0 t i2πνξ0 (t0 −t)
= e e sb(ξ)ei2πξ(t0 −t) dξ dρξ0 (ν, t)
K×[A,∞) R
Z
= e−αξ0 t ei2πνξ0 (t0 −t) s(t0 − t) dρξ0 (ν, t).
K×[A,∞)
We use the last expression as definition of our mathematical model

Z
def
Hs(t0 ) = e−αξ0 t (Tt Mνξ0 s)(t0 ) dρξ0 (ν, t). (3.9b)
K×[A,∞)
If s ∈ L2 (R), then by (3.9b) and the Minkowski integral inequality

°Z °
° °
kHsk2 = ° ° e−α ξ0 t
(T t M νξ 0 s)(·) dρξ0 (ν, t)°
°
K×[A,∞) 2
Z
≤ e−αξ0 t kTt Mνξ0 sk2 d |ρξ0 | (ν, t)
K×[A,∞)
≤ |ρξ0 | (K × [A, ∞))e−αξ0 A ksk2 .
Hence equation (3.9b) defines a bounded linear mapping H : L2 (R) → L2 (R). If,
in addition, ρξ0 is absolutely continuous with respect to the Lebesgue measure,
then we can write dρξ0 (ν, t) = rξ0 (ν, t) d(ν, t) where rξ0 is the function in (3.6),
which equals the inner integral in our physical model (3.5). In practice, rξ0 will be
bounded and thus also in L2 (R). Application of this to (3.9b) and the substitution
ν 0 = νξ0 gives
Z
Hs(·) = SH (ν 0 , t)(Tt Mν 0 s)(·) d(ν 0 , t), K 0 ⊆ (−ξ0 , ξ0 ) compact,
0
K ×[A,∞)
µ 0 ¶
0 def 1 ν
SH (ν , t) = rξ0 , t e−αξ0 t χK 0 ×[A,∞) (ν 0 , t) and rξ0 ∈ L1 ∩ L2 (R × R).
ξ0 ξ0
(3.9c)
Alternatively we can use (3.9b) as the definition of a larger space of opera-

tors on a smaller space of functions by allowing a larger subclass of the complex
Borel measures, such as, the class of all complex Borel measures ρξ0 for which the
mapping Z
ϕ 7→ ϕ(ν, t)e−αξ0 t dρξ0 (ν, t)
K×[A,∞)
16
defines a linear bounded functional SH on S0 (K × [A, ∞)). Then for all functions
s for which the mapping
(νξ0 , t) 7→ (Tt Mνξ0 s)(t0 ) ∈ S0 (K 0 × [A, ∞)) for all t0 , (3.10)
is well-defined, it follows that (3.9b) is pointwise well-defined for all t0 . Conse-

quently (3.9b) can be interpreted as the application of a functional SH ∈ S00 to
the test function (3.10), or, with the usual formal integral notation,
Z
Hs(·) = SH (ν 0 , t)(Tt Mν 0 s)(·) d(ν 0 , t), SH ∈ S00 (K 0 × [A, ∞)).
K 0 ×[A,∞)
Since the space S00 (R × R) includes Dirac delta distributions, this model includes
important idealized borderline cases such as the following:
Line-of-sight transmission: SH = aδν0 ,t0 , a Dirac distribution at (ν0 , t0 ) represent-
ing a time- and Doppler-shift with attenuation a.
Time-invariant systems: h(x, t) = h(0, t) and SH (ν, t) = h(t)δ0 (ν).

Moreover, S00 excludes derivatives of Dirac distributions, which can be used
to avoid complex-valued Hs with no physical meaning [Ric03, Sec. 3.1.1]. Fur-
ther, S0 is the smallest Banach space of test functions with some useful properties
like invariance under time-frequency shifts [Grö00, p.253], thus allowing for time-
frequency analysis on its dual S00 which is, in that particular sense, the largest
possible Banach space of tempered distributions that is useful for time-frequency
analysis. One more motivation for considering spreading functions in S00 is that
Hilbert–Schmidt operators are compact, hence, they exclude invertible operators,
such as the example SH = aδν0 ,t0 above, and small perturbations of invertible
operators, which are useful in the theory of radar identification and in some mo-
bile communication applications. For results using a Banach space setup, see for
example [PW06, Str06].
Hence it may come as a surprise that we will show in Section 3.3 that the
Hilbert–Schmidt model (3.9c) is a natural choice for wireless communications chan-
nels. This also justifies our use of this model in the remaining paper.
3.3 Finite lifelength channels

We shall show that wireless communications channels can be modelled well by
well-localized C (∞) -spreading functions. This fact allows for a minimal use of dis-
tribution theory in our analysis and adds simplicity to proofs, results and software
development both in this report and for mathematical and numerical analysis of
wireless communications channels in general.
In the following, we will assume that the bifrequency function
BH (ν, ·) is compactly supported. (3.11)
17
Â F t→ξ / σ (t , ξ) C (∞) 3 h(t0 , t) Â
F t→ξ
/ σH (t0 , ξ)
C (∞) 3 h(t0 , t) H 0
_ _ _ _
F t0 →ν F t0 →ν F t0 →ν F t0 →ν
² ² ² ²
F t→ξ F t→ξ
e−i2πhν,ti SH (ν, t) Â / BH (ν, ξ) C (∞) 3 e−i2πhν,ti SH (ν, t) Â / BH (ν, ξ) ∈ C (∞)
(a) (b)
Figure 3: (a) Bandlimiting properties of the physical channel provides us with

compactly supported BH . Thus h ∈ C (∞) , but SH may be a tempered distribution,
for example, if h(t0 , t) = h(0, t). (b) Our restriction to finite lifelength channels
gives a well localized SH ∈ C (∞) . In (a) and (b) we denote with underlining
subexponential decay and compact support.
This is not strictly true in general, but for communications applications we only
need to model the channel response to a family of signals s that all are bandlimited
to some common frequency band Ω1 . Hence it is clear from (2.11e) that Hs will
only depend on the restriction of BH to R × Ω1 , so that we are free to set BH
equal to zero outside R × Ω2 for an arbitrary and compact set Ω2 ⊇ Ω1 . Moreover,
recall from (3.9c) that SH (·, t) has compact support as well. Hence, combining
this with (2.11a) and (3.11), we see that the distribution BH (ν, ξ) is compactly
supported. Consequently, since the bifrequency function satisfies BH = b h, we
have that h ∈ C (∞) (R2 ). We summarize a straightforward generalization of these
properties to functions on Rd in Figure 3 (a).
Since we are modelling the short time input-output relationships of a channel,
any useful communications system must be constructed to be independent of the
properties of h outside some compact set Kh . Thus, we are free to multiply h
with a compactly supported function w ∈ C (∞) (R2d ) such that w = 1 on Kh
and w b is subexponentially decaying (as described in Section 2.2). It is easy to
check that it follows from this and from the compact support of b h that also the
c
convolution wh = w b
b ∗ h is subexponentially decaying. Now let H1 be the Hilbert–
Schmidt operator with time-varying impulse response wh ∈ C (∞) . From the fact
that the space of Schwartz functions is invariant under partial Fourier transforms,
it follows that also SH1 ∈ C (∞) . This gives an operator with system function
properties that we summarize in the multivariate case in Figure 3 (b). We will
also assume that w is chosen to be “wide and smooth enough” so that the smooth
cut-off of SH (ν, ·) only deletes a very small-amplitude and negligible part of its
exponential tail, and so that also the “blurring out” of the compact support of
SH (·, t) to subexponential decay has a rather small impact on the shape of SH ,
which therefore can be expected to resemble those given in the figures of Section 6.
The following sections are devoted to Hilbert–Schmidt operators having exactly
18
the properties that are summarized in Figure 3 (b). All results apply directly to
the communication channels described in this section as long as the narrowband
assumption of Section 3.2 holds for the entire frequency band of the transmitted
signals. We describe in a remark on page 29 how refined versions of our results
can be applied to wideband signals as long as the attenuation factor Aξ of (3.3)
is roughly frequency independent within the transmission frequency band.
Remark. The exponential decay of SH (ν, ·) only affects the just mentioned
shape of SH1 . In general, the reasoning in this and the following sections hold
also with exponential and subexponential decays replaced by other forms of rapid
decay.
4 Discretization of the channel model

In this section we derive finite sum formulas for the computation of the matrix rep-
resentation of the coefficient mapping G in (2.7) for classes of multivariate Hilbert–
Schmidt operators H that satisfy the properties summarized in Figure 3 (b).
For the series expansions Pin Gabor frames (gq,r ) and (γq ,r ) in Section 2.3, H
0 0
2
maps any L -function s = q,r∈Zd hs, geq,r i gq,r to
* +
X X X
Hs = hHs, γq0 ,r0 i γ
eq0 ,r0 = H hs, geq,r i gq,r , γq0 ,r0 γeq0 ,r0
q 0 ,r 0 ∈Zd q 0 ,r 0 ∈Zd q,r∈Zd
  (4.1)
X X
=  hs, geq,r i hHgq,r , γq0 ,r0 i γ
eq0 ,r0 ,
q 0 ,r 0 ∈Zd q,r∈Zd
where convergence in L2 (Rd ) follows from the continuity of H and of the inner
product. Thus the coefficient mapping G in (2.7) is
 
X
G : (hs, geq,r i)q,r 7→  hs, geq,r i hHgq,r , γq0 ,r0 i
q,r∈Zd
q 0 ,r 0
with biinfinite matrix representation
Gq0 ,r0 ;q,r = hHgq,r , γq0 ,r0 i , (4.2)
and with 2d-dimensional indices (q 0 , r 0 ) and (q, r) for rows and columns respec-
tively.
For communications applications with Q carrier frequencies, at least Q samples
of every received symbol are needed in the receiver. Thus a hasty and naive
approach to computing the matrix elements could start with a Q × Q matrix
representation of H for computing the samples of Hgq,r . If every R neighbouring
symbols have overlapping ²-essential support, then we need to compute (RQ)2
19
matrix elements hHgq,r , γq0 ,r0 i, which, with this approach, would require R2 ·O(Q5 )
arithmetic operations with Q typically being at least of the size 256–1024 in radio
communications, and with R = 4 for ² = 10−6 and the optimally well-localized
Gaussian windows that we shall use for our example applications in Section 6 (see
Figure 12). This is a quite demanding task, so there is a clear need for the more
efficient formulas and algorithms that we shall derive in the following sections.
We begin in Section 4.1 with motivating the use of well time-frequency localized
Gabor bases. Then, in Section 4.2 we introduce the tools and notation that we find
most suitable for discretization of Hilbert–Schmidt operators with the properties
summarized in Figure 3 (b). Finally, in Section 4.3 we use these tools for deriving
more efficient formulas for computing the matrix elements under a minimum of
further assumptions, simplifications or approximations.
Remark. If (gq,r ) is a frame for its closed linear span but not a Riesz basis,
then the coefficient subspace Cge of (2.7) will be a proper subspace of l2 . Thus an
orthogonal projection to Cge in the receiver would add some stability against noise
and perturbations, but at the cost of an additional requirement on the transmitter
to only use coefficient sequences in Cge.
For systems in use today, it is typical to instead use Gabor Riesz bases and arbi-
trary coefficient sequences with coefficients chosen from some standard coefficient
constellation, such as quadrature amplitude modulation (QAM), with maximum
coefficient amplitude given by transmission power regulations.
4.1 Gabor bases for near-diagonalization

Although (4.1) holds for any pair of frames (gi ) and (γj ), Gabor Riesz bases
are traditionally used for communications applications. Gabor systems give good
diagonalization of G, especially if all basis elements are narrowband (see below
and [KPZ02]).
We call an operator H time-invariant if it commutes with the translation
operator Tt0 for any t0 , that is, if Tt0 H = HTt0 . Consequently, the impulse
response h of H satisfy h(t0 , t) = h(0, t) and (2.11c) becomes a convolution h(0, ·)∗
s. Then the family of complex exponentials ei2πhξ,·i are “eigenfunctions” in the
sense that for inputs of the form s(·) = ei2πhξ,·i χ[0,L] (·), there is some complex
scalar λξ such that if supp h = [0, Lh ], then Hs = λξ s in the interval [Lh , L],
Figure 4: Convolution of a compactly supported complex exponential s(t) =

χ[0,L] (t)ei2πξt with a function h = hχ[0,Lh ] will on the interval [Lh , L] equal the
same complex exponential multiplied with a complex scalar.
20
as illustrated in Figure 4. Thus G is easily diagonalized by using Gabor windows
g = χ[0,L] , γ = χ[Lh ,L] and lattice constants such that the resulting Gabor systems
(gk,l ) and (γk,l ) are biorthogonal bases for their respective span. This trick is
used in wireline communications, where the smaller support of γ is obtained by
removing a guard interval (often called cyclic prefix) from g. See, for example,
[Gri02, Section 2.3] for more details and further references.
In wireless communications, we can at most hope for approximate diagonal-
ization of the channel operator due to its time varying nature. In general, two
different time-varying operators do not commute, so both cannot be diagonalized
with the same choice of bases. Thus, diagonalization is usually only possible in the
following sense: Typically, (Hgq,r ) is a finite³ and linearly
´ independent sequence,
and thus a Riesz basis with some dual basis Hg ^ q,r . In fact, since the computa-
tions in (4.1) are not restricted to Gabor frames, we can set γq0 ,r0 = Hg ^ q,r , which
^
gives true diagonalization of (4.2), but in general, Hg q,r will not be a Gabor frame
or have any other simple structure that enables efficient computation of all Hg ^
D E q,r
and all the inner products Hgq,r , Hg ^ q 0 ,r 0 .

Thus, for computational complexity to meet practical restrictions we will use
“almost dual” Gabor bases (gq,r ) and (γq0 ,r0 ), such as the Gabor bases proposed
in [MSG+ 05]. We are primarily interested in bases that are good candidates for
providing low intersymbol and interchannel interference (ISI and ICI). As proposed
in [MSG+ 05], we expect excellent joint time-frequency concentration of g and γ
to be the most important requirement for achieving that goal.
We will strive to customize decisions like the choice of discretization methods
for channels with the system function properties summarized in Figure 3 (b) and
for windows g and γ that are bandlimited and decay “fast enough” to be repre-
sented by a finite (and reasonably small) number of Nyquist frequency samples.
This is the topic of sections 4.2 and 4.3.
4.2 Discretization tools

Recall first from Section 3.3 and Figure 3 (b) that it is justified to work with
spreading functions that can be truncated to compact support with negligible
truncation errors.
Then Proposition 2 below shows that the functions Hgq,r and Hg \ q,r are well
localized around the lattice points of a lattice with the same lattice constants as
the lattice of the Gabor basis (gq,r ). This suggests to choose (γq0 ,r0 ) to be a Gabor
basis with the same lattice constants as (gq,r ). For this scenario, propositions 3–5
in Section 4.3 provide us with formulas that allow an efficient computation of the
matrix elements hHgq,r , γq0 ,r0 i.
Proposition 2. Suppose that H is a Hilbert–Schmidt operator on L2 (Rd ),
supp SH ⊆ Iωc ,ω × [a, b] and g ∈ L2 (Rd ). (4.3)
21
(a) If supp f ⊆ [A, B], then supp Hf ⊆ [A + a, B + b].
(b) If supp fb ⊆ IC f ,B f , then supp Hg

c ⊆ IC +ωc ,B +ω .
f f
Proof. (a) follows immediately from (2.11b). (b): For SH ∈ L2 (Rd × Rd ) and
f ∈ L2 (Rd ), a repeated use of the Hölder inequality gives that
µZ Z ¶2
2
|(Hf )(t0 )| ≤ |SH (ν, t)| dν |f (t0 − t)| dt
Rd Rd
Z ÃZ !2
≤ kf k2L2 (Rd ) 1 · |SH (ν, t)| dν dt
Rd Iωc ,ω
≤ kf k2L2 (Rd ) |ω| · kSH k2L2 (Rd ×Rd ) < ∞.
Similarly, for continuous f with compact support, Fubini–Tonelli gives that

Z Z Z
d
Hf (ξ) = SH (ν, t)ei2πhν,x−ti dν f (t0 − t) dt e−i2πhx,ξi dt0
R d [α,β] Iωc ,ω
Z Z Z
= f (t0 − t)e−i2πhξ−ν,x−ti dt0 SH (ν, t)e−i2πht,ξi dt dν
I [α,β] R d
Z ωc ,ω Z (4.4)
= fb(ξ − ν)SH (ν, t)e−i2πht,ξi dt dν
I [α,β]
Z ωc ,ω Z
= b \
f (ξ − ν)SH (ν, ·)(ξ) dν = fb(ξ − ν)BH (ν, ξ) dν.
Iωc ,ω Iωc ,ω
By the Hölder inequality, the L2 -norm squared of the righthand side equals
Z ¯¯Z ¯2
¯ ° °2 Z Z
¯ ¯ ° °
¯ fb(ξ − ν)BH (ν, ξ) dν ¯ dξ ≤ °fb° |BH (ν, ξ)|2 dν dξ.
Rd ¯ Iωc ,ω ¯ 2 Rd Iω ,ω
c
Since H is bounded and BH ∈ L2 (Rd × Rd ), both the lefthand and the righthand
side of (4.4) depend continuously on f , so that (4.4) holds for all f ∈ L2 (Rd )
by standard density arguments. Thus (4.3) implies that if supp fb ⊆ IC f ,B f , then
d ⊆ IC +ωc ,B +ω .
supp Hf f f
We will repeatedly apply the following special case of the Poisson Summation
formula [Grö00, Kat04] to the functions gq,r , Hgq,r and γq0 ,r0 :
Theorem 2 (Sampling theorems). Suppose that f ∈ L2 (Rd ) has a Fourier trans-
form with support supp fb ⊆ ICf ,Bf . Then fb ∈ L1 (Rd ) ∩ L2 (Rd ), f is continuous
(modulo modifications on a set of measure 0) and for all ξ ∈ R cd ,
X 1
fb(ξ) = |Tf | χICf ,Bf (ξ) f (kTf )e−i2πhξ,kTf i , Tf = . (4.5a)
k∈Zd
Bf
22
Equivalently (via (2.1)), (4.5a) is also known as the Nyquist sampling theorem
X
f (t) = |Tf | f (kTf )ei2πhCf ,t−kTf i sincBf (t − kTf ). (4.5b)
k∈Zd
For a ∈ Zd , some b, Ω ∈ Rd+ and f equal to gq,r or γq0 ,r0 , we will consider
def
time-frequency shifts fq,r = TraTg MqbΩ f , for which
def 1
Tg = , supp fd
q,r , ⊆ ICf +qbΩ,Bf , (4.6a)
Ω
and the Nyquist frequency samples are
fq,r (kTf ) = f (kTf − raTg )ei2πhkTf −raTg ,qbΩi . (4.6b)
If Tg = Tf , then (4.6b) becomes
fq,r (kTg ) = f ((k − ra)Tg )ei2πhk−ra,qbi .
If Tg 6= Tf , then we can instead compute the samples f (kTf − raTg ) by apply-

ing (4.5b), which will be a finite sum formula for those f that we will consider.
Proposition 2, Theorem 2, (4.6a) and (4.6b) are the basic tools that we will use
to compute the matrix elements hHgq,r , γq0 ,r0 i. We will apply these to bandlimited
g and γ that have only a finite number of nonzero samples in the Nyquist recon-
struction formula (4.5b). These samples should be chosen so that g and γ decay
fast enough to allow truncation to compact support with both the maximum and
the L2 -norm of the truncation error being less than some ² well below the overall
noise level of the application at hand. In Figure 7 on page 36, we show by example
how such window functions can be constructed. For the scenario that we summa-
rized in Figure 3 (b), we can, in practice, regard SH to be compactly supported
and conclude from Proposition 2 that Hgq,r and Hg \ q,r inherits the good localiza-
tion properties of gq,r , gd
q,r and SH . Hence, with negligible truncation errors, we
can assume also Hgq,r to be bandlimited, to be fully described by a finite number
of Nyquist frequency samples and to have ²-essential time-frequency support given
by Proposition 2. (In Section 5 we choose ² = 10−6 .) Throughout the paper, we
shall use the following notation for these compact supports and finite index sets.
def 1 def
1
supp gb ⊆IΩc ,Ω , Tg = Ω
, Tγ = Ω+ω , (4.7a)
supp SH ⊆Iωc ,ω × IC,L , c ⊆ supp γ
supp Hg b ⊆ IΩc +ωc ,Ω+ω , (4.7b)
d
K, M ⊂Z , |K| < ∞, |M| < ∞ and (4.7c)
g(mTg ) =γ(kTγ ) = (Hg)(kTγ ) = 0 for k ∈ Zd \ K and m ∈ Zd \ M. (4.7d)
The corresponding supports and index sets for gq,r , Hgq,r and γq,r are easily
computed from (4.6) above.
23
4.3 Computing the channel matrix
For spreading functions and window functions having the properties, supports and
index sets given in (4.7) and Figure 3 (b), the following proposition provides us
with a finite sum formula for computing the matrix elements hHgq,r , γq0 ,r0 i from
samples of Hgq,r and γq0 ,r0 :
Proposition 3. Let u, v ∈ L2 (Rd ) be functions with compact, overlapping fre-

quency supports satisfying
def
supp u
b ⊆ ICu ,B , supp vb ⊆ ICv ,B and ICuv ,Buv = ICu ,B ∩ ICv ,B 6= ∅.
For T = B1 and k ∈ Zd , suppose that u(kT ) 6= 0 and v(kT ) 6= 0 only for k in

some finite index sets Iu and Iv , respectively. Denote with vbpf the inverse Fourier
transform of the restriction (bandpass filtering) of vb to the support of u
b, that is,
def
vd
bpf (ξ) = v
b(ξ)χICuv ,Buv (ξ).
Then X
hu, viL2 (Rd ) = |T | u(kT )vbpf (kT ). (4.8a)
k∈Iu
Further, vbpf = v if Cu = Cv and, otherwise,

X 0
vbpf (kT ) = |T | v(k0 T )ei2πhCuv ,(k−k )T i sincBuv ((k − k0 )T ). (4.8b)
k0 ∈Iv
Proof. (4.8a) follows from the Plancherel theorem and Parseval’s identity:
X
hu, viL2 (Rd ) = hb
u, vbiL2 (Rd ) = hb
u, vd
bpf iL2 (IC ,B ) = |T | u(kT )vbpf (kT ).
v
k∈Iu
Moreover, the sampling theorem (4.5a) and (2.1) implies that

Z
vbpf (kT ) = vb(ξ)χICuv ,Buv (ξ)ei2πhξ,kT i dξ
Rd
X Z
0 0
= |T | v(k T ) ei2πhξ,(k−k )T i dξ
k0 ∈Iv ICuv ,Buv
X 0
= |T | v(k0 T )ei2πhCuv ,(k−k )T i sincBuv ((k − k0 )T ).
k0 ∈Iv
Using (4.7) and (4.8a), we get, for example, that

X
hHg, γiL2 (Rd ) = |Tγ | (Hg)(kTγ )γ(kTγ )
k∈K
24
and it remains only to derive efficient formulas for the computation of the samples
(Hg)(kTγ ). For this, we derive the finite sum formula (4.12) below with nonzero
terms only for m ∈ M when f = g (and similarly for f = gq,r )5 . It is clear that
the following proof holds for arbitrary samples Hf (t0 ), but by setting t0 = kTγ
we get our results in exactly the form that we shall need later on.
Proposition 4. Suppose that Cf , Ω ∈ Rd , f ∈ L2 (Rd ), supp fb ⊆ ICf ,Ω and
def 1
(f (mTf ))m∈Zd ∈ l1 (Zd ), with Tf = . (4.9)
Ω
Let H be a Hilbert–Schmidt operator with spreading function
SH ∈ C (∞) ∩ L2 (Rd ) and supp SH ⊆ Iωc ,ω × IC,L . (4.10a)

C ,Ω
Define SH f to be the convolution
³ ´
SH f (ν, t0 ) = SH (ν, ·) ∗ (ei2πhCf +ν,·i sincΩ (·)) (t0 ).
C ,Ω def
(4.10b)
Then for all k ∈ Zd and Tγ ∈ Rd+ , we have
(Hf )(kTγ ) =
Z X
C ,Ω
= |Tf | SH f (rTf + k(Tγ − Tf ), ν)(TrTf Mν f )(kTf ) dν (4.11)
Rd
r∈Zd
X ³ ´
C ,Ω b
= |Tf | f (mTf ) SH f (·, kTγ − mTf ) (−mTf ). (4.12)
m∈Zd
Proof. From (2.11b) and the Fubini–Tonelli theorem,

Z Z
(Hf )(t0 ) = SH (ν, t)f (t0 − t)ei2πhν,t0 −ti dt dν
d d
ZR R Z
¡ i2πhν,·i
¢ def
= SH (ν, ·) ∗ (f (·)e ) (t0 ) dν = fν (t0 ) dν, (4.13)
Rd Rd
where the Nyquist sampling theorem (4.5b) and then (2.1) gives that
X ³ ´
fν (t0 ) = |Tf | f (mTf ) SH (ν, ·) ∗ (ei2πhCf ,·−mTf i sincΩ (· − mTf )ei2πhν,·i ) (t0 )
m∈Zd
X ³ ´
= |Tf | f (mTf )ei2πhν,mTf i SH (ν, ·) ∗ (ei2πhCf +ν,·i sincΩ (·)) (t0 − mTf ).
m∈Zd
5
In (4.12) it is necessary for f to provide the finite summation index, since there can be
C ,Ω
a finite number of nonzero SH f (·, kTγ − mTg ) and (Hf )(kTγ ) (without alias effects due to
T 1
undersampling) only if Tγg ∈ Zd and Tγ = Ω+ω , that is, only if ω = 0.
25
Insertion of this in (4.13) gives
(Hf )(kTγ ) =
Z X
f (mTf )ei2πhν,mTf i SH f (ν, kTγ − mTf ) dν
C ,Ω
= |Tf | (4.14)
Rd
m∈Zd
Z X
f ((k − r)Tf )ei2πhν,(k−r)Tf i SH f (ν, rTf + k(Tγ − Tf )) dν.
C ,Ω
= |Tf |
Rd
r∈Zd
C ,Ω
This proves (4.11). By (4.10a) and (4.10b), SH f is bounded. Hence (4.9), (4.14)
and the Fubini–Tonelli theorem provide
X Z
ei2πhν,mTf i SH f (ν, kTγ − mTf ) dν
C ,Ω
(Hf )(kTγ ) = |Tf | f (mTf )
Iωc ,ω
m∈Zd
X ³ ´
C ,Ω b
= |Tf | f (mTf ) SH f (·, kTγ − mTf ) (−mTf ),
m∈Zd
which proves (4.12) and concludes the proof.

Next, we shall derive the finite sum formula (4.16) for computing the samples
Ωc +qbΩ,Ω
of (SH (·, t))b that are needed in (4.12) to compute (Hgq,r )(kTγ ) for a
Gabor system
def
gq,r = TraTg MqbΩ g, (q, r) ∈ Q × R
for a finite index set Q × R ⊆ Z2d . Equations (4.10b), (4.7a) and (4.7b) imply
Ωc +qbΩ,Ω
that SH (ν, ·) is obtained from bandpass filtering SH (ν, ·) to the frequency
def def
Figure 5: For Gabor frames gq,r = TraTg MqbΩ g with q ∈ Q = [q0 , q1 ] ∩ Zd , the
function |d gq,r | depends only on q. Proposition 4 shows how to compute all Hgq,r
only from the restriction of S\ H (ν, ·) to any interval IΩ0c ,Ω0 containing the supports
0
of all gd
q,r (· − ν) for all ν ∈ Iω c ,ω , as is illustrated here for minimal Ω and d = 1.
26
band
Ωc +qbΩ,Ω
supp gd
q,r (· − ν) = IΩc +qbΩ+ν,Ω , where ν ∈ supp SH (·, t) ⊆ Iωc ,ω .
Hence, we can compute all SH Ωc +qbΩ,Ω

from the restriction of S\
H (ν, ·) to the small-
est interval IΩ0c ,Ω0 containing the union of the supports of all gd q,r (· − ν) for all
ν ∈ Iωc ,ω , as illustrated for d = 1 in Figure 5. This means that if Q = [q0 , q1 ]∩Zd ,
then
def q0 + q1 def
Ω0c = Ωc + bΩ + ω c and Ω0 = (q1 − q0 )bΩ + Ω + ω. (4.15a)
2
This observation allows for a two-step discretization based on the system func-
tion properties summarized in Figure 3 (b), namely:
First, recall from Figure 3 (b) that h has compact support which is contained
in some “box” IC0 ,L0 × IC,L , where IC,L contains the support of SH (ν, ·) for all ν.
Recall also that due to the smooth truncation in Section 3.3, h(·, t) is supported
in an interval IC0 ,L0 , which is large enough not to affect the channel output in
the time interval under consideration. This together with the impulse response
integral representation (2.11c) of H gives that supp h(·, t) ⊆ IC0 ,L0 must be large
enough to contain the support of Hgq,r for all r ∈ R. Using element-wise addition
of sets, this condition has the form
supp g + RaTg + IC,L ⊆ IC0 ,L0 . (4.15b)
It follows that
© ª
supp S\
H (·, t)(t0 ) = supp h(t−t0 , t) = (t0 , t) ∈ R
2d
: t ∈ IC,L and t0 ∈ It−C0 ,L0 .
Due to this compact support set and the subexponential decay of SH (·, t), we can
apply and truncate the Nyquist sampling theorem (4.5b) to
X
SH (ν, t) = |ω 0 | SH (nω 0 , t)ei2πht−C0 ,ν−nω0 i sincL0 (ν − nω 0 ),
n∈N
(4.15c)
1
def
|N | <∞ and ω0 = .
L0
Second, we shall use that only the restriction of S\H (ν, ·) = BH (ν, · − ν) to
IΩc ,Ω is of importance for our computations. Since SH (ν, ·) has compact support,
0 0
we know from Section 2.2 that there is a smooth truncation of S\ H (ν, ·) to a C

(∞)
function S Ω\
00
Ωc 00 (ν, ·) such that
\
Ω
SΩ
00
\
00 (ν, ·) = SH (ν, ·) in IΩ0c ,Ω0 ,
\
Ω
supp SΩ
00
00 (ν, ·) ⊆ IΩ00 00 (4.15d)
c c c ,Ω
with IΩ00c ,Ω00 compact and such that

Ω 00
SΩ00 (ν, ·) decays subexponentially.
c
(4.15e)
27
Hence we can apply the sampling theorem (4.5a) to SΩ \
Ω 00
00 (ν, ·) and truncate it with
c
no or negligible error to a finite sum
X 00
Ω
SΩ\00 00
00 (ν, ·)(ξ) = |T | χI 00 00 (ξ) S Ω 00 −i2πhξ,pT 00 i
00 (ν, pT )e ,
c Ω ,Ω c
Ωc
p∈P
(4.15f)
1 00 def
|P| <∞ and T = 00 .
Ω
We are now ready to state our final proposition.
Proposition 5. For the functions and supports defined in (4.7), (4.10b) and (4.15),
def
set Ωc,q = Ωc + qbΩ for all q ∈ Q. Then
Ω\ ,Ω
X 00
SH c,q (·, t)(t0 ) = |ω 0 T 00 | χIC0 ,L0 (t − t0 ) ei2πhΩc,q ,t−pT i sincΩ (t − pT 00 )×
p∈P
X
Ω00 00 i2πht−t0 −pT 00 ,nω 0 i
× SΩ00 (nω 0 , pT )e
c
. (4.16)
n∈N
Proof. Note first that if ξ ∈ IΩc,q +ν,Ω , then for all nω 0 ∈ Iωc ,ω ,
ξ − (ν − nω 0 ) ∈ IΩc,q +nω0 ,Ω ⊆ IΩ0c ,Ω0 ,
so that (4.15d) implies that
S\ \
Ω00
H (ν, ·)(ξ − (ν − nω 0 ))χIΩc,q +ν,Ω (ξ) = SΩ00 (ν, ·)(ξ − (ν − nω 0 ))χIΩc,q +ν,Ω (ξ).
c
Hence, it follows from (4.15), the definition (4.10b) and (2.1) that for all ν ∈ Iωc ,ω ,
Ω\
SH c,q (ν, ·)(ξ) = S\
,Ω
H (ν, ·)(ξ)χIΩc,q +ν,Ω (ξ)
X
= |ω 0 | SH\
(nω 0 , ·)(ξ − (ν − nω 0 ))χIΩc,q +ν,Ω (ξ)e−i2πhC0 ,ν−nω0 i sincL0 (ν − nω 0 ) =
n∈N
XX 00 00
= |ω 0 T 00 | Ω
SΩ 00 −i2πhξ−(ν−nω 0 ),pT i
00 (nω 0 , pT )e
c
χIΩc,q +ν,Ω (ξ)×
n∈N p∈P
× e−i2πhC0 ,ν−nω0 i sincL0 (ν − nω 0 )
Inverse Fourier transformation and (2.1) again gives

X X 00 Z
Ωc,q ,Ω 00 Ω 00 00 i
SH (ν, t) = |ω 0 T | SΩ00c (nω 0 , pT ) ei2πhξ,t−pT dξ×
n∈N p∈P IΩc,q +ν,Ω
00
× ei2πhpT −C0 ,ν−nω0 i sincL0 (ν − nω 0 )
X X 00
00 i2πhΩc,q +ν,t−pT 00 i
= |ω 0 T 00 | Ω
SΩ 00 (nω 0 , pT )e
c
×
n∈N p∈P
00 −C 00 −C
× sincΩ (t − pT 00 )e−i2πhpT 0 ,nω 0 i
ei2πhpT 0 ,νi
sincL0 (ν − nω 0 ).
28
Since sincL0 (· − nω 0 ) is the inverse Fourier transform of χI0,L0 (·)e−i2πh·,nω0 i ,
Ω\ ,Ω
SH c,q (·, t)(t0 )
X X 00
00 i2πhΩc,q ,t−pT 00 i
= |ω 0 T 00 | Ω
SΩ00 (nω 0 , pT )e
c
sincΩ (t − pT 00 )×
n∈N p∈P
Z
00 −C
× e−i2πhν,t0 +C0 −ti sincL0 (ν − nω 0 ) dνe−i2πhpT 0 ,nω 0 i
Rd
XX 00 00 00 −C
= |ω 0 T 00 | Ω
SΩ 00 i2πhΩc,q ,t−pT i
00 (nω 0 , pT )e
c
sincΩ (t − pT 00 )e−i2πhpT 0 ,nω 0 i
×
n∈N p∈P
× χI0,L0 (t0 + C0 − t)e−i2πht0 +C0 −t,nω0 i .

This proves (4.16).
Remark. For physical wireless communications channels, we deduced a spread-

ing function integral representation (3.9c) of the channel which is valid for func-
tions with frequency localization “near ξ0 ”. For the OFDM and Satellite commu-
nications examples in Section 6, this assumption holds for the entire frequency
band IΩ0c ,Ω0 , since Ω0 /Ω0c is of the size 3 · 10−3 and 2 · 10−3 , respectively. Then,
Ω00 00
the coefficients SΩ 00 (nω 0 , pT ) above characterize the channel throughout this
c
frequency band. In the underwater communications example, however, this is not
the case, but the above theory can still be applied to each computed Hgq,r as
long as the relative bandwidth Ω of gq,r is much smaller than the carrier frequency
def
Ωc,q = Ωc +qbΩ. ¯ Normally
¯ we also want Ω to be larger than the maximum Doppler
¯ VP ¯
shift |νP ξ| = ¯ Vw ξ ¯ of (3.2), so as a rule of thumb |VP /Vw | should be very small,
so that a larger but still small Ω/ξ can be chosen. In this case, and provided
that the attenuation factor Aξ (P) in (3.3) is slowly varying with ξ, the following
modifications allow for a “wideband use” of the propositions in this section:
1. Equation (3.9c) splits into separate operators Hq for basis functions gq,r
with center frequency Ωc,q . For some compact sets Kq ⊆ (−Ωc,q , Ωc,q ), these
operators have integral representation
Z
Hq gq,r (·) = SHq (ν 0 , t)(Tt Mν 0 gq,r )(·) d(ν 0 , t),
Kq ×[A,∞)
µ 0 ¶ µ 0 ¶
0 def 1 ν −αξ0 t ν
SHq (ν , t) = rξ0 ,t e χK×[A,∞) ,t .
Ωc,q Ωc,q Ωc,q
2. A straightforward multivariate extension of (3.9c) gives an operator spread-
ing function SH for signals with frequency localization near some fixed center
frequency ξ 0 . Then the operator’s action on narrowband functions gq,r with
def
center frequency Ωc,q = Ωc + qbΩ is given by a spreading function
¯ ¯ µ ¶
def ¯¯ ξ 0 ¯¯ ξ0
SHq (ν, t) = ¯ SH ν, t (4.17)
Ω ¯ Ωc,q
c,q
29
3. With ξ 0 chosen so that all [0, Ωc,q ] ⊆ [0, ξ 0 ], we keep the notation (4.7),
now with Iωc ,ω × IC,L being the smallest interval containing supp SHq for all
q ∈ Q.
4. These changes do not affect Propositions 2–4, but if S\ H (·, t) has the support
given by (4.15b), then the dilation with a factor Ωc,q in (5) gives that S\
ξ0
Hq (·, t)
has its support in
def
IC0,q ,L0,q = I ξ0 C0 , ξ0 L0 .
Ωc,q Ωc,q
def 1
5. Hence we also know from (4.17) that if ω 0,q = L0,q
, then
¯ ¯ µ ¶
¯ ξ0 ¯ ξ0 def
¯
|ω 0,q | SHq (nω 0,q , t) = ¯ ¯
ω 0,q ¯ SH n ω 0,q , t = |ω 0 | SH (nω 0 , t)
Ω
c,q Ω c,q
so that (4.15c) takes the form

X
SHq (ν, t) = |ω 0 | SH (nω 0 , t)ei2πht−C0,q ,ν−nω0,q i sincL0,q (ν − nω 0,q ).
n∈N
Insertion of this in the proof of Proposition 5 then changes (4.16) into the
Ω\ ,Ω
following formula for computing all SHqc,q (·, t) from the same coefficients
Ω00 00
SΩ 00 (nω 0 , pT ):
c
Ω\ ,Ω
X 00
SHqc,q (·, t)(t0 ) = |ω 0 T 00 | χC0,q ,L0,q (t−t0 ) ei2πhΩc,q ,t−pT i sincΩ (t − pT 00 )×
p∈P
X
Ω00 00 ,nω
× S Ω00
c
(nω 0 , pT 00 )ei2πht−t0 −pT 0,q i
. (4.18)
n∈N
Remark. The described procedure can also be used for the analysis of channel
operators applied to signals not generated by a Gabor frame with narrowband
window function, as well as for the analysis of Hilbert–Schmidt operators satisfying
the properties in Figure 3 (b) in general. Then, the narrowband windows g and γ
are chosen to investigate properties of the operator. Diagonalization properties are
still useful, but now in the sense that they give a “simple” discretized descriptions
of the operator.
Remark. The derivation of equation (4.2) does not require (gq,r ) or (γq0 ,r0 ) to be
Gabor frames. Thus one possible future variation of the results of this paper is to
use, for example, compactly supported wavelet bases, such as B-spline, Daubechies
or Morlet wavelets, which also give synthesis and analysis of signals that can be
reconstructed from a finite number of sample values with bandlimiting conditions
replaced by projections on certain shift-invariant wavelet subspaces [EG05, EG04].
30
Gabor bases are a natural first choice for the OFDM applications that we have
described. Wavelet frame modifications of our algorithm might be more interesting
for a wideband communications scenario since the “frequency-dependent modula-
tion” in (3.2a) is actually a dilation that can only be reduced to a modulation in
the narrowband scenario described in Section 3.2.
5 The algorithm and its implementation

In Section 5.1 we summarize the results of the last section in an algorithm for
the coefficient operator matrix computation. Then we propose some further re-
finements for a fast implementation and compute its complexity in Section 5.2,
followed by suggestions for how to choose windows and parameters in Section 5.3.
5.1 The algorithm

The results of Section 4 lead to the following procedure for computing the co-
efficient operator matrix of a Hilbert–Schmidt operator satisfying the conditions
outlined in Figure 3(b).
Ω 00 00
1. Choose the spreading function coefficients SΩ 00 (nω 0 , pT ), the Gabor win-
c
dows g, γ and set all parameters to values typical for the application at hand.
We give suggestions for how to do this in Section 5.3 and 6.
2. For all q, q 0 ∈ Q and r, r 0 ∈ R, compute the matrix element hHgq,r , γq0 ,r0 i
in the following way:
(a) Compute the samples, index sets and supports of gq,r and γq0 ,r0 as
described in (4.6) and (4.7).
(b) Compute the matrix elements hHgq,r , γq0 ,r0 i by applying (4.8) to Hgq,r
and γq0 ,r0 . For this we need the samples (Hgq,r )(kTγ ), which we obtain
by setting f = gq,r in the finite sum formula (4.12), in which we get
the samples
³ ´
Ωc +QbΩ,Ω b
SH (·, kTγ − mTg ) (−mTg )
from the finite sum formula (4.16) or (4.18) for the more wideband
underwater communications example below.
See Appendix B for a detailed description of a Matlab implementation of this

procedure in the univariate case.
31
5.2 Refinements and complexity
For an efficient implementation, we also recommend the following two refinements
of the above algorithm:
1. In step 2 (b), the sum (B.9) can be obtained from a simple modulation of a
small number of sample values which should be computed in advance.
In fact, for the setup given in (4.7a) let Iq and Ibr be the smallest intervals
such that Hgq,r ⊆ Ir and Hg \ b
q,r ⊆ Iq for all q ∈ Q and r ∈ R. We shall
see below that there is only a small number of q 0 ∈ Q and r, r 0 ∈ R such
that the overlaps Ir ∩ supp γq0 ,r0 and supp Ib0 ∩ supp γd r 0 ,q 0 are nonempty.
With (B.9) computed for the above overlaps, we only need simple modulation
to also compute (B.9) for u = Hgq,r , v = γq0 ,r0 , q, q 0 ∈ Q and r, r 0 ∈ R in
Proposition 3.
2. Recall from (4.16) and Figure 5 that |N | · |P| is the number of Fourier
Ω00
coefficients needed to describe the smooth truncation SΩ00 of BH (ν, ·) to the
c
frequency interval IΩ00c ,Ω00 .
It is clear from Proposition 5 that |P| is proportional to the total bandwidth
IΩ0c ,Ω0 occupied by the Gabor basis (gq,r ) and thus also proportional to |Q|.
This dependence on |Q| comes from our choice to compute all Hgq,r from
Ω00
the same SΩ 00 . This simplified the presentation and is also memory-efficient,
c
since it minimizes the total bandwidth IΩ00c ,Ω00 \ IΩ0c ,Ω0 added for the smooth
truncation that is necessary to obtain a finite index set P finite. Thus it also
minimizes the total number of coefficients that are needed to describe the
channel behaviour in the entire frequency band IΩ0c ,Ω0 . The resulting algo-
rithm is also fast enough for the 2048 × 2048-matrix examples of Section 6.
For more computational efficiency when |Q| is large, however, it is favourable
Ω00
q
to do a separate smooth cut-off SΩ00c,q of BH (ν, ·) for every q to an inter-
Ω 00 00
val IΩ00c,q ,Ω00q ⊇ IΩc +qbΩ,Ω . This results in the coefficients SΩ00 (nω 0 , pT )
c
of (4.16) being replaced with a larger total number of coefficients, but with
index sets |Pq | proportional to |Ω|.
With these refinements, the complexity of the above procedure is that of one
nested sum over the index sets K, M, N and P for each q, q 0 ∈ Q and r, r 0 ∈ R.
Altogether, this requires O(|K| · |M| · |N | · |P| · (|Q| · |R|)2 ) arithmetic operations.
The following are typical index set sizes for the example applications of Section 6.
• |R| is the number of symbols for which the ISI shall be computed. For some
example applications with optimally well time-frequency localized Gaussian
window g, Figure 12 below shows that |R| of size 4–5 is enough to cover
a decay of the average ISI/ICI to 10−6 times the average of the diagonal
entries. For other windows g, we expect a need for larger |R|.
32
• |Q| is the number of carrier frequencies, which normally equals the number
of samples per received symbol. In radio communications |Q| is typically in
the range 128–1024. For the inherently more wideband underwater example
in Section 6, much smaller |Q| is possible. With the window and lattice
matching described in Section 5.3, |Q| = 47 in our underwater example
plots.
• |K| and |M| depend on the time-frequency localization of γ and g, respec-

tively. For the Gaussian windows used below, |K| · |M| = 28 · 19 = 532 for
the underwater channel and |K| · |M| = 19 · 19 = 361 for the other two.
• |N | · |P| is a constant that depends on |Ω| and is proportional to the area

hω, Li of the spreading function support. Below |N |·|P| equals 13·27 = 351,
13 · 58 = 754, and 59 · 237 = 13983 in the OFDM, satellite and underwater
example, respectively.
Hence, if g and γ have roughly M = |M| nonzero Nyquist frequency samples,

and if the received symbols have Q = |Q| nonzero Nyquist frequency samples
and if R = |R|, then our refined algorithm requires R2 · O(M 2 · Q2 ) arithmetic
operations. This can be compared to the R2 · O(Q5 ) operations of the more naive
and straightforward matrix computation approach discussed in section 4, which is
clearly slower when the number of carrier frequencies Q is larger than M .
Example 1. We show in Figure 7 that the Gaussian window approximation used

in Section 6 has M = 19 Nyquist frequency samples. We can compare this with
def def
the B-spline windows B0 = I0,1 and Bn = Bn−1 ∗ I0,1 , which also were proposed
in [MSG+ 05]. Since B cn (ξ) = sinc1 (ξ)n+1 , its ²-essential support is contained in
the interval I0,Ω for which 1/(πΩ/2)n+1 = ², that is, Ω = 2²−1/(n+1) /π. Since the
length of the support of Bn is n + 1, M is now the smallest integer larger than
2²−1/(n+1) (n + 1)/π + 1, which we plot for ² = 10−6 and some n in Figure 6. The
smallest M is 25, which we get for B12 , B13 and B14 .
3
10
M
2
10
20 40 60 80 100
n
Figure 6: Number of nonzero Nyquist frequency samples for some B-splines Bn .
33
5.3 Parameters and window functions
In the operator discretizations (4.1), the shape of the Gabor window g can be
optimized in different ways depending on the application. For wireless multicarrier
systems, the pulse shaping of g is an active research topic in its own, which we
do not pursue further here (see, for example, [MSG+ 05] for more details and
references). We have already motivated in Section 4 that both bases should be
Gabor bases with the same lattice parameters. For reasons described below, we
choose to set the lattice “time” parameter to aTg with a ∈ Z+ . Similarly, for
some b ∈ Rd+ we set the lattice frequency parameter to bΩ, so that the unit cell
“area” aTg bΩ = ab only depends on a and b. In the examples of Section 6, we
choose parameters so that ab > 1, which gives undersampling (see [Grö00]) and
a Gabor system that is a Riesz basis for its span. In general, for finite index sets
Q × R, we will consider Gabor Riesz bases
gq,r = TraTg MqbΩ g and γq,r = TraTg MqbΩ g, (q, r) ∈ Q × R. (5.1a)
We are primarily interested in windows with very good joint time-frequency
localization, which is of utmost importance for low ISI and ICI. Good frequency
localization also allows for high transmission power and, therefore, large signal-
to-noise ratio in the transmission band without exceeding power leakage bounds
for other frequency bands. Such leakage is strictly regulated for radio commu-
nications, but not (yet) for underwater sonar communications, where, however,
frequency bands already in use by, for example, the human ear or dolphins should
be respected.
In the following, we shall use standard Gaussian windows
s
π d −π2 h αν ,ν i
g(x) = e−hαx,xi with Fourier transform gb(ν) = e . (5.1b)
|α|
Gaussian windows have optimal time-frequency localization [Grö00, Section 2.2],

which results in very low ISI/ICI. This also has the advantage of making the effects
of truncation to ²-essential support clearly visible in the plots of Section 6.
The window γ is set to γ = Hg in the example applications of Section 6. Lower
ISI/ICI can be obtained with other choices described in [MSG+ 05].
Other application-dependent parameters we choose in the following way:
1. Index sets: Choose the Gabor basis index set R, the number of carrier
frequencies |Q| and the desired minimum size of the index sets N and P.
Larger |N | and |P| can be used for obtaining better control of the smoothness
Ω00
and decay of SΩ 00 .
c
2. Channel-dependent parameters: For the channel at hand, choose the total

allowed frequency band [Ω0 , Ω1 ], a “rectangle” Iωc ,ω × IC,L containing the
²-essential support of SH and the type of time-decay of SH . (See Section 6
for examples.)
34
3. Choose the lattice constants a and b, as well as the parameter α in (5.1b),
which uniquely determines the ²-essential support Ω = T1g of gb. For non-
Gaussian g, this step corresponds to the choice of a dilation parameter α in
def
an equation of the kind g(x) = g0 (αx).
We have done the following choice, inspired by a similar design criterion
in [KM98, Koz97, Mat00] For the smallest “rectangle” IAc ,A × IB c ,B that
contains the ²-essential support of the short-time Fourier transform of g
with window g (defined in (2.10)), choose a, b and α such that
ω B bΩ
= = (5.2)
L A aTg
and such that (gq,r ) is a Riesz basis for its span (obtained by setting ab > 1
in the results presented in Section 6). See Section B.13 for details.
4. Choose Q so that the supports of all gd

q,r are in the allowed frequency band
[Ω0 , Ω1 ]. Set Ωc and Ω according to (4.15a). Set Ω00c = Ω0c and choose Ω00
0 0
to be large enough to obtain the desired minimum size of P and also large
enough to for the desired smooth truncation (4.15d) and (4.15e) (Ω00 ≥ 1.5Ω0
1
in our implementation). Set Tγ = Ω+ω and T 00 = Ω100 .
5. Choose L0 and C0 so that (4.15b) holds and such that L0 is large enough
to obtain the desired minimum size N = ωω0 = ωL0 . Also set ω 0 = L10 .
6. Compute the Nyquist frequency samples of the restriction of g to its ²-

essential support (as illustrated in Figure 7 (a)).
Ω 00 00
7. Compute the samples SΩ 00 (nω 0 , pT ) that appear as coefficients in (4.16).
c
There are basically two ways to obtain these coefficients:
(a) Compute them from measurements or from a detailed model of a real

channel.
(b) Generate a random choice of the samples that satisfies the decay and
support properties described in Section 3 and Section 6. It seems
reasonable to assume that then there is always some constellation of
antennas and reflecting objects that correspond to this particular shape
of the spreading function.
For the results in Section 6 we have chosen approach (b).
6 Applications
We shall now describe parameter values and show sample plots for three channels
with different spreading function support areas.
35
6.1 Mobile phone communications
For a mobile phone moving at 100 km/h, equations (3.2) imply that the maximum
Doppler spread is VVwP ξ = 100/(299792.458 · 3600)ξ ≈ 10−7 ξ. For the GSM mobile
phone frequency bands ξ is ranging from 450.4–1990 MHz, so that frequency shifts
caused by the Doppler effect is in the interval [−184, 184] Hz. Typical time delays
are of the order 10−6 s, so the spreading function support area is of size 4 · 10−4 .
This support would exceed the critical value 1 (see page 11) if the highest channel
frequencies in use exceed 5 THz. These are infrared light frequencies, for which
the water in the Earth’s atmosphere absorbs too strongly for us to expect these
frequencies to be useful for wireless communications. Thus, with the terminology
of Section 2.4, mobile phone channels are inherently underspread.
With other parameters chosen as described in Section 5.3, we obtain a Gaussian
window g0 with Fourier transform ²-essential support Ω. Hence we can approx-
imate g0 by reconstructing it from the Nyquist sampling theorem with sample
interval Tg = 1/Ω and samples outside the ²-essential support of g0 discarded.
This construction gives a very small upper bound for both the supremum and the
L2 -norm of the resulting truncation error when only a small number of nonzero
samples are kept, as can be seen in the plots of g0 , g, gb0 and gb in Figure 7. In
practice, ² should be chosen small enough for the resulting errors to be dominated
by the overall noise level of the application at hand. For example, note in Fig-
ure 7 (a) how the truncation to a truly bandlimited g with a finite number of
0
10
−5
²g (0) 10
−5
10
²ĝ(0)
−10 ∧
10 |g(x)| |g (ξ)|
g(mT0) −10 g∧0 (ξ)
10
−5 0 5 −2 −Ω2 −1 0 1 Ω
2 2
t (s) x 10
−4 ξ (Hz) x 10
4
(a) (b)
Figure 7: We use the Gabor window/pulse shape g that we obtain from a Gaussian
2
g0 (x) = e−αx by assuming gb0 to have bandwidth Ω given by its ²-essential support
and truncating the reconstruction of g0 from its Nyquist frequency samples to
samples in the ²-essential support of g0 . We do this for ² = 10−6 . Plots (a) and
(b) show a comparison of the resulting g and gb with g0 and gb0 , respectively.
36
nonzero samples gives “Gaussian decay” down to amplitudes below ² and then
a slowly decaying tail with negligible amplitude and L2 -norm. In this and the
following examples, ² = 10−6 .
We show some example plots in Figure 8 for a system with 128 carrier fre-
quencies and 16 OFDM symbols. In (a) we show the ²-essential support of the
short-time Fourier transform Vgq,r gq,r (defined in (2.10)) for some neighbouring
basis functions gq,r , from which we se that we can expect nonzero ISI and ICI at
least for basis functions at distance
def
p
k(q, r) − (q 0 , r 0 )kl2 = (q − q 0 )2 + (r − r 0 )2 ≤ 4
in the Gabor lattice. In (b) we plot the “bandpass filtered spreading function”
9
x 10
1.87
Frequency (Hz)
1.87
1.87
1.87
1.87
−3 −2 −1 0 1 2 3
Time (s) −4
x 10
(a) (b)
−5
Average amplitude
10
−10
10
−15
10
0 5 10 15
||(r,q)−(r′,q′)||
l
2
(c) (d)
Figure 8: OFDM channel example. (a) ²-essential support of the short-time

Fourier transform Vgq,r gq,r for some neighbouring basis functions gq,r . (b) The
spreading function. (c) The 10-logarithm of |hHgq,r , γq0 ,r0 i| with index ordering
nq,r defined in (6.1). (d) Off-diagonal decay.
37
Ω00
SΩ 00 (ν, ·) computed from its samples by using (4.15c) and (4.15f) in a way similar
c
to the proof of Proposition 5 (see Appendix B.1 for details). In this plot we have
assumed an environment with a very large amount of small scatterers adding up
to a white Gaussian noise distribution of the coefficient values (see Section 6.4 for
examples with a “smother” spreading function with more correlated coefficients).
For plotting the coefficient operator matrix, we need to define a linear ordering
of the index sets
Q×R = {q0 , q0 + 1, q0 + 2, . . . , q0 + (|Q| − 1)}×{r0 , r0 + 1, r0 + 2, . . . , r0 + (|R| − 1)} .
We have chosen to group together indices belonging to the same OFDM symbol
by using the order
def
nq,r = (r − r0 ) · |Q| + q − q0 + 1, (6.1)
for which we have plotted the 10-logarithm of the matrix element amplitudes in
Figure 8 (c). With this ordering, the matrix is divided into 16 × 16 submatrices,
such that in each submatrix, r and r0 are fixed. The submatrices for which r = r0
show the ICI of symbol number r. Submatrices for which r 6= r0 show the ISI
between OFDM symbols number r and r0 .
Due to the matching of the Gabor lattice and the shape of the frequency lo-
calization of g to the shape of the spreading function support, which we explained
in Section 5.3, the size of |hHgq,r , γq0 ,r0 i| should mainly depend on the distance
k(q, r) − (q 0 , r 0 )kl2 between the time-frequency support centerpoints in the Gabor
lattice. For making this off-diagonal decay more visible, we have grouped together
matrix elements for which this distance is the same and plotted the average ampli-
tude in Figure 8 (d). Note the clearly visible effect that occurs at distances more
than four, which is caused by the truncations of windows and spreading functions
to their ²-essential support.
6.2 Satellite communications

The speed of a communications satellite in geostationary orbit is about 3 km per
second. Thus the maximum Doppler shift is VVwP ξ = (3/299792.458)ξ ≈ 10−5 ξ.
With typical transmission frequencies 1-30 GHz [MB02], we can expect Doppler
shifts up to some 104 Hz. Here, we will use an example from [MB02, p. 47] with
transmission frequency 6 GHz and Doppler shift 18 kHz. Again, we assume the
maximum time-spread to be some 10−6 s. Then the spreading function support
area is less than 0.036, so this is an underspread channel as well.
We show the same plots as for the OFDM example in Figure 9. Note that as a
result of the window and lattice constant matching, the plots (a) in Figure 8 an 9
are largely the same up to scaling.
38
6.3 Underwater sonar communications
For a vehicle travelling at 30 knot in sea water and using sonar communications, we
have VVwP ξ ≈ 30·0.51444
1531
ξ ≈ 10−2 ξ. We will use parameters typical for some medium
range systems described in [Sto99] with maximum time spread around 0.01 s and a
typical frequency band 20-35 kHz, so that the maximum Doppler shift is about 350
Hz. (More examples can be found, for example, in [LO97, Mid87, Sto96, ZK00,
ZT02].) These settings give spreading function support area 7 and an overspread
channel, which is typical for underwater sonar communications channels in general.
We show the same plots as for the previous two examples in Figure 10.
9
x 10
6.0005
6.0004
Frequency (Hz)
6.0003 10000
′′(ν,t)|
c
′′
6.0002
|SΩ
Ω
5000
6.0001
0
6 1
−6 0 1
x 10 0 4
−2 −1 0 1 2 −1 −1 x 10
Time (s) −5 t (s) ν (Hz)
x 10
(a) (b)
0
10
Average amplitude
−10
10
0 5 10 15
||(r,q)−(r′,q′)||
l
2
(c) (d)
Figure 9: Satellite channel example. (a) ²-essential support of the short-time

39
6.4 Further spreading function examples
In the above examples we used independent Gaussian distributions to obtain the
spreading function coefficients. For a “nicer” environment, one can expect more
correlation between the samples. We show examples of such spreading functions
in Figure 11. In Figure 12 we compare these two different spreading functions
with all other channel parameters being identical. The plots show that these two
different correlations of the spreading function samples do not affect the speed of
the off-diagonal decay significantly.
4
x 10
2.12
2.1
Frequency (Hz)
2.08
2.06
2.04
2.02
1.98
−0.01 −0.005 0 0.005 0.01
Time (s)
(a) (b)
−5
Average amplitude
10
−10
10
−15
10
0 5 10 15 20
||(r,q)−(r′,q′)||
l
2
(c) (d)
Figure 10: Underwater channel example. (a) ²-essential support of the short-time
40
6000
′′(ν,t)|
40
|SΩ′′(ν,t)|
4000
c
′′
|SΩ
Ω
c
′′
Ω
2000 20
0 0
1 1
−6 0 100 −6 0 1
x 10 0 x 10 0 4
−1 −100 −1 −1 x 10
t (s) ν (Hz) t (s) ν (Hz)
(a) (b)
(c)
Ω 00 00
Figure 11: Spreading functions with more correlated coefficients SΩ00 (nω 0 , pT ).
c
(a)–(c) shows the spreading functions for the OFDM, satellite and underwater
example, respectively. Figure 12 shows a comparison of the resulting off-diagonal
decays with those in figures 8–10.
41
7 Conclusions
Using a refinement of the standard multipath propagation model for the short
time behaviour of narrowband wireless channels, we have derived a spreading
function integral representation of such channels with a C (∞) spreading function
with subexponential decay.
This, together with a channel discretization using well time-frequency localized
Gabor bases, allowed us to derive formulas and an algorithms for the efficient
computation of certain matrix representations of communication channels. The
elements of this matrix describe the intersymbol and intercarrier interference for
the transmitted signal. We derived the algorithm, as well as some refinements of
it, under a minimum of assumptions or simplifications beyond the channel and
signal properties that are known from our channel and signal model.
Next, we discussed parameter choices in general and for three different chan-
nels. For these channels we used a Matlab implementation of our algorithm
to compute example plots showing the time-frequency localization of the Gabor
0
10
−5
10
Average amplitude
−10
10
OFDM, Gauss
−15 Satellite, Gauss
10 Underwater, Gauss
OFDM smooth,
Satellite smooth
Underwater smooth
0 1 2 3 4 5 6 7 8 9 10
||(r,q)−(r′,q′)||
l
2
Figure 12: Off-diagonal decays with spreading function coefficients normalized so

that the average of kHgq,r k2 equals one.
42
basis functions, the spreading functions, the coefficient operator matrix, and its
off-diagonal decay.
Our implementation is fast enough for at least 2048 × 2048 matrices and con-
siderably faster than a simpler and more straightforward approach to computing
the matrix elements.
Due to bandwidth and delay restrictions, multicarrier communications must
use bandlimited basis functions defined by a finite number of nonzero Nyquist
frequency samples. Our plots show clearly how such restrictions affect the off-
diagonal decay of the coefficient operator matrix. Thus the algorithm and software
can be useful for the numerical comparisons of the off-diagonal decay for different
pulse shapes and parameter settings.
Moreover, although we primarily consider communications applications in this
report, we derived our algorithm in a more general multivariate setting, as an
analysis tool for certain classes of Hilbert Schmidt operators with potential other
theoretical and practical applications.
Acknowledgements: The authors would like to thank Dr. Werner Kozek
at Siemens AG in Vienna, Austria, and Professor Harald Haas at Jacobs Uni-
versity Bremen, Germany, for insightful discussions about the physical properties
of wireless channels. We are also thankful to Professor Karlheinz Gröchenig at
the Numerical Harmonic Analysis Group in Vienna, Austria, for providing the
references on which Section 2.2 is based.
43
Appendices
A Fourier transform decay vs differentiability
We claimed at page 7 that
“. . . the function g(x) = (1 + cos(πx))4 χ[−1,1] (x) has Fourier transform

decay gb(ξ) = O((1 + |ξ|)−α ) for α = 9 but not for α > 9”.
This is not entirely obvious, but follows from Corollary 2 below.
Theorem 3 (Decay vs. Fourier transform smoothness).
(a) For compactly supported f ∈ C (n) (R), it holds that fb(ξ) = O(ξ −n ).
(b) If f is the inverse Fourier transform of some fb ∈ L1 (R) and ξ n fb(ξ) ∈ L1 (R),
then f is n times continuously differentiable and
fd
(n) (ξ) = (i2πξ)n fb(ξ).
Proof. (a) For compactly supported f ∈ C (n) (R) it follows that f (k) ∈ L1 (R)
for k = 0, 1, . . . , n, so that fd(n) (ξ) is continuous, fd(n) (ξ) → 0 as |ξ| → ∞ and
integration by parts gives that

Z ∞ Z ∞
d(n)
f (ξ) = (n)
f (x)e −i2πξx
dx = (i2πξ) n
f (x)e−i2πξx dx = (i2πξ)n fb(ξ).
−∞ −∞
Thus fb(ξ) = O(ξ −n ).

(b) is proved for n = 1 in [Kat04, Theorem 1.6, p. 136]. Hence the result follows
also for larger n, since if fb ∈ L1 (R) and ξ n fb(ξ) ∈ L1 (R), then ξ k fb(ξ) ∈ L1 (R) for
k = 1, 2, . . . , n − 1.
def
Corollary 2. The function g(x) = (1 + cos(πx))n χ[−1,1] (x) has Fourier transform
decay gb(ξ) = O((1 + |ξ|)−α ) for α = 2n + 1 but not for α > 2n + 1.
Proof. Let γ be the continuous function

def
γ(x) = (1 + cos(πx))χ[−1,1] (x).
Then
γ 0 (x) = − π sin(πx)χ[−1,1] (x) (continuous) and

γ 00 (x) = − π 2 cos(πx)χ[−1,1] (x) (discontinuity at ±1, γ 00 (±1) = π 2 ).
44
Set
( )
def
X
K = (k1 , k2 , . . . , kn ) ∈ N. : km = 2n .
m
Now if we do an iterative application of the product rule for differentiation, then

the first discontinuous term that appears is the term c(2,2,...,2) γ 00 (x)n in the 2nth
derivative
X
g (2n) (x) = ck γ (k1 ) (x)γ (k2 ) (x) · · · γ (kn ) (x)
k=(k1 ,k2 ,...,kn )∈K
X
=c(2,2,...,2) γ 00 (x)n + ck γ (k1 ) (cx)γ (k2 ) (x) · · · γ (kn ) (x)
k∈K\{(2,2,...,2)}
with
µ ¶µ ¶ µ ¶
k1
def k2 kn
ck = ··· .
2n 2n − k1 2n − k1 − k2 − · · · − kn−1
Hence g is 2n − 1 times continuously differentiable and 2n times differentiable.

Similarly to the proof of Theorem 3 (a), integration by parts gives that
Z Z 1
\ (2n+1) −i2πξx
£ (2n) ¤
−i2πξx 1−ε
g (2n+1) (ξ) = g (x)e dx = lim g (x)e x=−1+ε
+ g (2n) (x)e−i2πξx dx
R ε&0 −1
i2πξ Z 1
2n e − e−i2πξ
= − 2ic(2,2,...,2) π + i2πξ g (2n) (x)e−i2πξx dx
2i −1
Z 1
2n 2n+1
= − 2ic(2,2,...,2) π sin(2πξ) + (i2πξ) g(x)e−i2πξx dx
−1
2n 2n+1
= − 2ic(2,2,...,2) π sin(2πξ) + (i2πξ) gb(x).
From this and the fact that g, g (2n+1) ∈ L1 , it follows that gb, g\
(2n+1) ∈ C and
0
gb(ξ) =O(ξ −(2n+1) )
For the converse statement, note that we cannot have gb(ξ) = O(ξ −(2n+1+ε) ) for
any ε > 0, because then gb(ξ), ξ 2n gb(ξ) ∈ L1 and Theorem 3 would imply that g is
2n times continuously differentiable.
45
B Matlab implementation: overview and func-
tion reference
This appendix includes and explains the Matlab source code used for computing
the channel matrices and producing all the plots in this report.
Figure 13 shows the structure of our Matlab code for computing the channel
matrix. The main function ChannelSetupAndAnalysis (see page 50) begins with
setting a few basic parameters, such as the type of channel, matrix sizes and the
level of detail in some of the plots. Then the function SetParameters (page 87)
is called for computing the system dependent parameters and window functions
that were described in sections 5.3 and 6. Finally, the function CoeffOpMat is
called for computing and plotting the matrix elements. If the variable VERBOSE is
set to true, then ChannelSetupAndAnalysis will also call some other functions,
described in further subsections, that produce some of the other plots in this
report.
For computing the matrix elements, the function CoeffOpMat (see page 53)
first computes the Gabor basis elements γq0 ,r0 and gq,r by calling the functions
TFshiftAndResample (page 107), resample (page 85) and TFshift (page 105)
for the necessary time-frequency shifts . The application of H to each gq,r is
computed by the function H (page 76), which in turn calls the function FTBPFS
for computing S Ω\
H
c ,Ω
(·, t) (page 69). Finally the matrix element inner product
hHgq,r , γq0 ,r0 i is computed by calling the function l2Innerprod (page 78).
ChannelSetupAndAnalysis
SetParameters CoeffOpMatrix
TFshiftAndResample TFshift H l2InnerProd
resample FTBPFS
Figure 13: Program structure of the Matlab code for computing the channel
matrix.
46
00
Ω
B.1 BPFS: Compute SΩ00 (ν, t)
c
File path: Matlab/ChannelMatrix/BPFS.m

Ω00
Short description: Compute SΩ 00 (ν, t).
c
Ω00
Additional explanation: This function computes the samples of SΩ 00 (ν, t) that
c
we plotted, for example, in Figure 8 (b). This is done as follows:
Ω00
In the special case when SH = SΩ 00 , (4.15c) gives that
c
00
X 00
Ω Ω i2πht−C0 ,ν−nω 0 i
SΩ00 (ν, t) = |ω 0 |
c
SΩ00 (nω 0 , t)e
c
sincL0 (ν − nω 0 ), |N | < ∞.
n∈N
(B.1)
Moreover, by rewriting (4.15f) in the Nyquist sampling theorem form (4.5b) we
get
X 00
Ω00 00 Ω 00 i2πhΩ00 00
c ,t−pT i sinc 00 (t − pT 00 )
SΩ00 (nω 0 , t) = |T |
c
SΩ 00 (nω 0 , pT )e
c
Ω
p∈P
Insertion in (B.1) gives

00
X 00 00
Ω 00
SΩ00 (ν, t) = |ω 0 T |
c
ei2πhΩc ,t−pT i sincΩ00 (t − pT 00 )×
p∈P
X 00
Ω 00 i2πht−C0 ,ν−nω 0 i
× SΩ00 (nω 0 , pT )e
c
sincL0 (ν − nω 0 )
n∈N
The bandpass filtered spreading function samples are provided in the input
Ω 00 00
S np(n, p) = SΩ00 (N(n)ω 0 , P(p)T ).
c
(B.2)
Source code
function y=BPFS(N,P,S_np,Omega,omega_c,omega,Omega_cBis,OmegaBis,...
L_0,C_0,t_vec,nu_vec,VERBOSE)
% USAGE: y=BPFS(N,P,S_np,Omega,omega_c,omega,Omega_cBis,OmegaBis,...
% L_0,C_0,t_vec,nu_vec,VERBOSE)
%
% Computes spreading function at time(s) t_vec and fequency/ies xi.
%
% INPUT:
% N = Index set for n -> S_np(n,p)
% P = Index set for p -> S_np(n,p)
% S_np = Bandpass filtered spreading function samples
% Omega = Bandwidth of the analyzed Gabor window.
% omega = Bandwidth of the the spreading function eta.
47
% Omega_c = Center frequency of the analyzed Gabor window.
% omega_c = Center frequency of the spreading function eta.
% Omega_cBis = Center frequency of interval containing all basis
% function frequency supports (see documentation).
% OmegaBis = Length of interval containing all basis
% L_0 and C_0 = Length and centerpoint of the shortest interval
% containing the support of x -> h(x,t)
% t_vec = scalar/vector
% nu_vec = scalar/vector
%
% OUTPUT:
% y = matrix of size length(t_vec)xlength(nu_vec) with
% y(r,c) = spreading function value at time t_vec(r)
% and frequency nu=nu_vec(c).
omega_0 = 1/L_0;
T_0=1/Omega;
TBis=1./OmegaBis;
P=P(:).’; % Make row vector
N=N(:); % Make column vector
t_vec = t_vec(:).’; % Make row vector
nu_vec = nu_vec(:).’; % Make row vector
nu_vec_len = length(nu_vec);
t_vec_len = length(t_vec);
y=zeros( length(t_vec), length(nu_vec));
CompStartTime=toc;
r=0;
for t = t_vec
r=r+1;
c=0;
for nu = nu_vec
c=c+1;
if abs(nu-omega_c)<omega/2
for P_ind=1:length(P)
pTBis=P(P_ind)*TBis;
y(r,c)= y(r,c) ...
+ exp(i*2*pi*Omega_cBis*(t-pTBis))...
*SincOmega(OmegaBis,t-pTBis)...
*sum(S_np(:,P_ind).*exp(i*2*pi*(t-C_0).*(nu-N.*omega_0)) ...
.*SincOmega(L_0,nu-N.*omega_0));
end
48
end
end
if VERBOSE
disp([ int2str(r*nu_vec_len) ’ function values (’ ...
num2str(r/t_vec_len*100) ’ %) in ’ ...
num2str((toc-CompStartTime)/60) ’ minutes.’])
end
end
y=y.*omega_0.*TBis;
49
B.2 ChannelSetupAndAnalysis: The main program
File path: Matlab/ChannelMatrix/ChannelSetupAndAnalysis.m
Short description: The main program.
Additional explanation: See description on page 46.
Source code
function ChannelSetupAndAnalysis()
VERBOSE=false; % If true, then some additional output will be generated

close all hidden; % Close figures
SetPlot(1); % defined in the very end of this file
%% ================ CHOOSE CHANNEL AND LATTICE PARAMETERS ================
%% Choose channel type

%channel = ’OFDM’;
channel = ’Satellite’;
%channel = ’Underwater’;
%% Choose box-shaped or exponentially decaying envelope function for the

% spreading function time dependence
%SprdFctTime =’box’;
SprdFctTime =’exp’;
%% Choose how c_np depends on n

%SprdFctCoeffType =’Kronecker delta’; % Gives box-shaped spreading function
%SprdFctCoeffType =’pseudo random’;
%SprdFctCoeffType =’uniformly distributed random’;
SprdFctCoeffType =’Normally distributed random’;
%SprdFctCoeffType =’smooth’;
%% Pick size of Gabor lattice

%sys_size=’small’;
%sys_size=’medium’;
sys_size=’large’;
%% Choose minimum number of coefficients in the index set P

%P_min_size=0;%25;
P_min_size=25;
%% Choose fast or detailed plots when in VERBOSE mode:

DraftOrFinal=’draft’;
50
%DraftOrFinal=’medium’;
%DraftOrFinal=’final’;
%% Choose a desired TF-lattice unit cell area a*b

% (the real one will be equal to or smaller than the desired one,
% due to rounding off b to integer value.
D=1;
%% Choose amplitude threshold for "essential supports"

epsilon=1e-6;
%% =============== COMPUTE AND PLOT SYSTEM MATRIX ELEMENTS ===============
tic
SPstart=toc;
[Omega_c,Omega,omega_c,omega,Omega_cBis,OmegaBis,...
B_0,C_0,T_0,T,g,a,b,M,N,P,Q,R,S_np]=...
SetParameters(channel,SprdFctTime,SprdFctCoeffType,sys_size,...
P_min_size,’Gaussian’,D,epsilon,VERBOSE);
SPtime=toc-SPstart;
MatStart=toc;
K = CoeffOpMat(Omega_c,Omega,omega_c,omega,Omega_cBis,OmegaBis, ...
B_0,C_0,a,b,M,N,P,Q,R,S_np,g,epsilon,...
SprdFctCoeffType,channel,DraftOrFinal,VERBOSE);
disp(’Summary:’)
disp([’Channel: ’ channel]);
disp([’Omega_c = ’ num2str(Omega_c) ’ Hz.’])
disp([’Omega = ’ num2str(Omega) ’ Hz.’])
disp([’omega_c = ’ num2str(omega_c) ’ Hz.’])
disp([’omega = ’ num2str(omega) ’ Hz.’])
disp([’Omega_c" = ’ num2str(Omega_cBis) ’ Hz.’])
disp([’Omega" = ’ num2str(OmegaBis) ’ Hz.’])
disp([’B_0 = ’ num2str(B_0) ’ s.’])
disp([’C_0 = ’ num2str(C_0) ’ s.’])
disp([’a = ’ num2str(a)])
disp([’b = ’ num2str(b)])
disp([’a*b = ’ num2str(a*b)])
disp([’K = [ ’ int2str(K(1)) ’ : ’ int2str(K(end)) ’ ].’])
disp([’M = [ ’ int2str(M(1)) ’ : ’ int2str(M(end)) ’ ].’])
disp([’N = [ ’ int2str(N(1)) ’ : ’ int2str(N(end)) ’ ].’])
51
disp([’P = [ ’ int2str(P(1)) ’ : ’ int2str(P(end)) ’ ].’])
disp([’Q = [ ’ int2str(Q(1)) ’ : ’ int2str(Q(end)) ’ ].’])
disp([’R = [ ’ int2str(R(1)) ’ : ’ int2str(R(end)) ’ ].’])
disp([’epsilon = ’ num2str(epsilon)])
disp([’SprdFctCoeffType = ’ SprdFctCoeffType ])
disp(’’)
disp(’’)
disp([’Parameters set in ’ num2str(SPtime/60) ’ minutes.’])
disp([’Matrix computed in ’ num2str((toc-MatStart)/60) ’ minutes.’])
disp([’’])
%% ========================== INTERNAL FUNCTIONS ==========================
function [Mzp,Mindzp]=ZeropaddRows(M,Mind);
[r,c]=size(M);
if length(Mind)~=r
error(’Nonononononono!!!’)
end
Mzp=[zeros(r,2) M zeros(r,2)];
BorderWidth=0.2*(Mind(end)-Mind(1));
D=0.2;
Mindzp=[Mind(1)-BorderWidth 2*Mind(1)-Mind(2) Mind ...
2*Mind(end)-Mind(end-1) Mind(end)+BorderWidth];
%% ========================= That’s all, folks... =========================
52
B.3 CoeffOpMat: Compute the Coefficient Operator Matrix
File path: Matlab/ChannelMatrix/CoeffOpMat.m
Short description: Compute the Coefficient Operator Matrix.
Additional explanation: Computes the matrix elements Gq0 ,r0 ;q,r = hHgq,r , γq0 ,r0 i
of (4.2) for q, q 0 ∈ Q and r, r 0 ∈ R and with gq,r and γq0 ,r0 defined in (5.1a). Com-
puted here from the bandpass filtered spreading function samples provided in the
input S np defined in (B.2).
The computations are done using the algorithm that is described in Section 5.1
with the modifications of the remark on page 29 implemented, but with refinement
2 on page 32 not implemented.
Source code
function Hg_support ....
=CoeffOpMat(Omega_c,Omega,omega_c,omega,Omega_cBis, ...
OmegaBis,L_0,C_0,a,b,M,N,P,Q,R,S_np,g,epsilon,...
SprdFctCoeffType,channel,DraftOrFinal,VERBOSE)
% USAGE Hg_support=CoeffOpMatrix(Omega_c,Omega,omega_c,omega, ...
% Omega_cBis,OmegaBis, L_0,C_0,a,b,M,N,P,Q,R,...
% S_np,g,epsilon,SprdFctCoeffType,channel,...
% DraftOrFinal,VERBOSE)
% INPUTS
% Omega_c = Center frequency (Hz) of the synthesis Gabor
% window g
% Omega = Bandwidth (Hz) of synthesis Gabor window g
% omega_c = Center frequency (Hz) of spreading function eta
% omega = Bandwidth (Hz) of spreading function eta
% Omega_cBis = Centerpoint of interval containing all basis
% C_0 and L_0 = Centerpoint and length of the shortest interval
% containing the support of t_0 -> h(t_0,t)
% a and b : the Gabor system TF-lattice parameters are
% a/Omega and b*Omega.
% M = Vector containing all m for which g(mT_0) is
% nonzero.
% Q = We will consider modulations by Q(q)*b*Omega.
% R = We will consider shifts by R(r)*a/Omega.
% g = Gabor window.
% epsilon = Amplitude threshold used when restricting
% functions to their "essential support".
53
% SprFctCoeffSupport= Support of spreading function mapping
% t -> eta(t,nu)
% (before bandpass filtering with respect to t).
% SprdFctType = Type of random distribution for elements in S_np.
% Text string used for pathnames of saved plots.
% channel = Type of channel (text string for plot pathname
% generation).
% S_np = Spreading function coefficients
% (length(N)xlength(P)-matrix).
% DraftOrFinal = ’draft’, ’mediuum’ or ’final’. Decides how fast
% or detailed some plots shall be that are made
% if VERBOSE==true
% VERBOSE = true or false. Tells whether additional plots
% and text output shall be generated
% warning(’S_{np} Temporary change here!!!’)

% [Nlen,Plen]=size(S_np);
% NhalfLen=floor(Nlen/2);
% PhalfLen=floor(Plen/2);
% S_np=zeros(Nlen,Plen);
% %S_np(round(Nlen/2),round(Plen/2))=1;
% for n=1:Nlen
% for p=1:Plen
% S_np(n,p)=exp( log(epsilon)*( ((n-round(Nlen/2))/NhalfLen)^2 ...
% +((p-round(Plen/2))/PhalfLen)^2 ) );
% end
% end
% %S_np(round(Nlen/2),round(Plen/2))=1;
% Tbis=1./OmegaBis
% omega_0 = 1/L_0
% omega=omega
% T_0=1./Omega
% T =1./(Omega+omega)
% Nlen=length(N)
% Plen=length(P)
% aT0=a*T_0
% bOmega=b*Omega
%
% % End of warning =====================================================
omega_0 = 1/L_0;
tic
54
T_0=1./Omega;
T =1./(Omega+omega);
TBis=1./OmegaBis;
%% Choose support
% Support of g : in M*T_0
% Support of Hg: in M*T_0+P*TBis
suppHgstart = M( 1 )*T_0+P( 1 )*TBis;
suppHgend = M(end)*T_0+P(end)*TBis;
Kstart = floor(suppHgstart/T);
Kend = ceil( suppHgend /T);
K=[Kstart:Kend].’;
Klen=length(K); Mlen=length(M); Nlen=length(N);

Plen=length(P); Qlen=length(Q); Rlen=length(R);
% Choose Gabor windows

% ====================
% Contrary to PerMat Deluxe (or however we will name it),
% we do here keep track of sample values and sample points by
% storing them in two separate vectors. For example, if
% g_cont is defined on R and if g(m*T_0) is nonzero only for
% m= 3, 4 ,5, then
%
% M = [ 3 4 5 ] and
% g = [g_cont(3*T_0) g_cont(4*T_0) g_cont(5*T_0)]
% We choose some Gabor window g with support [-1/2,1/2]:

% g = ones(Mlen,1); % = Nonzero sample values of g.
StartComputeHgTime=toc;
WinDescr=’Gaussian’;
%gamma = ones(Klen,1); % = Nonzero sample values of gamma.
[Hg Hg_support] ... % Set gamma = Hg
= H(g,epsilon,L_0,C_0,K,M,N,P,S_np, ...
Omega_c,omega_c,Omega,omega,OmegaBis);
% Normalize S_nk so that Hg has L2-norm approximately = 1:

PrelL2Norm=sqrt(sum(abs(Hg(:)).^2)*T);
S_np=S_np./PrelL2Norm;
Hg=Hg/PrelL2Norm;
TypeOfDual = 1;
55
if TypeOfDual == 1
%Choose gamma = Hg ( = Cruuuuuuuude first choice.)
gamma = Hg;
gamma_support=Hg_support;
elseif TypeOfDual == 2
% Choose gamma to give "approximately dual" Gabor basis
[gamma,gamma_support] ...
= GaborDual(Hg,Hg_support, ...
Omega_c+omega_c,Omega+omega,a*T_0,b*Omega,Q,R,epsilon);
elseif TypeOfDual == 3
gamma = g;
gamma_support = K;
gamma=resample(Omega_c,Omega,M,g,K*T)
end
if VERBOSE
figure
OverSamplingFactor=10;
Hg_t=[Hg_support(1):1/OverSamplingFactor:Hg_support(end)].’.*T;
Hg_r=resample(Omega_c+omega_c,Omega+omega,Hg_support,Hg,Hg_t);
gamma_t ...
= [gamma_support(1):1/OverSamplingFactor:gamma_support(end)].’.*T;
gamma_r=resample(Omega_c+omega_c,Omega+omega,gamma_support,gamma, ...
gamma_t);
realArea=sum(abs(real(gamma_t)));
imagArea=sum(abs(imag(gamma_t)));
if (imagArea/realArea>epsilon)
realArea=sum(abs(real(gamma_t)))
imagArea=sum(abs(imag(gamma_t)))
figure
plot(gamma_t,real(gamma_r),’b-’,gamma_t,imag(gamma_r),’r:’);
legend(’real(\gamma(t))’,’imag(\gamma(t))’)
xlabel(’t’)
error(’Where did that imaginary part come from!? (Se Figure).’);
else
gamma =real(gamma );
gamma_t=real(gamma_t);
end
plot(gamma_t,gamma_r,’b-’,Hg_t,Hg_r,’r:’);
axis tight
legend(’\gamma(t)’,’Hg(t)’)
xlabel(’t (s)’)
PathExcludingExtension = ...
56
[’Figures\’ channel ’\GaborSystem\gammaAndHg’ ...
num2str((toc-StartComputeHgTime)/60) ’minutes’];
multiplot(PathExcludingExtension,’landscape’)
end
%% Print and plot some input parameters and functions
if VERBOSE
disp([’Index vector lengts: Klen = ’ int2str(Klen) ...
’,Mlen = ’ int2str(Mlen) ’,Nlen = ’ int2str(Nlen) ...
’,Plen = ’ int2str(Plen) ’,Qlen = ’ int2str(Qlen) ...
’,Rlen = ’ int2str(Rlen)’.’]);
disp([’Klen*Mlen*Nlen*Plen*Qlen*Rlen = ’ ...
int2str(Klen*Mlen*Nlen*Plen*Qlen*Rlen) ’.’])
disp([’Preparations finished in ’ num2str(toc) ’ seconds.’])
disp(’Choosing basis functions.’)
close all hidden % ...to avoid a sometimes occurring Matlab-bug

SetPlot(2) % Set large fonts and not-so-thin-lines
figure
PlotSpreadingFunction(L_0,C_0,N,P,S_np,SprdFctCoeffType,channel,...
omega_c,Omega,omega,Omega_cBis,...
OmegaBis,TBis,omega_0,DraftOrFinal,VERBOSE)
end
%% MAIN PROGRAM: Compute matrix elements ==========================

if VERBOSE
disp(’Starts computing <Hg_{q,r},gamma_{q’’,r’’}>...’)
end
%close all hidden % Save some memory
MatrixCompTimeStart=toc;
PlotTimeStart=toc;
GaborSystemSize=Qlen*Rlen;
% Operator matrix (G) and distances (D) for off-diagonal decay plot:
G=zeros(GaborSystemSize,GaborSystemSize);
D=zeros(GaborSystemSize,GaborSystemSize);
row=0;
ElementNr=1;
Omega_cqmax=Omega_c+Q(end)*b*Omega
for r=1:Rlen % raT_0
57
for q=1:Qlen %qbOmega
row=row+1;
col=0;
Omega_cq =Omega_c+Q(q)*b*Omega;
omega_q = Omega_cq/Omega_cqmax*omega;
C_0q=Omega_cqmax/Omega_cq*C_0;
L_0q=Omega_cqmax/Omega_cq*L_0;
[g_qr NewM g_NewOmega_c] ...

= TFshift(g ,M, Omega_c, Omega, R(r)*a, Q(q)*b);
Kshifted = [ floor((suppHgstart + R(r)*a*T_0)/T): ...

ceil( (suppHgend + R(r)*a*T_0)/T) ].’;
[Hg_qr Hg_qr_support] ...
= H(g_qr,epsilon,L_0q,C_0q,Kshifted,NewM,N,P,S_np, ...
g_NewOmega_c,omega_c,Omega,omega,OmegaBis);
if sum(abs(Hg_qr))==0
error(’xxx’)
end
for rprime=1:Rlen
for qprime=1:Qlen
ElementNr=ElementNr+1;
% Apply TF-shifts to gamma:
col=col+1;
[gamma_qr gamma_qr_support gamma_NewOmega_c] ...
= TFshiftAndResample(gamma,gamma_support, ...
Omega_c+omega_c,Omega+omega, ...
R(rprime)*a*T_0,Q(qprime)*b*Omega);
G(row,col)= l2InnerProd(Hg_qr, Hg_qr_support, gamma_qr,...
gamma_qr_support, ...
g_NewOmega_c+omega_c, ...
gamma_NewOmega_c,Omega+omega);
D(row,col)=sqrt( (Q(q)-Q(qprime))^2 + (R(r)-R(rprime))^2);
end
end
if r==Rlen && q==Qlen
figure(42)
subplot(1,4,1)
plot(Hg_qr_support*T,real(Hg_qr),’b’,gamma_qr_support*T,...
real(gamma_qr),’r:’);
xlabel(’t (s)’)
axis tight
legend(’real(Hg_{q,r})’,’real(\gamma_{q,r})’)
58
subplot(1,4,2)
plot(Hg_qr_support*T,imag(Hg_qr),’b’,gamma_qr_support*T,...
imag(gamma_qr),’r:’);
xlabel(’t (s)’)
axis tight
legend(’imag(Hg_{q,r})’,’imag(\gamma_{q,r})’)
subplot(1,4,3)
plot(Hg_qr_support*T,abs(Hg_qr),’b’,gamma_qr_support*T,...
abs(gamma_qr),’r:’);
xlabel(’t (s)’)
axis tight
legend(’abs(Hg_{q,r})’,’abs(\gamma_{q,r})’)
subplot(1,4,4)
plot(Hg_qr_support*T,angle(Hg_qr),’b’,gamma_qr_support*T,...
angle(gamma_qr),’r:’);
xlabel(’t (s)’)
axis tight
legend(’angle(Hg_{q,r})’,’angle(\gamma_{q,r})’)
grid on
axis tight
end
if VERBOSE
disp([Int2str(ElementNr) ’ matrix elements (=’ ...
num2str(ElementNr/(GaborSystemSize.^2)*100) ...
’% of the matrix) computed in ’ ...
num2str((toc-MatrixCompTimeStart)/60) ’ minutes.’]);
end
end
end
close all hidden
%% Plot the channel operator matrix =================================

SetPlot(2); % Set larger fontsize and wider lines in plots
figure
PcolorFull(abs(G));
figsize(’a4landscape’)
axis square
colorbar
[r c]=size(G);
xlabel(’\itn_{q{\prime},r{\prime}}’);
ylabel(’\itn_{q,r}’);
59
zlabel(’|G_{\itq,r;q{\prime},r{\prime}}|’)
[’Figures\’ channel ’\ChannelOpMatrix\matrix’ ...
num2str(toc/60) ’minutes’];
close all hidden
figure
PcolorFull(log10(abs(G)));
axis square
colorbar
[r c]=size(G);
xlabel(’\itn_{q{\prime},r{\prime}}’);
ylabel(’\itn_{q,r}’);
zlabel(’|G_{\itq,r;q{\prime},r{\prime}}|’)
[’Figures\’ channel ’\ChannelOpMatrix\log10-matrix’ ...
% %% Plot l^2-distances ================================================

% figure
% PcolorFull(D);
% figsize(’a4landscape’)
% axis square
% colorbar
% %surf(D);
% %[C,h]=contour(D);
% %set(h,’ShowText’,’on’)
% axis tight
% xlabel(’\itn_{q{\prime},r{\prime}}’);
% ylabel(’\itn_{q,r}’);
% zlabel(’||\it(r,q)-(r{\prime},q{\prime})||_{\itl_2}’)
% % title([’System matrix for ’ WinDescr ’ Gabor window and ’ ...
% % SprdFctCoeffType ’ spreading function coeffients (’ ...
% % num2str((toc-PlotTimeStart)/60) ’ minutes)’])
% PathExcludingExtension = ...
% [’Figures\’ channel ’\ChannelOpMatrix\L2-distances’ ...
% num2str(toc/60) ’minutes’];
% multiplot(PathExcludingExtension,’landscape’)
60
close all hidden
%% Prepare for off-diagonal decay plots =============================
G=G(:); % Matrix entries

D=D(:); % Distances
[D,ind]=sort(D); % Sort in order of increasing distance
G=G(ind); % Same sorting
NrOfBins = 100;
G_average=zeros(NrOfBins,1);
BinWalls=linspace(D(1),D(end),NrOfBins+1);
D_remaining=D;
G_remaining=abs(G);
for BinNr = 1:NrOfBins
inds=find(D_remaining<=BinWalls(BinNr+1));
if length(inds)>0;
ind=max(inds);
G_average(BinNr)=sum(abs(G_remaining(1:ind)))./ind;
if (BinNr <NrOfBins)
D_remaining=D_remaining(ind+1:end);
G_remaining=G_remaining(ind+1:end);
end % if
end % if
end % for
if (ind < length(D_remaining))
error(’This should not have happened...’)
end % if
ind=1;
tmp=0;
NrOfTerms=0;
for n=1:length(G)
if abs(D(n)-D(ind)) <= 10*eps
tmp=tmp+abs(G(n)); % TEST!!!
NrOfTerms=NrOfTerms+1;
else
61
G(ind)=tmp/NrOfTerms; % TEST!!!
ind=ind+1;
NrOfTerms=1;
D(ind)=D(n);
tmp=abs(G(n)); % TEST!!!
end
end
G(ind)=tmp/NrOfTerms; % TEST!!!
D(ind)=D(end);
D=D(1:ind);
G=G(1:ind);
[’Figures\’ channel ’\ChannelOpMatrix\Off-diagDecay’ ...
save([PathExcludingExtension ’.mat’], ’D’, ’G’)
%% Off-diagonal decay, binned linear plot... ========================

figure
x=(BinWalls(1:end-1)+BinWalls(2:end))/2;
bar(x,G_average)
xlabel(’||\it(r,q)-(r{\prime},q{\prime})||_{\itl_2}’)
ylabel(’Average amplitude’)
axis tight
[’Figures\’ channel ’\ChannelOpMatrix\Off-diagDecayBins’ ...
%% Off-diagonal decay, linear plot... ===============================

figure
plot(D,G,’r.:’)
axis tight
[’Figures\’ channel ’\ChannelOpMatrix\Off-diagDecayLin’ ...
if VERBOSE
disp([’...and those computations took ’ num2str(toc/60) ’ minutes.’])
disp(’Thank you for your patience. :-)’)
end
62
%% Off-diagonal decay, logarithmic plot... ==========================
figure
ind = find (abs(G)>0)
semilogy(D(ind),G(ind),’b.:’)
axis tight
[’Figures\’ channel ’\ChannelOpMatrix\Off-diagDecayLog’ ...
D=D(ind);
G=G(ind);
%% ======================= INTERNAL FUNCTIONS =======================
function y=PcolorFull(varargin)
% Use: y=pcolor2(C)
% y=pcolor2(X,Y,C)
%
% Description: Modified version of the Matlab command pcolor.
% The difference is that PcolorFull does not remove the last row and
% column of the matrix before plotting.
if (nargin==1)
C=varargin{1};
[m,n]=size(C);
y=pcolor([C ones(m,1).*C(1,1)
ones(1,n+1).*C(1,1)]);
elseif(nargin==3)
X=varargin{1};
Y=varargin{2};
C=varargin{3};
X=X(:); % Make column vector

Y=Y(:); % Make column vector
[m,n]=size(C);
X=[X ; X(end)+(X(end)-X(end-1))];
Y=[Y ; Y(end)+(Y(end)-Y(end-1))];
y=pcolor(X,Y,[C ones(m,1).*C(1,1)
63
ones(1,n+1).*C(1,1)]);
else
error(’WrongINcorrect number of inputs.’)
end;
shading flat
64
B.4 CutToEssSupp: Truncate a function to its essential sup-
port
File path: Matlab/ChannelMatrix/CutToEssSupp.m
Short description: Truncate a function to its essential support.
Additional explanation: Truncate a function to its essential support, as defined
on page 5.
Source code
function [gc_sample_ind,gc_samples] ...
= CutToEssSupp(epsilon,OversampFact,g_centerfreq, ...
g_bandwidth,g_sample_ind,g_samples)
% Reduces the support of a function g sampled at Nyquist frequency
% by applying an amplitude threshold epsilon on g sampled at
% OverSamplingFact times the Nyquist frequency
%
% USAGE: [gc_sample_ind,gc_samples] ...
% = CutToEssSupp(epsilon,OversampFact,g_centerfreq, ...
% g_bandwidth,g_sample_ind,g_samples)
% INPUT
% epsilon = Cutting threshold
% OversampFact = Oversampling factor
% g_centerfreq = Center frequency of g.
% g_bandwidth = Bandwidth of g.
% g_sample_ind = Integer sample indices.
% g_samples = Samples of g such that the mth sample
% is sampled at "time" g_sample_ind(m)*T
% with T = 1/g_bandwidth
%
% OUTPUT
% g = Number, vector or matrix of the same size as t, containing the
% computed new samples
T = 1/g_bandwidth;
t=[g_sample_ind(1)*T:T/OversampFact:g_sample_ind(end)*T].’;
g=resample(g_centerfreq,g_bandwidth,g_sample_ind,g_samples,t);
absg=abs(g);
ind = find( absg./max(absg(:))>=epsilon );
tstart=t(min(ind));
tend =t(max(ind));
ind=find( ((g_sample_ind*T)>=tstart) & ((g_sample_ind*T)<=tend));
65
gc_sample_ind=g_sample_ind(min(ind):max(ind));
gc_samples =g_samples(min(ind):max(ind));
if length(gc_samples) == length(g_samples)
warning([’Sorry,no cutting was possible. Smallest relative amplitude was ’...
num2str(min(absg(:)./max(absg(:)))./epsilon) ’\epsilon.’])
end
66
B.5 figsize: Set the orientation and size of the current
figure
File path: Matlab/ChannelMatrix/figsize.m
Short description: Set the orientation and size of the current figure.
Additional explanation:
Source code
function figsize(size)
% Sets the size of the current figure
%
% Use: figsize(size)
if strcmp(size,’a4’)
%orient portrait
set(gcf,’PaperType’,’a4’);
wh=get(gcf,’PaperSize’);
width=wh(1);
height=wh(2);
set(gcf,’PaperPosition’,[0.25,0.25,width,height]);
elseif strcmp(size,’a4landscape’)
%orient landscape
width=wh(2);
height=wh(1);
elseif strcmp(size,’a4HalfLandscape’)
width=wh(2);
height=wh(1)/2;
elseif strcmp(size,’subplot13’)
width=wh(1);
height=wh(2)/6;
width=wh(1);
67
height=wh(2)/4;
width=wh(1);
height=wh(2)/3;
width=wh(1);
height=wh(2)/2;
else
error([’Plotsize ’’’ size ’’’ not known.’]);
end
68
Ω\
c ,Ω
B.6 FTBPFS: Compute SH (·, t)(t0 )
File path: Matlab/ChannelMatrix/FTBPFS.m
Short description: Compute S Ω\
c ,Ω
H(·, t)(t ).
0
Ω\c ,Ω
Additional explanation: Computes SH (·, t)(t0 )
using the formula (4.16).
Computed here with coefficients given as the input S np defined in (B.2), page (B.2).
Source code
function y=FTBPFS(Omega_c,omega_c,Omega,omega,OmegaBis, ...
B_0,C_0,N,P,S_np,t,t_0)
% y=FTBPFS(Omega_c,omega_c,Omega,omega,OmegaBis, ...
% B_0,C_0,N,P,S_np,t,t_0)
%
% Computes Fourier transform of bandpass filtered spreading function
% at time t and fequency xi
%
% USAGE:
%
% INPUT:
% B_0 and C_0 = Length and centerpoint of the shortest interval
% t = real number
% t_0 = real number
T_0=1/Omega;
TBis=1./OmegaBis;
omega_0=1/B_0;
P=P(:).’; % Make row vector
N=N(:); % Make column vector
y=0;
tmint0 = t-t_0; % Compute only once to save time
69
if ( (C_0-B_0/2<=tmint0) && (tmint0<=C_0+B_0/2) )
for P_ind=1:length(P)
p=P(P_ind);
pTBis = p*TBis;
y = y +...
exp(i*2*pi*Omega_c*(t-pTBis))...
*SincOmega(Omega,t-pTBis)...
*sum(S_np(:,P_ind).*exp(i*2*pi*(tmint0-pTBis).*N.*omega_0));
end
y=y*omega_0*TBis;
else
warning(’Oh dearie me, this should really never happen...’)
end
70
B.7 GaborDual: Compute the dual of a Gabor Riesz basis
File path: Matlab/ChannelMatrix/GaborDual.m
Short description: Compute the dual of a Gabor Riesz basis.
Additional explanation: Not carefully tested!!! This function was at first
considered for use but then never needed in this report due to the good diagonaliza-
tion properties of Gaussian windows in use. It might, however, be modified, tested
and included in future version of the software for use with other windows g.
Consider a Gabor Riesz basis for some subspace of L2 (Rd ) with the elements
def
gq,r = Tra Mqb g, Q, R ⊆ Zd
Fix some (q0 , r0 ) ∈ Z2 . For γL2 (Rd ) it follows that γ is dual to

def P
-series expansion γ = q,r∈Zd cq,r gq,r . recall from page 9 that γ is the canon-
ical dual window if and only if the coefficient sequence (cq,r )q,r∈Zd is the unique
minimum l2 -norm solution of
* +
X
δ(q0 ,r0 ) (q 0 , r 0 ) = hγq0 ,r0 , gq0 ,r0 i = cq,r gq+q0 ,r+r0 , gq0 ,r0
q,r∈Zd
X
= cq,r hgq+q0 ,r+r0 , gq0 ,r0 i ,
q,r∈Zd
(q’,r’;q,r)=(Q(1),R(1);Q(1),R(1)) (q’,r’;q,r)=(Q(1),R(1);Q(end),R(end))
G(1,1) G(1,NrOfBasisFct)
G(NrOfBasisFct,1) G(1,NrOfBasisFct)
(q’,r’;q,r)=(Q(end),R(end);Q(1),R(1)) (q’,r’;q,r)=(Q(end),R(end);Q(end),R(end))
Figure 14: Indexing of G and the corresponding values of (q 0 , r 0 ; q, r).
71
or in biinfinite vector notation with indices in Zd × Zd ,
Gc = δ (0,0) , where (G)q0 ,r0 ;q,r = hgq,r , gq0 ,r0 i = hgq−q0 ,r−r0 , g0,0 i
and (
def 1 if (q 0 , r 0 ) = (0, 0),
(δ (0,0) )
(q 0 ,r 0 ) =
0 otherwise.
G is known as the Gram matrix. Figure (14) together with the code shows how G
is indexed.
Source code
function [g_dual_samples,g_dual_sample_ind] ...
= GaborDual(g_samples,g_sample_ind, ...
g_center_freq,g_bandwidth,a,b,Q,R,epsilon)
% Computes the canonical dual of a Gabor window g
% that is completely determined by a finite number of
% input Nyquist frequency samples.
% QUICK & DIRTY implementation for finite index sets.
% E.g. Fourier series approach should be much faster!
%
% USAGE: [g_dual_samples,g_dual_sample_ind]
% = GaborDual(g_samples,g_sample_ind, ...
% g_center_freq,g_bandwidth,a,b,Q,R)
%
% INPUTS
% with T = 1/g_bandwidth, where ...
% g_center_freq = Center frequency of g.
% a = TF-lattice time-constant.
% b = TF-lattice frequency-constant.
% Q = Integer frequency-shifts in TF-lattice.
% R = Integer time-shifts in TF-lattice.
% epsilon = Amplitude threshold for truncation to
% epsilon-essential support.
%
% OUTPUT
% g_dual_samples = Sample values of the dual window.
% g_dual_sample_ind = Sample indices of the dual window.
T=1/g_bandwidth;
Qlen=length(Q);
72
Rlen=length(R);
q0=Q(round(Qlen/2));
r0=R(round(Rlen/2));
NrOfBasisFcts=Qlen*Rlen;
G=zeros(NrOfBasisFcts,NrOfBasisFcts);
%% Compute row 1 of the Gram matrix =================================

qprime=Q(1); rprime=R(1);
row=1;
col=0;
[g_qprimerprime_samples,...
g_qprimerprime_sample_ind,...
g_qprimerprime_center_freq] ...
= TFshiftAndResample(g_samples,g_sample_ind, ...
g_center_freq,g_bandwidth,rprime*a,qprime*b);
for r=R(:).’;
for q=Q(:).’;
col=col+1;
[g_qr_samples,g_qr_sample_ind,g_qr_center_freq] ...
g_center_freq,g_bandwidth,(r+r0)*a,(q+q0)*b);
G(row,col)=l2InnerProd(g_qr_samples,g_qr_sample_ind, ...
g_qprimerprime_samples,...
g_qprimerprime_sample_ind,T);
end
end
%% Compute Gram matrix elements to the right of the main diagonal ====
for row = 2:NrOfBasisFcts-1
G(row,row+1:end)=G(1,2:end-(row-1));
end
%% Compute the last row of the Gram matrix ===========================

qprime=Q(end); rprime=R(end);
row=NrOfBasisFcts;
col=0;
[g_qprimerprime_samples,...
g_qprimerprime_sample_ind,...
g_qprimerprime_center_freq] ...
g_center_freq,g_bandwidth,rprime*a,qprime*b);
73
for r=R(:).’;
for q=Q(:).’;
col=col+1;
if ( (r==r0) && (q==q0))
zero_ind= col;
end
g_center_freq,g_bandwidth,(r+r0)*a,(q+q0)*b);
G(row,col)=l2InnerProd(g_qr_samples,g_qr_sample_ind, ...
g_qprimerprime_samples,...
g_qprimerprime_sample_ind,T);
end
end
%% Compute Gram matrix elements to the left of the main diagonal =====
for row = 2:NrOfBasisFcts-1
G(row,1:row-1)=G(end,end-(row-1):end-1);
end
%% Compute coefficients ==============================================
pinvG=pinv(G); % Invertiing the Gram matrix

clear G;
c=pinvG(:,zero_ind); % = pinvG*(Kronecker Delta)
clear pinvG;
%% Compute the dual window from the coefficients =====================

g_dual_ind_start = floor(R( 1 )/T)+g_sample_ind( 1 );
g_dual_ind_end = floor(R(end)/T)+g_sample_ind(end);
g_dual_sample_ind = [g_dual_ind_start:g_dual_ind_end].’;
g_dual_samples = zeros(size(g_dual_sample_ind));
c_ind=0;
for r=R(:).’;
for q=Q(:).’;
c_ind=c_ind+1;
74
g_center_freq,g_bandwidth,r*a,q*b);
indstart = find(g_dual_sample_ind == g_qr_sample_ind( 1 ));
indend = find(g_dual_sample_ind == g_qr_sample_ind(end));
g_dual_samples(indstart:indend) ...
= g_dual_samples(indstart:indend) + c(c_ind)*g_qr_samples;
end
end
ind=find(abs(g_dual_samples)>=epsilon);
g_dual_sample_ind = g_dual_sample_ind(ind(1):ind(end));
g_dual_samples = g_dual_samples(ind(1):ind(end));
75
B.8 H: Compute samples (Hg)(kTγ ) using (4.12)
File path: Matlab/ChannelMatrix/H.m
Short description: Compute samples (Hg)(kTγ ) using (4.12).
Additional explanation: Function for computing the samples
X ³ ´
Ωc ,Ω b
(Hg)(kTγ ) = |Tg | g(mTg ) SH (·, kTγ − mTg ) (−mTg ).
m∈M
Source code
function [Hg Hg_support]=H(g,epsilon,L_0,C_0,K,M,N,P,S_np, ...
Omega_c,omega_c,Omega,omega,OmegaBis)
% USAGE: [Hg Hg_support]=H(g,epsilon,L_0,C_0,K,M,N,P,S_np, ...
% Omega_c,omega_c,Omega,omega,OmegaBis)
% Function for computating samples of Hg
%
% INPUT:
% g = Input sample values of some function g_cont:
% g(m)=g_cont(M(m)*T_g)
% epsilon = Amplitude threshold used for restricting the output to
% its "essental support"
% C_0 & L_0 = Centerpoint ansd length of the shortest interval
% K = Index set for which Hg shall be computed.
% M = Index set for g.
% N = Index set for spreading function trig poly coefficients.
% P = Index set for nonzero sample velues of spreading function
% S_np = Spreading function coefficients
% (length(N)xlength(P)-matrix).
% OmegaBis = Length of interval containing all basis function
% frequency supports (see the software documentation).
%
% OUTPUT:
% Hg = Vector of the same dimension and length as K with
% Hg(k) = sample value of corresponding continuous
% function Hg_cont. Hg(k)=Hg_cont(K(k)*T_gamma).
% Hg_support = index such that Hg(j)=0 for j not in Hg_support.
76
T_g=1/Omega;
T_gamma=1/(Omega+omega);
Klen=length(K);
Mlen=length(M);
Hg=zeros(size(K));
K=K(:); % Make column vector
g=g(:); % Make column vector
for k= 1:Klen
for m= 1:Mlen
FBS=FTBPFS(Omega_c,omega_c,Omega,omega,OmegaBis, ...
L_0,C_0,N,P,S_np,K(k)*T_gamma-M(m)*T_g,-M(m)*T_g);
if FBS==0
ind = find( (K*T_gamma<C_0-L_0/2) | (K*T_gamma>C_0+L_0/2) )
warning([’Element(s) nr ’ num2str(ind) ’ of ’ ...
int2str(length(K)) ...
’in K outside expecteded range! (’ ...
’Perhaps L_0 is too small?)’]);
end
Hg(k)=Hg(k)+g(m)*FBS;
end
end
Hg=T_g*Hg;
% Cut g and Hg to the essential support of Hg :
OversamplFact=2;
[Hg_support,Hg] = ...
CutToEssSupp(epsilon,OversamplFact,Omega_c+omega_c,Omega+omega,K,Hg);
77
B.9 l2Innerprod: Computes the l2 inner product (4.8a)
File path: Matlab/ChannelMatrix/l2Innerprod.m
Short description: Computes the l2 inner product (4.8a).
Additional explanation: This function computes the inner products (4.8), that
is, the inner products
X
hu, viL2 (Rd ) = |T | u(kT )vbpf (kT ),
k∈Iu
with vbpf = v if Cu = Cv and

X 0
vbpf (kT ) = |T | v(k0 T )ei2πhCuv ,(k−k )T i sincBuv ((k − k0 )T ).
k0 ∈Iv
otherwise.
Source code
function uv = l2InnerProd(u,u_ind,v,v_ind,C_u,C_v,B)
% USAGE: uv = l2InnerProd(u,u_ind,v,v_ind,C_u,C_v,B)
%
% INPUT:
% u = Vector of sample values of some function u_cont, ...
% u_ind = ... such that u(k) = u_cont(u_ind*T)
% v = Vector of sample values of some function v_cont, ...
% v_ind = ... such that v(k) = v_cont(v_ind*T)
% C_u = Center frequency of u.
% C_v = Center frequency of v.
% B = Bandwidth of u = Bandwidth of v.
%
% OUTPUT
% y = L_2 inner product of u and v
%
% NOTE: u_ind and v_ind MUST be row- or column vectors of the type
% [a a+1 ... b-1 b]
% but they do not need to be of equal length or even overlap.
% (Yes, it’s a waste of space to save more than the starting
% point of each index sets, but this is just a quick & dirty
% test-implementation without intentions to do any clever
% programming. =) )
T=1/B;
u_ind=u_ind(:); % Make column vector
78
u=u(:); % Make column vector
v_bpf=zeros(size(u_ind));
v_ind_len=length(v_ind);
JointFreqSuppStart = max(C_u-B/2,C_v-B/2);
JointFreqSuppEnd = min(C_u+B/2,C_v+B/2);
B_uv=JointFreqSuppEnd-JointFreqSuppStart;
C_uv=(JointFreqSuppEnd+JointFreqSuppStart)/2;
if B_uv <= 0
uv=0; % No overlap of frequency support
else
for kprime_ind = 1:v_ind_len
kprime=v_ind(kprime_ind);
k_minus_kprime_T=(u_ind-kprime).*T;
v_bpf = v_bpf ...
+ v(kprime_ind)*exp(i*2*pi*C_uv.*k_minus_kprime_T) ...
.*SincOmega(B_uv,k_minus_kprime_T);
end
v_bpf=T*v_bpf;
uv=T*sum(u.*conj(v_bpf));
end
79
B.10 multiplot: Command for saving the current figure in
.eps and .emf format
File path: Matlab/ChannelMatrix/multiplot.m
Short description: Command for saving the current figure in .eps and .emf
format.
Source code
function multiplot(PathExcludingExtension,orientation)
% Saves current figure in sizes and (colour or black and white) file
% formats suitable for importing into LaTeX and Powerpoint documents.
if nargin==1 || strcmp(orientation,’portrait’)
orient portrait
figsize(’a4’)
elseif strcmp(orientation,’a4HalfLandscape’)
orient portrait
figsize(’a4HalfLandscape’)
else
orient portrait
end
%print(’-depsc’, [PathExcludingExtension ’.eps’])
print(’-depsc2’,[PathExcludingExtension ’2.eps’])
colormap(gray)
%print(’-deps’ ,[PathExcludingExtension ’BW.eps’])
print(’-deps2’,[PathExcludingExtension ’BW2.eps’])
%print(’-r300’,’-dmeta’,[PathExcludingExtension ’.emf’])
print(’-dmeta’,[PathExcludingExtension ’.emf’])
if nargin==1 || strcmp(orientation,’landscape’)
%orient landscape
orient portrait
elseif strcmp(orientation,’a4HalfLandscape’)
orient portrait
figsize(’a4HalfLandscape’)
else
orient portrait
80
figsize(’a4’)
end
%print(’-r300’,’-dmeta’,[PathExcludingExtension ’.emf’])
print(’-dmeta’,[PathExcludingExtension ’BW.emf’])
81
B.11 PlotSpreadingFunction: Plots bandpass filtered spread-
ing function
File path: Matlab/ChannelMatrix/PlotSpreadingFunction.m
Short description: Plots bandpass filtered spreading function.
Ω00
Additional explanation: Plots the bandpass filtered spreading function SΩ 00 ,
c
which is computed here by calling the function BPFS (see Appendix B.1 with input
coefficients S np defined in (B.2), page 47).
Source code
function PlotSpreadingFunction(...
L_0,C_0,N,P,S_np,SprdFctCoeffType,channel,omega_c,Omega, ...
omega,Omega_cBis,OmegaBis,TBis,omega_0,DraftOrFinal,VERBOSE)
% USAGE: PlotSpreadingFunction(...
% L_0,C_0,N,P,S_np,SprdFctCoeffType,channel,omega_c,Omega, ...
% omega,Omega_cBis,OmegaBis,TBis,omega_0,DraftOrFinal,VERBOSE)
%
% Plots the bandpass filtered spreading function.
%
% INPUTS:
% L_0 and C_0 = Length and centerpoint of the shortest interval
% SprdFctCoeffType = Type of random distribution for elements in S_np.
% Text string used for pathnames of saved plots.
% channel = Type of channel
% (text string for plot pathname generation).
% Omega_cBis = Center frequency of interval containing all basis
% TBis = 1/OmegaBis. (Hmmmm... why was this imported?)
% omega_0 = 1/L_0. (Hmmmm... why was this imported?)
% DraftOrFinal = true or false, with true giving more detailed
% plots.
% VERBOSE = true or false, with true giving more textual
% output and more plots.
PlotTimeStart=toc;
82
if strcmp(DraftOrFinal,’draft’)
elseif strcmp(DraftOrFinal,’medium’)
elseif strcmp(DraftOrFinal,’final’)
else
error(’Ho hum...’);
end;
%% Compute plotting intervals ---------------------------------------

% For both nu and t, we plot the bandpass filtered spreading function
% in an interval that is the smallest interval I containing all
% samples extended at each endpoint with an extension of length
% PlotIntervalExtension=0;*|I|:
SetPlot(2); % Set larger fontsize and wider lines in plots
PlotIntervalExtension=0;
t_start=P(1)*TBis;
t_end=P(end)*TBis;
t_ext=PlotIntervalExtension*(t_end-t_start)/2;
t_start=TBis*floor((t_start-t_ext)/TBis);
t_end=TBis*floor((t_end+t_ext)/TBis);
t=[t_start:TBis/OverSamplingFactor:t_end];
t=[t_start:TBis:t_end]; warning(’Temporary removal of time oversampling!’)
NrOfSamples_t=length(t);
nu_start=N(1)*omega_0;
nu_end=N(end)*omega_0;
nu_ext=PlotIntervalExtension*(nu_end-nu_start)/2;
nu_start=omega_0*floor((nu_start-nu_ext)/omega_0);
nu_end=omega_0*floor((nu_end+nu_ext)/omega_0);
nu=[nu_start:omega_0/OverSamplingFactor:nu_end];
NrOfSamples_nu=length(nu);
if VERBOSE
disp([’Computing and plotting spreading function on a ’ ...
int2str(NrOfSamples_t) ’x’ int2str(NrOfSamples_nu) ’ grid’ ...
’(oversampling factor ’ int2str(OverSamplingFactor) ’)...’])
end
83
S_H=BPFS(N,P,S_np,Omega,omega_c,omega,Omega_cBis,OmegaBis,...
L_0,C_0,t,nu,VERBOSE);
%surf(nu,t,abs(S_H));
mesh(nu,t,abs(S_H));
%set(gca,’YTick’,[t(1) P(1)*TBis P(end)*TBis t(end)]);
%set(gca,’YTickLabel’,{num2str(t(1)), num2str(P(1)*TBis), ...
% num2str(P(end)*TBis), num2str(t(end))});
xlabel(’\nu (Hz)’)
ylabel(’t (s)’)
zlabel(’|S_{\Omega_c^{\prime\prime}}^{\Omega^{\prime\prime}}(\nu,t)|’);
% title([’Spreading function |(S_h(\nu,t))| computed from ’ ...
% SprdFctCoeffType ’ coefficients.’]);
axis tight
[’Figures\’ channel ’\SpreadingFunction\’ ...
84
B.12 resample: Compute new samples from Nyquist fre-
quency samples
File path: Matlab/ChannelMatrix/resample.m
Short description: Compute new samples from Nyquist frequency samples.
Additional explanation: The resampling done below is a more or less direct
application of the Nyquist sampling theorem (4.5b), with the input parameters
appearing as follows:
X
g(t) = |T | g(mT )ei2πhg centerfreq,t−mT i sincg bandwidth (t − mT ).
m∈g sample ind
Source code
function g=resample(g_centerfreq,g_bandwidth,g_sample_ind,g_samples,t)
% USAGE: g=resample(g_centerfreq,g_bandwidth,g_sample_ind,g_samples,t)
%
% Computes new samples of a function g that is completely determined
% by a finite number of input Nyquist frequency samples.
%
% INPUT
% g_centerfreq = Center frequency of g.
% with T = 1/g_bandwidth
% t = Number, vector or matrix containing the point at which
% g shall be computed.
%
% OUTPUT
% g = Number, vector or matrix of the same size as t, containing the
% computed new samples
g = zeros(size(t));
T = 1/g_bandwidth;
NrOfInputSamples = length(g_samples);
for m_ind=1:NrOfInputSamples
m = g_sample_ind(m_ind);
t_min_mT = t-m*T;
g = g + g_samples(m_ind)*exp(i*2*pi*g_centerfreq.*t_min_mT) ...
.*SincOmega(g_bandwidth,t_min_mT);
85
% % Correction for real-valued signals:
% g = g + g_samples(m_ind)*cos(2*pi*g_centerfreq.*t_min_mT) ...
% .*SincOmega(g_bandwidth,t_min_mT);
end
g=g*T;
86
B.13 SetParameters: Set system parameters for different
applications
File path: Matlab/ChannelMatrix/SetParameters.m
Short description: Set system parameters for different applications.
Additional explanation: Here the system parameters are set with channel-
dependent parameters chosen as described in Section 6 and other parameters
chosen according to the following algorithm. The steps of this algorithm are
enumerated with the same numbers as the corresponding parts of the source code.
1. Choose the Gabor system index set R the size of the index sets Q, N and
the minimum desired size of the index set P.
2. Choose a spreading function according to the channel-dependent character-

istics described in Section 6 by setting corresponding values on the Doppler
spread parameters ω c and ω as well as the total bandwidth
¡ [Ω0 , Ω1 ] and¢
a vector t exp decay par containing the parameters A ² t0 t1 t2 t3
that describe the spreading function time decay envelope function a in (B.8)
below.
3. Compute the bandwidth Ω and sampling interval Tg of g according to (B.6)

below. This is the bandwidth that we get for a Gaussian window
s
π d −π2 h αν ,ν i
g(x) = e−hαx,xi with Fourier transform gb(ν) = e .
|α|
(B.3)
and with lattice and window parameters chosen to match the shape of the
spreading function support as proposed in [KM98, Koz97] (but here with
an “²-essential modification” (see (B.4) below) that seemed reasonable but
not has been theoretically justified). We do this as follows: Let IC,L × Iωc ,ω
be the smallest “rectangle” in Rd × Rd that contains the ²-essential support
of SH . Similarly, let IAc ,A × IB c ,B be the smallest rectangle containing the
²-essential support of the short-time Fourier transform Vg g of g with window
g (as defined in (2.10)), which for our choice of g is
Z
Vg g(t, ν) = e−hαx,xi e−hα(x−t),x−ti e−i2πhν,x−ti dx
Rd
Z
e−2hα(x− 2 ),x− 2 i− 2 hαt,ti e−i2πhν,x−ti dx
t t 1
=
Rd
Z
i2π hν, 2t i− 12 hαt,ti
=e e−h2αy,yi e−i2πhν,yi dy
R d
s s
1
d
π −π2 h 2α ,ν iν π d − 12 (hαt,ti+π2 h αν ,ν i)
|Vg g(t, ν)| =e− 2 hαt,ti e = e .
|2α| |2α|
87
Our “²-essential” version of the design criterion in [KM98, Koz97] is
ω B bΩ
= = . (B.4)
L A aTg
The boundary of the ²-essential support of Vg g is the set of all points (t, ν)
such that ¿ À
2
ν
hαt, ti + π , ν = −2 ln(²).
α
1
The smallest rectangular box 2
[−(A, B), (A, B)] that contains all such
points (t, ν) has sides
− 21
p α2 p
1
A=α −2 ln(²) and B = −2 ln(²)

π
Hence from (B.4) we get
ω
α=π (B.5)
L
Similarly, if we define I0,Ω to be the smallest rectangle containing the ²-
®
essential support of gb, then by (B.3) αν , ν = − ln(²)
π2
, so that
2α 2 p
1
1
Ω= − ln(²) and Tg = . (B.6)
π Ω
4. Compute the lattice constants a and b using equations (B.7) below. Recall
from (5.1a) that we use Gabor frames gq,r = TraTg MqbΩ g with lattice con-
stants bΩ and aTg . We have b ∈ Rd+ but choose to require that a ∈ Zd+ ,
because then any element in the frame is easily computed from only a few
sample values g(mTg ) as described in Theorem 2, page 22.
Hence, given an input target sample density D = ab = (aTg )(bΩ), insertion
in (B.4) gives that ωL
bΩ
= aT g
= a2DT 2 . For communications applications we
g
want to have D ∈ R+ \ [0, 1] (for undersampling and unique Gabor basis
series expansion coefficients, see, e.g., [Grö00]). Thus we choose
Ã ! 12  Ã ! 12 
DL L
a= 2 rounded off to the nearest element in Zd+ \ 0, .
Tg ω Tg2 ω
(B.7a)
After rounding off a we recompute D and get
ω D
D = a2 Tg2 and b= . (B.7b)
L a
88
5. Check that Q not is too large for the Gabor system to satisfy the bandwidth
requirements [Ω0 , Ω1 ]. Set Ωc = 0 (which is required for real-valued win-
dows) and choose R so that it satisfies the bandwidth requirements [Ω0 , Ω1 ].
Compute Ω0c and Ω0 using (4.15a). Set T 0 = Ω10 and Ω00c = Ω0c . Choose Ω00
to be large enough to obtain a desired minimum size (P min size) of P and
1
such that Ω00 ≥ 1.5Ω0 . Set Tγ = Ω+ω and T 00 = Ω100 .
6. Compute the samples of the restriction of g to its ²-essential support (as

plotted in Figure 7 (a)). We get this support with the same arguments as
in step 3 by noting that g(x) = ² if hαx, xi = − ln(²).
7. Choose L0 and C0 according to (4.15b) and large enough for the number
|N | of samples nω 0 to be at least N desired length.
Ω 00 00
8. Compute the samples SΩ 00 (nω 0 , pT ) that appear as coefficients in (4.16).
c
There are basically two ways to obtain these coefficients:
(a) Compute them from measurements on a real channel.

(b) Generate a random choice of the samples that satisfy the decay and
support properties describes in sections 3 and 6, based on the underly-
ing assumption that then, there is always some constellation of anten-
nas and reflecting objects that correspond to this particular spreading
function.
For the results in Section 5 we have chosen approach 8b. The exponential
decay is computed as follows, where the input parameters are chosen in
step 2 above.
 ¡ ¢
ktk2 −t0

 ² · A² if t0 < ktk2 ≤ t1 ,


A if t1 < ktk2 ≤ t2 ,
a(t) = ¡ ² ¢ktk2 −t2 (B.8)

 A · if t < ktk ≤ t ,

 A 2 2 3
0 otherwise.
The subexponential decay is computed in the same way, since the differ-
ence between the exponential and subexponential decay is entirely in the
negligible (and neglected) tail.
Source code
function [Omega_c,Omega,omega_c,omega,Omega_cBis,OmegaBis,...
B_0,C_0,T_g,T_gamma,g,a,b,M,N,P,Q,R,S_np]=...
SetParameters(channel,SprdFctTime,SprdFctType,sys_size,...
P_min_size,window,D,epsilon,VERBOSE);
% USAGE:
% [Omega_c,Omega,omega_c,omega,Omega_cBis,OmegaBis,...
89
% B_0,C_0,T_g,T_gamma,g,a,b,M,N,P,Q,R,S_np] = ...
% SetParameters(channel,SprdFctTime,SprdFctType,sys_size,...
% P_min_size,window,D,epsilon,VERBOSE);
%
% The following is a very brief description of the input and output
% parameters. For a more comprehensive description, see the software
% documentation.
%
% INPUT:
% channel = ’OFDM’, ’satellite’ or ’underwater’
% SprdFctTime = ’box’ or ’exp’. Tells wether output parameter
% p_envelope shall be constant or exponential decaying
% SprdFctType = ’Kronecker delta’, ’pseudo random’, ...
% ’uniformly distributed random’ or ...
% ’Normally distributed random’.
% sys_size = ’small’, medium or ’large’. Tells whether to do
% quick computations for a small system, half-slow for
% a larger system orreally big ones.
% P_min_size = Minimum number of elements desired for the index set P.
% window = ’Gaussian’ (or all you’ll get back is an error message =) )
% D = TF-lattice unit cell area.
% D<1: oversampling
% D=1: critical sampling
% D>1: undersampling
% epsilon = amplitude threshold for epsilon-essential supports.
% VERBOSE = true or false. Tells whether additional plots
% and text output shall be generated.
%
% OUTPUT:
% Omega_c = Center frequency (Hz) of synthesis Gabor window g.
% Omega = Bandwidth (Hz) of synthesis Gabor window g.
% omega_c = Center frequency (Hz) of spreading function eta.
% omega = Bandwidth (Hz) of spreading function eta.
% Chosen so that the length of P is equal to or
% larger than the input parameter P_min_size.
% T_g = 1/Omega. (See the software documentation.)
% T_gamma = 1/(Omega+omega). (See the software documentation.)
% g = Gabor window.
% a and b : the Gabor system TF-lattice parameters are
% a*T_g and b*Omega.
% Q = We will consider modulations by Q(q)*b*Omega.
90
% R = We will consider shifts by R(r)*a/Omega.
% M = Vector containing all m for which g(mT_g)
% is nonzero.
% N = Index set for the function n -> S_np(n,p).
% S_np = Bandpass filtered spreading function samples.
% P = Index set for the function p -> S_np(n,p).
%% 1: Choose Gabor System size and desired length of index set N =====
% For a transmitter Gabor wondow (OFDM pulse shape) g with bandwidth
% Omega, we will consider modulations by Q(j)*b*Omega and .
% and shifts by R(l)*a*(sampling frequency of g or gamma).
%NlenDesired=0;%12;
NlenDesired=12;
% if strcmp(sys_size,’small’)
% Q = [-2:2].’;
% R = [-2:2].’;
% elseif strcmp(sys_size,’medium’)
% Q = [-8:7].’;
% R = [-8:7].’;
% elseif strcmp(sys_size,’large’)
% Q = [-16:15].’;
% R = [-16:15].’;
% end
% Set nr of frequecies/carriers (Qlen) and the
if strcmp(sys_size,’small’)
Qlen = 16; % 0.19 minutes for OFDM
R = [0:5].’;
if strcmp(channel,’Underwater’)
Qlen=47;
R = [0:2].’;
end
elseif strcmp(sys_size,’medium’)
Qlen = 128; %
Qlen=66;
%R = [0:16].’;
end
R = [0:5].’;
elseif strcmp(sys_size,’large’)
%Qlen=47;
%R = [0:20].’;
91
Qlen=47;
R = [0:2].’;
else
%Qlen = 128;
%R = [0:15].’;
Qlen = 512;
R = [0:2].’;
end
end
% if ( (sum( Q(:) == 0 )~=1) || (sum( R(:) == 0 )~=1))
% warning(’Both Q and R must contain zero for GaborDual to work.’)
% end
%% 2: Set spreading function parameters depending on channel type ====

if strcmp(channel(1:4),’OFDM’)
% OFDM channel with exponential time decay
omega_c = 0; % Hz. Center frequency of doppler spread interval.
omega = 184*2; % Hz. Length of doppler spread interval.
% This is, roughly, the

% http://www.gsmworld.com/technology/spectrum/frequencies.shtml :
Omega_0=1870e6;
Omega_1=1990e6;
if strcmp(SprdFctTime,’exp’)
t_exp_decay_par=[1 epsilon 1e-6*[0.8 1.3 2.5 3.5]];
%t_exp_decay_par=[1 epsilon 1e-6*[1.5 1.5 1.7 1.7]]; % TEMP!
elseif strcmp(SprdFctTime,’box’)
else
error([’’’’ SprdFctTime ...
’’’ is not a valid value of input parameter SprdFctTime.’])
end
elseif strcmp(channel,’Satellite’)
% Example from p.47 in satellite communications book:
omega = 18000*2; % Hz. Length of doppler spread interval.
Omega_0=6e9;
Omega_1=30e9;
t_exp_decay_par=[1 epsilon 1e-6*[0.25 1 2.5 3.5]];
else
92
end
elseif strcmp(channel,’Underwater’)
% Example from p.47 in satellite communications book:
%
% omega = 100; % Hz. Length of doppler spread interval.
% Omega_0=8000; % Hz
% Omega_1=11000; % Hz
Omega_0=20000; % Hz
Omega_1=35000; % Hz
omega = 30*0.51444/1531*Omega_1*2; % Hz.
% Doppler spread interval length.
%t_exp_decay_par=[1 epsilon [0 0.4 0.6 1]];
t_exp_decay_par=[1 epsilon [0 0.3 0.8 1]*0.01];warning(’Test...’);
t_exp_decay_par=[1 epsilon 1e-6*[0 0 1 1]];
else
end
warning(’underwater channel parameters not carefully chosen yet’)
else
error([’’’’ channel
’’’ is not a valid value of input parameter ’’channel’’.’]);
end
% Check later that the code is independent of time-shifts of
% the support (i.e. add +/- the below shift to Q’ instead):
t_exp_decay_par(3:6)=t_exp_decay_par(3:6) ...
-(t_exp_decay_par(3)+t_exp_decay_par(6))/2;
%% 3: Compute bandwidth and sampling density for g ===================
tau=t_exp_decay_par(6)-t_exp_decay_par(3);
alpha_par=pi*omega/tau % =alpha in the documentation,
% but that name is occupied
Omega = 2*sqrt(-alpha_par*log(epsilon))/pi
T_g=1/Omega;
%% 4: Compute lattice constants a and b ===============================

% such that a is a positive integer
% We choose to round off towards zero whenever possible...
93
a=round(sqrt((D*tau)/(T_g^2*omega)));
if a < sqrt(( tau)/(T_g^2*omega)) % would give oversampling
a=a+1;
end
if a==0;
%... and set a=1 otherwise:
a=1;
end
D=a^2*T_g^2*omega/tau;
b=D/a;
%% 5: Compute Q, sampling period times and choose OmegaBis ===========

% so that Tbis is small enough to fit at least P_min_size
% sample values with sampling period Tbis in an interval
% with length equal to the spreading function time support
% % Let’s choose the nr of samples such that ...

% RelAmplitude=0.4;
% Qlen=floor(2*sqrt(-log(RelAmplitude)/alpha_par)/T_g+1)
Omega_c = 0;
Q1 = ceil( (Omega_0+Omega/2)/(b*Omega) );
Q = [Q1:Q1+(Qlen-1)].’;
if ( Omega_1 < Omega_c+Q(end)*b*Omega -Omega/2 )
error([’Your Gabor basis exceeds the allowed frequency band. ’ ...
’Reducing |Q| from ’ int2str(length(Q)) ’ to ’ ...
int2str(floor((Omega_1-Omega_0)/(b*Omega))) ...
’ would solve the problem.’])
end
OmegaPrime =(length(Q)-1)*b*Omega+Omega+omega+omega_c;
Omega_cPrime=Omega_c+(Q(1)+Q(end))/2*b*Omega;
Tprime=1/OmegaPrime;
Omega_cBis=Omega_cPrime;
spr_fct_time_support_length=t_exp_decay_par(6)-t_exp_decay_par(3);
OversamplingFactor=1.5;
OmegaBis=max(OmegaPrime*OversamplingFactor,...
ceil(P_min_size/spr_fct_time_support_length));
T_gamma=1/(Omega+omega);
Tbis=1/OmegaBis;
%% 6: Compute g and plot Gabor frame STFT supports ===================
94
if strcmp(window,’Gaussian’)
M_max=ceil(sqrt(-log(epsilon)/alpha_par)/T_g);
M=[-M_max:M_max].’; % Vector containing all m for which
% |g(mT_g)|>=epsilon
g=exp(-alpha_par.*((M.*T_g).^2));
g=g./sqrt(sum(abs(g(:)).^2)*T_g); % Normalize to L^2-norm=1
if VERBOSE
plot_Gauss_frame(alpha_par,epsilon,Q,R,a,b,Omega_c,Omega, ...
g,M,channel);
end
else
error(’Only Gaussian window implemented, so far’)
end
%% 7: Choose omega_0 and B_0 =========================================

% Choose the support of t_0 -> h(t_0,t) so that it includes all
% points in time for which we possibly could be interested in
% computing H(g_{q,r})
t_start = floor(((M( 1 )+a*R( 1 ))*T_g+t_exp_decay_par(3)) ...
/T_gamma-1)*T_gamma;
t_end = ceil( ((M(end)+a*R(end))*T_g+t_exp_decay_par(6)) ...
/T_gamma+1)*T_gamma;
C_0 = (t_start + t_end )/2;

B_0 = t_end - t_start;
Nlen = ceil(max(NlenDesired,omega*B_0)); % since B_0=1/omega_0
B_0=Nlen/omega;
omega_0=1/B_0;
%% 8: Spreading function coefficients ================================

% Compute envelope:
[p_envelope P] = ...
spr_fct_time_decay(t_exp_decay_par,Tbis,VERBOSE,epsilon,...
SprdFctType,channel);
[n_envelope N] = ...
spr_fct_nu_decay(omega_c,omega,omega_0,VERBOSE,epsilon, ...
SprdFctType,channel);
S_np = zeros(length(N),length(P));
for p=1:length(P)
S_np(:,p)=p_envelope(p).*n_envelope;
95
end
% Normalize to L^2-norm 1
S_np=S_np/sqrt(sum(abs(S_np(:)).^2).*(omega_0*Tbis));
%surf(S_np)
if strcmp(SprdFctType,’pseudo random’);
% Used for getting "seemingly random" coefficients
% but in a way that can be exactly reproduced in any other
% (for example C++) implementation.
for P=1:length(P)
t=p*TBis;
S_np(:,p)=S_np(:,p).*cos(1000*10^6*N*t)+i*sin(12001*10^6*N*t);
end
elseif strcmp(SprdFctType,’smooth’);
% Coefficients uniformly random distributedin [-1-,1]x[-i,i]:
Plen=length(P)
for n=1:Nlen
for p=1:Plen
S_np(n,p)=2*S_np(n,p).*(5+exp(i*2*pi*(n/Nlen+p/Plen)));
end
end
elseif strcmp(SprdFctType,’uniformly distributed random’);
% Coefficients uniformly random distributedin [-1-,1]x[-i,i]:
S_np=2*S_np.*((rand(size(S_np))-0.5)+i*(rand(size(S_np))-0.5));
elseif strcmp(SprdFctType,’Normally distributed random’);
% Normally distributed random coefficients
S_np=S_np.*randn(size(S_np));
else
error(’Input parameter SprdFctType has an incorrect value’)
end
%% ====================================================================
%% ======================== INTERNAL FUNCTIONS ========================
%% ====================================================================
%% Envelope function for Spr.Fct.exponential decay in t=p*Tbis =======

function [p_envelope P] = ...
spr_fct_time_decay(t_exp_decay_par,Tbis,VERBOSE,epsilon, ...
SprdFctType,channel)
% Function for computing the envelope function governing the time-decay
% of the spreading function trigonometric polynomial coefficients:
A =t_exp_decay_par(1);
96
epsilon=t_exp_decay_par(2);
t0 =t_exp_decay_par(3);
Pstart=floor(t0/Tbis);
Pend =ceil( t3/Tbis);
P=[Pstart:Pend].’;
p_envelope=zeros(size(P));
for p_ind=1:length(P)
t=P(p_ind)*Tbis;
if t <= t0
p_envelope(p_ind)=0;
elseif t <= t1
p_envelope(p_ind)=epsilon*(A/epsilon)^((t-t0)/(t1-t0));
elseif t <= t2
p_envelope(p_ind)=A;
elseif t <= t3
p_envelope(p_ind)=A*(epsilon/A)^((t-t2)/(t3-t2));
else
p_envelope(p_ind)=0;
end
end
% We do not want Hf or system matrix entries to be overly small
% (in worst case, close to eps), so we normalize to
% coefficients of "L^1(R)-norm 1:
p_envelope = p_envelope./(Tbis*sum(abs(p_envelope)));
if VERBOSE
figure
stem(P*Tbis,p_envelope);
axis tight
xlabel(’p{\cdot}T’’’’ (s)’);
ylabel(’a(p{\cdot}T’’’’)’);
title([’Spreading function {\itp} with amplitude threshold ’ ...
’{\in}=’ num2str(epsilon) ’.’])
PathExcludingExtension = [’Figures\’ channel ’\CoeffEnvelopes\p’];
end
%% Envelope function for Spr.Fct. exponential decay in nu=n*omega_0 ==

function [n_envelope N] = ...
97
spr_fct_nu_decay(omega_c,omega,omega_0,VERBOSE,epsilon, ...
SprdFctType,channel)
% Here we are cheating slightly by making exponential rather
% than subexponential decay, but the difference between the two is
% all in the negligible (and here neglected) tail
Nstart=floor((omega_c-omega/2)/omega_0);
Nend =ceil( (omega_c+omega/2)/omega_0);
N=[Nstart:Nend].’;
n_envelope=zeros(size(N));
DecayingPartRelLen=0.7;
A =1;
nu0 =omega_c-omega/2;
nu1 =omega_c-omega/2*(1-DecayingPartRelLen);
nu2 =omega_c+omega/2*(1-DecayingPartRelLen);
nu3 =omega_c+omega/2;
for n_ind=1:length(N)
nu=N(n_ind)*omega_0;
if nu <= nu0
n_envelope(n_ind)=0;
elseif nu <= nu1
n_envelope(n_ind)=epsilon*(A/epsilon)^((nu-nu0)/(nu1-nu0));
elseif nu <= nu2
n_envelope(n_ind)=A;
elseif nu <= nu3
n_envelope(n_ind)=A*(epsilon/A)^((nu-nu2)/(nu3-nu2));
else
n_envelope(n_ind)=0;
end
end
%% We do not want Hf or system matrix entries to be overly small
% (in worst case, close to eps), so we normalize to
% coefficients of "L^1(R)-norm 1:
n_envelope = n_envelope./(omega_0*sum(abs(n_envelope)));
if VERBOSE
figure
stem(N*omega_0,n_envelope);
axis tight
xlabel(’n{\cdot}\omega_0 (s)’);
ylabel(’a(p{\cdot}T’’’’)’);
title([’Spreading function {\nu} with amplitude threshold ’ ...
’{\in}=’ num2str(epsilon) ’.’])
98
PathExcludingExtension = [’Figures\’ channel ’\CoeffEnvelopes\nu’];
end
%% Function for plotting TF-localization of the chosen frame =========

function plot_Gauss_frame(alpha_par,epsilon,Q,R,a,b,Omega_c,...
Omega,g_sampled,M,channel)
Omega_epsilon=Omega; % TMP!!!
T_g=1/Omega;
Q=Q(:).’; % Make row vector
R=R(:).’; % Make row vector
%% Plot TF-localization of the chosen Gaussian window Gabor frame:

SetPlot(2)
figure
line_colours = {’k’,’r’,’b’,’g’,’m’,’c’};
line_styles = {’-’,’:’,’-.’,’--’,’:.’};
NrOfStyles=length(line_styles);
NrOfColours=length(line_colours);
dot_colours ={’k.’};
nr_of_colours=length(line_colours);
phi=linspace(0,2*pi,100);
time_axis=sqrt(-2*log(epsilon)/alpha_par);
freq_axis=sqrt(-2*alpha_par*log(epsilon))/pi;
g_ellipse_t=time_axis*cos(phi);
g_ellipse_f=freq_axis*sin(phi);
for r=R;
for q=Q;
line_style=[line_colours{mod(q-Q(1),NrOfColours)+1} ...
line_styles{mod(r-R(1),NrOfStyles)+1}];
%col_nr=mod(q+r-1,nr_of_colours)+1;
%col_nr=mod(q+r-1,nr_of_colours)+1;
plot(r*a*T_g,Omega_c+q*b*Omega,dot_colours{1});
hold on;
plot(g_ellipse_t+r*a*T_g,g_ellipse_f+Omega_c+q*b*Omega,line_style);
end
end
hold off
% Zoom in one black ellipse
r=R(1);
q=Q(1);
99
axis([r*a*T_g-time_axis r*a*T_g+time_axis ...
q*b*Omega+Omega_c-freq_axis q*b*Omega+Omega_c+freq_axis]);
xlabel(’Time (s)’);
ylabel(’Frequency (Hz)’);
% title([num2str(length(R)) ’x’ num2str(length(Q)) ...
% ’ Gaussian Gabor system essential supports. ’ ...
% ’Lattice unit cell area a\cdot b=’ num2str(a*b) ’.’]);
multiplot([’Figures/’ channel ’/GaborSystem/EssentialSupports’], ...
’landscape’)
%% Plot the window and its sample values

figure
SetPlot(1)
subplot(1,2,1)
N=ceil(time_axis/T_g);
Tcont=T_g/1000;
t=[-2*M(end)*T_g:Tcont:2*M(end)*T_g].’;
g=resample(Omega_c,Omega,M,g_sampled,t);
g_max=max(g_sampled);
% CutToEssSupp(epsilon,OversampFact,g_centerfreq, ...
% g_bandwidth,g_sample_ind,g_samples)
semilogy(t,abs(g),’b-’,M.*T_g,g_sampled,’ro’, ...
[t(1);t(end)],g_max*[epsilon;epsilon],’r:’);
axis([t(1) t(end) 1e-13 g_max]);
legend(’|g(x)|’,’{g}({\itmT_g})’,’Location’,’South’)
%%
text(-0.27*10^(-4),10^(-3.5),’$\epsilon g(0)$’,’Interpreter’,’latex’)
%title([’Transmitter window {\it g} and its nonzero Nyquist freq. ’...
% ’samples. Amplitude threshold \in=’ num2str(epsilon) ’.’])
xlabel(’t (s)’);
subplot(1,2,2)
Nstep=20;
step=Omega/2/Nstep;
xi=[-1.2*Nstep:1.2*Nstep].’.*step;
Fg=zeros(size(xi));
Fgsupport=find(abs(xi)<=Omega/2);
for mind=1:length(M)
m=M(mind);
Fg(Fgsupport)=Fg(Fgsupport) ...
+g_sampled(mind)*exp(-i*2*pi*xi(Fgsupport)*m*T_g);
end
Fg=real(Fg)*T_g;
100
gauss=sqrt(pi/alpha_par)*exp(-pi^2*xi.^2./alpha_par)*g_max;
y_min=min(gauss);
y_max=max([gauss(:);Fg(:)]);
semilogy(xi,abs(Fg),’b’,xi,gauss,’r:.’,...
[xi(1);xi(end)], ...
sqrt(pi/alpha_par)*g_max*[epsilon;epsilon],’r-.’, ...
[-Omega_epsilon/2;-Omega_epsilon/2],[y_min y_max],’b:’,...
[ Omega_epsilon/2; Omega_epsilon/2],[y_min y_max],’b:’ )
xlabel(’{\xi} (Hz)’);
axis tight
%set(’Interpreter’,’latex’);
legend(’|{\itg}^{\wedge}(\xi)|’, ...
’{\itg}_0^{\wedge}(\xi)’, ...
’Location’,’South’)
text(-0.2*10^4,10^(-7.5),’$\epsilon\hat{g}(0)$’,’Interpreter’,’latex’)
text(-1.9*10^4,10^(-11),’$-\frac{\Omega}{2}$’,’Interpreter’,’latex’)
text( 1.65*10^4,10^(-11),’$\frac{\Omega}{2}$’,’Interpreter’,’latex’)
%title([’Comparison of {\itg}^{\wedge}(\xi) with the ’ ...
% ’approximated Gaussian.’])
multiplot([’Figures/’ channel ...
’/GaborSystem/GaborWindowAndCriticalSampleValues’], ...
’a4HalfLandscape’)
101
B.14 SetPlot: Changes of line thickness, font size etc in
plots
File path: Matlab/ChannelMatrix/SetPlot.m
Short description: Changes of line thickness, font size etc in plots.
Additional explanation: Thicker lines, smaller fonts etc is good for smaller
plots and for slides.
Source code
function SetPlot(mode);
switch (mode)
case 0
disp(’ Thin-line plots selected’);
pos = [678 36 600 400];
font_size = 11;
line_width = 1;
case 1
disp(’ Thick-line plots selected’);
% pos = [400 36 800 600];
pos = [200 36 1000 700];
font_size = 16;
line_width = 2;
case 2
disp(’ Thick-line plots with even larger font selected’);
% pos = [400 36 800 600];
pos = [200 36 1000 700];
font_size = 30;
line_width = 2;
end
disp(’ ’);
set(0,’defaultaxesfontsize’,font_size);
set(0,’defaulttextfontsize’,font_size);
set(0,’DefaultAxesLineWidth’,line_width);
set(0,’DefaultlineLineWidth’,line_width);
set(0,’DefaultLineMarkerSize’,12);
102
% MATLAB defines the figure Position property as a vector.

%
% [left bottom width height]
set(0,’defaultfigurepos’,pos);
103
B.15 SincOmega: Computes the function sincω (x)
File path: Matlab/ChannelMatrix/SincOmega.m
Short description: Computes the function sincω (x) .
Additional explanation: Computes the function sincω (x), which is defined
before (2.1).
Source code
function y=SincOmega(omega,x)
% USAGE: y=SincOmega(omega,x)
%
% INPUT:
% omega = vector of length d
% x = vector of length d
%
% OUTPUT
% y=sin(pi*omega*x)./(pi*x)
y=omega.*sinc(omega.*x);
104
B.16 TFshift: Time-frequency shifts (integer time-shifts,
fast)
File path: Matlab/ChannelMatrix/TFshift.m
Short description: Time-frequency shifts (integer time-shifts, fast).
Additional explanation: Compute the time-frequency shifts (4.6b) in the sim-
plest case when Tg = Tf , that is
fq,r (kTg ) = f ((k − ra)Tg )ei2πhk−ra,qbi .
Source code
function [TFg_samples TFg_sample_ind TFg_center_freq ] ...
= TFshift(g_samples,g_sample_ind,g_center_freq, ...
g_bandwidth,t_shift_ind,nu_shift_ind)
% Applies a time-frequency shift to a function g that is completely
% determined by a finite number of Nyquist frequency input samples.
%
% USAGE: [TFg_samples TFg_sample_ind TFg_center_freq ] ...
% = TFshift(g_samples,g_sample_ind,g_center_freq, ...
% g_bandwidth,t_shift_ind,nu_shift_ind)
%
% INPUT:
% t_shift_ind = INTEGER!
% The time-shift to compute is t_shift_ind*T.
% nu_shift_ind = The freq-shift to compute is nu_shift_ind*g_bandwidth.
%
% OUTPUT
% TFg_samples = Sample values of TF-shifted g
% TFg_sample_ind = Corresponding new sample indices.
% TFg_center_freq = Corresponding new sample indices.
%
% SE ALSO: TFshiftAndResample for noninteger t_shift_ind.
if isintegers(t_shift_ind)
T=1/g_bandwidth;
TFg_sample_ind=g_sample_ind+t_shift_ind;
TFg_samples=g_samples.*exp(i*2*pi*g_sample_ind.*nu_shift_ind);
105
TFg_center_freq=g_center_freq+nu_shift_ind*g_bandwidth;
else
error([’Fifth input ’ num2str(t_shift_ind) ’is no integer. ’ ...
’Use TFshiftAndResample instead.’])
end
%% Internal function ==================================================

function b=isintegers(N)
% INPUT
% N = number, vector or matrix
%
% OUTPUT
% b = 1 if all elements in N are integers
% b = 0 otherwise.
b = min(mod(N(:),1) == 0);
106
B.17 TFshiftAndResample: Time-frequency shifts
File path: Matlab/ChannelMatrix/TFshiftAndResample.m
Short description: Time-frequency shifts.
Additional explanation: Compute the time-frequency shifts (4.6b) in the when
possibly Tg 6= Tf , that is
fq,r (kTf ) = f (kTf − raTg )ei2πhkTf −raTg ,qbΩi .
with the samples f (kTf − raTg ) computed by applying (4.5b), which will be a
finite sum formula for those f that we will consider.
Source code
function [TFg_samples TFg_sample_ind TFg_center_freq ] ...
g_center_freq,g_bandwidth,t_shift,nu_shift)
% USAGE:
% [TFg_samples TFg_sample_ind TFg_center_freq ] ...
% = TFshiftAndResample(g_samples,g_sample_ind, ...
% g_center_freq,g_bandwidth,t_shift,nu_shift)
%
% Applies a time-frequency shift to a function g that is completely
% determined by a finite number of input Nyquist frequency samples.
%
% INPUT:
% t_shift = Desired time-shift
% nu_shift = Desired frequency-shift
%
% OUTPUT
% TFg_samples = Sample values of TF-shifted g
% TFg_sample_ind = Corresponding new sample indices.
% TFg_center_freq = Corresponding new sample indices.
%
% SE ALSO: TFshift (does the same faster when t_shift/T is an integer).
T = 1/g_bandwidth;
TFg_sample_ind_start=g_sample_ind( 1 ) + floor(t_shift/T);
TFg_sample_ind_end =g_sample_ind(end) + ceil( t_shift/T);
TFg_sample_ind = [TFg_sample_ind_start:TFg_sample_ind_end].’;
107
TFg_center_freq = g_center_freq + nu_shift;
t = TFg_sample_ind.*T;
TFg_samples = exp(i*2*pi*nu_shift.*(t-t_shift)) ...
.*resample(g_center_freq,g_bandwidth,g_sample_ind, ...
g_samples,t-t_shift);
108
C Functions for two of the plots
This appendix contains two additional Matlab functions, that were used to make
some other plots in this report and not are related to the functions presented in
Appendix B.
C.1 FourierTranformDecay: Matlab-file for producing Fig-

ure 1, p. 8
File path: Matlab/FourierTransformDecays/FourierTranformDecay.m
Short description: Matlab-file for producing Figure 1, p. 8.
Source code
function F = FourierTranformDecay(xi)
global xi_global cos_degree_global
tic
if nargin == 0
xi_min=0;
xi_max=50;
xi = [xi_min:0.001:xi_max].’;
else
xi_min=min(xi(:));
xi_max=max(xi(:));
end
y_min=1e-16; % = where numerical precision seem to start suffering
SetPlot(2)
%NumericalPrecision = 1e-11; % Precision in integration

%NumericalPrecision = 1e-1; % Precision in integration
NumericalPrecision = 5e-12; % Precision in integration
x=[-1.5:0.01:1.5].’;
f_exp=f(x);
aint=ApproxIntegral(f_exp,x);
f_exp=f_exp./aint;
109
[r,c] = size(xi);
F=zeros(r,c);
for n = 1:length(xi(:))
xi_global=xi(n);
F(n)=2*quadl(@integrand,0,1,NumericalPrecision);
end
absF=abs(F)./aint;
figure(1);
plot(x,f_exp);
axis tight
xlabel(’{\itx}’);
ylabel(’{\itf}({\itx})’);
print(’-depsc2’,’..\..\Figures\DecayComparison\BumpFunction’)
figure(2);
semilogy(xi,absF);
axis([xi_min xi_max y_min max(absF)])
xlabel(’{\it\xi}’);
ylabel(’{}_{ \wedge}\newline|{\itf} ({\it\xi})|’);
grid on
print(’-depsc2’, ...
’..\..\Figures\DecayComparison\FourierTransformOfBumpFunction’)
for cos_power=1:5
cos_degree_global=cos_power;
%% Old slow version kept for error-checking
%
%for n = 1:length(xi(:))
% xi_global=xi(n);
% F(n)=2*quadl(@integrand_RaisedCos,0,1,NumericalPrecision);
%end
fRC=f_RaisedCos(x);
aint=ApproxIntegral(fRC,x);
fRC=fRC./aint;
absF = abs(Four_raised_cos(xi,cos_power))./aint;;
figure(2+2*cos_power-1);
110
plot(x,fRC);
axis tight
xlabel(’{\itx}’);
ylabel([’{\itf}({\itx})’]);
print(’-depsc2’,[’..\..\Figures\DecayComparison\RaisedCos’ ...
int2str(cos_power)])
figure(2+2*cos_power);
semilogy(xi,absF);
axis([xi_min xi_max y_min max(absF)])
%axis tight
xlabel(’{\it\xi}’);
ylabel(’{}_{ \wedge}\newline|{\itg} ({\it\xi})|’);
grid on
print(’-depsc2’, ...
[’..\..\Figures\DecayComparison\FourierTransformOfRaisedCos’...
int2str(cos_power)])
end
disp([’Execution finished in ’ num2str(toc/3600) ’hours, hurrah! :-D’])
%% -------------------------------------------------------------------------
function aif=ApproxIntegral(f,x);
dx=x(2)-x(1);
aif=sum(f)*dx;
%% Candidate function 1 -----------------------------------------------------

function y = f(x)
[r,c] = size(x);
y=zeros(r,c);
pos_out_ind=find( (x>-1) & (x<1) ); % Indices giving positive output
y(pos_out_ind)=exp(-1./(1-x(pos_out_ind).^2));
%% Candidate function 2 -----------------------------------------------------

function y = f_RaisedCos(x)
global cos_degree_global
[r,c] = size(x);
y=zeros(r,c);
111
pos_out_ind=find( (x>-1) & (x<1) ); % Indices giving positive output
y(pos_out_ind)=(1+cos(pi*x(pos_out_ind))).^cos_degree_global;
%% Integrand 1 ------------------------------------------------------------
function y = integrand(x)
global xi_global
y=f(x).*cos(2*pi*xi_global.*x);
%% Integrand 2 ------------------------------------------------------------
function y = integrand_RaisedCos(x)
global xi_global
y=f_RaisedCos(x).*cos(2*pi*xi_global.*x);
%% Fourier transform of cos-window raised to power n ----------------------
function y = Four_raised_cos(xi,n);
n=round(n); % First make sure that we rase to an integer power
if n==0
y=2*sinc(2*xi); % Four. tr. of "box with support in [-1,1]";
else
y = Four_raised_cos(xi, n-1) ...
+( Four_raised_cos(xi+1/2,n-1) ...
+Four_raised_cos(xi-1/2,n-1))./2;
end
%% ------------------------------------------------------------------------
function SetPlot(mode);
switch (mode)
case 0
disp(’ Thin-line plots selected’);
pos = [678 36 600 400];
font_size = 11;
line_width = 1;
112
case 1
disp(’ Thick-line plots selected’);
% pos = [400 36 800 600];
pos = [200 36 1000 700];
font_size = 16;
line_width = 2;
case 2
disp(’ Thick-line plots with even larger font selected’);
% pos = [400 36 800 600];
pos = [200 36 1000 700];
font_size = 30;
line_width = 2;
end
disp(’ ’);
set(0,’defaultaxesfontsize’,font_size);
set(0,’defaulttextfontsize’,font_size);
set(0,’DefaultAxesLineWidth’,line_width);
set(0,’DefaultlineLineWidth’,line_width);
% MATLAB defines the figure Position property as a vector.

%
% [left bottom width height]
set(0,’defaultfigurepos’,pos);
113
C.2 NrOfNonzeroNyquistFreqSamples: Produce Figure 6, p. 33
File path: Matlab/ChannelMatrix/NrOfNonzeroNyquistFreqSamples.m
Short description: Produce Figure 6, p. 33.
Source code
function NrOfNonzeroNyquistFreqSamples()
Mmax=100;
close all hidden
setplot(2)
figure(1)
M=zeros(Mmax,1);
n=[1:Mmax].’;
M=ceil(2*10.^(6./(n+1)).*(n+1)/pi)+1;
semilogy(n,M);
xlabel(’\itn’);
ylabel(’\itM’)
axis tight
minM=min(M)
m=find(M==minM)
grid on
[’..\..\Figures\OffDiagDecay-plots\SplineM’];
114
References
[BAB+ 01] Ernst Bonek, henrik Asplund, Conor Brennan, Christian Bergljung,
Peter Cullen, Dirk Didascalou, Patrick C.F. Eggers, José Fernandes,
Cristophe Grangeat, Ralf Heddergott, Peter Karlsson, Ralf Katten-
bach, Mikael B. Knudsen, Preben E. Mogensen, Andreas F. Molisch,
Bo Olsson, Jörg Pamp, Gert Frølund Pedersen, I. Pedersen, Martin
Steinbauer, martin Weckerle, and Thomas Zwick. Antennas and Prop-
agation, chapter 3, pages 77–306. John Wiley & Sons, 2001.
[Bel63] Philip A. Bello. Characterization of randomly

time-variant linear channels. IEEE T. Com-
mun. Syst., 11:360–393, December 1963. WWW:
http://ieeexplore.ieee.org/xpl/abs free.jsp?arNumber=1088793.
[Bel64] P. Bello. Time-frequency duality. IEEE

Trans.
Inform. Theory, 10(1):18–33, 1964. WWW:
[Beu38] A. Beurling. Sur les intégrales de Fourier absolument convergentes

et leur application à une transformation fonctionnelle. In Neuvième
congrès des mathématiciens Scandinaves, Helsingfors, 1938. Reprinted
as [Beu89, p. 39–60].
[Beu89] Arne Beurling. The Collected Works of Arne Beurling. Volume 2,

Harmonic Analysis, volume 2. Birkhäuser, Boston ; Basel, 1989.
[BLM04] Imad Barhumi, Geert Leus, and Marc Moonen. Time-domain and
frequency-domain per-tone equalization for OFDM over doubly se-
lective channels. Signal Process., 84(11):2055–2066, November 2004.
WWW: http://www.sciencedirect.com/science/journal/01651684.
[BLM05] Imad Barhumi, Geert Leus, and Marc Moonen. Time-

varying FIR equalization for doubly selective channels. IEEE
T. Wirel. Commun., 4(1):202–214, January 2005. WWW:
http://ieeexplore.ieee.org/xpl/abs free.jsp?arnumber=1381438.
[BS01] Scott Beaver and Thomas Strohmer. Optimal OFDM pulse and lat-
tice design for doubly dispersive channels. In Conference Record of
the Thirty-Fifth Asilomar Conference on Signals, Systems and Com-
puters, pages 1763–1767, Pacific Grove, CA, USA, 2001. WWW:
http://www.math.ucdavis.edu/∼strohmer.
[Chr02] Ole Christensen. An Introduction to Frames and Riesz Bases.

Birkhäuser, 2002.
[Con00] John B. Conway. A course in operator theory. 2000.
115
[DH98] Jacek Dziubański and Eugenio Hernández. Band-
limited wavelets with subexponential decay. Canad.
Math. Bull., 41(4):398–403, 1998. WWW:
http://www.journals.cms.math.ca/CMB/v41n4/index.en.html.
[Dom56] Y. Domar. Harmonic analysis based on certain commutative Banach

algebras. Acta Math., 96:1–66, 1956.
[EG04] Stefan Ericsson and Niklas Grip. Efficient wavelet pre-

filters with optimal time-shifts. IEEE Trans. Sig-
nal Process., 53(7):2451–2461, July 2004. WWW:
http://www.sm.luth.se/∼grip/Research/publications.php.
[EG05] Stefan Ericsson and Niklas Grip. An analysis method for

sampling in shift-invariant spaces. Int. J. Wavelets Multires-
olut. Inf. Process., 3(3):301–319, September 2005. WWW:
[FK98] Hans G. Feichtinger and Werner Kozek. Quantization of TF lattice–

invariant operators on elementary LCA groups. In Hans G. Feichtinger
and Thomas Strohmer, editors, Gabor Analysis and Algorithms, chap-
ter 7, pages 233–266. Birkhäuser, Boston, MA, USA, 1998.
[Fol95] Gerald B Folland. A Course in Abstract Harmonic Analysis. Studies

in advanced mathematics. CRC Press, Boca Raton, FL, 1995.
[Fol99] Gerald B. Folland. Real analysis. Modern techniques and their appli-
cations. Pure and Applied Mathematics. A Wiley-Interscience Series
of Texts, Monographs, and Tracts., NY, second edition, 1999.
[FZ98] Hans G. Feichtinger and Georg Zimmermann. A Banach space of test

functions for Gabor analysis. In Hans G. Feichtinger and Thomas
Strohmer, editors, Gabor Analysis and Algorithms, chapter 3, pages
123–170. Birkhäuser, Boston, MA, USA, 1998.
[GHS+ 01] Calle Gustavsson, Christian Henriksson, Mårten Sander, Pe-

ter Sidén, Roland Standert, and Andreas Vedin. Full du-
plex OFDM modem over a frequency selective channel.
Project course report, Department of Signals, Sensors and
Systems, KTH, Stockholm, Sweden, May 2001. WWW:
http://www.s3.kth.se/signal/edu/projekt/students/01/lightbrown/www/final report.PDF.
[GP07] Niklas Grip and Götz Pfander. A discrete model for the efficient
analysis of time-varying narrowband communication channels. Mul-
tidim. Syst. Sign. Process., To appear:38 pages, 2007. WWW:
http://www.sm.luth.se/∼grip/Research/publications.php#GrPf07a.
116
[Gri02] Niklas Grip. Wavelet and Gabor Frames and Bases: Approx-
imation, Sampling and Applications. Doctoral thesis 2002:49,
Luleå University of Technology, SE-971 87 Luleå, 2002. WWW:
[Grö00] Karlheinz Gröchenig. Foundations of Time-Frequency Analysis.

Birkhäuser, 2000.
[Hör03] Lars Hörmander. The Analysis of Linear Partial Differential Opera-

tors I; Distribution Theory and Fourier Analysis. Classics in mathe-
matics. Springer, second edition, 2003.
[Jan98] A. J. E. M. Janssen. The duality condition for Weyl-Heisenberg

frames. In Hans G. Feichtinger and Thomas Strohmer, editors, Gabor
Analysis and Algorithms, chapter 1, pages 33–84. Birkhäuser, Boston,
MA, USA, 1998.
[JYGR74] W. C. Jakes, Y. S. Yeh, M. J. Gans, and D. O. Reudnik. Large-

scale variations of the average signal. In William C. Jakes, editor,
Microwave Mobile Communications, chapter 5, pages 309–387. John
Wiley & Sons, NY, 1974. IEEE 1994 reprint with corrections.
[Kat04] Yitzhak Katznelson. An Introduction to Harmonic Analysis. Cam-

bridge University Press, third corrected edition, 2004.
[KM98] Werner Kozek and Andreas F. Molisch. Nonorthogonal pulse shapes

for multicarrier communications in doubly dispersive channels. IEEE
J. Sel. Area. Comm., 16(8):1579–1589, October 1998. WWW:
http://ieeexplore.ieee.org/iel4/49/15739/00730463.pdf.
[Koz97] Werner Kozek. Matched Weyl-Heisenberg Expansions

of Nonstationary Environments. Ph.D. thesis, Tech-
nical University of Vienna, March 1997. WWW:
http://www.mat.univie.ac.at/ nuhag/papers/1997/koz0897.html.
[KP06] W. Kozek and G.E. Pfander. Identification of operators with bandlim-

ited symbols. SIAM J. Math. Anal., 37(3):867–888, 2006 2006.
[KPZ02] Werner Kozek, Götz Pfander, and Georg Zimmermann. Perturba-

tion stability of coherent Riesz systems under convolution operators.
Appl. Comput. Harmon. Anal., 12(3):286–308, May 2002. WWW:
http://www.math.iu-bremen.de/pfander/publications.php.
[LM04] Geert Leus and Marc Moonen. Equalization techniques for fading
channels. In Mohamed Ibnkahla, editor, Handbook on Signal Process-
ing for Communications, chapter 16, pages 16.1–16.30. CRC Press,
2004.
117
[LO97] W. K. Lam and R. F. Ormondroyd. A coherent COFDM modula-
tion system for a time-varying frequency-selective underwater acous-
tic channel. In Seventh International Conference on Electronic En-
gineering in Oceanography : technology transfer from research to in-
dustry, pages 198–203, Southampton, UK, June 1997. IEEE. WWW:
[LPW05] James Lawrence, Götz Pfander, and David Walnut. Linear in-
dependence of Gabor systems in finite dimensional vector spaces.
J. Fourier Anal. Appl., 11(6):715–726, 2005. To appear. WWW:
http://www.math.iu-bremen.de/pfander/publications.php.
[LZGk03] Geert Leus, Shengli Zhou, and Georgios B. Giannakis. Orthogonal

multiple access over time- and frequency-selective channels. IEEE
Trans. Inform. Theory, 49(8):1942–1950, August 2003. WWW:
[Mat00] Gerald Matz. A time-frequency calculus for time-varying systems

and nonstationary processes with applications. Ph.d. thesis, E
389, Vienna University of Technology, November 2000. WWW:
http://www.lss.supelec.fr/perso/matz gerald/GM other.html.
[MB02] Gérard Maral and Michel Bousquet. Satellite Communications Sys-

tems: Systems, Techniques and Technology. John Wiley & Sons,
fourth edition, 2002.
[MGk02] Xiaoli Ma and Georgios B. Giannakis. Maximum-diversity

transmissions over time-selective wireless channels. In Wire-
less Communications and Networking Conference, 2002.
(WCNC2002), volume 1, pages 497–501, March 2002. WWW:
[MGk03a] Xiaoli Ma and Georgios B. Giannakis. Full-diversity full-

rate complex-field space-time coding. IEEE Trans. Sig-
nal Process., 51(11):2917–2930, November 2003. WWW:
[MGk03b] Xiaoli Ma and Georgios B. Giannakis. Maximum-diversity

transmissions over doubly selective wireless channels. IEEE
Trans. Inform. Theory, 49(7):1832–1840, July 2003. WWW:
[Mid87] D. Middleton. Channel modeling and threshold signal pro-

cessing in underwater acoustics: An analytical overview.
IEEE J. Oceanic Eng., 12(1):4–28, January 1987. WWW:
118
[MLGk05] Xiaoli Ma, Geert Leus, and Georgios B. Giannakis. Space-time-
Doppler block coding for correlated time-selective fading channels.
IEEE Trans. Signal Process., 53(6):2167–2181, June 2005. WWW:
[MMH+ 02] Gerald Matz, Andreas F. Molisch, Franz Hlawatsch, Mar-

tin Steinbauer, and Ingo Gaspard. On the systematic mea-
surement errors of correlative mobile radio channel sounders.
IEEE T. Commun., 50(5):808–821, May 2002. WWW:
http://www.lss.supelec.fr/perso/matz gerald/GM jrnl.html.
[MSG+ 05] G. Matz, D. Schafhuber, K. Gröchenig, M. Hartmann, and

F. Hlawatsch. Analysis, optimization, and implementation of low-
interference wireless multicarrier systems. preprint (submitted to IEEE
Transactions on Wireless Communications), 2005.
[PW06] Götz E. Pfander and David F. Walnut. Measure-

ment of time-variant linear channels. IEEE Trans. In-
form. Theory, 52(11):4808–4820, November 2006. WWW:
http://ieeexplore.ieee.org/iel5/18/36107/01715527.pdf?isnumber=36107&arnumber=1715527.
[Rap02] Theodore S. Rappaport. Wireless communications: principles and

practice. Prentice Hall PTR, second edition, 2002.
[Reu74] D. O. Reudnik. Large-scale variations of the average signal. In

William C. Jakes, editor, Microwave Mobile Communications, chap-
ter 2, pages 79–131. John Wiley & Sons, NY, 1974. IEEE 1994 reprint
with corrections.
[Ric03] Scott Rickard. Time-frequency and time-scale representations of dou-

bly spread channels. Ph.D. dissertation, Princeton University, Novem-
ber 2003.
[RS80] Michael Reed and Barry Simon. Methods of modern mathematical

physics; I: Functional analysis. Academic Press, NY, revised and
enlarged edition, 1980.
[RS00] Hans Reiter and Jan D. Stegeman. Classical Harmonic Analysis and
Locally Compact Groups. Number 22 in London Mathematical Soci-
ety monographs; New Series. Oxford University Press, second edition,
2000.
[Rud87] Walter Rudin. Real and Complex Analysis. McGraw-Hill, third edi-
tion, 1987.
[SA99] Akbar M. Sayeed and Behnaam Aazhang. Joint multipath-

Doppler diversity in mobile wireless communications. IEEE
119
T. Commun., 47(1):123–132, January 1999. WWW:
[Sto96] Milica Stojanovic. Recent advances in high-speed

underwater acoustic communications. IEEE J.
Oceanic Eng., 21(2):125–136, April 1996. WWW:
http://intl.ieeexplore.ieee.org/xpl/tocresult.jsp?isNumber=27330&puNumber=26.
[Sto99] Milica Stojanovic. Underwater Acoustic Communications, vol-

ume 22, pages 688–698. John Wiley & Sons, 1999. WWW:
http://www.mit.edu/people/millitsa/publications.html.
[Str06] Thomas Strohmer. Pseudodifferential operators and Ba-

nach algebras in mobile communications. Appl. Com-
put. Harmon. Anal., 20(2):237–249, March 2006. WWW:
http://www.math.ucdavis.edu/∼strohmer.
[TL04] Jitendra K. Tugnait and Weilin Luo. Blind identification of

time-varying channels using multistep linear predictors. IEEE
Trans. Signal Process., 52(6):1739–1749, June 2004. WWW:
[Zad50] Lofti A. Zadeh. Frequency analysis of variable networks. Proc. IRE,

38(3):291–299, March 1950.
[ZK00] Y.V. Zakharov and V.P. Kodanev. Multipath-Doppler diversity of

OFDM signals in an underwater acoustic channel. In Proceedings of
ICASSP’2000, volume 5, pages 2941–2944. IEEE, June 2000. WWW:
[ZT02] M. Zatman and B. Tracey. Underwater acoustic mimo

channel capacity. In Conference Record of the Thirty-Sixth
Asilomar Conference on Signals, Systems and Computers, vol-
ume 2, pages 1364–1368. IEEE, November 2002. WWW:
http://intl.ieeexplore.ieee.org/xpl/abs free.jsp?arNumber=1197002.
Niklas Grip received his Ph.D. degree in Mathematics from the Luleå University
of Technology, in Sweden 1997, where he now has a junior research fellow position.
His main scientific interests are in the fields of wavelets, Fourier and Gabor analysis
with applications to signal processing and multicarrier communications.
Götz E. Pfander received his Ph.D. degree in Mathematics from the University
of Maryland in 1999. Since 2002 he is an Assistant Professor of Mathematics
at Jacobs University Bremen, Germany. His interests include harmonic analysis,
wavelet and Gabor theory, and signal processing.
120

Time Varying Narrowband Communications Channels: Analysis and Implementation

Uploaded by

Copyright:

Available Formats

Time Varying Narrowband Communications Channels: Analysis and Implementation

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Time Varying Narrowband Communications Channels: Analysis and Implementation

Uploaded by

Copyright:

Available Formats

Niklas Grip and Götz E.

Time Varying Narrowband Communications

Technical Report No. 12

School of Engineering and Science

Niklas Grip Götz E. Pfander

Department of Mathematics School of Engineering and Science

E-Mail: [email protected], [email protected],

3 The channel model 11

4 Discretization of the channel model 19

5 The algorithm and its implementation 31

A Fourier transform decay vs differentiability 44

B Matlab implementation: overview and function reference 46

C Functions for two of the plots 109

1. We review the most important physical properties of wireless channels and

2. We employ the channel model described above to derive an efficient algo-

Here, sincω is extended continuously to Rd and we shall frequently use that

For ² > 0 we define the ²-essential support of a bounded function f : Rd → C to

w(x) ≥ 1 and w(x + y) ≤ w(x)w(y). (2.2)

2.2 Frequency localization of compactly supported func-

2.3 Gabor analysis

The corresponding synthesis or reconstruction operator

Rg∗ : L2 → l2 , Rg∗ f = (hf, gq,r i)q,r∈Zd ,

We will provide an explicit expression for G in Section 4.

2.4 Hilbert–Schmidt operators

Definition 1 (Hilbert–Schmidt operator). Let H be a separable Hilbert

Proposition 1 ([RS80, Theorem VI.23]). Let {M, µ} be a measure space and

Within this report, we shall only consider Hilbert–Schmidt operators on the

Corollary 1. A linear bounded operator H : L2 (Rd ) → L2 (Rd ) is Hilbert–Schmidt

The integral in (2.9) is defined in a weak sense. In fact, for s, g ∈ L2 (Rd ),

on L2 (Rd × Rd ). Hence, Hs is defined to be the unique L2 (Rd )-function with

hHs, giL2 (Rd ) = hSH , Vs giL2 (Rd ×Rd ) .

3 The channel model

3.1 Single path frequency response

lP (t) = LP + VP · (t − T0 ) with |VP · (t − T0 )| ¿ LP . (3.1)

Vw > Vmax ≥ |VP | for all paths P.

This mapping from (VP , LP ) to (νP , tP ) is invertible with inverse

VP = −νP Vw and LP = Vw (tP (1 + νP ) − νP T0 ) . (3.2c)

1. A multiplication by a transmitter amplitude gain GT (P).

2. The time-frequency shift by (νP ξ, tP ) that is given in (3.2a).

Thus, in this mathematical model, (3.5) takes the form

3.2 Narrowband signals

νξ ≈ νξ0 , e−αξ t ≈ e−αξ0 t and ρξ ≈ ρξ0 (3.8)

Then, at least for bandlimited and thus continuous narrowband L1 -signals s, we

We use the last expression as definition of our mathematical model

If s ∈ L2 (R), then by (3.9b) and the Minkowski integral inequality

≤ |ρξ0 | (K × [A, ∞))e−αξ0 A ksk2 .

Alternatively we can use (3.9b) as the definition of a larger space of opera-

(νξ0 , t) 7→ (Tt Mνξ0 s)(t0 ) ∈ S0 (K 0 × [A, ∞)) for all t0 , (3.10)

is well-defined, it follows that (3.9b) is pointwise well-defined for all t0 . Conse-

Time-invariant systems: h(x, t) = h(0, t) and SH (ν, t) = h(t)δ0 (ν).

3.3 Finite lifelength channels

BH (ν, ·) is compactly supported. (3.11)

Figure 3: (a) Bandlimiting properties of the physical channel provides us with

4 Discretization of the channel model

with biinfinite matrix representation

Gq0 ,r0 ;q,r = hHgq,r , γq0 ,r0 i , (4.2)

4.1 Gabor bases for near-diagonalization

Figure 4: Convolution of a compactly supported complex exponential s(t) =

and all the inner products Hgq,r , Hg ^ q 0 ,r 0 .