About this Journal Submit a Manuscript Table of Contents 
Mathematical Problems in Engineering
Volume 2012 (2012), Article ID 478295, 27 pages
doi:10.1155/2012/478295
Research Article

Multidimensional Wave Field Signal Theory: Transfer Function Relationships

Natalie Baddour

Department of Mechanical Engineering, University of Ottawa, Ottawa, ON, K1N 6N5, Canada

Received 29 August 2011; Accepted 20 September 2011

Academic Editor: Carlo Cattani

Copyright © 2012 Natalie Baddour. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The transmission of information by propagating or diffusive waves is common to many fields of engineering and physics. Such physical phenomena are governed by a Helmholtz (real wavenumber) or pseudo-Helmholtz (complex wavenumber) equation. Since these equations are linear, it would be useful to be able to use tools from signal theory in solving related problems. The aim of this paper is to derive multidimensional input/output transfer function relationships in the spatial domain for these equations in order to permit such a signal theoretic approach to problem solving. This paper presents such transfer function relationships for the spatial (not Fourier) domain within appropriate coordinate systems. It is shown that the relationships assume particularly simple and computationally useful forms once the appropriate curvilinear version of a multidimensional spatial Fourier transform is used. These results are shown for both real and complex wavenumbers. Fourier inversion of these formulas would have applications for tomographic problems in various modalities. In the case of real wavenumbers, these inversion formulas are presented in closed form, whereby an input can be calculated from a given or measured wavefield.

1. Introduction

The transmission of information over space and time is often governed by the theory of waves. Many important physical phenomena are described by a Helmholtz equation, for example, in the fields of electromagnetism, acoustics, and optics [1–4]. Other physical phenomena are well described by the use of damped or diffusion waves, such as the propagation of heat or photonic waves. For example, ultrasound has long been used for medical imaging [3] while the use of optical radiation via diffuse photon density waves for imaging inhomogeneities in turbid media is a newer development [5–8]. Similarly, photothermal tomographic imaging methods have also been used for nondestructive evaluation [9–12], and more recently for biomedical imaging [13, 14]. Furthermore, a combination of ultrasonic and optical techniques known as photoacoustics has also gained prominence in recent years [4, 15].

Modelled as a linear process, the propagation and scattering of waves obey the principles of superposition and homogeneity. One common method of studying linear processes is to view them as linear systems with an input and output. The input/output relationships of the system can then be characterized via the system transfer function, which implies interpreting the problem as a signal theory problem and thus enabling the application of powerful system theoretic results and concepts. This type of approach applied to propagating or diffusion wave problems has already yielded promising results [16–18] and the goal of this paper is to expand the scope of this approach.

The seeds of such an approach to the solution of tomographic-type problems can be seen in [1, 8, 10, 19, 20]; however, in those papers there is no clear reference to any such system theoretic concepts. As such, a unifying framework for a signal theoretic approach to wavefields does not exist. The goal of this work is to propose a clear and cohesive signal theoretic framework for working with acoustic, thermal, photonic, or other wavefield-related problems. This would be particularly powerful, enabling the application of a vast array of results from a signal processing point of view and in particular, recent novel results from algebraic signal processing [21–27] should find ground-breaking applications.

With these ideas in mind, the goal of this exposition is to derive multidimensional spatial transfer function relationships for the Helmholtz and pseudo-Helmholtz equation so that the tools of linear system theory can be used to solve related problems.

We consider the input quantity to be the time- and space-dependent inhomogeneity term (the “forcing” term in the relevant partial differential equation). The output quantity is then taken as the resulting field present at any point in space. The goal is to develop transfer function relationships in spatial multidimensions so that general results of the theory of signals and systems may be used. It is known that the relationship between the input (inhomogeneity term) and the output (resulting wavefield) is one of convolution with the impulse response (Green's function in classical terms). In the Fourier domain, the equivalent relationship is a transfer function relationship, that is to say one of multiplication. In this paper, key transfer function relationships are shown to hold for the spatial (not Fourier) domain with the use of appropriate coordinate systems. That is to say, the relationship between input and output is shown to be one of multiplication (not convolution) even in the spatial domain, hence yielding the signal theoretic construction that is the aim of the paper.

Rather than focusing on a specific application area (acoustics, medical imaging, seismic, etc.), the end goal of this paper is a unifying framework. The backbone of this framework is Fourier transforms in multidimensions which are also transcribed into curvilinear coordinate systems. The representations in this paper are embedded in a time-dependent, spatially one- to three-dimensional description.

The outline of this paper is as follows. Sections 2 to 7 establish the background preliminaries, definitions, and sign conventions that are important to the developments in the rest of the paper. The definition of multidimensional Fourier transforms is given in Section 2. Section 3 presents the mathematical forms of the types of signal to which this development applies, namely, any signal governed by a Helmholtz or pseudo-Helmholtz equation. Section 4 introduces the Green's and transfer functions for the (generalized) Helmholtz equation and points out that the spatial transfer function assumes a particularly simple form in curvilinear coordinates. Section 5 explains the sign conventions that will be employed in this paper to ensure consistency and simplicity of results. Sections 6 and 7 develop multidimensional Fourier transforms in curvilinear coordinates. Sections 8, 9, and 10 present the relevant transfer function relationships in curvilinear coordinates in 3D, 2D, and 1D, respectively, and form the core results of this paper. Section 11 summarizes the results of the previous three sections in a tabular format that is easy to refer to. Sections 12 and 13 present some applications while Section 14 concludes the paper.

2. Fourier Transforms in 𝑛 -Dimensions

The theory of Fourier transforms can be extended to ℝ 𝑛 in a way that is completely analogous to the treatment of one-dimensional Fourier transforms. Let elements of ℝ 𝑛 be denoted by ⃗ 𝑥 = ( 𝑥 1 , 𝑥 2 , … , 𝑥 𝑛 ) or more generally ⃗ 𝑟 , and elements of the corresponding Fourier dual space Ω 𝑛 be denoted by ⃗ 𝜔 = ( 𝜔 1 , 𝜔 2 , … , 𝜔 𝑛 ) . Under the suitability of integration of 𝑓 , then the Fourier transform in ℝ 𝑛 is defined for ⃗ 𝜔 in Ω 𝑛 as 𝐹 = ⃗ 𝜔 ℝ 𝑛 𝑓 𝑒 ⃗ 𝑥 − 𝑗 ⃗ 𝜔 ⋅ ⃗ 𝑥 𝑑 ⃗ 𝑥 . ( 2 . 1 ) Under suitable conditions, the function 𝑓 can be recovered from the inverse transform through 𝑓 = 1 ⃗ 𝑥 ( 2 𝜋 ) 𝑛 Ω 𝑛 𝐹 𝑒 ⃗ 𝜔 𝑗 ⃗ 𝜔 ⋅ ⃗ 𝑥 𝑑 ⃗ 𝜔 . ( 2 . 2 ) Various other conventions are possible regarding the location of the positive and negative signs and also the factors of 2 𝜋 in (2.1) and (2.2).

3. The Helmholtz and Pseudo-Helmholtz Equation

All wave fields governed by the wave equation (such as acoustic or electromagnetic waves) lead to the Helmholtz equation once a Fourier transform is used to transform the time domain to the frequency domain: ∇ 2 𝑢 ⃗ 𝑟 , 𝜔 + 𝑘 2 𝑠 𝑢 ⃗ 𝑟 , 𝜔 = − 𝑠 ⃗ 𝑟 , 𝜔 , ( 3 . 1 ) where 𝑘 2 𝑠 = 𝜔 2 / 𝑐 2 𝑠 . Here, 𝑠 ( ⃗ 𝑟 , 𝜔 ) is the temporal Fourier transform of the inhomogeneous time- and space-dependent source term for the wave equation. From a signal theory perspective, this is considered to be the input to the system. The variable 𝑢 ( ⃗ 𝑟 , 𝜔 ) represents a physical variable that is governed by the wave equation, for example, acoustic pressure, and is considered to be the output from a signal theory point of view. Both 𝑠 ( ⃗ 𝑟 , 𝜔 ) and 𝑢 ( ⃗ 𝑟 , 𝜔 ) are functions of position, ⃗ 𝑟 , and (temporal) frequency, 𝜔 . The variable 𝑐 𝑠 represents the speed of the wave, which for an acoustic wave would be the speed of sound and for an electromagnetic wave would be the speed of light. For wavefields governed by the wave equation, 𝑘 2 𝑠 is a real (and positive) quantity.

Other physical phenomena lead to a “pseudo-” Helmholtz equation upon Fourier transformation of the time variable to a (temporal) frequency variable. For example, the equation for a diffuse photon density wave (DPDW) which describes the photon density 𝑢 ( ⃗ 𝑟 , 𝑡 ) in a solid due to incident energy intensity 𝑠 ( ⃗ 𝑟 , 𝑡 ) (optical source function) is given in the time domain by [2] D ∇ 2 𝑢 − 1 ⃗ 𝑟 , 𝑡 𝑐 𝜕 𝑢 𝜕 𝑡 ⃗ 𝑟 , 𝑡 − 𝜇 𝑎 𝑢 ⃗ 𝑟 , 𝑡 = − 𝑠 ⃗ 𝑟 , 𝑡 . ( 3 . 2 )

In the above equation, 𝜇 𝑎 is the optical absorption coefficient (m−1), 𝑐 is the speed of light in the turbid medium (m/s), and D is the optical diffusion coefficient (m). Taking the temporal Fourier transform (FT) of (3.2) to transform time to frequency gives the following: ∇ 2 𝑢 ⃗ 𝑟 , 𝜔 + 𝑘 2 𝑝 𝑢 𝑠 ⃗ 𝑟 , 𝜔 = − ⃗ 𝑟 , 𝜔 D , ( 3 . 3 ) where 𝑘 2 𝑝 = − ( 𝜇 𝑎 / D ) − ( 𝑖 𝜔 / 𝑐 D ) is the complex wavenumber. Equation (3.3) is a pseudo-Helmholtz equation. The term “pseudo” is used because although (3.3) has the form of a Helmholtz equation similar to (3.1), the wavenumber 𝑘 2 𝑝 is a complex variable, with the imaginary part indicating a decaying or damped wave. In this development the terms “Helmholtz” and “pseudo-Helmholtz” will generally be used interchangeably unless otherwise noted. Similarly, the standard heat equation is given by ∇ 2 𝑢 − 1 ⃗ 𝑟 , 𝑡 𝛼 𝜕 𝑢 𝜕 𝑡 ⃗ 𝑟 , 𝑡 = − 𝑠 ⃗ 𝑟 , 𝑡 , ( 3 . 4 ) where 𝑠 ( ⃗ 𝑟 , 𝑡 ) is the time- and space-dependent heat source, 𝑢 ( ⃗ 𝑟 , 𝑡 ) describes the temperature in the material as a function of time and space, and 𝛼 is the thermal diffusivity of a material. As done previously, we take the temporal FT of (3.4) to obtain ∇ 2 𝑢 ⃗ 𝑟 , 𝜔 + 𝑘 2 𝑡 𝑢 ⃗ 𝑟 , 𝜔 = − 𝑠 ⃗ 𝑟 , 𝜔 , ( 3 . 5 ) where 𝑘 2 𝑡 = − 𝑖 𝜔 / 𝛼 is the complex wavenumber. Once again (3.5) has the form of a Helmholtz equation, although it is a pseudo-Helmholtz equation to be exact since the wavenumber is complex and indicates a damped wave.

Thus, we see that the general Helmholtz form of (3.1) can be used to describe several different physical phenomena, from the propagation of light or acoustic waves to the heavily damped nature of photonic or thermal waves. The exact form of the wavenumber in each case is the best indicator as to the propagation characteristics of a wave with a real 𝑘 2 indicating a propagating wave and a complex 𝑘 2 indicating a damped wave.

4. Green's Function and Transfer Function for the Helmholtz Equation

Taking the full spatial Fourier transform of the Helmholtz (or pseudo-Helmholtz) equation and rearranging yields 𝑈 𝜔 ⃗ 𝜔 , 𝜔 = 𝑈 𝑥 , 𝜔 𝑦 , 𝜔 𝑧 = 1 , 𝜔 𝜔 2 𝑥 + 𝜔 2 𝑦 + 𝜔 2 𝑧 − 𝑘 2 𝑆 ⃗ 𝜔 , 𝜔 , ( 4 . 1 ) where we have used the shorthand notation ( ⃗ 𝜔 ) = ( 𝜔 𝑥 , 𝜔 𝑦 , 𝜔 𝑧 ) to denote a point in 3D spatial Fourier frequency space. The wavenumber may be real or complex. A capital letter is used to denote the full 3D spatial and temporal Fourier transform of a function, although it should be clear from the arguments which is being indicated. In cases where it may not be clear, a tilde ( ∼ ) will be used to denote the function in spatial Fourier space. For shorthand notation, let 𝜔 2 𝑘 = 𝜔 2 𝑥 + 𝜔 2 𝑦 + 𝜔 2 𝑧 , so that 𝜔 𝑘 is the length of the spatial Fourier vector. The Fourier space has the same spatial frequency dimension as the spatial dimension of the problem. So in 2D space, a 𝜔 𝑧 spatial frequency variable would not be required and so forth.

By inverse spatial Fourier transformation of (4.1), the wavefield in 𝑛 dimensions is given by 𝑢 = 1 ⃗ 𝑟 , 𝜔 ( 2 𝜋 ) 𝑛 ∞ − ∞ 𝑆 ⃗ 𝜔 , 𝜔 𝜔 2 𝑘 − 𝑘 2 𝑒 𝑖 ⃗ 𝜔 ⋅ ⃗ 𝑟 𝑑 ⃗ 𝜔 . ( 4 . 2 ) Using the definition of the Fourier transform of 𝑠 ( ⃗ 𝑟 , 𝜔 ) , the above equation can be rewritten as 𝑢 = 1 ⃗ 𝑟 , 𝜔 ( 2 𝜋 ) 𝑛 ∞ − ∞ 𝑒 𝑖 ⃗ 𝜔 ⋅ ⃗ 𝑟 𝜔 2 𝑘 − 𝑘 2 ∞ − ∞ 𝑠 𝑒 ⃗ 𝑥 , 𝜔 − 𝑖 ⃗ 𝜔 ⋅ ⃗ 𝑥 𝑑 ⃗ 𝑥 𝑑 ⃗ 𝜔 . ( 4 . 3 ) Let us define the spatial Green’s function in 𝑛 -dimensional space as 𝑔 = 1 ⃗ 𝑟 ∣ ⃗ 𝑥 , 𝜔 = 𝑔 ⃗ 𝑟 − ⃗ 𝑥 , 𝜔 ( 2 𝜋 ) 𝑛 ∞ − ∞ 𝑒 𝑖 ⃗ 𝜔 ⋅ ( ⃗ 𝑟 − ⃗ 𝑥 ) 𝜔 2 𝑘 − 𝑘 2 𝑑 ⃗ 𝜔 . ( 4 . 4 ) This is the Green's function for the (pseudo-) Helmholtz equation. By switching the order of integration, (4.3) can be rewritten such that it can be clearly interpreted as a spatial convolution of Green's function with the input source: 𝑢 = ⃗ 𝑟 , 𝜔 ∞ − ∞ 𝑠 𝑔 ⃗ 𝑥 , 𝜔 ⃗ 𝑟 − ⃗ 𝑥 , 𝜔 𝑑 ⃗ 𝑥 = 𝑠 ⃗ 𝑟 , 𝜔 ∗ ⃗ 𝑟 𝑔 ⃗ 𝑟 , 𝜔 . ( 4 . 5 ) The notation ∗ ⃗ 𝑟 has been used to denote a (multidimensional) space-only convolution. The frequency domain equivalent to (4.5) is (4.1), which can be interpreted as a multiplication in terms of the spatial transfer function as 𝑈 𝑆 ⃗ 𝜔 , 𝜔 = 𝐺 ⃗ 𝜔 , 𝜔 ⃗ 𝜔 , 𝜔 , ( 4 . 6 ) where the Fourier transform of the spatial Green's function is the spatial transfer function: 𝐺 𝑔 = 1 ⃗ 𝜔 , 𝜔 = 𝔽 ⃗ 𝑟 , 𝜔 𝜔 2 𝑘 − 𝑘 2 . ( 4 . 7 ) The dependence of the spatial Green's function on temporal frequency is via the wavenumber 𝑘 of the wavefield (which may be real or complex) and this is emphasized by writing 𝐺 ( ⃗ 𝜔 , 𝜔 ) as a function of spatial and temporal frequencies.

4.1. Spatial Transfer Function in Curvilinear Coordinates

Equations (4.4) and (4.7) hold in any 𝑛 -dimensional space with Green's function of (4.4) being the 𝑛 -dimensional spatial inverse Fourier transform of the spatial transfer function of (4.7). More specifically, in 1D, 𝜔 2 𝑘 = 𝜔 2 𝑥 , and in 2D 𝜔 2 𝑘 = 𝜔 2 𝑥 + 𝜔 2 𝑦 . Transforming to polar coordinates in the 2D spatial frequency domain via the transformation 𝜔 𝑥 = 𝜌 c o s 𝜓 , 𝜔 𝑦 = 𝜌 s i n 𝜓 gives 𝜔 2 𝑘 = 𝜔 2 𝑥 + 𝜔 2 𝑦 = 𝜌 2 where 𝜌 is the radial frequency variable in polar coordinates. Similarly, in 3D, 𝜔 2 𝑘 = 𝜔 2 𝑥 + 𝜔 2 𝑦 + 𝜔 2 𝑧 becomes 𝜔 2 𝑘 = 𝜔 2 𝑥 + 𝜔 2 𝑦 + 𝜔 2 𝑧 = 𝜌 2 where 𝜌 is the spherical radial frequency variable in spherical polar coordinates in spatial frequency space via the coordinate transformation given by 𝜔 𝑥 = 𝜌 s i n 𝜓 𝜔 c o s 𝜃 𝜔 , 𝜔 𝑦 = 𝜌 s i n 𝜓 𝜔 s i n 𝜃 𝜔 and 𝜔 𝑧 = 𝜌 c o s 𝜓 𝜔 . Thus, the spatial transfer functions assume particularly simple, radially symmetric forms in polar coordinates for 2D space and spherical polar coordinates for 3D space. These particularly simple forms for the 2D and 3D transfer functions in polar coordinates motivates the development of the general 2D and 3D Fourier transforms in curvilinear coordinates [28, 29] so that this simple form of the transfer function may be exploited.

5. Notation and Sign Conventions

Wavenumbers are defined by the squares of their quantities and arise as a result of taking the Fourier transform of the corresponding propagation equation (be it acoustic, thermal, or otherwise), which in turn leads to a Helmholtz or pseudo-Helmholtz equation. These wavenumbers are defined so that the wave propagation equation transformed to the temporal frequency domain all have the form of a pseudo-Helmholz equation as given by (3.1) or (3.5) with a complex or real wavenumber. For the rest of the paper, we will consider that 𝑘 𝑡 represents a generic complex wavenumber while 𝑘 𝑠 represents a generic real wavenumber.

Rather than the squared wavenumber, the quantity of interest will prove to be the wavenumber itself, namely, 𝑘 𝑠 and 𝑘 𝑡 , which are the square roots of the given squared wavenumber in the Helmholtz equation. Each 𝑘 can be considered as the sum of a real and an imaginary part, so that 𝑘 𝑡 = 𝑘 𝑡 𝑟 + 𝑖 𝑘 𝑡 𝑖 with 𝑘 𝑡 𝑟 denoting the real part of 𝑘 𝑡 and 𝑘 𝑡 𝑖 denoting the imaginary part. Since the square root of any 𝑘 2 can be ± 𝑘 , we will use the convention that for a given complex 𝑘 , the required square root of the corresponding 𝑘 2 is defined such that the imaginary part of 𝑘 is negative. Hence, 𝑘 𝑡 is chosen as the square root of 𝑘 2 𝑡 such that 𝑘 𝑡 𝑖 < 0 and so forth. If this sign convention is adopted, then it was shown in [36] that the many results for complex wavenumber use the same mathematical form of travelling wave solution as for the real wavenumbers.This makes the notation and book-keeping considerably simpler. It is noted that this sign convention is the opposite of what this author adopts in [30].

5.1. Sommerfeld Radiation Condition

To aid in the selection of a causal solution, the Sommerfeld radiation condition is required. The Sommerfeld radiation condition states that the sources in the field must be sources not sinks of energy. Therefore, energy radiated from sources must scatter to infinity and cannot radiate from infinity into the field. Mathematically, a solution 𝑢 ( 𝑥 ) , where 𝑥 is the spatial variable, to the Helmholtz equation is considered to be radiating if it satisfies l i m | 𝑥 | → ∞ | 𝑥 | ( 𝑛 − 1 ) / 2 𝜕 𝜕 | 𝑥 | + 𝑖 𝑘 𝑢 ( 𝑥 ) = 0 , ( 5 . 1 ) where 𝑛 is the dimension of the space and 𝑘 is the wavenumber in the Helmholtz equation. This is the radiation equation based on an implied time variation of 𝑒 𝑖 𝜔 𝑡 which is implicit in our chosen definition of the Fourier transform, (2.1). For example, with the standard definition of the Fourier transform, the transform of 𝑓 ( 𝑡 ) is 𝑖 𝜔 𝐹 ( 𝜔 ) , clearly implying the 𝑒 𝑖 𝜔 𝑡 dependence. Had a different definition of the Fourier transform been used, the implied time variation would have been 𝑒 − 𝑖 𝜔 𝑡 , in which case the sign of 𝑖 in (5.1) would be reversed.

6. 2D Fourier Transforms in Terms of Polar Coordinates

Given the desire to exploit the simple form of the spatial transfer function in curvilinear coordinates, we consider 2D and 3D Fourier transform in terms of curvilinear coordinates. Let us first consider the Fourier transforms in 2D in polar coordinates. In order to define this, some preliminary definitions of Hankel transforms are required first.

6.1. Hankel Transforms

The Hankel transform of order 𝑛 is defined by the integral [31] ⌢ 𝐹 𝑛 ( 𝜌 ) = ℍ 𝑛 ( 𝑓 ( 𝑟 ) ) = ∞ 0 𝑓 ( 𝑟 ) 𝐽 𝑛 ( 𝜌 𝑟 ) 𝑟 𝑑 𝑟 , ( 6 . 1 ) where 𝐽 𝑛 ( 𝑧 ) is the 𝑛 th order Bessel function. If 𝑛 > − 1 / 2 , the transform is self-reciprocating and the inversion formula is given by 𝑓 ( 𝑟 ) = ℍ 𝑛 − 1 ⌢ 𝐹 𝑛 = ( 𝜌 ) ∞ 0 ⌢ 𝐹 𝑛 ( 𝜌 ) 𝐽 𝑛 ( 𝜌 𝑟 ) 𝜌 𝑑 𝜌 . ( 6 . 2 )

6.2. Connection between the 2D Fourier Transform and Hankel Transform

To develop Fourier transforms in 2D in polar coordinates, both the function 𝑓 ( ⃗ 𝑟 ) and its Fourier transform 𝐹 ( ⃗ 𝜔 ) are expressed in polar coordinates. In general 𝑓 ( ⃗ 𝑟 ) = 𝑓 ( 𝑟 , 𝜃 ) is not radially symmetric and is a function of both 𝑟 and 𝜃 so that the 𝜃 dependence can be expanded into a Fourier series due to the 2 𝜋 periodicity of the function in 𝜃 : 𝑓 ( 𝑟 , 𝜃 ) = ∞ 𝑛 = − ∞ 𝑓 𝑛 ( 𝑟 ) 𝑒 𝑗 𝑛 𝜃 , ( 6 . 3 ) where the Fourier coefficients 𝑓 𝑛 ( 𝑟 ) can be found from 𝑓 𝑛 1 ( 𝑟 ) = 2 𝜋 0 2 𝜋 𝑓 ( 𝑟 , 𝜃 ) 𝑒 − 𝑗 𝑛 𝜃 𝑑 𝜃 . ( 6 . 4 ) Similarly, the 2D Fourier transform 𝐹 ( ⃗ 𝜔 ) = 𝐹 ( 𝜌 , 𝜓 ) can also be expanded into its own Fourier series so that 𝐹 ( 𝜌 , 𝜓 ) = ∞ 𝑛 = − ∞ 𝐹 𝑛 ( 𝜌 ) 𝑒 𝑗 𝑛 𝜓 , ( 6 . 5 ) and where those Fourier coefficients are similarly found from 𝐹 𝑛 1 ( 𝜌 ) = 2 𝜋 0 2 𝜋 𝐹 ( 𝜌 , 𝜓 ) 𝑒 − 𝑗 𝑛 𝜓 𝑑 𝜓 , ( 6 . 6 )

The relationship between the Fourier coordinates in normal space 𝑓 𝑛 ( 𝑟 ) in (6.3) and the spatial frequency space Fourier coordinates 𝐹 𝑛 ( 𝜌 ) in (6.5) is desired. The details of this are given in [28], and the results are summarized here. It is emphasized that the relationship between 𝑓 𝑛 ( 𝑟 ) and 𝐹 𝑛 ( 𝜌 ) is not a Fourier transform. The relationship between them is given by 𝐹 𝑛 ( 𝜌 ) = 2 𝜋 𝑖 − 𝑛 ∞ 0 𝑓 𝑛 ( 𝑟 ) 𝐽 𝑛 ( 𝜌 𝑟 ) 𝑟 𝑑 𝑟 = 2 𝜋 𝑖 − 𝑛 ℍ 𝑛 𝑓 𝑛 ( 𝑟 ) , ( 6 . 7 ) where ℍ 𝑛 is the 𝑛 th order Hankel transform. The reverse relationship is given by 𝑓 𝑛 𝑖 ( 𝑟 ) = 𝑛 2 𝜋 ∞ 0 𝐹 𝑛 ( 𝜌 ) 𝐽 𝑛 𝑖 ( 𝜌 𝑟 ) 𝜌 𝑑 𝜌 = 𝑛 ℍ 2 𝜋 𝑛 𝐹 𝑛 ( 𝜌 ) . ( 6 . 8 ) The 𝑛 th term in the Fourier series for the original function will Hankel transform into the 𝑛 th term of the Fourier series of the Fourier transform. However, it is an 𝑛 th order Hankel transform for the 𝑛 th term, namely, all the terms are not equivalently transformed. The mapping from 𝑓 𝑛 ( 𝑟 ) to 𝐹 𝑛 ( 𝜌 ) is one of 𝑛 th-order Hankel transform, which in general is not a 2D Fourier transform.

The operation of taking the 2D Fourier transform of a function is thus equivalent to (1) first finding its Fourier series expansion in the angular variable and (2) then finding the 𝑛 th-order Hankel transform (of the radial variable to the spatial radial variable) of the 𝑛 th coefficient in the Fourier series. Since each of these operations involves integration over one variable only with the others being considered parameters vis-à-vis the integration, the order in which these operations are performed is interchangeable.

7. Spherical Hankel Transform

We introduce the spherical Hankel transform, which will form part of the 3D Fourier transform in spherical coordinates.

7.1. Definition of the Spherical Hankel Transform

The spherical Hankel transform can then be defined as [32, 33] ⌢ 𝐹 𝑛 ( 𝜌 ) = 𝕊 𝑛 { 𝕗 𝑓 ( 𝑟 ) } = ∞ 0 𝑓 ( 𝑟 ) 𝑗 𝑛 ( 𝜌 𝑟 ) 𝑟 2 𝑑 𝑟 . ( 7 . 1 ) 𝕊 𝑛 is used to specifically denote the spherical Hankel transform of order 𝑛 . The inverse transform is given by [29] 2 𝑓 ( 𝑟 ) = 𝜋 ∞ 0 ⌢ 𝐹 𝑛 ( 𝜌 ) 𝑗 𝑛 ( 𝜌 𝑟 ) 𝜌 2 𝑑 𝜌 . ( 7 . 2 ) The spherical Hankel transform is particularly useful for problems involving spherical symmetry.

7.2. Spherical Harmonics

The spherical harmonics are the solution to the angular portion of Laplace's equation in spherical polar coordinates and can be shown to be orthogonal. These spherical harmonics are given by [32] 𝑌 𝑚 𝑙 ( 𝜓 , 𝜃 ) = ( 2 𝑙 + 1 ) ( 𝑙 − 𝑚 ) ! 𝑃 4 𝜋 ( 𝑙 + 𝑚 ) ! 𝑚 𝑙 ( c o s 𝜓 ) 𝑒 𝑖 𝑚 𝜃 , ( 7 . 3 ) where 𝑌 𝑚 𝑙 is called a spherical harmonic function of degree 𝑙 and order 𝑚 , 𝑃 𝑚 𝑙 is an associated Legendre function, 0 ≤ 𝜓 ≤ 𝜋 represents the colatitude and 0 ≤ 𝜃 ≤ 2 𝜋 represents the longitude. With the normalization of the spherical harmonics as given in (7.3), the spherical harmonics are orthonormal so that 0 2 𝜋 𝜋 0 𝑌 𝑚 𝑙 𝑌 𝑚 𝑙 s i n 𝜓 𝑑 𝜓 𝑑 𝜃 = 𝛿 𝑙 𝑙 𝛿 𝑚 𝑚 . ( 7 . 4 ) Here 𝛿 𝑖 𝑗 is the kronecker delta and the overbar indicates the complex conjugate. It is important to note that there are several different normalizations of the spherical harmonics that are possible, that will differ from (7.3) and thus lead to a different version of (7.4). Several of these are nicely catalogued in [34], including which disciplines tend to use which normalization. This is important to note since any result that uses orthogonality will differ slightly depending on the choice of normalization.

The spherical harmonics form a complete set of orthonormal functions and thus form a vector space. When restricted to the surface of a sphere, functions may be expanded on the sphere into a series approximation much like a Fourier series. This is in fact a spherical harmonic series. Any square-integrable function may be expanded as 𝑓 ( 𝜓 , 𝜃 ) = ∞ 𝑙 𝑙 = 0 𝑚 = − 𝑙 𝑓 𝑚 𝑙 𝑌 𝑚 𝑙 ( 𝜓 , 𝜃 ) , ( 7 . 5 ) where 𝑓 𝑚 𝑙 = 0 2 𝜋 𝜋 0 𝑓 ( 𝜓 , 𝜃 ) 𝑌 𝑚 𝑙 ( 𝜓 , 𝜃 ) s i n 𝜓 𝑑 𝜓 𝑑 𝜃 . ( 7 . 6 ) The coefficients 𝑓 𝑚 𝑙 are sometimes referred to as the spherical Fourier transform of 𝑓 ( 𝜓 , 𝜃 ) [35].

7.3. 3D Fourier Transforms in Spherical Polar Coordinates

To find 3D Fourier transforms in spherical polar coordinates, both the function and its 3D Fourier transform are written in terms of polar coordinates in the spatial and spatial frequency domains. That is, a function in 3D space is expressed as 𝑓 ( ⃗ 𝑟 ) = 𝑓 ( 𝑟 , 𝜓 𝑟 , 𝜃 𝑟 ) as a function of spherical polar coordinates and its 3D Fourier transform is written as 𝐹 ( ⃗ 𝜔 ) = 𝐹 ( 𝜌 , 𝜓 𝜔 , 𝜃 𝜔 ) in frequency spherical polar coordinates. Certain relationships can be shown to hold. The relationship between the function and its transform are summarized here, with the relevant details omitted for brevity. The function itself can be expanded as a series in terms of the spherical harmonics as 𝑓 ⃗ 𝑟 = 𝑓 𝑟 , 𝜓 𝑟 , 𝜃 𝑟 = ∞ 𝑙 𝑙 = 0 𝑘 = − 𝑙 𝑓 𝑘 𝑙 ( 𝑟 ) 𝑌 𝑘 𝑙 𝜓 𝑟 , 𝜃 𝑟 , ( 7 . 7 ) where 𝑌 𝑘 𝑙 ( 𝜓 𝑟 , 𝜃 𝑟 ) are the spherical harmonics and the Fourier coefficients are given by 𝑓 𝑘 𝑙 ( 𝑟 ) = 0 2 𝜋 𝜋 0 𝑓 𝑟 , 𝜓 𝑟 , 𝜃 𝑟 𝑌 𝑘 𝑙 𝜓 𝑟 , 𝜃 𝑟 s i n 𝜓 𝑟 𝑑 𝜓 𝑟 𝑑 𝜃 𝑟 . ( 7 . 8 ) The 3D Fourier transform of 𝑓 can also be written in Fourier space in spherical polar coordinates as 𝐹 ⃗ 𝜔 = 𝐹 𝜌 , 𝜓 𝜔 , 𝜃 𝜔 = ∞ 𝑙 𝑙 = 0 𝑘 = − 𝑙 𝐹 𝑘 𝑙 ( 𝜌 ) 𝑌 𝑘 𝑙 𝜓 𝜔 , 𝜃 𝜔 , ( 7 . 9 ) where 𝐹 𝑘 𝑙 ( 𝜌 ) = 0 2 𝜋 𝜋 0 𝐹 𝜌 , 𝜓 𝜔 , 𝜃 𝜔 𝑌 𝑘 𝑙 𝜓 𝜔 , 𝜃 𝜔 s i n 𝜓 𝜔 𝑑 𝜓 𝜔 𝑑 𝜃 𝜔 . ( 7 . 1 0 ) We emphasize again that the relationship between 𝑓 𝑘 𝑙 ( 𝑟 ) and 𝐹 𝑘 𝑙 ( 𝜌 ) is not that of a Fourier transform. In fact, this relationship is given by [29] 𝐹 𝑘 𝑙 ( 𝜌 ) = 4 𝜋 ( − 𝑖 ) 𝑙 ∞ 0 𝑓 𝑘 𝑙 ( 𝑟 ) 𝑗 𝑙 ( 𝜌 𝑟 ) 𝑟 2 𝑑 𝑟 = 4 𝜋 ( − 𝑖 ) 𝑙 𝕊 𝑙 𝑓 𝑘 𝑙 , ( 𝑟 ) ( 7 . 1 1 ) where 𝕊 𝑙 denotes a spherical Hankel transform of order 𝑙 . The inverse relationship is given by 𝑓 𝑘 𝑙 1 ( 𝑟 ) = 4 𝜋 ( 𝑖 ) 𝑙 2 𝜋 ∞ 0 𝐹 𝑘 𝑙 ( 𝜌 ) 𝑗 𝑙 ( 𝜌 𝑟 ) 𝜌 2 1 𝑑 𝜌 = 4 𝜋 ( 𝑖 ) 𝑙 𝕊 𝑙 − 1 𝐹 𝑘 𝑙 ( 𝜌 ) . ( 7 . 1 2 )

8. Spatial Transfer Function Relationship in 3D

With this long set of preliminaries, definitions, and conventions established, we can now begin the derivation of the 3D spatial transfer function relationship between input 𝑠 ( ⃗ 𝑟 , 𝜔 ) and output 𝑢 ( ⃗ 𝑟 , 𝜔 ) . The foregoing section makes use of several key results in [36] which will be presented here without proof. The results in [36] address general cases, not all of which are relevant to the work here. The versions of the theorems required herein are stated based on the sign conventions presented in Section 5 with regards to the sign conventions for the complex wavenumbers and also for the chosen definition of the Fourier transform. These affect the implied time dependence, the Sommerfeld radiation condition, and thus the physically correct solution that is chosen from several mathematical possibilities. The versions of the theorems presented here ensure that the chosen waves are outwardly propagating (Sommerfeld radiation condition) and also decay to zero at infinity for a damped wavefield, thus ensuring bounded and physically meaningful solutions.

Theorem 8.1. It is shown in [36] that the following result holds true: 𝐼 = ∞ 0 𝜙 ( 𝑥 ) 𝑥 2 − 𝑘 2 𝑗 𝑛 ( 𝑥 𝑟 ) 𝑥 2 𝑘 𝑑 𝑥 = − 𝜋 𝑖 2 ℎ 𝑛 ( 2 ) ( 𝑘 𝑟 ) 𝜙 ( 𝑘 ) . ( 8 . 1 ) Here, 𝜙 is any analytic function defined on the positive real line that remains bounded as 𝑥 goes to infinity, 𝑗 𝑛 ( 𝑥 ) is a spherical Bessel function of order 𝑛 , ℎ 𝑛 ( 2 ) ( 𝑥 ) is a spherical Hankel function of order 𝑛 , and 𝑘 is a wavenumber which may be real or complex, chosen such that the imaginary part of 𝑘 is negative. Given the definition of the Fourier transform that is being currently used, the presented result satisfies the Sommerfeld radiation condition, ensuring an outwardly propagating wave.

8.1. 3D Transfer Function Relationship in Spherical Polar Coordinates

From (4.7) and using the conversion to spherical polar coordinates in 3D, the transfer function for the Helmholtz equation in the Fourier domain can be written as 𝐺 = 1 ⃗ 𝜔 , 𝜔 𝜌 2 − 𝑘 2 , ( 8 . 2 ) which depends on the frequency spherical radial variable 𝜌 only.

In general, as the imaginary part of the wavenumber 𝑘 gets smaller, the transfer function becomes more frequency selective, in the sense of strongly passing a smaller bandwidth. As the imaginary part of the wavenumber 𝑘 gets larger, the transfer function becomes more low-pass in nature, passing lower frequencies more strongly than the higher frequencies, with a corresponding larger bandwidth of frequencies being passed. This is physically meaningful as the imaginary part of 𝑘 represents the damping inherent in this wave modality. Wavenumbers with larger imaginary parts are more heavily damped and that damping tends to affect the higher frequencies more—that is, the higher frequencies are attenuated and the lower ones are passed through, implying a low-pass nature to the physical system.

Recall that 𝜌 is the magnitude of the spatial frequency variables. The case for the imaginary part of 𝑘 equal to zero (in other words, a real wavenumber) is easy to visualize. For the case of a real 𝑘 , the function is discontinuous at 𝜌 = 𝑘 , meaning that only those spatial frequencies where exactly 𝜌 = 𝑘 are passed and the rest are attenuated. In other words, the resulting wave must have 𝜌 2 = 𝜔 2 𝑥 + 𝜔 2 𝑦 + 𝜔 2 𝑧 = 𝑘 2 . While these statements are fairly mathematical in nature, they explain the nature of the resolution of various imaging modalities. For example, it is far more difficult to achieve the resolution of acoustic imaging with thermal imaging and one of the explanations for this is the “low-pass” blurring nature of the complex wavenumber of a thermal wave versus the highly selective band-pass nature of the acoustic wavenumber. This is discussed further in other papers, for example, in [18].

8.2. 3D Green Function Coefficients

Let us define a set of functions that will be referred to as the Green function coefficients. The actual Green function for the system is the full 3D inverse Fourier transform of the transfer function, which for a spherically symmetric function is equivalent to a spherical Hankel transform of order zero only [33]. With the help of Theorem 8.1, we define these Green function coefficients as 𝑔 3 D 𝑛 ( 𝑟 , 𝑘 ) = 𝕊 𝑛 − 1 1 𝜌 2 − 𝑘 2 = 2 𝜋 ∞ 0 1 𝜌 2 − 𝑘 2 𝑗 𝑛 ( 𝜌 𝑟 ) 𝜌 2 𝑑 𝜌 = − 𝑖 𝑘 ℎ 𝑛 ( 2 ) ( 𝑘 𝑟 ) , ( 8 . 3 ) where 𝑔 3 D 𝑛 ( 𝑟 , 𝑘 ) has been used to denote the 𝑛 th order Green function coefficient, with the subscript indicating the order of the coefficient and the wavenumber 𝑘 included as a parameter of the function.

8.3. Transfer Function Relationship in Spherical Polar Coordinates

The 3D Fourier transform of the input (source) function is written in polar coordinates and expanded in terms of a spherical harmonic series as 𝑆 = ⃗ 𝜔 , 𝜔 ∞ 𝑙 𝑙 = 0 𝑚 = − 𝑙 𝑆 𝑚 𝑙 ( 𝜌 , 𝜔 ) 𝑌 𝑚 𝑙 𝜓 𝜔 , 𝜃 𝜔 . ( 8 . 4 )

Hence, (4.1) for the output wavefield becomes 𝑈 = ⃗ 𝜔 , 𝜔 ∞ 𝑙 𝑙 = 0 𝑚 = − 𝑙 𝑆 𝑚 𝑙 ( 𝜌 , 𝜔 ) 𝜌 2 − 𝑘 2 𝑌 𝑚 𝑙 𝜓 𝜔 , 𝜃 𝜔 . ( 8 . 5 ) The symbol 𝑘 is used to denote the wavenumber. The spatial inverse Fourier transform and thus the output in spatial coordinates is given by 𝑢 𝑟 , 𝜓 𝑟 , 𝜃 𝑟 = , 𝜔 ∞ 𝑙 𝑙 = 0 𝑚 = − 𝑙 ( 𝑖 ) 𝑙 2 4 𝜋 𝜋 ∞ 0 𝑆 𝑚 𝑙 ( 𝜌 , 𝜔 ) 𝑗 𝑙 ( 𝜌 𝑟 ) 𝜌 2 − 𝑘 2 𝜌 2 𝑌 𝑑 𝜌 𝑚 𝑙 𝜓 𝑟 , 𝜃 𝑟 . ( 8 . 6 ) The temporal frequency 𝜔 is a constant as far as the spatial inverse Fourier transformation is concerned. We recognize that the quantity within the curly brackets can be evaluated with the help of Theorem 8.1 and using the definition of the Green function coefficients in (8.3) to obtain 2 𝜋 ∞ 0 𝑆 𝑚 𝑙 ( 𝜌 , 𝜔 ) 𝑗 𝑙 ( 𝜌 𝑟 ) 𝜌 2 − 𝑘 2 𝜌 2 𝑑 𝜌 = 𝑔 3 D 𝑙 ( 𝑟 , 𝑘 ) 𝑆 𝑚 𝑙 ( 𝑘 , 𝜔 ) . ( 8 . 7 ) It now follows that the wavefield expression becomes 𝑢 = ⃗ 𝑟 , 𝜔 ∞ 𝑙 𝑙 = 0 𝑚 = − 𝑙 𝑖 𝑙 𝑔 4 𝜋 3 D 𝑙 ( 𝑟 , 𝑘 ) 𝑆 𝑚 𝑙 ( 𝑘 , 𝜔 ) 𝑌 𝑚 𝑙 𝜓 𝑟 , 𝜃 𝑟 . ( 8 . 8 ) If the measured wavefield itself, namely, the left-hand side of (8.8), is expanded as in a spherical harmonic series so that 𝑢 ⃗ 𝑟 , 𝜔 = 𝑢 𝑟 , 𝜓 𝑟 , 𝜃 𝑟 = , 𝜔 ∞ 𝑙 𝑙 = 0 𝑚 = − 𝑙 𝑢 𝑚 𝑙 ( 𝑟 , 𝜔 ) 𝑌 𝑚 𝑙 𝜓 𝑟 , 𝜃 𝑟 , ( 8 . 9 ) then (8.8) gives us the simple input-output transfer function relationship we seek: 𝑢 𝑚 𝑙 𝑖 ( 𝑟 , 𝜔 ) = 𝑙 𝑔 4 𝜋 3 D 𝑙 ( 𝑟 , 𝑘 ) 𝑆 𝑚 𝑙 ( 𝑘 , 𝜔 ) . ( 8 . 1 0 ) Note how this relationship is in the spatial domain (not the frequency spatial domain) and gives a multiplicative transfer function relationship between the input coefficients 𝑆 𝑚 𝑙 ( 𝑘 , 𝜔 ) and the resulting wavefield 𝑢 𝑚 𝑙 ( 𝑟 , 𝜔 ) .

Using the relationships proposed between the coefficients of the function in the spatial domain and the coefficients in the Fourier domain as given in (7.11), the general relationship between Hankel and Fourier transforms is given by 𝑆 𝑚 𝑙 ( 𝜌 , 𝜔 ) = 4 𝜋 ( − 𝑖 ) 𝑚 𝕊 𝑙 𝑠 𝑚 𝑙 ( 𝑟 , 𝜔 ) = 4 𝜋 ( − 𝑖 ) 𝑙 ∞ 0 𝑠 𝑚 𝑙 ( 𝑟 , 𝜔 ) 𝑗 𝑙 ( 𝜌 𝑟 ) 𝑟 2 𝑑 𝑟 . ( 8 . 1 1 )

Hence (8.10) can be written as 𝑢 𝑚 𝑙 ( 𝑟 , 𝜔 ) = 𝑔 3 D 𝑙 ( 𝑟 , 𝑘 ) ∞ 0 𝑠 𝑚 𝑙 ( 𝑥 , 𝜔 ) 𝑗 𝑙 ( 𝑘 𝑥 ) 𝑥 2 𝑑 𝑥 , ( 8 . 1 2 ) so that the value of the measured wavefield is related to the spherical Hankel transform of the input function, evaluated at the wavenumber of the wavefield in question. Using the definition of the spherical Hankel transforms, this becomes a very compact yet powerful expression: 𝑢 𝑚 𝑙 ( 𝑟 , 𝜔 ) = 𝑔 3 D 𝑙 ( 𝑟 , 𝑘 ) 𝕊 𝑙 𝑠 𝑚 𝑙 ( 𝑘 , 𝜔 ) . ( 8 . 1 3 ) The value of defining the Green function coefficients now becomes apparent in that it permits the relationship between 𝑢 𝑚 𝑙 ( 𝑟 , 𝜔 ) and 𝑠 𝑚 𝑙 ( 𝑘 , 𝜔 ) to have a simple multiplicative transfer function relationship. In other words, the relationship between output wave coefficients and input is one of multiplication with the transfer function coefficients instead of a complicated convolution-type relationship.

If we also assume that the input function is separable in time and space (a fairly general assumption) so that 𝑠 ( ⃗ 𝑟 , 𝜔 ) = 𝑞 ( ⃗ 𝑟 ) 𝜂 ( 𝜔 ) . This would be the case for a spatial inhomogeneity and a temporal input. In this case, then 𝑠 𝑚 𝑙 ( 𝑟 , 𝜔 ) = 𝑞 𝑚 𝑙 ( 𝑟 ) 𝜂 ( 𝜔 ) and the formula for the forward problem is given by 𝑢 𝑚 𝑙 ( 𝑟 , 𝜔 ) = 𝑔 3 D 𝑙 ( 𝑟 , 𝑘 ) 𝜂 ( 𝜔 ) ∞ 0 𝑞 𝑚 𝑙 ( 𝑟 ) 𝑗 𝑙 ( 𝑘 𝑟 ) 𝑟 2 𝑑 𝑟 = 𝑔 3 D 𝑙 ( 𝑟 , 𝑘 ) 𝜂 ( 𝜔 ) 𝕊 𝑙 𝑞 𝑚 𝑙 ( 𝑘 ) . ( 8 . 1 4 )

8.4. Discussion and Relationship with the Fourier Diffraction Theorem

Equation (8.13) is the key input-output relationship that is sought. It gives the relationship between input and output and the relationship is a transfer function (multiplication) type of relationship, even in the spatial domain where the normal relationship to be expected is one of convolution. This relationship is for the spherical harmonic expansion coefficients, not between the full functions themselves. However, because it is a direct, proportional type of relationship between the ( 𝑙 , 𝑚 ) th term of the input and the ( 𝑙 , 𝑚 ) th term of the output, it is particularly useful and simple to apply.

In particular, (8.13) states that the wavefield (output) coefficients are directly proportional to the transfer function coefficients and the proportionality factor between them is the spherical Hankel transform of the input coefficients, evaluated on the sphere 𝜌 = 𝑘 . Note that ( 𝑙 , 𝑚 ) th order output coefficients are related to ( 𝑙 , 𝑚 ) th order transfer function coefficients in proportion to the 𝑙 th order spherical Hankel transform of the ( 𝑙 , 𝑚 ) th order input coefficient. This is in fact a generalization of the Fourier diffraction theorem of tomography [3] which loosely states that the output wavefield is proportional to the Fourier transform of the input inhomogeneity evaluated somewhere on the 𝜌 = 𝑘 sphere in Fourier space. This is still true but we have expressed a more precise version of this theorem in the sense that we have removed any ambiguity regarding the “somewhere” on the 𝜌 = 𝑘 sphere and replaced it with a relationship (8.14) that does not depend on angular location. This was enabled by the use of the Fourier transform in curvilinear coordinates so that a full Fourier transform requires an ( 𝑙 , 𝑚 ) set of Fourier coefficients and spherical Hankel transforms. Furthermore, the idea of proportionality between Fourier transform of the inhomogeneity and output wavefield is also made more precise by stating that the proportionality between them is actually the coefficients of the Green function itself. In essence, the output wavefield can be seen as being the Green function coefficients, with each term weighted by the spherical Hankel transform of the input inhomogeneity.

9. 2D Transfer Function Relationship in 2D Polar Coordinates

We proceed to develop the transfer function expressions for the 2D case in polar coordinates. As for the 3D case, we begin with a necessary theorem which is stated without proof. As previously mentioned, the results in [36] address general cases, and the versions of these theorems required herein are stated here based on the sign conventions presented in Section 5 for the complex wavenumbers and also for the Fourier transform which affects the implied time dependence and thus the Sommerfeld radiation condition.

Theorem 9.1. It is shown in [36] that the following result holds true: 𝐼 = ∞ 0 𝜙 ( 𝜌 ) 𝜌 2 − 𝑘 2 𝐽 𝑛 1 ( 𝜌 𝑟 ) 𝜌 𝑑 𝜌 = − 𝜋 𝑖 2 𝐻 𝑛 ( 2 ) ( 𝑘 𝑟 ) 𝜙 ( 𝑘 ) . ( 9 . 1 ) Here, 𝜙 is an analytic function defined on the positive real line that remains bounded as 𝑥 goes to infinity, 𝐽 𝑛 ( 𝑥 ) is a Bessel function of order 𝑛 , 𝐻 𝑛 ( 2 ) ( 𝑥 ) is a Hankel function of order 𝑛 , and 𝑘 is a wavenumber which may be real or complex, chosen such that the imaginary part of 𝑘 is negative. Given the definition of the Fourier transform that is being currently used, the presented result satisfies the Sommerfeld radiation condition, ensuring an outwardly propagating wave.

9.1. 2D Green Function Coefficients

The overall system transfer function in 2D is given by 𝐺 = 1 ⃗ 𝜔 , 𝜔 𝜌 2 − 𝑘 2 , ( 9 . 2 ) which depends on the frequency radial variable only.

As for the 3D case, we define the 2D Green function coefficients as the 𝑛 th-order inverse Hankel transform of the overall system spatial transfer function. These will be shown to be the required transfer functions for working in the 2D polar formulation. From the definition of the inverse Hankel transform, these are given by 𝑔 2 D 𝑛 ( 𝑟 , 𝑘 ) = ℍ 𝑛 − 1 1 𝜌 2 − 𝑘 2 = ∞ 0 1 𝜌 2 − 𝑘 2 𝐽 𝑛 ( 𝜌 𝑟 ) 𝜌 𝑑 𝜌 = − 𝜋 𝑖 2 𝐻 𝑛 ( 2 ) ( 𝑘 𝑟 ) . ( 9 . 3 ) Recall that the actual (full) Green function for the system is the full 2D inverse Fourier transform of the transfer function, which is an inverse Hankel transform of order zero only. In (9.3), 𝑔 2 D 𝑛 ( 𝑟 , 𝑘 ) has been used to denote the nth order Green function coefficients with the subscript indicating the order.

9.2. Fourier Theorem for 2D Fourier Transforms in Polar Coordinates

The 2D Fourier transform of the input function is written in polar coordinates as 𝑆 ( 𝜌 , 𝜓 , 𝜔 ) = ∞ 𝑛 = − ∞ 𝑆 𝑛 ( 𝜌 , 𝜔 ) 𝑒 𝑗 𝑛 𝜓 , ( 9 . 4 ) where 𝑆 𝑛 ( 𝜌 , 𝜔 ) can be found from 𝑆 𝑛 1 ( 𝜌 , 𝜔 ) = 2 𝜋 0 2 𝜋 𝑆 ( 𝜌 , 𝜓 , 𝜔 ) 𝑒 − 𝑗 𝑛 𝜓 𝑑 𝜓 . ( 9 . 5 ) The output wavefield 𝑢 ( ⃗ 𝑟 , 𝜔 ) is similarly expanded as a series so that in the spatial domain we can write 𝑢 ⃗ 𝑟 , 𝜔 = 𝑢 ( 𝑟 , 𝜃 , 𝜔 ) = ∞ 𝑛 = − ∞ 𝑢 𝑛 ( 𝑟 ) 𝑒 𝑗 𝑛 𝜃 , ( 9 . 6 ) while in the spatial Fourier domain this is 𝑈 ⃗ 𝜔 , 𝜔 = 𝑈 ( 𝜌 , 𝜓 , 𝜔 ) = ∞ 𝑛 = − ∞ 𝑈 𝑛 ( 𝜌 , 𝜔 ) 𝑒 𝑗 𝑛 𝜓 . ( 9 . 7 ) The output wavefield is given by 𝑈 ( ⃗ 𝜔 , 𝜔 ) = 𝐺 ( ⃗ 𝜔 , 𝜔 ) 𝑆 ( ⃗ 𝜔 , 𝜔 ) and since 𝐺 ( ⃗ 𝜔 , 𝜔 ) is radially symmetric and does not require a full series, the output wavefield in the Fourier domain is given by 𝑈 ( 𝜌 , 𝜓 , 𝜔 ) = ∞ 𝑛 = − ∞ 𝑈 𝑛 ( 𝜌 , 𝜔 ) 𝑒 𝑗 𝑛 𝜓 = ∞ 𝑛 = − ∞ 𝑆 𝑛 ( 𝜌 , 𝜔 ) 𝜌 2 − 𝑘 2 𝑒 𝑗 𝑛 𝜓 , ( 9 . 8 ) or equivalently as 𝑈 𝑛 𝑆 ( 𝜌 , 𝜔 ) = 𝑛 ( 𝜌 , 𝜔 ) 𝜌 2 − 𝑘 2 . ( 9 . 9 ) The equivalent expression to (9.8) in the spatial domain is given by 𝑢 = ⃗ 𝑟 , 𝜔 ∞ 𝑛 = − ∞ 𝑢 𝑛 ( 𝑟 ) 𝑒 𝑗 𝑛 𝜃 = ∞ 𝑛 = − ∞ 𝑖 𝑛 2 𝜋 ∞ 0 𝑆 𝑛 ( 𝜌 , 𝜔 ) 𝜌 2 − 𝑘 2 𝐽 𝑛 𝑒 ( 𝜌 𝑟 ) 𝜌 𝑑 𝜌 𝑗 𝑛 𝜃 , ( 9 . 1 0 ) or 𝑢 𝑛 𝑖 ( 𝑟 ) = 𝑛 2 𝜋 ∞ 0 𝑆 𝑛 ( 𝜌 , 𝜔 ) 𝜌 2 − 𝑘 2 𝐽 𝑛 ( 𝜌 𝑟 ) 𝜌 𝑑 𝜌 . ( 9 . 1 1 ) The temporal frequency 𝜔 is a constant as far as the spatial inverse Fourier transformation is concerned. The quantity within the curly brackets can be evaluated with the help of Theorem 9.1 as ∞ 0 𝑆 𝑛 ( 𝜌 , 𝜔 ) 𝜌 2 − 𝑘 2 𝐽 𝑛 ( 𝜌 𝑟 ) 𝜌 𝑑 𝜌 = − 𝜋 𝑖 2 𝐻 𝑛 ( 2 ) ( 𝑘 𝑟 ) 𝑆 𝑛 ( 𝑘 , 𝜔 ) = 𝑔 2 D 𝑛 ( 𝑟 , 𝑘 ) 𝑆 𝑛 ( 𝑘 , 𝜔 ) . ( 9 . 1 2 ) Hence, the output wavefield expression is given by 𝑢 = ⃗ 𝑟 , 𝜔 ∞ 𝑛 = − ∞ 𝑢 𝑛 ( 𝑟 ) 𝑒 𝑗 𝑛 𝜃 = ∞ 𝑛 = − ∞ 𝑖 𝑛 𝑔 2 𝜋 2 D 𝑛 ( 𝑟 , 𝑘 ) 𝑆 𝑛 ( 𝑘 , 𝜔 ) 𝑒 𝑗 𝑛 𝜃 . ( 9 . 1 3 ) The definition of forward and inverse Hankel transforms, the general relationship between Hankel and Fourier transforms is given by 𝑆 𝑛 ( 𝜌 ) = 2 𝜋 𝑖 − 𝑛 ℍ 𝑛 𝑠 𝑛 ( 𝑟 ) = 2 𝜋 𝑖 − 𝑛 ∞ 0 𝑠 𝑛 ( 𝑟 ) 𝐽 𝑛 ( 𝜌 𝑟 ) 𝑟 𝑑 𝑟 , ( 9 . 1 4 ) and may be used to simplify (9.13) to yield 𝑢 = ⃗ 𝑟 , 𝜔 ∞ 𝑛 = − ∞ 𝑢 𝑛 ( 𝑟 ) 𝑒 𝑗 𝑛 𝜃 = ∞ 𝑛 = − ∞ 𝑔 2 D 𝑛 ( 𝑟 , 𝑘 ) ∞ 0 𝑠 𝑛 ( 𝑥 , 𝜔 ) 𝐽 𝑛 ( 𝑘 𝑥 ) 𝑥 𝑑 𝑥 𝑒 𝑗 𝑛 𝜃 , ( 9 . 1 5 ) or more compactly as 𝑢 𝑛 ( 𝑟 , 𝜔 ) = 𝑔 2 D 𝑛 ( 𝑟 , 𝑘 ) ∞ 0 𝑠 𝑛 ( 𝑥 , 𝜔 ) 𝐽 𝑛 ( 𝑘 𝑥 ) 𝑥 𝑑 𝑥 . ( 9 . 1 6 ) It is in (9.15) that the interpretation of the evaluation of the 2D Fourier transform becomes apparent, namely through the evaluation of the Bessel function at the (real or complex) wavenumbers. Using the definition of the Hankel transform, (9.16) can be written in the compact and powerful formulation of 𝑢 𝑛 ( 𝑟 , 𝜔 ) = 𝑔 2 D 𝑛 ( 𝑟 , 𝑘 ) ℍ 𝑛 𝑠 𝑛 ( 𝑘 , 𝜔 ) . ( 9 . 1 7 ) If it is also further assumed that the inhomogeneity function is separable in time and space (a fairly general assumption) so that 𝑠 ( ⃗ 𝑟 , 𝜔 ) = 𝑞 ( ⃗ 𝑟 ) 𝜂 ( 𝜔 ) , and ∞ 𝑛 = − ∞ 𝑠 𝑛 ( 𝑟 , 𝜔 ) 𝑒 𝑗 𝑛 𝜃 = 𝜂 ( 𝜔 ) ∞ 𝑛 = − ∞ 𝑞 𝑛 ( 𝑟 ) 𝑒 𝑗 𝑛 𝜃 , ( 9 . 1 8 ) then (9.16) can be further reduced to 𝑢 𝑛 ( 𝑟 , 𝜔 ) = 𝑔 2 D 𝑛 ( 𝑟 , 𝑘 ) 𝜂 ( 𝜔 ) ℍ 𝑛 𝑞 𝑛 ( 𝑘 ) . ( 9 . 1 9 ) It is noted that the value of the measured wavefield at any position 𝑟 is related to the Hankel transform of the heterogeneity function evaluated at the wavenumber of the wavefield in question, with the Green function coefficient acting as the proportionality term.

Equation (9.19) is the 2D equivalent of (8.14). The same comments can be made as those made in the discussion after the 3D version of this problem, as the relationship is identical, save for translating 3D tools such as the spherical harmonic expansions and spherical Hankel transforms into 2D tools involving Fourier series and Hankel transforms.

10. Transfer Function Relationship in 1D

We proceed to develop the equivalent relationship for the 1D case. We have seen from (4.7) that the 1D transfer function in the Fourier domain is given by 𝐺 𝜔 𝑥 = 1 , 𝜔 𝜔 2 𝑥 − 𝑘 2 , ( 1 0 . 1 ) which depends on only a single spatial frequency. Equation (4.2) in 1D then becomes 1 𝑢 ( 𝑟 , 𝜔 ) = ( 2 𝜋 ) ∞ − ∞ 𝑆 𝜔 𝑥 , 𝜔 𝜔 2 𝑥 − 𝑘 2 𝑒 𝑖 𝜔 𝑥 𝑟 𝑑 𝜔 𝑥 . ( 1 0 . 2 ) The temporal frequency 𝜔 is a constant as far as the spatial inverse Fourier transformation is concerned. We are thus interested in calculating integrals of the form 1 𝐼 = 2 𝜋 ∞ 0 𝜙 ( 𝑥 ) 𝑥 2 − 𝑘 2 𝑒 𝑖 𝑥 𝑟 𝑑 𝑥 , ( 1 0 . 3 ) where 𝜙 is an analytic function defined on the positive real line that approaches zero as 𝑥 goes to infinity. The required theorem is given in [36] and the relevant result is presented below.

Theorem 10.1. It is shown in [36] that the following is true 1 2 𝜋 ∞ − ∞ 𝜙 ( 𝜌 ) 𝜌 2 − 𝑘 2 𝑒 𝑖 𝜌 𝑟 ⎧ ⎪ ⎨ ⎪ ⎩ 1 𝑑 𝜌 = 2 𝑖 𝑘 𝜙 ( − 𝑘 ) 𝑒 − 𝑖 𝑘 𝑟 1 , 𝑟 > 0 , 2 𝑖 𝑘 𝜙 ( 𝑘 ) 𝑒 𝑖 𝑘 𝑟 , 𝑟 < 0 . ( 1 0 . 4 ) Here, 𝜙 is an analytic function defined on the positive real line that remains bounded as 𝑥 goes to infinity and 𝑘 is a wavenumber which may be real or complex, chosen such that the imaginary part of 𝑘 is negative. Given the definition of the Fourier transform that is being currently used, the presented result satisfies the Sommerfeld radiation condition, ensuring an outwardly propagating wave.

10.1. 1D Green's Functions via Inverse Fourier Transformation of the Spatial Transfer Functions

The Green's function for the Helmholtz equation is 1D and is given by inverse Fourier transform of the spatial transfer function as 𝑔 1 D 1 ( 𝑟 , 𝑘 ) = 2 𝜋 ∞ − ∞ 𝑒 𝑖 𝜔 𝑥 𝑟 𝜔 2 𝑥 − 𝑘 2 𝑑 𝜔 𝑥 , ( 1 0 . 5 ) and can be evaluated with the help of Theorem 10.1 to give 𝑔 1 D ( ⎧ ⎪ ⎨ ⎪ ⎩ 1 𝑟 , 𝑘 ) = 𝑒 2 𝑖 𝑘 − 𝑖 𝑘 𝑟 1 𝑟 > 0 𝑒 2 𝑖 𝑘 𝑖 𝑘 𝑟 ⎫ ⎪ ⎬ ⎪ ⎭ = 1 𝑟 < 0 𝑒 2 𝑖 𝑘 − 𝑖 𝑘 | 𝑟 | . ( 1 0 . 6 )

10.2. Transfer Function Relationship in 1D

Equation (10.2) can now be evaluated with the help of Theorem 10.1 as well as the result in (10.6) as 1 𝑢 ( 𝑟 , 𝜔 ) = ( 2 𝜋 ) ∞ − ∞ 𝑆 𝜔 𝑥 , 𝜔 𝜔 2 𝑥 − 𝑘 2 𝑒 𝑖 𝜔 𝑥 𝑟 𝑑 𝜔 𝑥 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 2 𝑖 𝑘 𝑆 ( − 𝑘 , 𝜔 ) 𝑒 − 𝑖 𝑘 𝑟 1 𝑟 > 0 2 𝑖 𝑘 𝑆 ( 𝑘 , 𝜔 ) 𝑒 𝑖 𝑘 𝑟 𝑟 < 0 = 𝑔 1 D ( 𝑟 , 𝑘 ) 𝑆 ( − s g n ( 𝑟 ) 𝑘 , 𝜔 ) , ( 1 0 . 7 ) where 𝑔 1 𝐷 ( 𝑟 , 𝑘 ) = ( 1 / 2 𝑖 𝑘 ) 𝑒 − 𝑖 𝑘 | 𝑟 | . This is the 1D version of the transfer function relationship, with the wavefield being directly proportional to the Fourier transform of the object function evaluated at 𝑘 . This is similar to the results for the 2D and 3D cases. For a real wavenumber, a “sphere” in 1D becomes the two points on the real line at ± 𝑘 and the proportionality term is the Green's function for the space in question.

As before, we assume that the inhomogeneity function is separable in time and space (a fairly general assumption) so that 𝑠 ( 𝑟 , 𝜔 ) = 𝑞 ( 𝑟 ) 𝜂 ( 𝜔 ) → 𝑆 ( 𝜔 𝑥 , 𝜔 ) = 𝑄 ( 𝜔 𝑥 ) 𝜂 ( 𝜔 ) , where 𝑄 ( 𝜔 𝑥 ) is the 1D Fourier transform of 𝑞 ( 𝑟 ) . The relevant result now reads 𝑢 ( 𝑟 , 𝜔 ) = 𝑔 1 D ( 𝑟 , 𝑘 ) 𝜂 ( 𝜔 ) 𝑄 ( − s g n ( 𝑟 ) 𝑘 ) . ( 1 0 . 8 ) This result is directly applicable for computational purposes.

11. Summary of Results

The relationships in the previous sections are summarized in Table 1.

tab1
Table 1

12. Applications of the Transfer Function Relationships to the Wave Equation

The transfer function relationships can be used to find time domain Green's functions where the temporal Fourier integral can be easily inverted. In particular, we consider the case where the input to the standard wave equation is a Dirac-delta function at the origin ∇ 2 𝑢 − 1 ⃗ 𝑟 , 𝑡 𝑐 2 𝜕 2 𝜕 𝑡 2 𝑢 ⃗ 𝑟 , 𝑡 = − 𝛿 ⃗ 𝑟 𝛿 ( 𝑡 ) . ( 1 2 . 1 )

This corresponds to the Helmholtz equation above with 𝑘 = 𝜔 / 𝑐 . The temporal and spatial Fourier transform of 𝛿 ( ⃗ 𝑟 ) 𝛿 ( 𝑡 ) is 1, which is spherically symmetric so that the transfer function relationships given above only need the zeroth-order component and simplify considerably. In 3D, the transfer function relationship is given by 𝑢 ( 𝑟 , 𝜔 ) = 𝑔 3 D 0 ( 𝑟 , 𝑘 ) ⋅ 1 = − 𝑖 𝑘 ℎ 0 ( 2 ) ( 𝑘 𝑟 ) . ( 1 2 . 2 ) In 2D, the relationship is given by 𝑢 ( 𝑟 , 𝜔 ) = 𝑔 2 D 0 ( 𝑟 , 𝑘 ) ⋅ 1 = − 𝜋 𝑖 2 𝐻 0 ( 𝑘 𝑟 ) . ( 1 2 . 3 ) In 1D, the relationship is 𝑢 ( 𝑟 , 𝜔 ) = 𝑔 1 D 1 ( 𝑟 , 𝑘 ) ⋅ 1 = 𝑒 2 𝑖 𝑘 − 𝑖 𝑘 | 𝑟 | . ( 1 2 . 4 ) These three relationships can now be inverse Fourier transformed in time.

For the 3D case, we use the fact that ℎ 0 ( 2 ) ( 𝑘 𝑟 ) = 𝑗 0 ( 𝑘 𝑟 ) − 𝑖 𝑦 0 ( 𝑘 𝑟 ) = 𝑖 e x p ( − 𝑖 𝑘 𝑟 ) / 𝑘 𝑟 so that the inverse Fourier transform of (12.2) gives 1 𝑢 ( 𝑟 , 𝑡 ) = 2 𝜋 ∞ − ∞ 1 𝑟 𝑒 − 𝑖 𝜔 𝑟 / 𝑐 𝑒 𝑖 𝜔 𝑡 = 1 𝑑 𝜔 𝑟 𝛿 𝑟 𝑡 − 𝑐 . ( 1 2 . 5 ) The 1D time response can be found from the inverse Fourier transform of  (12.4) 1 𝑢 ( 𝑟 , 𝑡 ) = 2 𝜋 ∞ − ∞ 𝑐 𝑒 2 𝑖 𝜔 − 𝑖 𝜔 | 𝑟 | / 𝑐 𝑒 𝑖 𝜔 𝑡 = 𝑐 𝑑 𝜔 4 2 𝐻 𝑡 − | 𝑟 | 𝑐 = 𝑐 − 1 4 s g n 𝑡 − | 𝑟 | 𝑐 , ( 1 2 . 6 ) where 𝐻 ( 𝑥 ) is the Heaviside unit step function.

The 2D case is a little more complicated but can nevertheless be evaluated in closed form by an inverse Fourier transform of (12.3). To do this, we write the zeroth-order Hankel function in its integral form as [37] 𝐻 0 ( 2 ) ( 𝑥 ) = 𝐽 0 ( 𝑥 ) − 𝑖 𝑌 0 ( 𝑥 ) = 2 𝑖 𝜋 ∞ 1 𝑒 − 𝑖 𝑥 𝜏 √ 𝜏 2 − 1 𝑑 𝜏 . ( 1 2 . 7 ) Hence from (12.7) and (12.3), the inverse Fourier transform of (12.3) gives 1 𝑢 ( 𝑟 , 𝑡 ) = 2 𝜋 ∞ − ∞ ∞ 1 𝑒 − 𝑖 𝜔 𝑟 𝜏 / 𝑐 √ 𝜏 2 − 1 𝑑 𝜏 𝑒 𝑖 𝜔 𝑡 𝑑 𝜔 . ( 1 2 . 8 ) Changing the order of integration gives 𝑢 ( 𝑟 , 𝑡 ) = ∞ 1 1 √ 𝜏 2 1 − 1 2 𝜋 ∞ − ∞ 𝑒 𝑖 𝜔 [ 𝑡 − 𝑟 𝜏 / 𝑐 ] 𝑑 𝜔 𝑑 𝜏 . ( 1 2 . 9 ) But 1 2 𝜋 ∞ − ∞ 𝑒 𝑖 𝜔 [ 𝑡 − 𝑟 𝜏 / 𝑐 ] 𝑑 𝜔 = 𝛿 𝑡 − 𝑟 𝜏 𝑐 , ( 1 2 . 1 0 ) so that (12.9) becomes 𝑢 ( 𝑟 , 𝑡 ) = ∞ 1 1 √ 𝜏 2 𝛿 − 1 𝑡 − 𝑟 𝜏 𝑐 𝑑 𝜏 = ∞ − ∞ 1 √ 𝜏 2 − 1 𝐻 ( 𝜏 − 1 ) 𝛿 𝑡 − 𝑟 𝜏 𝑐 𝑑 𝜏 . ( 1 2 . 1 1 ) Changing variables so that 𝑥 = 𝑟 𝜏 / 𝑐 gives 𝑐 𝑢 ( 𝑟 , 𝑡 ) = 𝑟 ∞ − ∞ 𝐻 𝑐 𝑥 𝑟 1 − 1 √ ( 𝑐 𝑥 / 𝑟 ) 2 𝑟 − 1 𝛿 ( 𝑡 − 𝑥 ) 𝑑 𝑥 = 𝐻 𝑡 − 𝑐 1 √ 𝑡 2 − 𝑟 2 / 𝑐 2 , ( 1 2 . 1 2 ) which yields the desired result in closed form.

13. Application: Inversion Formulas for Tomographic Applications with Real Wavenumbers

In the case of a real wavenumber 𝑘 , an immediate application of the above formulas is for closed-form inversion formulas which are useful for tomographic applications. These inversion formulas are applicable to applications where the wavefield (e.g., acoustic field) is measured and the goal is to reconstruct the source (input) that led to that measured wavefield (output). As we have seen, formulas above give the output (wavefield) in terms of the transform of the input evaluated at 𝑘 . Thus, if the output wavefield is measured, this implies knowledge of the transform of the input. An inverse transform then leads to the input itself.

13.1. Inversion Formula in 3D

Equation (8.14) admits an inversion useful for tomographic applications. It can be written as 𝑖 𝑢 𝑘 𝑙 𝑟 0 , 𝜔 𝑘 𝑠 ℎ 𝑙 ( 2 ) 𝑘 𝑠 𝑟 0 = 𝜂 ( 𝜔 ) ∞ 0 𝜙 𝑘 𝑙 ( 𝑥 ) 𝑗 𝑙 𝑘 𝑠 𝑥 𝑥 2 𝑑 𝑥 , ( 1 3 . 1 ) where 𝑥 has been used as a dummy integration variable in order to avoid possible confusion and 𝑟 0 has been used as the radial variable and indicates the position where a measurement of the wavefield is made. Multiplying both sides by 𝑗 𝑙 ( 𝑘 𝑠 𝑟 ) 𝑘 2 𝑠 and integrating over 𝑘 𝑠 gives ∞ 0 𝑖 𝑢 𝑘 𝑙 𝑟 0 , 𝜔 𝑘 𝑠 ℎ 𝑙 ( 2 ) 𝑘 𝑠 𝑟 0 𝑗 𝜂 ( 𝜔 ) 𝑙 𝑘 𝑠 𝑟 𝑘 2 𝑠 𝑑 𝑘 𝑠 = ∞ 0 ∞ 0 𝜙 𝑘 𝑙 ( 𝑥 ) 𝑗 𝑙 𝑘 𝑠 𝑥 𝑥 2 𝑑 𝑥 𝑗 𝑙 𝑘 𝑠 𝑟 𝑘 2 𝑠 𝑑 𝑘 𝑠 = ∞ 0 ∞ 0 𝜙 𝑘 𝑙 ( 𝑥 ) 𝑗 𝑙 𝑘 𝑠 𝑥 𝑗 𝑙 𝑘 𝑠 𝑟 𝑘 2 𝑠 𝑑 𝑘 𝑠 𝑥 2 𝑑 𝑥 . ( 1 3 . 2 )

Using the orthogonality of the spherical Bessel functions, it follows that (13.2) becomes ∞ 0 𝑖 𝑢 𝑘 𝑙 𝑟 0 , 𝜔 𝑘 𝑠 ℎ 𝑙 ( 2 ) 𝑘 𝑠 𝑟 0 𝑗 𝜂 ( 𝜔 ) 𝑙 𝑘 𝑠 𝑟 𝑘 2 𝑠 𝑑 𝑘 𝑠 = ∞ 0 𝜙 𝑘 𝑙 𝜋 ( 𝑥 ) 2 𝑥 2 𝛿 ( 𝑥 − 𝑟 ) 𝑥 2 𝜋 𝑑 𝑥 = 2 𝜙 𝑘 𝑙 ( 𝑟 ) . ( 1 3 . 3 ) Recalling that 𝑘 𝑠 = 𝜔 / 𝑐 𝑠 , the inversion formula for the spatial source becomes 𝜙 𝑘 𝑙 ( 𝑟 ) = 2 𝑖 𝜋 𝑐 2 𝑠 ∞ 0 𝑢 𝑘 𝑙 𝑟 0 , 𝜔 ℎ 𝑙 ( 2 ) 𝜔 𝑟 0 / 𝑐 𝑠 𝑗 𝜂 ( 𝜔 ) 𝑙 𝜔 𝑟 𝑐 𝑠 𝜔 𝑑 𝜔 . ( 1 3 . 4 ) If the measured wavefield is also separable in space and time so that 𝑢 𝑘 𝑙 ( 𝑟 0 , 𝜔 ) = 𝜒 𝑘 𝑙 ( 𝑟 0 ) 𝑇 ( 𝜔 ) , then (13.4) can be further simplified to 𝜙 𝑘 𝑙 ( 𝑟 ) = 2 𝑖 𝜒 𝑘 𝑙 𝑟 0 𝜋 𝑐 2 𝑠 ∞ 0 𝑇 ( 𝜔 ) ℎ 𝑙 ( 2 ) 𝜔 𝑟 0 / 𝑐 𝑠 𝑗 𝜂 ( 𝜔 ) 𝑙 𝜔 𝑟 𝑐 𝑠 𝜔 𝑑 𝜔 . ( 1 3 . 5 ) Equation (13.5) gives a simple form for finding the spatial source function at any position 𝑟 once having measured the temporal response of the wavefield at a position 𝑟 0 . However, since the spherical harmonic transform 𝑢 𝑘 𝑙 ( 𝑟 0 , 𝜔 ) of the wavefield is required, this implies that sufficient angular information about the wavefield 𝑢 must be obtained in order to permit this spherical harmonic calculation.

13.2. Inversion Formula in 2D

Equation (9.19) admits an inversion leading to the input function. We assume that the source function is separable in time and space (a fairly general assumption) so that 𝑠 ( ⃗ 𝑟 , 𝜔 ) = 𝜙 ( ⃗ 𝑟 ) 𝜂 ( 𝜔 ) . In this case, then 𝑠 𝑘 𝑙 ( 𝑟 , 𝜔 ) = 𝜙 𝑘 𝑙 ( 𝑟 ) 𝜂 ( 𝜔 ) and (8.12) becomes 2 𝑖 𝑢 𝑛 𝑟 0 , 𝜔 𝜋 𝐻 𝑛 ( 2 ) 𝑘 𝑠 𝑟 0 = 𝜂 ( 𝜔 ) ∞ 0 𝜙 𝑛 ( 𝑥 ) 𝐽 𝑛 𝑘 𝑠 𝑥 𝑥 𝑑 𝑥 , ( 1 3 . 6 ) where 𝑥 has been used as a dummy integration variable in order to avoid possible confusion and 𝑟 0 indicates the radial position where the wavefield is measured. Multiplying both sides of (13.6) by 𝐽 𝑛 ( 𝑘 𝑠 𝑟 ) 𝑘 𝑠 and integrating over 𝑘 𝑠 gives ∞ 0 2 𝑖 𝑢 𝑛 𝑟 0 , 𝜔 𝜋 𝐻 𝑛 ( 2 ) 𝑘 𝑠 𝑟 0 𝐽 𝜂 ( 𝜔 ) 𝑛 𝑘 𝑠 𝑟 𝑘 𝑠 𝑑 𝑘 𝑠 = ∞ 0 ∞ 0 𝜙 𝑛 ( 𝑥 ) 𝐽 𝑛 𝑘 𝑠 𝑥 𝑥 𝑑 𝑥 𝐽 𝑛 𝑘 𝑠 𝑟 𝑘 𝑠 𝑑 𝑘 𝑠 , = ∞ 0 𝜙 𝑛 ( 𝑥 ) ∞ 0 𝐽 𝑛 𝑘 𝑠 𝑥 𝐽 𝑛 𝑘 𝑠 𝑟 𝑘 𝑠 𝑑 𝑘 𝑠 𝑥 𝑑 𝑥 . ( 1 3 . 7 ) Using the orthogonality of the spherical Bessel functions, it follows that (13.7) becomes ∞ 0 2 𝑖 𝑢 𝑛 𝑟 0 , 𝜔 𝜋 𝐻 𝑛 ( 2 ) 𝑘 𝑠 𝑟 0 𝐽 𝜂 ( 𝜔 ) 𝑛 𝑘 𝑠 𝑟 𝑘 𝑠 𝑑 𝑘 𝑠 = ∞ 0 𝜙 𝑛 1 ( 𝑥 ) 𝑥 𝛿 ( 𝑥 − 𝑟 ) 𝑥 𝑑 𝑥 = 𝜙 𝑛 ( 𝑟 ) . ( 1 3 . 8 ) Recalling that 𝑘 𝑠 = 𝜔 / 𝑐 𝑠 , the inversion formula for the inhomogeneity becomes 𝜙 𝑛 ( 𝑟 ) = 2 𝑖 𝜋 𝑐 2 𝑠 ∞ 0 𝑢 𝑛 𝑟 0 , 𝜔 𝐻 𝑛 ( 2 ) 𝜔 𝑟 0 / 𝑐 𝑠 𝐽 𝜂 ( 𝜔 ) 𝑛 𝜔 𝑟 𝑐 𝑠 𝜔 𝑑 𝜔 . ( 1 3 . 9 ) If the measured wavefield is also separable in space and time so that 𝑢 𝑛 ( 𝑟 0 , 𝜔 ) = 𝜒 𝑛 ( 𝑟 0 ) 𝑇 ( 𝜔 ) , then (13.4) can be further simplified to 𝜙 𝑛 ( 𝑟 ) = 2 𝑖 𝜒 𝑛 𝑟 0 𝜋 𝑐 2 𝑠 ∞ 0 𝑇 ( 𝜔 ) 𝐻 𝑛 ( 2 ) 𝜔 𝑟 0 / 𝑐 𝑠 𝐽 𝜂 ( 𝜔 ) 𝑛 𝜔 𝑟 𝑐 𝑠 𝜔 𝑑 𝜔 . ( 1 3 . 1 0 ) Equation (13.10) gives a simple form for finding the spatial source function at any position 𝑟 once having measured the wavefield response at a radial position 𝑟 0 . Comparing (13.10) and (13.5), it is noted that they have identical forms with the exception of the replacement of Hankel and Bessel functions for the 2D case with spherical Hankel and Bessel functions for the 3D case.

13.3. Inversion Formula in 1D

We can also further assume that the wavefield is also separable in time and space so that 𝑢 ( 𝑟 , 𝜔 ) = 𝜒 ( 𝑟 ) 𝑇 ( 𝜔 ) , which leads to to 𝜙 − 𝑘 𝑠 = 𝜒 ( 𝑟 ) 𝑇 ( 𝜔 ) 𝜂 ( 𝜔 ) 2 𝑖 𝑘 𝑠 𝑒 𝑖 𝑘 𝑠 𝑟 = ∞ − ∞ 𝜙 ( 𝑦 ) 𝑒 𝑖 𝑦 𝑘 𝑠 𝜙 𝑘 𝑑 𝑦 𝑟 > 0 , 𝑠 = 𝜒 ( 𝑟 ) 𝑇 ( 𝜔 ) 𝜂 ( 𝜔 ) 2 𝑖 𝑘 𝑠 𝑒 − 𝑖 𝑘 𝑠 𝑟 = ∞ − ∞ 𝜙 ( 𝑦 ) 𝑒 − 𝑖 𝑦 𝑘 𝑠 𝑑 𝑦 𝑟 < 0 . ( 1 3 . 1 1 ) Using the orthogonality of the Fourier kernel, both sides are multiplied by 𝑒 ± 𝑖 𝑥 𝑘 𝑠 and integrated over all 𝑘 𝑠 to give 𝜒 ( 𝑟 ) ∞ − ∞ 𝑇 ( 𝜔 ) 𝜂 ( 𝜔 ) 2 𝑖 𝑘 𝑠 𝑒 𝑖 𝑘 𝑠 𝑟 𝑒 − 𝑖 𝑥 𝑘 𝑠 𝑑 𝑘 𝑠 = ∞ − ∞ 𝜙 ( 𝑦 ) ∞ − ∞ 𝑒 𝑖 𝑦 𝑘 𝑠 𝑒 − 𝑖 𝑥 𝑘 𝑠 𝑑 𝑘 𝑠 𝑑 𝑦 = 2 𝜋 ∞ − ∞ 𝜙 ( 𝑦 ) 𝛿 ( 𝑦 − 𝑥 ) 𝑑 𝑦 𝑟 > 0 ⟹ 𝜙 ( 𝑥 ) = 𝑖 𝜒 ( 𝑟 ) 𝜋 𝑐 2 𝑠 ∞ − ∞ 𝑇 ( 𝜔 ) 𝑒 𝜂 ( 𝜔 ) 𝑖 𝜔 ( ( 𝑟 − 𝑥 ) / 𝑐 𝑠 ) 𝜔 𝑑 𝜔 𝑟 > 0 , ( 1 3 . 1 2 ) and similarly 𝜒 ( 𝑟 ) ∞ − ∞ 𝑇 ( 𝜔 ) 𝜂 ( 𝜔 ) 2 𝑖 𝑘 𝑠 𝑒 − 𝑖 𝑘 𝑠 𝑟 𝑒 𝑖 𝑥 𝑘 𝑠 𝑑 𝑘 𝑠 = ∞ − ∞ 𝜙 ( 𝑦 ) ∞ − ∞ 𝑒 − 𝑖 𝑦 𝑘 𝑠 𝑒 𝑖 𝑥 𝑘 𝑠 𝑑 𝑘 𝑠 𝑑 𝑦 = 2 𝜋 ∞ − ∞ 𝜙 ( 𝑦 ) 𝛿 ( 𝑦 − 𝑥 ) 𝑑 𝑦 𝑟 < 0 ⟹ 𝜙 ( 𝑥 ) = 𝑖 𝜒 ( 𝑟 ) 𝜋 𝑐 2 𝑠 ∞ − ∞ 𝑇 ( 𝜔 ) 𝑒 𝜂 ( 𝜔 ) 𝑖 𝜔 ( ( 𝑥 − 𝑟 ) / 𝑐 𝑠 ) 𝜔 𝑑 𝜔 𝑟 < 0 . ( 1 3 . 1 3 ) Several points need to be made regarding (13.12) and (13.13). First, they both give the source function as a function of position 𝑥 along the real line. The interpretation of the variable 𝑟 is that of the position at which the measurement is made and is considered to be a fixed quantity. Both (13.12) and (13.13) are inverse Fourier transforms of ( 𝑇 ( 𝜔 ) / 𝜂 ( 𝜔 ) ) 𝜔 , evaluated at ( 𝑥 − 𝑟 ) / 𝑐 𝑠 or ( 𝑟 − 𝑥 ) / 𝑐 𝑠 , depending on whether measurements are made in transmission or reflection. Clearly | 𝑥 − 𝑟 | / 𝑐 𝑠 is the time taken for a wave to travel the distance | 𝑥 − 𝑟 | .

The 2D and 3D equivalent to (13.12) and (13.13) are (13.5) and (13.10) which also similarly involve an inverse transform of ( 𝑇 ( 𝜔 ) / 𝜂 ( 𝜔 ) ) 𝜔 but have a different kernel for the integration. The reason for the differences in the nature of the kernels from the 2D/3D cases to the 1D case is that in the 2D and 3D cases the kernels used for the Fourier transforms are the Bessel and spherical Bessel functions, respectively, which are the standing wave solutions in those dimensions. However, Green's functions in those dimensions are the Hankel and spherical Hankel functions which are the travelling wave solutions. We note that the kernels for the inversion formulas of (13.5) and (13.10) involve inverse Hankel and spherical Hankel functions, which represent the inverse of Green's function. In contrast, the travelling wave solutions in 1D are the complex exponential and those are also the kernel of 1D Fourier transform. The ratio of the two complex exponentials (one represents the Fourier kernel, the other represents the inverse of 1D Green's function) finally give another complex exponential which represents the shift | 𝑥 − 𝑟 | / 𝑐 𝑠 .

14. Summary and Conclusions

This paper presents the derivation of spatial transfer function relationships for the Helmholtz equation. The focus has been on deriving forward transfer function relationships so that once given a particular input function, the resulting output wavefield can be calculated at any point in space. These are termed “transfer function” relationships because the relationship between the input and output quantities is one of multiplication and no convolution is involved, even in the spatial domain where normally a convolution would be required. Interestingly, instead of Green's function itself, Green's function coefficients for the space are required. These Green function coefficients form the kernel of the output response. Each element of the kernel is then weighted (multiplicatively) by the relevant Fourier/Hankel/Spherical Hankel transform of the input, evaluated on the 𝑛 -dimensional sphere of radius 𝑘 , where 𝑘 is the wavenumber and 𝑛 is the dimension of the space. Combined together, these form the final output response. This is true for dimensions 𝑛 = 1 , 2, 3. These simple but powerful results serve to cast the entire problem as an input-output problem with a transfer function relating input to output. In this view of the problem, the input (inhomogeneity) and output (resulting wavefield) in space are related by a simple multiplicative transfer function relationship and not via a convolution. Some applications of these results were shown in the manuscript.

References

  1. A. Mandelis, “Theory of photothermal-wave diffraction and interference in condensed media,” Journal of the Optical Society of America A, vol. 6, no. 2, pp. 298–308, 1989. View at Publisher · View at Google Scholar
  2. A. Mandelis, Diffusion-Wave Fields. Mathematical methods and Green Function, Springer-Verlag, New York, NY, USA, 2001. View at Zentralblatt MATH
  3. M. Slaney and A. C. Kak, Principles of Computerized Tomographic Imaging, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 1988. View at Zentralblatt MATH
  4. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Review of Scientific Instruments, vol. 77, no. 4, pp. 041101.1–041101.22, 2006. View at Publisher · View at Google Scholar
  5. A. Mandelis and C. Feng, “Frequency-domain theory of laser infrared photothermal radiometric detection of thermal waves generated by diffuse-photon-density wave fields in turbid media,” Physical Review E—Statistical Nonlinear, and Soft Matter Physics, vol. 65, no. 2, pp. 021909/1–021909/19, 2002.
  6. A. Mandelis, Y. Fan, G. Spirou, and A. Vitkin, “Development of a photothermoacoustic frequency swept system: theory and experiment,” Journal De Physique. IV : JP, vol. 125, pp. 643–647, 2005. View at Publisher · View at Google Scholar
  7. A. Mandelis, R. J. Jeon, S. Telenkov, Y. Fan, and A. Matvienko, “Trends in biothermophotonics and bioacoustophotonics of tissues,” in Proceedings of the International Society for Optics and Photonics (SPIE '05), vol. 5953, pp. 1–15, The International Society for Optical Engineering, Warsaw, Poland, 2005.
  8. C. L. Matson, “A diffraction tomographic model of the forward problem using diffuse photon density waves,” Optics Express, vol. 1, no. 1, pp. 6–11, 1997.
  9. L. Nicolaides and A. Mandelis, “Image-enhanced thermal-wave slice diffraction tomography with numerically simulated reconstructions,” Inverse Problems, vol. 13, no. 5, pp. 1393–1412, 1997. View at Publisher · View at Google Scholar
  10. L. Nicolaides, M. Munidasa, and A. Mandelis, “Thermal-wave infrared radiometric slice diffraction tomography with back-scattering and transmission reconstructions: experimental,” Inverse Problems, vol. 13, no. 5, pp. 1413–1425, 1997.
  11. O. Pade and A. Mandelis, “Computational thermal-wave slice tomography with backpropagation and transmission reconstructions,” Review of Scientific Instruments, vol. 64, no. 12, pp. 3548–3562, 1993. View at Publisher · View at Google Scholar
  12. O. Pade and A. Mandelis, “Thermal-wave slice tomography using wave field reconstruction,” in Proceedings of the 8th International Topical Meeting on Photoacoustic and Photothermal Phenomena, vol. 4, 1994.
  13. S. A. Telenkov, G. Vargas, J. S. Nelson, and T. E. Milner, “Coherent thermal wave imaging of subsurface chromophores in biological materials,” Physics in Medicine and Biology, vol. 47, no. 4, pp. 657–671, 2002. View at Publisher · View at Google Scholar
  14. S. A. Telenkov, J. S. Nelson, and T. E. Milner, “Tomographic reconstruction of tissue chromophores using thermal wave imaging,” Laser Tissue Interaction XIII, vol. 4617, pp. 40–46, 2002. View at Publisher · View at Google Scholar
  15. S. A. Telenkov and A. Mandelis, “Fourier-domain biophotoacoustic subsurface depth selective amplitude and phase imaging of turbid phantoms and biological tissue,” Journal of Biomedical Optics, vol. 11, no. 4, 2006. View at Publisher · View at Google Scholar · View at PubMed
  16. N. Baddour, “Theory and analysis of frequency-domain photoacoustic tomography,” Journal of the Acoustical Society of America, vol. 123, no. 5, pp. 2577–2590, 2008. View at Publisher · View at Google Scholar · View at PubMed
  17. N. Baddour, “A multidimensional transfer function approach to photo-acoustic signal analysis,” Journal of the Franklin Institute, vol. 345, no. 7, pp. 792–818, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  18. N. Baddour, “Fourier diffraction theorem for diffusion-based thermal tomography,” Journal of Physics A, vol. 39, no. 46, pp. 14379–14395, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  19. A. Mandelis, “Theory of photothermal wave diffraction tomography via spatial Laplace spectral decomposition,” Journal of Physics A, vol. 24, no. 11, pp. 2485–2505, 1991. View at Publisher · View at Google Scholar
  20. A. Mandelis, “Green's functions in thermal-wave physics: cartesian coordinate representations,” Journal of Applied Physics, vol. 78, no. 2, pp. 647–655, 1995. View at Publisher · View at Google Scholar
  21. M. Puschel and J. M. F. Moura, “Algebraic signal processing theory: cooley-Tukey type algorithms for DCTs and DSTs,” IEEE Transactions on Signal Processing, vol. 56, no. 4, pp. 1502–1521, 2008. View at Publisher · View at Google Scholar
  22. M. Puschel, “Algebraic signal processing theory: An overview,” in Proceedings of the IEEE 12th Digital Signal Processing Workshop and 4th IEEE Signal Processing Education Workshop, pp. 386–391, Moose, Wyo, USA, 2006. View at Publisher · View at Google Scholar
  23. M. Puschel and J. Moura, “Algebraic signal processing theory: foundation and 1-D time,” IEEE Transactions on Signal Processing, vol. 56, no. 8, pp. 3572–3585, 2008. View at Publisher · View at Google Scholar
  24. M. Puschel and J. M. F. Moura, “Algebraic signal processing theory: 1-D space,” IEEE Transactions on Signal Processing, vol. 56, no. 8, pp. 3586–3599, 2008. View at Publisher · View at Google Scholar
  25. M. Puschel and J. M. F. Moura, “The algebraic structure in signal processing: time and space,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP '06), vol. 5, 2006. View at Publisher · View at Google Scholar
  26. M. Puschel and M. Rotteler, “Algebraic signal processing theory: 2-D spatial hexagonal lattice,” IEEE Transactions on Image Processing, vol. 16, no. 6, pp. 1506–1521, 2007. View at Publisher · View at Google Scholar
  27. M. Puschel and M. Rotteler, “Algebraic signal processing theory: cooley-Tukey type algorithms on the 2-D hexagonal spatial lattice,” Applicable Algebra in Engineering, Communication and Computing, vol. 19, no. 3, pp. 259–292, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  28. N. Baddour, “Operational and convolution properties of two-dimensional Fourier transforms in polar coordinates,” Journal of the Optical Society of America A, vol. 26, no. 8, pp. 1767–1778, 2009. View at Publisher · View at Google Scholar
  29. N. Baddour, “Operational and convolution properties of three-dimensional Fourier transforms in spherical polar coordinates,” Journal of the Optical Society of America, vol. 27, no. 10, pp. 2144–2155, 2010.
  30. N. Baddour, “Theory and analysis of frequency-domain photoacoustic tomography,” Journal of the Acoustical Society of America, vol. 123, no. 5, pp. 2577–2590, 2008.
  31. G. S. Chirikjian and A. B. Kyatkin, Engineering Applications of Noncommutative Harmonic Analysis, CRC Press, Boca Raton, Fla, USA, 2000.
  32. G. S. Chirikjian and A. Kyatkin, Engineering Applications of Noncommutative Harmonic Analysis: With Emphasis on Rotation and Motion Groups, Academic Press, New York, NY, USA, 2001.
  33. R. Piessens, “The Hankel transform,” in The Transforms and Applications Handbook, vol. 2, pp. 9.1–9.30, CRC Press, Boca Raton, Fla, USA, 2000.
  34. “Spherical harmonics—wikipedia, the free encyclopedia,” http://en.wikipedia.org/wiki/Spherical_harmonics, Accessed: 16-Nov-2009.
  35. J. R. Driscoll and D. M. Healy, “Computing Fourier transforms and convolutions on the 2-sphere,” Advances in Applied Mathematics, vol. 15, no. 2, pp. 202–250, 1994. View at Publisher · View at Google Scholar · View at MathSciNet
  36. N. Baddour, “Multidimensional wave field signal theory: mathematical foundations,” AIP Advances, vol. 1, no. 2, p. 022120, 2011. View at Publisher · View at Google Scholar
  37. M. Abramowitz and I. Stegun, Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover Publications, New York, NY, USA, 1964.