# Functional Programming in Learning Electromagnetic Theory

Scott N. Walck Department of Physics
Lebanon Valley College
Annville, Pennsylvania, USA<EMAIL_ADDRESS>
###### Abstract

Electromagnetic theory is central to physics. An undergraduate major in
physics typically takes a semester or a year of electromagnetic theory as a
junior or senior, and a graduate student in physics typically takes an
additional semester or year at a more advanced level. In fall 2023, the author
taught his undergraduate electricity and magnetism class using numerical
methods in Haskell in parallel with traditional analytical methods. This
article describes what functional programming has to offer to physics in
general, and electromagnetic theory in particular. We give examples from
vector calculus, the mathematical language in which electromagnetic theory is
expressed, and electromagnetic theory itself.

##  1 Introduction

James Clerk Maxwell put the finishing touches on modern electromagnetic theory
in 1865, publishing the famous equations that still bear his name. It is the
oldest piece of theoretical physics that is part of humanity’s current best
understanding of the physical world.

Electromagnetic (EM) theory describes one of the four fundamental forces of
nature. Today’s physicists believe that all known interactions can be
classified as one of four forces. These are the strong force, the
electromagnetic force, the weak force, and gravity. Electricity holds atoms
together, which seems enough to give it a claim to importance.

Electromagnetic theory serves as the model for modern field theories of
elementary particles. EM theory is the prototypical example of a _field
theory_ , that is, a theory whose important quantities depend on space or
spacetime. Modern theories of elementary particle physics are based on quantum
field theories, the quantum versions of theories like electromagnetic theory.

Surprisingly, Maxwell’s 1865 electromagnetic theory did not need to be
modified to incorporate Einstein’s 1905 relativity theory or the quantum
theory of 1925. In the case of relativity theory, electromagnetic waves began
to be interpreted differently, as physicists discarded the notion of an ether
in which EM waves traveled. In the case of quantum theory, mathematical
objects like electric field and magnetic field were interpreted differently,
as operators rather than numbers or vectors, but the Maxwell equations
remained intact. It seems that Faraday and Maxwell were onto something
fundamental about the universe.

Electromagnetic theory is the earliest theory that is still part of our
current best understanding of the universe. Newtonian mechanics is incredibly
useful, and beautiful, but it cannot be said to express our current best ideas
about the universe. The 20th-century ideas of quantum mechanics and relativity
have each led to newer theories that are slightly different from Newtonian
mechanics, although they give essentially the same answers for slow massive
things. Electromagnetic theory, on the other hand, passed unchanged through
the 20th-century quantum and relativity revolutions.

Electromagnetic theory unites electricity, magnetism, and light into a single
theory. Light is an electromagnetic wave, a wave of electric and magnetic
fields. Electromagnetic theory is three theories in one, describing
electricity, magnetism, and optics.

We take the attitude of Papert[4] and others[5, 6, 7, 8] that students are
aided in their learning by having building blocks with which to create
interesting structures, that such creative activity is a motivating and
effective way to learn, and that the feedback provided by computer-language-
based building blocks can expose our confusions and produce delight in our
achievements.

The author has written previously about the benefits of functional programming
to physics. In [7], we showed how Newtonian mechanics benefits from expression
in a functional programming language. In [8], we showed how a course in
quantum mechanics can benefit from functional programming. The textbook [9]
treats both undergraduate Newtonian mechanics and electromagnetic theory. It
can be used as a textbook for a course in computational physics or as a
companion text for courses in classical mechanics or electromagnetism. The
author of this paper uses [9] in both ways. The present paper focuses on
functional programming in electromagnetic theory, and in particular on the
author’s experience incorporating it into a traditional undergraduate
electromagnetic theory course.

The plan for the paper is as follows. In Section 2, we describe what
functional programming has to offer to physics in general and to
electromagnetic theory in particular. Section 3 describes field theory and how
it can be encoded in Haskell. In Section 4, we give examples of what vector
calculus can look like in a functional programming language. Section 5
describes an example of electromagnetic theory, namely how magnetic field is
produced by electric current.

##  2 What Does Functional Programming Offer to Physics Pedagogy?

In this section, we describe how programming in general, and functional
programming in particular, can help people learn physics, especially
electromagnetic theory. The first three items deal with programming in
general, while the last five focus on functional programming.

###  2.1 With a little programming, we can study physical situations that are
not exactly solvable.

Most physical situations are not exactly solvable. The traditional tools of
algebra, calculus, and differential equations only get us so far. They can’t
solve most problems. The theory tells you how to make (partial) differential
equations, but you can’t solve them exactly. We don’t know how to write down
functions that exactly satisfy the differential equations. Traditionally, you
spend a lot of time studying situations that _can_ be solved exactly, and then
come up with tricks to approximate the other situations that can’t be solved
exactly.

In EM theory, the electric field inside a capacitor with infinite plates is
exactly solvable, but the electric field inside a capacitor with finite plates
is not. The magnetic field produced by a circular loop carrying current is
exactly solvable along the axis that runs through the center of the circle
(perpendicular to the plane of the circle), but not exactly solvable anywhere
else.

There are some simple, standard approximations (numerical methods), that are
known to be reasonably good, and don’t involve tricks. The problem is that
they require an enormous number of exceedingly boring computations. And here
is where the computer comes to the rescue. Just by being able and willing to
do an enormous number of simple, boring computations, the computer can give us
good, approximate solutions to many problems.

Figure 1 shows the numerically calculated magnetic field of a circular current
loop, a problem which, surprisingly, is not analytically solvable.

Figure 1: Magnetic field produced by a current loop. The current loop lies in
the $xy$ plane. The figure shows the $yz$ plane. A darker arrow indicates a
stronger magnetic field. The current loop (not shown) pierces the $yz$ plane
near the darkest arrows. This situation, simple as it seems, is exactly
solvable only on the $z$ axis, running vertically through the center of the
figure.

###  2.2 Computer programming is a modern tool that can help solve problems.

Why not use all the tools at our disposal, and in particular the computer, to
help us solve our problems? Computers can do things that calculators can’t. I
don’t predict the death of the pocket calculator, but who knows? The slide
rule is, if not dead, an eclectic tool used by very few people.

Physics societies like the American Association of Physics Teachers (AAPT)
have been encouraging physics teachers to include more computer techniques in
their courses.

###  2.3 Programming is a valuable skill in itself.

Programming, like mathematics and like writing, is a valuable skill in itself.
Knowing something about programming is useful for getting a job, but it’s also
useful for organizing your thinking. (Mathematics and writing also help us
organize our thinking.) Writing code is not just about getting the computer to
do something. It’s also about expressing ideas using formal language in a way
that makes sense to people. And so, for the same reason that essay writing can
help you clarify your ideas about art or politics, code writing can help you
clarify your ideas about physics.

###  2.4 Writing in a language that will be read by a computer forces a
precision that human language and mathematical notation sometimes lack.

In the preface to their wonderful book _Structure and Interpretation of
Classical Mechanics_ (SICM)[5], Sussman and Wisdom note

> The traditional use of ambiguous notation is convenient in simple
> situations, but in more complicated situations it can be a serious handicap
> to clear reasoning. In order that the reasoning be clear and unambiguous, we
> have adopted a more precise mathematical notation. Our notation is
> functional and follows that of modern mathematical presentations.

and then

> We require that our mathematical notations be explicit and precise enough
> that they can be interpreted automatically, as by a computer.

Sussman was the first pioneer in the application of functional programming to
physics education. He showed the power of higher-order functions in physics
using _Scheme_ , a dialect of Lisp he co-invented with Steele, and used with
Abelson in their classic computer science text _Structure and Interpretation
of Computer Programs_[2].

The approach described in this paper is similar to SICM in that both are
interested in using functional programming as a language for physics. Our
approach differs from SICM in several ways. First, we use Haskell instead of
Scheme, a difference which is noticeable primarily in our reliance on types
and Haskell’s type system. Second, SICM is a graduate-level textbook while we
target undergraduates. Third, SICM restricts its attention to classical
mechanics, the theory of motion created by Newton and developed by Lagrange
and Hamilton, while our concern in this paper is with a different theory of
physics, namely electromagnetic theory.

###  2.5 Avoiding mutation makes a programming language closer to mathematics.

The core of a pure functional programming language, while providing many
interesting challenges for the compiler designer, presents a very simple model
of computation to the user of the language. Names stand for values, not memory
locations. Names refer to a single thing, and don’t change over time.

These are the principles upon which mathematical notation is based. Physics
already uses mathematical language, so it is a smaller step for the physics
student to learn the core of a pure functional programming language than it is
to learn an imperative language. Obviously, there are trade-offs; functional
programming’s affection for recursion can make teaching iteration more
difficult.

In my courses, I teach students how to use the Haskell functions map, iterate,
filter, take, and takeWhile as well as list comprehensions, but I do not teach
them how to write explicitly recursive functions. Teaching recursion is a
substantial endeavor, and I need to use the time to teach ideas of physics.

###  2.6 Types in a language help the reader understand and the writer clarify
and organize.

I have come to believe that a strong system of types is exceedingly helpful in
a language for physics. Algebraic data types are useful in a language for
physics for the same reason they are helpful in any domain-specific language,
namely that the key ideas of the domain can be encoded as types.

The code writer can experience the benefits of clarity and organization of
types even if she is not the person writing the data type definition.

An example of an algebraic data type in physics is _charge distribution_.
Physicists like to talk about point charges, but they also like to talk about
electric charge that is distributed along a curve, across a surface, or
throughout a volume. A _charge distribution_ is a specification of electric
charge that could be at a point, along a curve, across a surface, throughout a
volume, or some combination of these.

We use an algebraic data type to make charge distribution into a type.

    data ChargeDistribution
        = PointCharge   Charge      Position
        | LineCharge    ScalarField Curve
        | SurfaceCharge ScalarField Surface
        | VolumeCharge  ScalarField Volume
        | MultipleCharges [ChargeDistribution]

The definition tells us that a charge distribution could be a point charge, a
line charge (physicists use the term _line charge_ for any one-dimensional
distribution of charge; it need not occur along a straight line), a surface
charge, a volume charge, or a combination of these. The ChargeDistribution
data type is algebraic in having multiple constructors, and recursive in the
last constructor.

We can see from its definition that the ChargeDistribution data type is based
on other data types for physics like Charge, Position, ScalarField, Curve,
Surface, and Volume. Charge is just a type synonym for double-precision
floating point number. Position is a point in three-dimensional space. We’ll
discuss the Position data type in the next section. Curve, Surface, and Volume
are geometric objects that describe finite curves, surfaces, and volumes.
These are helpful because EM theory has a substantial geometric content.

I do not demand that students be able to write algebraic data type definitions
like this. Some students are interested in how this works, and I’m always
happy to share what I know, but my main interest is in having students use
this data type to describe a charge distribution. Using this data type brings
one face to face with precisely what a charge distribution is, how we talk
about it, and what’s it’s good for.

In EM theory, scalar and vector fields play a prominent role. A _field_ in
physics is a function of space or spacetime. In an undergraduate setting, it
is more common to treat fields as functions of space, and I follow that custom
with the following type synonyms for scalar field and vector field.

    type ScalarField = Position -> R
    type VectorField = Position -> Vec

In the first definition, R is a type synonym for Double, a double-precision
floating point number. I like to think of these numbers as approximations to
real numbers, and they occur so often in physics that it is convenient to give
them a shorter name that comes closer to how we think about them. In the
second definition, Vec is a data type for three-dimensional vectors. Such
vectors are ubiquitous in physics, especially at the undergraduate pre-
relativity level.

The central idea of electrostatics, the portion of EM theory focusing on
electrical phenomena when nothing is changing in time, is that a charge
distribution produces an electric field. We encode this central idea as a
function.

    eField :: ChargeDistribution -> VectorField

The type signature of eField sings the central idea of electrostatics. The
definition of eField takes some effort to develop, as we need various kinds of
integrals in order to deal with the continuous charge distributions. The type
signature, however, serves as a wonderful summary and reminder of the big
picture, and hopefully serves as a counterbalance to the possibility that
people might get lost in the details of the function definition.

In fact, students routinely get lost in the traditional analytical
presentation of electric field produced by a line charge, electric field
produced by a surface charge, and electric field produced by a volume charge.
There are so many details to worry about, and nothing in the traditional
mathematical notation is reminding us of the central idea. A language with a
strong system of types can help us keep the important ideas in mind.

###  2.7 Higher-order functions are a natural way to express the higher-order
ideas that physics trades in.

Take the derivative from calculus as a higher-order idea that physics uses.
The derivative of a function at a point expresses the rate at which the
function changes at that point. We can write a numerical derivative in
Haskell.

    type R = Double

    derivative :: R -> (R -> R) -> R -> R
    derivative dt f t = (f (t + dt/2) - f (t - dt/2)) / dt

In this definition, dt is a numerical step size for the numerical derivative
we are calculating. The second argument is the function f :: R -> R we want to
differentiate. The third argument t :: R is the value of the independent
variable where we want to estimate the slope.

Physicists like to think of the derivative as an “operator” that takes a
function as input and gives another function as output. Haskell’s currying
allows us to think just like that. Instead of thinking of the derivative as a
function that takes three inputs (step size, function, and independent
variable), we can think of derivative as a function that takes two inputs
(step size and function) and returns a function. Even better, we can come to
see that these two ways of thinking are the same for the computer because
there is a sense in which there is no difference. The two ways of thinking are
equivalent. Better yet, the function derivative 1e-6 really is the “operator”
that takes a function as input and gives a function as output.

Numeric integration is another example of a higher-order function that encodes
a basic idea of physics. And since calculus, vector calculus, and differential
equations are so ubiquitous in physics, there are a plethora of higher-order
functions serving calculus-related roles, some of which we will see later in
this article.

Electromagnetic theory is primarily about vector fields, namely the electric
and magnetic fields. We saw above that a vector field is a function, so any
function that deals with a vector field is a higher-order function.

In physics, ideas we regard as basic, like the derivative, the integral, and
the vector field, naturally invite us to use higher-order functions. Because
higher-order functions in Haskell and other functional programming languages
are easy to write and use, the language of the code can help you understand
the theory. When a language allows a succinct expression of an idea that a
domain regards as basic, using that language allows one to experiment and play
with the idea, which is a great way to learn. Reading and writing in a
language that makes it easy (or at least easier) to express the ideas in a
theory promotes understanding of the theory. I claim that typed functional
languages are especially suitable for physics, and I show examples below as
evidence of this claim. I have anecdotal evidence that students understand
analytical methods in electromagnetic theory better for having studied
numerical methods in a functional programming language.

###  2.8 Functional programming offers a vision of a language for both
calculating and proving, two activities that theoretical physics is full of.

Physics has some results that can be obtained from deductive reasoning. The
Maxwell equations imply conservation of electric charge, for example. In a
course in EM theory, we show how to derive the continuity equation, an
expression of local charge conservation, from the Maxwell equations.

This derivation takes a handful of steps, and although it doesn’t rise to the
stature of a theorem as appreciated by mathematicians, nevertheless the
language and pattern of a theorem are present. If the Maxwell equations are
true, then the continuity equation is true. Such relationships are helpful to
recognize, because a deep understanding of any theory requires going beyond
“what is true” into knowing when ideas are independent, and knowing when one
idea is a logical consequence of another.

Haskell, the language that I use in my EM theory class, is great for
calculation, but not really intended for theorem proving. However, a handful
of functional languages based on dependent types, such as Coq, Agda, Idris,
and Lean, are capable of theorem proving (meaning proof checking with some
automation) in addition to calculation.

It seems that a language that could calculate and prove would be a wonderful
language in which to express the principles of physics. While I have great
interest in this possibility, it is in no sense a reality in my teaching.

It seems to me that the possibility of using a dependently-typed functional
language for physics in the future is a reason to learn a functional language
like Haskell for physics now.

##  3 Field Theory in Haskell

Electromagnetic theory is a _field theory_ , where a _field_ in physics is a
quantity that depends on location in space or spacetime. (The definition of
_field_ in mathematics is different, and unrelated.) The most common
quantities that depend on space or spacetime are scalars (basically numbers)
and vectors. In undergraduate electromagnetic theory, it is common to view
fields as functions of space, and we will follow this practice. A scalar
quantity that depends on space is called a _scalar field_ , and a vector
quantity that depends on space is called a _vector field_. We gave type
definitions for ScalarField and VectorField in section 2.6.

EM theory can be expressed (quite beautifully, in fact) in relativistic
notation where the vectors are 4-vectors (elements of a 4-dimensional vector
space) and the fields are functions of spacetime, but in an undergraduate
setting, it is much more common to use a notation in which vectors are
3-vectors (members of a 3-dimensional vector space), and the fields are
functions of space. This does not preclude the possibility that fields can
change in time; it merely places space and time on an apparently different
footing.

Having decided to pursue an undergraduate notation in which 3-dimensional
space is central, we turn to the question of how to represent points in
3-dimensional space. Our first choice is Cartesian coordinates $x,y,z$, in
which the three coordinate axes are mutually perpendicular. However, several
important situations in EM theory have spherical or cylindrical symmetry,
which motivates the use of spherical and cylindrical coordinates. Even
situations which do not have full cylindrical or spherical symmetry can
benefit from these coordinates. In fact, there are many possible coordinate
systems for 3-dimensional space, but we will focus on the 3 most commonly
used, namely Cartesian, cylindrical, and spherical. Figure 2 shows how
cylindrical and spherical coordinates are defined.

$y$$z$$x$$s$$\phi$$z$$\hat{z}$$\hat{s}$$\hat{\phi}$

$y$$z$$x$$r$$\theta$$\phi$$\hat{r}$$\hat{\theta}$$\hat{\phi}$

Figure 2: Cylindrical coordinates (left) and spherical coordinates (right)

In cylindrical coordinates, we use $s$ to represent the distance from the $z$
axis, $\phi$ to represent the angle shown in Figure 2, and $z$ to mean the
same thing it means in Cartesian coordinates. The cylindrical coordinates of a
point in space are then given by the triple $(s,\phi,z)$. In spherical
coordinates, we use $r$ to represent the distance from the origin, $\theta$ to
represent the angle shown in Figure 2, and $\phi$ to mean the same thing it
means in cylindrical coordinates. The spherical coordinates of a point in
space are then given by the triple $(r,\theta,\phi)$. This notation for
cylindrical and spherical coordinates comes from Griffiths’ popular
textbook[3].

The purpose of the Position data type is to free our minds from the detail of
establishing a single, canonical way of referring to a point in space. We wish
to freely use any of the three coordinate systems in constructing a Position,
and any of the three in extracting coordinates. Here is the definition of the
Position data type.

    data Position = Cart R R R
                    deriving (Show)

We store coordinates in Cartesian form, and name the constructor Cart to
remind us of this fact, but users of Position need not be aware of the
specifics of how the coordinates are stored.

The following three functions construct a Position from coordinates.

    cartesian   :: (R,R,R) -> Position   -- (x,y,z)
    cylindrical :: (R,R,R) -> Position   -- (s,phi,z)
    spherical   :: (R,R,R) -> Position   -- (r,theta,phi)

For convenience, we define curried forms cart, cyl, and sph.

Finally, we have three functions to extract coordinates from a Position in any
of the three systems.

    cartesianCoordinates   :: Position -> (R,R,R)   -- (x,y,z)
    cylindricalCoordinates :: Position -> (R,R,R)   -- (s,phi,z)
    sphericalCoordinates   :: Position -> (R,R,R)   -- (r,theta,phi)

Since the VectorField is the type of the electric field and the magnetic
field, and since we can develop some geometric intuition by making pictures of
vectors fields, we spend some time doing this in my class. The function vfGrad
uses a technique I call _gradient visualization_ , in which the strength of
the field is denoted by the darkness of arrows. This is the kind of
visualization used in Figure 1 for the magnetic field produced by a circular
current loop.

    vfGrad :: (R -> R)
           -> ((R,R) -> Position)
           -> (Vec -> (R,R))
           -> FilePath
           -> Int    -- n for n x n
           -> VectorField
           -> IO ()

There are several parameters that the function vfGrad needs in order to
produce a picture like Figure 1. The first argument to vfGrad is a monotonic
function that allows the transition from light arrows to dark arrows to occur
in a nonlinear way; using this we can avoid pictures with a small number of
black arrows on an otherwise white background that would occur with a linear
scaling because the field is very strong at one location. The next two
arguments to vfGrad control the relationship between the 3-dimensional field
and the 2-dimensional picture; they specify which plane we are looking at and
how to convert 3-dimensional vectors into 2-dimensional arrows. The last three
arguments to vfGrad are a file name for the graphics file, an integer
describing how many arrows we want in each direction, and finally the vector
field itself. The definition of this function uses Brent Yorgey’s _diagrams_
package[10].

With types for scalar fields and vector fields, we are ready to explore the
mathematics used in electromagnetic theory, vector calculus.

##  4 Vector Calculus in Haskell

Vector calculus is a large subject. Undergraduate physics majors at Lebanon
Valley College take one semester of vector calculus, but that is not quite
enough for them to get to results like Stokes’ theorem, which we need in
electromagnetic theory. So, part of our electricity and magnetism course is a
review and deeper dive into vector calculus.

I have always included some vector calculus in my electricity and magnetism
course, but in fall 2023, I taught analytical methods for vector calculus in
parallel with numerical methods using Haskell. This section gives a sampling
of what vector calculus can look like in Haskell.

In order to do vector calculus, we are first going to need a data type for
vectors. The Vec type was seen in passing in Section 2.6, and is for
3-dimensional vectors. Details about the Vec type can be found in [7] and [9].

Electromagnetic theory is a geometric subject, and the vector calculus it uses
has a geometric flavor. To see why, consider the following form of the Maxwell
equations.

###  4.1 Maxwell Equations, Integral Form

There are several ways of writing the Maxwell equations. Below we give a
version known as the “integral form”. (There is also a “differential form”
that is probably seen more often, but hides the geometric character of the
theory.) For each of the four equations, I give a statement in words, followed
by an equation in mathematical notation that expresses the same thing.

  1. 1.

Ampere-Maxwell Law: The rate of change of electric flux through a surface is
the magnetic circulation around its boundary minus the current flowing through
the surface.

$\frac{d}{dt}\int_{S}\vec{E}\cdot
d\vec{a}=\frac{1}{\epsilon_{0}\mu_{0}}\int_{\partial S}\vec{B}\cdot
d\vec{\ell}-\frac{1}{\epsilon_{0}}I$

  2. 2.

Faraday’s Law: The rate of change of magnetic flux through a surface is the
opposite of the electric circulation around its boundary.

$\frac{d}{dt}\int_{S}\vec{B}\cdot d\vec{a}=-\int_{\partial S}\vec{E}\cdot
d\vec{\ell}$

  3. 3.

Gauss’s Law: The electric flux through a closed surface is the charge
contained inside.

$\int_{S}\vec{E}\cdot d\vec{a}=\frac{1}{\epsilon_{0}}Q$

  4. 4.

No Magnetic Monopoles Law: The magnetic flux through a closed surface is zero.

$\int_{S}\vec{B}\cdot d\vec{a}=0$

In these equations, $\vec{E}$ is the electric field and $\vec{B}$ is the
magnetic field. Electric and magnetic fields are the fundamental physical
quantities that electromagnetic theory deals with. They are created by charged
particles, and they exert forces on charged particles. Most (perhaps all, I
have no counterexample) physicists regard electric and magnetic fields to be
just as real as the particles they interact with. Particles and fields, then,
have a symbiotic relationship. Conceptually, they are equal partners in the
explanation of electromagnetic phenomena, both dynamic actors on the
electromagnetic stage. Most people, students included, have a much easier time
imagining a particle than imagining a field. Electric and magnetic fields have
type VectorField, which is a type synonym for Position -> Vec.

In the first two equations, $S$ is a surface and $\partial S$ is its boundary
curve. An electric current $I$ flows through surface $S$. In the third and
fourth equations, $S$ is a closed surface, containing some volume, and $Q$ is
the total electric charge in that volume. The constants $\epsilon_{0}$ and
$\mu_{0}$ appear in the equations because they use the International System
(SI) of units. These constants are elided from the English-language
descriptions.

Notice that surfaces and their boundary curves make an explicit appearance in
this form of the Maxwell equations. This is evidence of the geometrical
character of electromagnetic theory.

Beyond their role in the statement of the Maxwell equations, curves and
surfaces are used to describe continuous charge and current distributions. We
saw an example of this with the definition of the ChargeDistribution type in
Section 2.6. The ChargeDistribution data type makes reference to the types
Curve, Surface, and Volume. The simplest of these is the Curve data type,
defined as follows.

    data Curve = Curve { curveFunc          :: R -> Position
                       , startingCurveParam :: R  -- t_a
                       , endingCurveParam   :: R  -- t_b
                       }

A Curve is a map of one parameter into space, along with initial and final
values of that parameter. A Surface is a map of two parameters into space,
along with boundary data.

    data Surface = Surface ((R, R) -> Position) R R (R -> R) (R -> R)

The first argument to Surface is the parameterizing function that maps $(s,t)$
into a Position. The second and third arguments are lower and upper limits of
the parameter $s$. The fourth and fifth arguments are functions of $s$ that
give the lower and upper limits of $t$. Allowing the limits of the second
parameter to depend on the first parameter makes it easier to describe a
surface like a triangle. The boundary of a Surface is a Curve.

A Volume is a map of three parameters into space, along with boundary
information. The boundary of a Volume is a Surface.

Once we have data types for curves, surfaces, and volumes, vector calculus
wants us to perform integrals over them.

###  4.2 Flux and Circulation Integrals

For any vector field $\vec{F}$ and any surface $S$, the _flux_ of the field
through the surface is the integral

$\int_{S}\vec{F}\cdot d\vec{a}.$

The integral is calculated by approximating the surface $S$ by many small
triangles. Each triangle has a vector area whose magnitude is the area of the
triangle and whose direction is perpendicular to the triangle. At each
triangle, we form the inner (dot) product of the vector field $\vec{F}$ at the
location of the triangle and the vector area of the triangle; this inner
product is a number. We add up these numbers for all of the triangles and call
the result the flux of $\vec{F}$ through $S$. If the vector field is the
electric field, we call the flux electric flux. If the vector field is the
magnetic field, we call the flux magnetic flux. If $\vec{F}$ represented the
flow of stuff, then the flux of $\vec{F}$ through $S$ is how much stuff flows
through $S$.

Here is the Haskell code to approximate a flux integral.

    dottedSurfaceIntegral :: SurfaceApprox -> VectorField -> Surface -> R
    dottedSurfaceIntegral approx vF s
        = sum [vF r’ <.> da’ | (r’,da’) <- approx s]

The function dottedSurfaceIntegral needs an approximation method (type
SurfaceApprox) to divide the surface into triangles.

    type SurfaceApprox = Surface -> [(Position, Vec)]

This approximation method, called approx in the code, produces a list of pairs
of positions and vector areas. The left side of the list comprehension shows
that we take the dot product (using the operator <.>) of the vector field vF
evaluated at position r’ with the vector area da’. The numerical results of
the inner products are then added up with sum.

For any vector field $\vec{F}$ and any curve $C$, the _circulation_ of the
field along the curve is the integral

$\int_{C}\vec{F}\cdot d\vec{\ell}.$

The integral is calculated by approximating the surface $C$ by many small
segments. Each segment has a vector length whose magnitude is the length of
the segment and whose direction is along the segment. At each segment, we form
the inner (dot) product of the vector field $\vec{F}$ at the location of the
segment and the vector length of the segment; this inner product is a number.
We add up these numbers for all of the segments and call the result the
circulation of $\vec{F}$ along $C$. If the vector field is the electric field,
we call the circulation electric circulation. If the vector field is the
magnetic field, we call the circulation magnetic circulation.

Here is the Haskell code to approximate a circulation integral.

    dottedLineIntegral :: CurveApprox -> VectorField -> Curve -> R
    dottedLineIntegral approx f c = sum [f r’ <.> dl’ | (r’,dl’) <- approx c]

The function dottedLineIntegral needs an approximation method to divide the
curve into segments. This approximation method, called approx in the code,
produces a list of pairs of positions and vector lengths. The left side of the
list comprehension shows that we take the dot product (using the operator <.>)
of the vector field f evaluated at position r’ with the vector length dl’. The
numerical results of the inner products are then added up with sum.

There are quite a number of these line, surface, and volume integrals that are
used in electromagnetic theory. The table below shows the nine integrals that
EM theory needs from vector calculus.

Haskell function | mathematical notation | requires
---|---|---
scalarLineIntegral | $\int_{C}f\ d\ell$ | scalar field | curve
vectorLineIntegral | $\int_{C}\vec{F}\ d\ell$ | vector field | curve
dottedLineIntegral | $\int_{C}\vec{F}\cdot d\vec{\ell}$ | vector field | curve
crossedLineIntegral | $\int_{C}\vec{F}\times d\vec{\ell}$ | vector field | curve
scalarSurfaceIntegral | $\int_{S}f\ da$ | scalar field | surface
vectorSurfaceIntegral | $\int_{S}\vec{F}\ da$ | vector field | surface
dottedSurfaceIntegral | $\int_{S}\vec{F}\cdot d\vec{a}$ | vector field | surface
scalarVolumeIntegral | $\int_{V}f\ dv$ | scalar field | volume
vectorVolumeIntegral | $\int_{V}\vec{F}\ dv$ | vector field | volume

Vector calculus has its own versions of the fundamental theorem of calculus,
the theorem that makes precise the idea that derivatives and integrals are
inverse operations of each other. The fundamental theorems of vector calculus
are known as the gradient theorem, Stokes’ theorem, and the divergence
theorem. They can be expressed using the integrals in the table above.

###  4.3 A Homework Problem

After students understand the meaning of the integrals in the table above,
they are in a position to study the three fundamental theorems of vector
calculus. Since each theorem claims that an integral of a derivative over some
geometric object is equal to a different integral over the boundary of that
object, checking these fundamental theorems provides opportunities to practice
evaluating the integrals. Even better, the two sides, which are each a bit of
a challenge to set up, must give the same number in the end, which serves as
satisfying evidence that the integrals were done correctly.

Stokes’ theorem claims that for any vector field $\vec{F}$ and any surface
$S$, the following equality holds.

$\int_{S}(\vec{\nabla}\times\vec{F})\cdot d\vec{a}=\int_{\partial
S}\vec{F}\cdot d\vec{\ell}$

The expression $\vec{\nabla}\times\vec{F}$ is called the curl of $\vec{F}$.
The curl is a vector derivative that takes a VectorField as input and produces
a VectorField as output. The physical meaning of the curl is the extent to
which a vector field, if it were a velocity field of water, encourages
rotation of a small object placed in the water. A vector field with zero curl
is called conservative, meaning that a dotted line integral of the field
around any closed curve is zero.

The following homework problem, which was assigned in fall 2023, asks students
to check Stokes’ theorem for a particular vector field and a particular
surface.

> Homework Problem:  Check Stokes’ theorem using the vector field
>
> $\vec{F}(x,y,z)=-z\hat{y}+y\hat{z}$
>
> and the rectangular region with corners $(x,y,z)=(0,0,-4)$, $(0,2,-4)$,
> $(0,2,4)$, $(0,0,4)$, in which the orientation is in the positive $x$
> direction.
>
> Checking the fundamental theorem for curls means (a) finding the curl of
> your vector field, (b) finding the flux integral of the curl over the
> surface, (c) finding the line integral (the circulation) of your vector
> field over the boundary of the surface (which has 4 parts), and (d)
> confirming that the results from (b) and (c) are the same. You may do this
> problem analytically, numerically, or a combination of both.

To solve this homework problem numerically with Haskell, we first encode the
vector field $\vec{F}$ as vF and the rectangle as rect.

    vF :: VectorField
    vF r = let (x,y,z) = cartesianCoordinates r
           in (-z) *^ yHat r ^+^ y *^ zHat r

    rect :: Surface
    rect = Surface (\(y,z) -> cart 0 y z) 0 2 (const (-4)) (const 4)

To calculate the left-hand side of Stokes’ theorem, we need to take the curl
of the vector field. As with a numerical derivative, we need to give a step
size over which to evaluate the curl; we choose $10^{-6}$.

    curlF :: VectorField
    curlF = curl 1e-6 vF

Now we can calculate the left side of Stokes’s theorem.

    leftSide :: R
    leftSide = dottedSurfaceIntegral (surfaceSample 200) curlF rect

We get the result 32.00000000236293. For the right side, we need the boundary
curve of the reactangle.

    boundaryParam :: R -> Position
    boundaryParam t
        | t < 1      = cart 0 (2*t) (-4)
        | t < 2      = cart 0 2 (8*(t-1) - 4)
        | t < 3      = cart 0 (2 - 2 * (t-2)) 4
        | otherwise  = cart 0 0 (4 - 8 * (t-3))

    boundaryOfRect :: Curve
    boundaryOfRect = Curve boundaryParam 0 4

    rightSide :: R
    rightSide = dottedLineIntegral (curveSample 1000) vF boundaryOfRect

The right side gives 32.0, which is reasonably close to the left side, and can
be explained by numerical error.

In this homework problem, I give students the choice about whether to use
analytic (pencil and paper) methods, numerical methods, or a combination of
both. I do not always give students this choice. Sometimes I demand they use
analytic methods; sometimes I demand they use numerical methods. At least one
third of the time, however, I give students the choice, which they seem to
appreciate.

For this homework problem in fall 2023, 5 of the 10 students who submitted it
chose to do it analytically, and 5 chose to do it numerically using Haskell.

As a last example of the power of types and higher order functions in EM
theory, we turn to calculation of the magnetic field by the Biot-Savart law.

##  5 Magnetic Field Produced by a Current Distribution

Earlier, we saw that the central idea of electrostatics is that a charge
distribution produces an electric field. Analogously, the central idea of
magnetostatics is that a current distribution produces a magnetic field. One
can define a type CurrentDistribution, and a function

    bField :: CurrentDistribution -> VectorField

that calculates the magnetic field produced by any current distribution.[9]

In this section, we focus on a particular type of current distribution, the
one we probably have in mind when thinking of electric current, namely the
current flowing through a wire. The equation for calculating the magnetic
field produced by current flowing through a wire is called the Biot-Savart
law, and is shown below.

$\vec{B}(\vec{r})=-\frac{\mu_{0}I}{4\pi}\int_{C}\frac{\vec{r}-\vec{r}^{\prime}}{~{}\left|\vec{r}-\vec{r}^{\prime}\right|^{3}}\times
d\vec{\ell}^{\prime}$ (1)

In this equation, $C$ is a curve of arbitrary shape that describes the
configuration of the wire. The wire could run along a straight line, around in
a circle, in a helix, or any other configuration. The integral asks the
_source point_ $\vec{r}^{\prime}$ to take all possible values along the wire,
the current at each position along the wire being a source of magnetic field.
The _field point_ $\vec{r}$ is the place where we want to find the magnetic
field; it is a constant from the perspective of the integral. The current $I$
is just a number, usually measured in Amperes. The symbol
$d\vec{\ell}^{\prime}$ represents a small section of the wire. The
multiplication $\times$ is the vector cross product.

Here is the Biot-Savart law, translated into Haskell.

    bFieldFromLineCurrent
        :: Current      -- current (in Amps)
        -> Curve
        -> VectorField  -- magnetic field (in Tesla)
    bFieldFromLineCurrent i c r
        = let coeff = -mu0 * i / (4 * pi)  -- SI units
              integrand r’ = d ^/ magnitude d ** 3
                  where d = displacement r’ r
          in coeff *^ crossedLineIntegral (curveSample 1000) integrand c

The type Current is a synonym for R or Double; it is the number of Amperes of
current flowing through the wire. We define a local constant coeff to hold the
numerical value of $-\mu_{0}I/4\pi$ in SI units, and we define a local
function integrand to hold the integrand. We want to define a local variable d
for the displacement from r’ to r, but because r’ exists locally to the
function integrand, the definition for d must occur within the definition for
integrand and cannot be placed parallel to the definitions of coeff and
integrand. We use the crossedLineIntegral to do the integration; this function
appears in the table of vector calculus integrals in Section 4.2.

The type signature of bFieldFromLineCurrent makes clear, in a computer-checked
way, the two inputs required to find the magnetic field: the Current and the
Curve along which the current flows. For a reader of Haskell, the function
bFieldFromLineCurrent is a clearer description of what is going on than
Equation 1 since the latter does not make it terribly clear that the magnetic
field depends only on the curve and the current.

The function bFieldFromLineCurrent was used to produce Figure 1 in Section
2.1.

Next we describe a homework problem which requires the Biot-Savart law. In
Section 2.1, we mentioned that the magnetic field produced by a circular loop
carrying current is exactly solvable along the axis that runs through the
center of the circle (perpendicular to the plane of the circle), but not
exactly solvable anywhere else. This homework problem asks students to use
analytical methods to find the magnetic field on the axis, where a closed-form
algebraic result can be obtained, and then to use numerical methods to extend
their solution to any other point in space.

> Homework Problem:  Find the magnetic field on the $z$-axis produced by a
> circular current loop with radius $R$ lying in the $xy$-plane, centered at
> the origin, carrying current $I$. Do this analytically, then show how to use
> Haskell to find the magnetic field anywhere.

##  6 Conclusion

We’ve shown a number of features that functional programming languages have,
especially types, higher-order functions, and referential transparency, that
make them particularly appropriate for expressing the elegant ideas of
electromagnetic theory. We’ve seen how scalar fields and vector fields can be
encoded as types, and how these types make central ideas of EM theory clear.
Aspects of vector calculus, like line and surface integrals, can be used to
express the Maxwell equations. A Haskell version of the Biot-Savart law is
arguably more understandable than its traditional mathematical form.
Additional examples of electromagnetic theory in Haskell can be found in [9].
My anecdotal evidence that this is an effective way to teach the subject comes
from students this year performing better on the _analytical_ portions of
exams than in prior years. Although the number of students (10) is small, I
take this to mean that the process of programming in Haskell helped students
get a better handle on the basic ideas I’m trying to teach. The author
anticipates that dependently-typed functional languages could have even more
to offer to physics and EM theory. For those who want to express the ideas of
electromagnetic theory in beautiful code, functional languages are hard to
beat.

## References

  * [1]
  * [2] Harold Abelson, Gerald Jay Sussman & Julie Sussman (1996): _Structure and Interpretation of Computer Programs_ , second edition.  MIT Press.
  * [3] David J. Griffiths (2017): _Introduction to Electrodynamics_.  Cambridge University Press, 10.1017/9781108333511.
  * [4] Seymour A. Papert (1993): _Mindstorms: Children, Computers, And Powerful Ideas_ , 2 edition.  Basic Books.
  * [5] Gerald Jay Sussman & Jack Wisdom (2001): _Structure and Interpretation of Classical Mechanics_.  The MIT Press.
  * [6] Gerald Jay Sussman & Jack Wisdom (2013): _Functional Differential Geometry_.  The MIT Press.
  * [7] Scott N. Walck (2014): _Learn Physics by Programming in Haskell_.  In James Caldwell, Philip Hölzenspies & Peter Achten, editors: Proceedings 3rd International Workshop on Trends in Functional Programming in Education, Soesterberg, The Netherlands, 25th May 2014, Electronic Proceedings in Theoretical Computer Science 170, Open Publishing Association, pp. 67–77, 10.4204/EPTCS.170.5.
  * [8] Scott N. Walck (2016): _Learn Quantum Mechanics with Haskell_.  In Johan Jeuring & Jay McCarthy, editors: Proceedings of the 4th and 5th International Workshop on Trends in Functional Programming in Education, Sophia-Antipolis, France and University of Maryland College Park, USA, 2nd June 2015 and 7th June 2016, Electronic Proceedings in Theoretical Computer Science 230, Open Publishing Association, pp. 31–46, 10.4204/EPTCS.230.3.
  * [9] Scott N. Walck (2023): _Learn Physics with Functional Programming: A Hands-on Guide to Exploring Physics with Haskell_.  No Starch Press.
  * [10] Brent Yorgey (2011–2016): _The diagrams package_.  https://hackage.haskell.org/package/diagrams.


# Generalized Time Domain Velocity Vector

Srđan Kitić and Jérôme Daniel
Orange Labs Both authors have equally contributed to the submitted article.

###### Abstract

We introduce and analyze Generalized Time Domain Velocity Vector (GTVV), an
extension of the previously presented acoustic multipath footprint extracted
from the Ambisonic recordings. GTVV is better adapted to adverse acoustic
conditions, and enables efficient parameter estimation of multiple plane wave
components in the recorded multichannel mixture. Experiments on simulated data
confirm the predicted theoretical advantages of these new spatio-temporal
features.

###### Index Terms:

localization, echo, delay, relative transfer function, spherical harmonics

##  I Introduction

Estimating the direction of an active sound source is an important task for a
number of applications, including speech enhancement [1], robot navigation
[2], compression [3] or augmented reality [4]. On the one hand, localizing an
acoustic source is often made difficult by acoustic reflections and
reverberation that interfere with direct wave, which generally causes an
angular bias and uncertainty [5, 6]. On the other hand, the knowledge of
directions and delays of reflected waves, with respect to the direct one,
enables some innovative use-cases such as the inference of room geometry [7],
source separation aided by the acoustic echoes [8], or localization of a
source hidden by a soundproof obstacle [9].

In First Order Ambisonics (FOA), _i.e._ the first order spherical harmonic
(SH) representation of a sound field [10, 11], a common source localization
method is through the use of so-called pseudo-intensity [12], _i.e._ an
estimate of active sound intensity [13]. In favorable acoustic conditions
(free field propagation, relatively low noise), this vector quantity is
orthogonal to the propagating wavefront, and allows for the Direction-of-
Arrival (DoA) estimation at low computational cost [14, 12]. Unfortunately, in
the reverberant environments, the predictions produced by this method are
often biased. Some methods to make the intensity estimate more robust aim at
identifying the time-frequency bins dominated by the direct sound [15], or on
exploiting the higher order spherical harmonic coefficients (_a.k.a._ Higher
Order Ambisonics (HOA)), such as in [16, 17].

While the common wisdom in state-of-the-art was, thus, to suppress the effects
of acoustic reflections, in the earlier work [6], we have shown that it is
possible to improve the DoA estimate and estimate the direction and the delay
of a dominant horizontal acoustic reflection, at the same time. To achieve
this, we use two analysis tools: Frequency Domain Velocity Vector (FDVV), and
Time Domain Velocity Vector (TDVV). The former is a normalized version of
complex sound intensity [13, 14], aggregating both active intensity (in-phase
part), and reactive intensity (quadrature part), while the latter is its
temporal representation. We have demonstrated that TDVV has an interpretable
structure that enables straightforward estimation of the acoustic multipath
parameters: the propagating directions, relative delays and attenuation
factors. However, the TDVV analysis has been theoretically justified only
under a restrictive condition that the magnitude of the direct wave exceeds
the cumulative magnitude of the reflections.

In this article we generalize the TDVV framework in two ways: first, we
include in the representation the HOA coefficients of any order, and second,
the reference component is allowed to be a linear combination of all channels.
Both extensions bring their own benefits: while the addition of HOA improves
the discrimination of two (or more) wavefront directions, the appropriate
choice of the reference component enables the time domain analysis, even under
the adverse acoustic conditions. We remark that, when the reference component
is the omnidirectional (zero-order) Ambisonic channel, the frequency and time
domain formulations are equivalent to Relative Transfer Function (ReTF) [18,
12] (also known as “relative harmonic coefficients” in some recent works [19,
20]), and Relative Impulse Response (ReIR) [21], in spherical harmonics
domain, respectively. To the best knowledge of the authors, the only other
work that mentions non-conventional reference component is [22], which
discusses the ReTF/ReIR in spherical harmonic domain for a denoising
application, under the assumption that the directions of acoustic reflections
are already known.

The article is organised as follows: in Section II we introduce Generalized
Time domain Velocity Vector (GTVV) and discuss its properties. Next, in
Section III a robust estimator of GTVV is presented. The inference of acoustic
parameters from an estimated GTVV is discussed in Section IV. The proof-of-
concept experimental results are given in Section V. The Section VI concludes
the article.

##  II Generalized Time domain Velocity Vector

Consider a single, far field sound source emitting the signal $S(f)$ at
frequency $f$, in an reverberant environment. Our aim is to model the received
signal which includes the direct path propagation, as well as several dominant
acoustic reflections, which are the delayed and attenuated copies of $S(f)$.
Assuming that mode strength compensation [12] has been applied at the
recording microphone, the spherical harmonics expansion coefficient
$B^{(\mathsf{lm})}(f)$, of order $\mathsf{l}$ and degree $\mathsf{m}$, due to
$\mathsf{N}+1$ incoming acoustic waves (“early” reverberation), is modelled as
follows [12, 10]:

$B^{(\mathsf{lm})}(f)=S(f)\sum\limits_{\mathsf{n}=0}^{\mathsf{N}}a_{\mathsf{n}}(f)Y^{(\mathsf{lm})}(\Omega_{\mathsf{n}}).$
(1)

The frequency-dependent factor $a_{\mathsf{n}}(f)$ models the complex
amplitude $c_{\mathsf{n}}(f)$ and phase shift of the $\mathsf{n}$th wavefront,
_i.e._ $a_{\mathsf{n}}(f)=c_{\mathsf{n}}(f)e^{-j2\pi
f\bar{\tau}_{\mathsf{n}}}$, where $\bar{\tau}_{\mathsf{n}}$ is its Time-of-
Arrival (ToA). The scalar $Y^{(\mathsf{lm})}(\Omega_{\mathsf{n}})$ represents
the (order $\mathsf{l}$, degree $\mathsf{m}$) spherical harmonic function
evaluated at the wavefront direction
${\Omega_{\mathsf{n}}:=(\theta_{\mathsf{n}},\varphi_{\mathsf{n}})}$, with
$\theta_{\mathsf{n}}$ and $\varphi_{\mathsf{n}}$ being its azimuth and
elevation coordinates, respectively.

Let $\boldsymbol{\mathbf{{y}}}_{\mathsf{n}}$ designate the vector whose
entries are the SH coefficients (up to some maximal order
$\mathsf{l}=\mathsf{L}$), corresponding to the $\mathsf{n}$th wavefront only:

$\boldsymbol{\mathbf{{y}}}_{\mathsf{n}}:=\boldsymbol{\mathbf{{y}}}(\Omega_{\mathsf{n}})=[Y^{(00)}(\Omega_{\mathsf{n}})\;\boldsymbol{\mathbf{{y}}}^{(1*)}(\Omega_{\mathsf{n}})^{\mathsf{T}}\dots\;\boldsymbol{\mathbf{{y}}}^{(\mathsf{L}*)}(\Omega_{\mathsf{n}})^{\mathsf{T}}]^{\mathsf{T}},$
(2)

where each subvector
$\boldsymbol{\mathbf{{y}}}^{(\mathsf{l}*)}(\Omega_{\mathsf{n}})$ is composed
of $2\mathsf{l}+1$ coefficients of order $\mathsf{l}$, and
$(\cdot)^{\mathsf{T}}$ is the transpose. Then the observation vector
$\boldsymbol{\mathbf{{b}}}(f)$, consisting of stacked expansion coefficients
$\left\\{B^{(\mathsf{lm})}(f)\right\\}_{\mathsf{l}\in[0,\mathsf{L}],\mathsf{m}\in[-l,l]}$
can be represented as

$\boldsymbol{\mathbf{{b}}}(f)=S(f)\sum\limits_{\mathsf{n}=0}^{\mathsf{N}}a_{\mathsf{n}}(f)\boldsymbol{\mathbf{{y}}}_{\mathsf{n}}.$
(3)

We denote by $\boldsymbol{\mathbf{{v}}}(f)$ the Generalized Frequency domain
Velocity Vector (GFVV), defined as follows:

$\boldsymbol{\mathbf{{v}}}(f)=\frac{\boldsymbol{\mathbf{{b}}}(f)}{\boldsymbol{\mathbf{{w}}}(f)^{\mathsf{T}}\boldsymbol{\mathbf{{b}}}(f)}=\frac{\sum\limits_{\mathsf{n}=0}^{\mathsf{N}}a_{n}(f)\boldsymbol{\mathbf{{y}}}_{\mathsf{n}}}{\sum\limits_{\mathsf{n}=0}^{\mathsf{N}}a_{\mathsf{n}}(f)\beta_{\mathsf{n}}(f)},$
(4)

where $\boldsymbol{\mathbf{{w}}}(f)$ is a complex weight vector (_e.g._ , a
beamformer), and
$\beta_{\mathsf{n}}(f)=\boldsymbol{\mathbf{{w}}}(f)^{\mathsf{T}}\boldsymbol{\mathbf{{y}}}_{\mathsf{n}}$.
In the noiseless setting, provided that $S(f)$ is non-negligible, GFVV does
not depend on the content of the source signal. We remark that the frequency
domain velocity vector (FDVV), introduced in [6], is a particular case of
GFVV, for which $\boldsymbol{\mathbf{{b}}}$ is restricted to $4$ FOA channels,
and
$\boldsymbol{\mathbf{{w}}}=[1\;\boldsymbol{\mathbf{{0}}}^{\mathsf{T}}]^{\mathsf{T}}$,
_i.e._ the reference component is the omnidirectional channel $B^{(00)}(f)$
($\boldsymbol{\mathbf{{0}}}$ is the all-zero vector of an appropriate size).

Let ${a_{\mathsf{n}}(f)/a_{\mathsf{p}}(f)=g_{\mathsf{n}}(f)e^{-j2\pi
f(\bar{\tau}_{\mathsf{n}}-\bar{\tau}_{\mathsf{p}})}}$, _i.e._ we declare
$g_{\mathsf{n}}(f)=c_{\mathsf{n}}(f)/c_{\mathsf{p}}(f)$, which is a relative
gain between two wavefronts. Without loss of generality, let $\mathsf{p}=0$
and set $\beta_{0}=1$, which gives

$\boldsymbol{\mathbf{{v}}}(f)=\frac{\boldsymbol{\mathbf{{y}}}_{0}+\sum\limits_{\mathsf{n}=1}^{\mathsf{N}}\frac{\gamma_{\mathsf{n}}}{\beta_{\mathsf{n}}}\boldsymbol{\mathbf{{y}}}_{\mathsf{n}}}{1+\sum\limits_{\mathsf{n}=1}^{\mathsf{N}}\gamma_{\mathsf{n}}},$
(5)

where $\gamma_{\mathsf{n}}=g_{\mathsf{n}}(f)\beta_{\mathsf{n}}(f)e^{-j2\pi
f(\bar{\tau}_{\mathsf{n}}-\bar{\tau}_{0})}$ (thus, $\gamma_{0}=1$).

Assume now that $|\sum_{\mathsf{n}=1}^{\mathsf{N}}\gamma_{\mathsf{n}}(f)|<1$:
a sufficient condition for this inequality to hold would be
${\sum_{\mathsf{n}=1}^{\mathsf{N}}|g_{\mathsf{n}}(f)\beta_{\mathsf{n}}(f)|<1}$.
This is rarely the case for FDVV [6] or ReTF [12, 19, 20], for which the
reference component is omnidirectional ($\beta_{\mathsf{n}}=1$ for all
$\mathsf{n}$). However, a suitable choice of the filter
$\boldsymbol{\mathbf{{w}}}(f)$ could ensure that this condition is satisfied.
For example, if $\boldsymbol{\mathbf{{w}}}$ was a signal-independent
beamformer [12], it could be oriented towards $\Omega_{0}$, such that the
remaining directions are sufficiently attenuated. If this is indeed the case,
the denominator part of the eq. (5) admits a Taylor (specifically, geometric)
series expansion:

$\rho:=\frac{1}{1+\sum\limits_{\mathsf{\mathsf{n}}=1}^{\mathsf{N}}\gamma_{\mathsf{n}}}=\sum\limits_{\mathsf{k}=0}^{\infty}\left(-\sum\limits_{\mathsf{n}=1}^{\mathsf{N}}\gamma_{\mathsf{n}}\right)^{\mathsf{k}}=\sum\limits_{\mathsf{k}=0}^{\infty}\rho_{\mathsf{k}},$
(6)

where

$\rho_{\mathsf{k}}=\sum_{\mathsf{i}_{1}+\mathsf{i}_{2}+\ldots+\mathsf{i}_{N}=k}\frac{\mathsf{k}!}{\mathsf{i}_{1}!\mathsf{i}_{2}!\ldots\mathsf{i}_{\mathsf{N}}!}\prod\limits_{\mathsf{q}=1}^{\mathsf{N}}(\gamma_{\mathsf{n}})^{\mathsf{q}}.$

This expression is cumbersome, however, one may observe that each
$\rho_{\mathsf{k}}$ would have an $\mathsf{k}$th-order term that involves only
the parameters of a single wavefront, while the remainder contains “cross-
terms” between different wavesfronts. These cross-terms are grouped together
under the variable $\eta(f)$:

$\rho=1+\sum\limits_{\mathsf{k}=1}^{\infty}\sum\limits_{\mathsf{q}=1}^{\mathsf{N}}\left(-g_{\mathsf{q}}(f)\beta_{\mathsf{q}}(f)\right)^{\mathsf{k}}e^{-j2\pi
f\mathsf{k}(\bar{\tau}_{\mathsf{q}}-\bar{\tau}_{0})}+\eta(f).$ (7)

Since $|g_{\mathsf{q}}(f)\beta_{\mathsf{q}}(f)|<1$, the higher order
polynomial terms in the expression above have progressively decreasing
magnitudes (the same holds true for the $\eta(f)$ part).

At this point we make one further assumption: we consider the gains
$g_{\mathsf{n}}$ to be frequency-independent. In practice, the frequency
response of the gains depends on the environment: for instance, it may be
dictated by the impedance of the materials that generate acoustic reflections
[13]. Hence, this may not be a fully plausible assumption, but one could argue
that, in many practical settings, the impedance varies smoothly within certain
frequency range. Importantly, considering $g_{\mathsf{n}}$ to be a constant
reduces the number of unknowns in our model. Although it is a parameter we are
in control of, we also choose the frequency-independent spatial filter
($\boldsymbol{\mathbf{{w}}}(f)=\boldsymbol{\mathbf{{w}}}$, thus
$\beta_{\mathsf{n}}(f)=\beta_{\mathsf{n}}$ for all frequencies), for
mathematical convenience, since it simplifies the ensuing time domain
expressions.

Hence, the GFVV is expressed as

$\boldsymbol{\mathbf{{v}}}(f)=\left(\boldsymbol{\mathbf{{y}}}_{0}+\sum\limits_{\mathsf{n}=1}^{\mathsf{N}}g_{\mathsf{n}}\beta_{\mathsf{n}}e^{-j2\pi
f(\bar{\tau}_{\mathsf{n}}-\bar{\tau}_{0})}\boldsymbol{\mathbf{{y}}}_{\mathsf{n}}\right)\\\
\times\left(1+\sum\limits_{\mathsf{k}\geq
1}\sum\limits_{\mathsf{q}=1}^{\mathsf{N}}(-1)^{\mathsf{k}}g_{\mathsf{q}}^{\mathsf{k}}\beta_{\mathsf{q}}^{\mathsf{k}}e^{-j2\pi
f\mathsf{k}(\bar{\tau}_{\mathsf{q}}-\bar{\tau}_{0})}+\eta(f)\right).$ (8)

Applying the inverse Fourier transform to $\boldsymbol{\mathbf{{v}}}(f)$, and
manipulating the resulting expression, produces GTVV:

$\boldsymbol{\mathbf{{v}}}(t)=\delta(t)\boldsymbol{\mathbf{{y}}}_{0}+\\\
\sum\limits_{\mathsf{k}\geq
1}\sum\limits_{\mathsf{n}=1}^{\mathsf{N}}(-g_{\mathsf{n}}\beta_{\mathsf{n}})^{\mathsf{k}}\left(\boldsymbol{\mathbf{{y}}}_{0}-\frac{1}{\beta_{\mathsf{n}}}\boldsymbol{\mathbf{{y}}}_{\mathsf{n}}\right)\delta(t-\mathsf{k}(\bar{\tau}_{\mathsf{n}}-\bar{\tau}_{0}))+\tilde{\eta}(t),$
(9)

where again the terms depending on multiple wavefronts are isolated in the
single variable $\tilde{\eta}(t)$. Note that $\tilde{\eta}(t)$ has discrete
support on the temporal axis, which never includes $t\leq 0$ and, usually,
does not coincide with
$\delta(t-\mathsf{k}(\bar{\tau}_{\mathsf{n}}-\bar{\tau}_{0}))$ either. The
GTVV structure in (9) reveals very useful properties, telling us that a “well-
behaved” GTVV is necessarily causal (_i.e._ ,
${\boldsymbol{\mathbf{{v}}}(t<0)=\boldsymbol{\mathbf{{0}}}}$), and sparse. And
equally important, evaluating GTVV at $t=0$ provides an immediate estimate of
$\boldsymbol{\mathbf{{y}}}_{0}$.

##  III Estimation of GFVV and GTVV

The energy of speech signals varies over time and frequency, hence we cannot
directly use the expression (4) to estimate the GFVV. Particularly, as some
frequency bins have negative Signal-to-Noise Ratio (SNR), such an estimate
would be prone to errors, or even numerical instabilities. Instead, we adapt a
well-known ReTF estimator presented in [18], that exploits the
(non)stationarity properties of speech and noise signals [23]. We chose this
approach due to its simplicity and the fact that it does not require a prior
knowledge of voice activity, but other estimators could be used as well,
_e.g._ based on covariance subtraction or whitening [12].

For reasons that will become obvious, in this section we use a time-frequency
representation of the considered quantities, such as Short-Time Fourier
Transform (STFT). Let $N^{(\mathsf{lm})}(\mathsf{u},f)$ denote the additive
noise of the $(\mathsf{l,m})$th channel of
$\boldsymbol{\mathbf{{b}}}(\mathsf{u},f)$, assumed independent of the source
signal $S(\mathsf{u},f)$, at frame $\mathsf{u}$ and frequency $f$. The goal is
to estimate a ratio of source signal-excited parts of the given Ambisonic
channel and the reference, hence one can write:

$B^{(\mathsf{lm})}(\mathsf{u},f)=v^{(\mathsf{lm})}(f)\sum\limits_{\mathsf{l}^{\prime},\mathsf{m}^{\prime}}w^{(\mathsf{l}^{\prime}\mathsf{m}^{\prime})}B^{(\mathsf{l}^{\prime}\mathsf{m}^{\prime})}(\mathsf{u},f)+U^{(\mathsf{lm})}(\mathsf{u},f),$
(10)

with

$U^{(\mathsf{lm})}(\mathsf{u},f)=N^{(\mathsf{lm})}(\mathsf{u},f)-v^{(\mathsf{lm})}(f)\sum\limits_{\mathsf{l}^{\prime},\mathsf{m}^{\prime}}w^{(\mathsf{l}^{\prime}\mathsf{m}^{\prime})}N^{(\mathsf{l}^{\prime}\mathsf{m}^{\prime})}(\mathsf{u},f).$

The variables $v^{(\mathsf{lm})}(f)$ and $w^{(\mathsf{lm})}$ denote the
$(\mathsf{l,m})$th entries of $\boldsymbol{\mathbf{{v}}}(f)$ and
$\boldsymbol{\mathbf{{w}}}$, respectively. Note that $v^{(\mathsf{lm})}(f)$
does not depend on $\mathsf{u}$: this is valid when neither the source nor the
microphone are mobile during the considered temporal segment.

Figure 1: Magnitudes of H-TDVV (left) and GTVV (right) at different times $t$.

Since $B^{(\mathsf{lm})}(\mathsf{u},f)$ and $U^{(\mathsf{lm})}(\mathsf{u},f)$
are correlated, $v^{(\mathsf{lm})}(f)$ and $U^{(\mathsf{lm})}(\mathsf{u},f)$
should be estimated simultaneously [23]. Multiplying both sides of (10) by
$B^{(\mathsf{lm})}(f)^{*}$ and taking the expectation yields

$\phi_{B^{(\mathsf{lm})},B^{(\mathsf{lm})*}}(\mathsf{u},f)=\\\
v^{(\mathsf{lm})}(f)\sum\limits_{\mathsf{l}^{\prime},\mathsf{m}^{\prime}}w^{(\mathsf{l}^{\prime}\mathsf{m}^{\prime})}\phi_{B^{(\mathsf{l}^{\prime}\mathsf{m}^{\prime})},B^{(\mathsf{lm})*}}(\mathsf{u},f)+\phi_{U^{(\mathsf{lm})},B^{(\mathsf{lm})*}}(\mathsf{u},f),$
(11)

with $\phi_{XY}=\mathbb{E}[XY]$, and $\hat{\phi}_{XY}$ being its estimate.
Assuming that noise statistics varies slower than that of speech, one can
approximate that
$\phi_{U^{(\mathsf{lm})},B^{(\mathsf{lm})*}}(\mathsf{u},f)\approx\phi_{U^{(\mathsf{lm})},B^{(\mathsf{lm})*}}(\mathsf{u}^{\prime},f)$
holds within certain temporal segment
$\mathsf{u}^{\prime}\in[\mathsf{u}-\mathsf{u}_{0},\mathsf{u}]$. This gives
rise to the overdetermined system of equations

$\hat{\boldsymbol{\mathbf{{\phi}}}}_{B^{(\mathsf{lm})},B^{(\mathsf{lm})*}}=\left[\begin{matrix}\hat{\boldsymbol{\mathbf{\Phi}}}_{B^{(\mathsf{lm})}}\boldsymbol{\mathbf{{w}}}&\boldsymbol{\mathbf{{1}}}\end{matrix}\right]\left[\begin{matrix}v^{(\mathsf{lm})}(f)\\\
\phi_{U^{(\mathsf{lm})},B^{(\mathsf{lm})*}}(f)\end{matrix}\right]+\boldsymbol{\mathbf{{\varepsilon}}},$
(12)

where
$\hat{\boldsymbol{\mathbf{{\phi}}}}_{B^{(\mathsf{lm})},B^{(\mathsf{lm})*}}$
and $\boldsymbol{\mathbf{{\varepsilon}}}$ are the vectors of concatenated
estimates of
$\phi_{B^{(\mathsf{lm})},B^{(\mathsf{lm})*}}(\mathsf{u}^{\prime},f)$, and
errors
$\varepsilon(\mathsf{u}^{\prime},f)=\phi_{U^{(\mathsf{lm})},B^{(\mathsf{lm})*}}(f)-\hat{\phi}_{U^{(\mathsf{lm})},B^{(\mathsf{lm})*}}(\mathsf{u}^{\prime},f)$,
respectively, for
${\mathsf{u}^{\prime}\in[\mathsf{u}-\mathsf{u}_{0},\mathsf{u}]}$. Likewise,
$\hat{\boldsymbol{\mathbf{\Phi}}}_{B^{(\mathsf{lm})}}$ is a matrix whose
columns correspond to the estimates of cross-correlation between
$B^{(\mathsf{lm})*}(\mathsf{u}^{\prime},f)$ and
$B^{(\mathsf{l}^{\prime}\mathsf{m}^{\prime})}(\mathsf{u}^{\prime},f)$, for all
$(\mathsf{l}^{\prime},\mathsf{m}^{\prime})$, and the same temporal segment,
while $\boldsymbol{\mathbf{{1}}}$ is the all-one vector. Solving such a system
in the least-squares sense (with respect to
$\boldsymbol{\mathbf{{\varepsilon}}}$), for all $f$ and
$(\mathsf{l},\mathsf{m})$, and applying the inverse Fourier transform, yields
a robust estimate of $\boldsymbol{\mathbf{{v}}}(t)$, $t\in(-T/2,T/2]$, where
$T$ is the support size of the STFT window. Therefore, in practice, a GTVV
estimate is a _matrix_
$\boldsymbol{\mathbf{V}}\in\mathbb{R}^{\mathsf{\mathsf{(L+1)^{2}}\times\mathsf{T}}}$.

Algorithm 1 S-OMP for the inference of GTVV parameters

1:$\mathsf{i}=0,\;\boldsymbol{\mathbf{R}}^{(0)}=\boldsymbol{\mathbf{V}},\;\hat{\boldsymbol{\mathbf{\Lambda}}}=\\{\emptyset\\}$

2:repeat

3: $\mathsf{i}\leftarrow\mathsf{i}+1$

4:
$s^{(\mathsf{i})}=\operatorname*{argmax\,}_{s}\|{\boldsymbol{\mathbf{R}}^{(\mathsf{i}-1)}}^{\mathsf{T}}\boldsymbol{\mathbf{Y}}_{:,s}\|_{\infty},\;\;\boldsymbol{\mathbf{Y}}_{:,s}=\boldsymbol{\mathbf{{y}}}(\Omega_{s}^{\prime})$

5:
$\hat{\boldsymbol{\mathbf{\Lambda}}}\leftarrow\hat{\boldsymbol{\mathbf{\Lambda}}}\cup\Omega_{s^{(\mathsf{i})}}^{\prime}$

6:
$\tau^{(\mathsf{i})}:=\bar{\tau}^{(\mathsf{i})}-\bar{\tau}_{0}=\frac{1}{f_{s}}\operatorname*{argmax\,}_{\mathsf{q}}|\boldsymbol{\mathbf{{y}}}(\Omega_{s^{(\mathsf{i})}}^{\prime})^{\mathsf{T}}\boldsymbol{\mathbf{R}}_{:,\mathsf{q}}^{(\mathsf{i}-1)}|$

7:
$\boldsymbol{\mathbf{Z}}^{(\mathsf{i})}=\operatorname*{argmax\,}_{\boldsymbol{\mathbf{Z}}}\|\boldsymbol{\mathbf{Y}}_{\hat{\boldsymbol{\mathbf{\Lambda}}}}\boldsymbol{\mathbf{Z}}-\boldsymbol{\mathbf{V}}\|_{\text{F}}^{2}$

8:
$\boldsymbol{\mathbf{R}}^{(\mathsf{i})}=\boldsymbol{\mathbf{Y}}_{\hat{\boldsymbol{\mathbf{\Lambda}}}}\boldsymbol{\mathbf{Z}}^{(\mathsf{i})}-\boldsymbol{\mathbf{V}}$

9:until stopping criterion

10:return
$\hat{\boldsymbol{\mathbf{\Lambda}}},\\{\bar{\tau}^{(\mathsf{i})}\\}$

Examples of the estimated GTVVs for
$\boldsymbol{\mathbf{{w}}}=[1\;\boldsymbol{\mathbf{{0}}}^{\mathsf{T}}]^{\mathsf{T}}$
(termed here “H-TDVV” as an HOA extension of the TDVV model [6]), and
${\boldsymbol{\mathbf{{w}}}=\boldsymbol{\mathbf{{y}}}(\Omega_{0})}$ (the
maximum directivity beamformer [12], oriented toward DoA), are given in Fig.
1.

HOA order | 1 | 2 | 3 | 4
---|---|---|---|---
Low reverberation
H-TDVV | $12.8^{\circ}$ | $0.17$ |  $2\cdot 10^{-3}$ s | $15.0^{\circ}$ | $0.47$ |  $5.1\cdot 10^{-4}$ s | $10.9^{\circ}$ | $0.52$ |  $10^{-4}$ s | $6.1^{\circ}$ | $0.87$ |  $\mathbf{6.1\cdot 10^{-5}}$ s
GTVV | $\mathbf{12.7}^{\circ}$ | $\mathbf{0.20}$ |  $\mathbf{10^{-4}}$ s | $\mathbf{9.8}^{\circ}$ | $\mathbf{1.01}$ |  $\mathbf{9.9\cdot 10^{-5}}$ s | $\mathbf{8.3}^{\circ}$ | $\mathbf{1.70}$ |  $\mathbf{6.3\cdot 10^{-5}}$ s | $\mathbf{5.8}^{\circ}$ | $\mathbf{3.12}$ |  $6.7\cdot 10^{-5}$ s
High reverberation
H-TDVV | $14.6^{\circ}$ | $\mathbf{0.06}$ |  $1.7\cdot 10^{-2}$ s | $14.7^{\circ}$ | $0.18$ |  $3.8\cdot 10^{-3}$ s | $11.9^{\circ}$ | $0.22$ |  $1.9\cdot 10^{-4}$ s | $8.0^{\circ}$ | $0.12$ |  $\mathbf{2.6\cdot 10^{-5}}$ s
GTVV | $\mathbf{13.9}^{\circ}$ | $\mathbf{0.06}$ |  $\mathbf{4.7\cdot 10^{-4}}$ s | $\mathbf{12.6}^{\circ}$ | $\mathbf{0.77}$ |  $\mathbf{3.3\cdot 10^{-4}}$ s | $\mathbf{10.9}^{\circ}$ | $\mathbf{0.77}$ |  $\mathbf{10^{-4}}$ s | $\mathbf{6.1}^{\circ}$ | $\mathbf{2.21}$ |  $1.1\cdot 10^{-4}$ s
TABLE I: Reflection estimation performance (from left to right of each cell:
angular error, number of detections and the delay error).

##  IV GTVV parameter inference via S-OMP

Once the GTVV has been estimated, one can exploit its structure (9) to obtain
the acoustic parameters of individual wavefronts. For that purpose, we adopt
the Simultaneous Orthogonal Matching Pursuit (S-OMP) algorithm [24], a greedy
method for estimating the sparse support shared by a group of observation
vectors. The aim of using such a simple algorithm is to emphasize the
performance gains due to the presented GTVV model. Its pseudocode is given in
Algorithm 1, which we describe and motivate in the following.

Let
$\boldsymbol{\mathbf{Y}}\in\mathbb{R}^{\mathsf{\mathsf{(L+1)^{2}}\times\mathsf{Y}}}$
be a dictionary whose columns (“atoms”) are SH vectors111We use (without loss
of generality) real-valued SH functions.
$\\{\boldsymbol{\mathbf{{y}}}(\Omega_{\mathsf{j}}^{\prime})\\}_{\mathsf{j}\in[0,\mathsf{Y}-1]}$,
obtained by sampling $\mathsf{Y}$ directions
${\boldsymbol{\mathbf{\Omega}}^{\prime}=\\{\Omega_{j}^{\prime}\\}}_{\mathsf{j}\in[0,\mathsf{Y}-1]}$
quasi-uniformly on the unit sphere. The model (9) suggests that the GTVV
matrix $\boldsymbol{\mathbf{V}}$ could be approximated by a linear combination
of few atoms from $\boldsymbol{\mathbf{Y}}$, _i.e._ ,
$\boldsymbol{\mathbf{V}}\approx\boldsymbol{\mathbf{Y}}_{\boldsymbol{\mathbf{\Lambda}}}\boldsymbol{\mathbf{Z}}$,
where $\boldsymbol{\mathbf{Y}}_{\boldsymbol{\mathbf{\Lambda}}}$ is the set of
columns of $\boldsymbol{\mathbf{Y}}$ corresponding to the support
$\boldsymbol{\mathbf{\Lambda}}$, and
$\boldsymbol{\mathbf{Z}}\in\mathbb{R}^{\mathsf{\mathsf{\\#\boldsymbol{\mathbf{\Lambda}}}\times\mathsf{T}}}$
represents the contributions of different wavefronts. We assume that the
number $\mathsf{N}$ of strong reflections in (1) is small compared to the
number of dictionary atoms $\mathsf{Y}$, _i.e._ a sparse (synthesis) model.

The orthogonalization step 7 of Algorithm 1 projects the observed
$\boldsymbol{\mathbf{V}}$ onto the subspace spanned by
$\boldsymbol{\mathbf{Y}}_{\boldsymbol{\mathbf{\Lambda}}}$, where
$\hat{\boldsymbol{\mathbf{\Lambda}}}$ is the estimated support at the current
iteration. The new residual $\boldsymbol{\mathbf{R}}^{(\mathsf{i})}$ now
contains only the contributions of wavefronts from directions
$\boldsymbol{\mathbf{\Omega}}^{\prime}\setminus\hat{\boldsymbol{\mathbf{\Lambda}}}$.
In a favorable case, the selection rule in the step 4, based on the
$\ell_{\infty}$ norm, would choose an atom
$\boldsymbol{\mathbf{{y}}}(\Omega_{s^{(\mathsf{i})}}^{\prime})$ the most
aligned with the wavefront $\boldsymbol{\mathbf{{y}}}(\Omega_{\mathsf{n}})$,
bearing the largest magnitude in this residual mixture. The pair
$\Omega_{s^{(\mathsf{i})}}^{\prime}=(\theta_{s^{(\mathsf{i})}}^{\prime},\varphi_{s^{(\mathsf{i})}}^{\prime})$
is then added to the set of estimated directions. However, should a scalar
product with some of the “cross-terms” $\tilde{\eta}$ in (9) become larger
than this value, the selected atom will be erroneous. This becomes less likely
for higher Ambisonic orders, due to the completeness property of spherical
harmonics [12].

According to (9), the GTVV magnitudes associated with a SH vector
$\boldsymbol{\mathbf{{y}}}(\Omega_{\mathsf{n}>0})$ are strictly decaying (due
to ${|g_{n}\beta_{n}|<1}$). Hence, the largest correlation between
$\boldsymbol{\mathbf{{y}}}(\Omega_{\mathsf{n}})$ and the current residual
matrix $\boldsymbol{\mathbf{R}}^{(\mathsf{i}-1)}$ should be observed for a
column corresponding to
$t\approx\tau_{\mathsf{n}}:=\bar{\tau}_{\mathsf{n}}-\bar{\tau}_{0}$, which
allows for the (relative) delay estimation in the step 6 (with $f_{s}$ being
the sampling rate). Finally, under the assumption that the DoA component is
the dominant wavefront, the dictionary atom chosen in the first iteration
should coincide with $\Omega_{0}$ at $\bar{\tau}=0$.

##  V Computer experiments

HOA order | 1 | 2 | 3 | 4
---|---|---|---|---
Low / High reverberation
TRAMP |  $18.8^{\circ}$/$29.7^{\circ}$ |  $11.5^{\circ}$/$14.9^{\circ}$ |  $10.9^{\circ}$/$14.2^{\circ}$ |  $8.1^{\circ}$/$10.2^{\circ}$
H-TDVV |  $14.6^{\circ}$/$25.9^{\circ}$ |  $4.5^{\circ}$/$6.2^{\circ}$ |  $7.0^{\circ}$/$10.5^{\circ}$ |  $3.0^{\circ}$/$5.7^{\circ}$
GTVV |  $\mathbf{7.5}^{\circ}$/$\mathbf{18.9}^{\circ}$ |  $\mathbf{3.5}^{\circ}$/$\mathbf{5.3}^{\circ}$ |  $\mathbf{3.7}^{\circ}$/$\mathbf{7.4}^{\circ}$ |  $\mathbf{2.6}^{\circ}$/$\mathbf{4.8}^{\circ}$
TABLE II: DoA estimation error (best results are in bold).

To evaluate the presented GTVV model, we use the shoebox acoustic simulator
software MCRoomSim [25] to generate Room Impulse Responses (RIRs). This
provides us with the complete ground truth data, containing not only the DoA
information, but also the delays and directions of reflected components. The
shoebox room dimensions are $5\times 4\times 2.8$ m3, and the reverberation
time is about $0.16$ s and $0.44$ s, for the low- and high-reverberant
conditions, respectively. The received signals are obtained by convolving five
RIRs, for each reverberation time condition, with a dry speech source at the
sampling rate $f_{s}=16$ kHz. In order to simulate the diffuse interference,
white Gaussian noise (in the Ambisonic encoding format) is added to each
microphone channel, such that SNR$=20$ dB. The STFT frames (multiplied by the
Hamming window) are $0.064$ long, and overlap by $75\%$.

The baselines are the HOA version of the TRAMP algorithm [26] (an efficient
variant of the Ambisonic SRP-PHAT beamformer), and S-OMP applied to the
“omnidirectional” H-TDVV. Moreover, we use the latter’s DoA estimate
$\hat{\Omega}_{0}$ to steer the beamformer
${\boldsymbol{\mathbf{{w}}}=\boldsymbol{\mathbf{{y}}}(\hat{\Omega}_{0})}$ when
computing GTVV. For simplicity, we consider only the DoA and first order
reflections, despite the fact that S-OMP is oblivious to the reflection order.
Thus, we choose the number of iterations (_i.e._ , the number of estimated
directions) as a simple stopping criterion. Clearly, S-OMP cannot estimate
more directions than there are Ambisonic channels, hence we limit their number
to $4$ (in the case of FOA), and $7$ (otherwise). All methods use the same
dictionary $\boldsymbol{\mathbf{Y}}$, obtained by discretizing
$\mathsf{Y}=770$ directions on the spherical Lebedev grid [27].

The DoA estimation performance, in terms of average angular error per frame,
is presented in Table II. Likewise, the errors in estimating reflections are
provided in Table I. Since the algorithm does not discriminate first order
reflections, the results are restricted to those deviating no more than
$20^{\circ}$ from the ground truth. Thus, we additionally provide the average
detection results, as well as the delay estimation errors.

The presented results indicate that all methods, generally, tend to perform
better at higher Ambisonic orders, and worse when reverberation is significant
(this is especially the case for FOA). Due to spatial limitations we are
unable to show the results of the experiments with varying SNR levels. We
report, though, that lowering SNR mainly impacts the detection rate, while the
angular accuracy remains relatively stable. Overall, the results confirm the
predicted theoretical advantages of GTVV, as it largely outperforms the other
two methods in all tested scenarios.

##  VI Conclusion

We have presented a Generalized Time Domain Velocity Vector, an extension of
the TDVV representation introduced in the earlier work. GTVV is better
conditioned than TDVV, and more potent at estimating individual wavefronts of
the recorded Ambisonic mixture, as confirmed by computer experiments using the
S-OMP algorithm. The future investigations will focus on more elaborate
inference techniques, relaxation of some working assumptions, and the
quantification of uncertainty of estimated acoustic parameters.

## References

  * [1] A. Xenaki, J. Bünsow Boldt, and M. Græsbøll Christensen, “Sound source localization and speech enhancement with sparse bayesian learning beamforming,” _The Journal of the Acoustical Society of America_ , vol. 143, no. 6, pp. 3912–3921, 2018.
  * [2] C. Evers and P. A. Naylor, “Acoustic SLAM,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 26, no. 9, pp. 1484–1498, 2018\.
  * [3] V. Pulkki, “Spatial sound reproduction with directional audio coding,” _Journal of the Audio Engineering Society_ , vol. 55, no. 6, pp. 503–516, 2007.
  * [4] F. Ribeiro, D. Florencio, P. A. Chou, and Z. Zhang, “Auditory augmented reality: Object sonification for the visually impaired,” in _2012 IEEE 14th international workshop on multimedia signal processing (MMSP)_. IEEE, 2012, pp. 319–324.
  * [5] C. Blandin, A. Ozerov, and E. Vincent, “Multi-source TDOA estimation in reverberant audio using angular spectra and clustering,” _Signal Processing_ , vol. 92, no. 8, pp. 1950–1960, 2012.
  * [6] J. Daniel and S. Kitić, “Time domain velocity vector for retracing the multipath propagation,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 421–425.
  * [7] I. Dokmanić, R. Parhizkar, A. Walther, Y. M. Lu, and M. Vetterli, “Acoustic echoes reveal room shape,” _Proceedings of the National Academy of Sciences_ , vol. 110, no. 30, pp. 12 186–12 191, 2013.
  * [8] R. Scheibler, D. Di Carlo, A. Deleforge, and I. Dokmanic, “Separake: Source separation with a little help from echoes,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 6897–6901.
  * [9] S. Kitić, N. Bertin, and R. Gribonval, “Hearing behind walls: localizing sources in the room next door with cosparsity,” in _2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2014, pp. 3087–3091.
  * [10] F. Zotter and M. Frank, _Ambisonics: A practical 3D audio theory for recording, studio production, sound reinforcement, and virtual reality_. Springer Nature, 2019.
  * [11] E. Fernandez-Grande and A. Xenaki, “Compressive sensing with a spherical microphone array,” _The Journal of the Acoustical Society of America_ , vol. 139, no. 2, pp. EL45–EL49, 2016.
  * [12] D. P. Jarrett, E. A. Habets, and P. A. Naylor, _Theory and applications of spherical microphone array processing_. Springer, 2017, vol. 9.
  * [13] T. D. Rossing and T. D. Rossing, _Springer handbook of acoustics_. Springer, 2007, vol. 1.
  * [14] J. Merimaa, “Analysis, synthesis, and perception of spatial sound: binaural localization modeling and multichannel loudspeaker reproduction,” Ph.D. dissertation, 2006.
  * [15] A. Moore, C. Evers, P. A. Naylor, D. L. Alon, and B. Rafaely, “Direction of arrival estimation using pseudo-intensity vectors with direct-path dominance test,” in _2015 23rd European Signal Processing Conference (EUSIPCO)_. IEEE, 2015, pp. 2296–2300.
  * [16] A. H. Moore, C. Evers, and P. A. Naylor, “Direction of arrival estimation in the spherical harmonic domain using subspace pseudointensity vectors,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 25, no. 1, pp. 178–192, 2016.
  * [17] S. Hafezi, A. H. Moore, and P. A. Naylor, “Augmented intensity vectors for direction of arrival estimation in the spherical harmonic domain,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 25, no. 10, pp. 1956–1968, 2017.
  * [18] S. Gannot, D. Burshtein, and E. Weinstein, “Signal enhancement using beamforming and nonstationarity with applications to speech,” _IEEE Transactions on Signal Processing_ , vol. 49, no. 8, pp. 1614–1626, 2001.
  * [19] Y. Hu, P. N. Samarasinghe, T. D. Abhayapala, and S. Gannot, “Unsupervised multiple source localization using relative harmonic coefficients,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 571–575.
  * [20] Y. Hu, P. N. Samarasinghe, S. Gannot, and T. D. Abhayapala, “Semi-supervised multiple source localization using relative harmonic coefficients under noisy and reverberant environments,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 28, pp. 3108–3123, 2020.
  * [21] R. Talmon, I. Cohen, and S. Gannot, “Relative transfer function identification using convolutive transfer function approximation,” _IEEE Transactions on audio, speech, and language processing_ , vol. 17, no. 4, pp. 546–555, 2009\.
  * [22] Y. Biderman, B. Rafaely, S. Gannot, and S. Doclo, “Efficient relative transfer function estimation framework in the spherical harmonics domain,” in _2016 24th European Signal Processing Conference (EUSIPCO)_. IEEE, 2016, pp. 1658–1662.
  * [23] O. Shalvi and E. Weinstein, “System identification using nonstationary signals,” _IEEE transactions on signal processing_ , vol. 44, no. 8, pp. 2055–2063, 1996.
  * [24] J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Algorithms for simultaneous sparse approximation. Part I: Greedy pursuit,” _Signal processing_ , vol. 86, no. 3, pp. 572–588, 2006.
  * [25] A. Wabnitz, N. Epain, C. Jin, and A. Van Schaik, “Room acoustics simulation for multichannel microphone arrays,” in _Proceedings of the International Symposium on Room Acoustics_. Citeseer, 2010, pp. 1–6.
  * [26] S. Kitić and A. Guérin, “TRAMP: Tracking by a Real-time AMbisonic-based Particle filter,” in _IEEE-AASP Challenge on Acoustic Source Localization and Tracking-LOCATA_ , 2018.
  * [27] V. I. Lebedev and D. Laikov, “A quadrature formula for the sphere of the 131st algebraic order of accuracy,” in _Doklady Mathematics_ , vol. 59, no. 3. Pleiades Publishing, Ltd., 1999, pp. 477–481.