Ultrasound and Fourier optics

Propagation as a linear operation

As discussed in Section 1.2, linear propagation is governed by the linear wave equation. An important consequence of representing propagation of sound as a linear process is that we can apply the superposition principle.

The Huygens-Fresnel principle states that wavefronts can be decomposed into a collection of point sources, each the origin of a spherical, expanding, wave that can be represented as a free space Green's function. This concept underlies the derivation of an important tool, the Fraunhofer Approximation.

The significance of the Fraunhofer approximation

The Fraunhofer approximation (FA) plays a pivotal role in our exploration of ultrasound k-space. In a nutshell, this well-known expression from the optics literature states that the far-field complex amplitude pattern produced by a complex aperture amplitude function is approximately equal to 2-D Fourier transform of that function. Applied to ultrasound, this approximation states that the ultrasound beam's pressure amplitude pattern can be estimated by taking the 2-D Fourier transform of the transducer aperture. Naturally, this approximation is based on several assumptions and requirements that constrain its application to the far field of an unfocused transducer or the focal plane of a focused transducer. While the approximation is based on a single frequency, or monochromatic solution, we can extend its application to broadband cases through superposition. As we shall see, the pattern is parameterized by the frequency of the wave emitted at the aperture.

As long as we don't violate the assumptions made in formulating the FA, this powerful approximation allows us to extend our intuition regarding linear systems to the study of ultrasound beamforming. Restricting the discussion for a moment to the lateral and axial dimensions, the most obvious example of an aperture function is the rectangular aperture of a linear array lying along the lateral coordinate axis, emitting a single frequency of sound. Representing this aperture function as a rect and applying the FA, we can guess that the ultrasound pressure amplitude pattern in the far field of this aperture is a sinc function.

The development of the Fraunhofer approximation

The complete theoretical development of the FA, from which these excerpts are drawn, can be found in Goodman's Introduction to Fourier Optics[12].

Figure 2.1: Coordinate system used in the derivation of the Fraunhofer approximation.

The coordinate system used is this discussion is shown schematically in Figure 2.1. We consider an arbitrary aperture $\Sigma$ in the aperture plane $(x_{1},y_{1})$ emitting monochromatic radiation of wavenumber $k$, spatially weighted with a complex amplitude function $U(x_{1},y_{1})$. The object is to determine the complex field pattern at an observation plane $U(x_{0},y_{0})$. The geometric path length between any two points $P_{0}$ and $P_{1}$ is simply:

$$r_{01}=\sqrt{z^2+(x_{1}-x_{0})^2+(y_{1}-y_{0})^2}$$ (2.1)

Consider the application of Huygens' principle to the aperture, i.e. that this arbitrary, spatially extended source is represented by an infinite number of point sources lying across $\Sigma$, each with a particular amplitude and phase defined by $U(x_{1},y_{1})$. The amplitude $U(P_{0})$ observed at point $P_{0}$ can be determined by adding up the contributions of all of these point sources after taking into account the geometric path lengths from all of these sources to $P_{0}$. This concept is mathematically expressed in a Huygens-Fresnel integral:

$$\begin{split}U(P_{0})&=\iint_{\Sigma}h(P_{0},P_{1})U(P_{1}) d... ...1}}\cos(\overrightarrow{n},\overrightarrow{r}_{01}) \end{split}$$ (2.2)

The function $h(P_{0},P_{1})$ combines a Green's function $\exp(jkr_{01}) / r_{01}$ representing a spherically expanding wave with what is known as an obliquity term, $\cos(\overrightarrow{n},\overrightarrow{r}_{01})$ (the cosine of the angle between $\overrightarrow{n}$ and $\overrightarrow{r}_{01}$). This Huygens-Fresnel integral sums the contributions from an infinite number of sources across the aperture surface $\Sigma$. One can also think of the propagation of waves from these sources over the distance between $U(x_{1},y_{1})$ and $U(x_{0},y_{0})$ as a kind of filter that has an impulse response $h(P_{0},P_{1})$.

Elaborating on Eq. 2.2:

$$\begin{split}U(x_{0},y_{0})&=\iint_{\Sigma}h(x_{0},y_{0};x_{1... ...x_{1},y_{1})&=0 \text{ outside the aperture }\Sigma \end{split}$$ (2.3)

We now introduce the first assumption, known as the small angle approximation:

$$\begin{split}cos(\overrightarrow{n},\overrightarrow{r}_{01})&... ...},y_{1})&\approx \frac{1}{j\lambda z}\exp(jkr_{01}) \end{split}$$ (2.4)

The Fresnel approximation

This is an intermediate step on the way to the Fraunhofer approximation. Given the expression for path length,

$$\begin{split}r_{01}&=\sqrt{z^2+(x_{1}-x_{0})^2+(y_{1}-y_{0})^... ...igg{)}^2+\bigg{(}\frac{(y_{1}-y_{0})}{z}\bigg{)}^2} \end{split}$$ (2.5)

we introduce a truncated binomial expansion:

$$\begin{split}\sqrt{1+b}&=1 + \frac{1}{2}b - \frac{1}{8}b^2 + ... ...2}\bigg{(}\frac{(y_{1}-y_{0})}{z}\bigg{)}^2\Bigg{]} \end{split}$$ (2.6)

Therefore:

$$h(x_{0},y_{0};x_{1},y_{1})\approx \frac{\exp(jkz)}{j\lambda z} \exp\Bigg{[} \frac{jk}{2z} \bigg{[}(x_{1}-x_{0})^2+(y_{1}-y_{0})^2\bigg{]}\Bigg{]}$$ (2.7)

The Fraunhofer approximation

After rearranging terms:

$$\begin{split}h(x_{0},y_{0};x_{1},y_{1})\approx \frac{\exp(jkz... ...{z} \bigg{[}(x_{0}x_{1}+y_{0}y_{1})\bigg{]}\Bigg{]} \end{split}$$ (2.8)

The Fraunhofer assumption is that the distance between the aperture plane and the observation plane is much greater than the dimensions of the aperture function, i.e.:

$$z»\frac{k(x_{1}^2+y_{1}^2)_{\text{max}}}{2}$$ (2.9)

The distance at which this approximation becomes reasonable is known as the Fraunhofer distance. Thus this quadratic phase term in Eq. 2.8 is assumed to be $\approx 1$. Finally:

$$U(x_{0},y_{0}) \approx \frac {\exp(jkz)\exp\Big{[}\frac{jk}{2z}(x... ...xp\bigg{[} -\frac{j2\pi}{\lambda z}(x_{0}x_{1}+y_{0}y_{1})\bigg{]}dx_{1} dy_{1}$$ (2.10)

Aside from the quadratic phase term outside the integral, this result is the Fourier transform of $U(x_{1},y_{1})$ evaluated at spatial frequencies

$$f_{x}=\frac{x_{0}}{\lambda z} f_{y}=\frac{y_{0}}{\lambda z}.$$ (2.11)

The two-dimensional Fourier transform

The Fraunhofer approximation is a powerful tool for understanding ultrasound beamforming. Its intuitive application requires knowledge of 2-D Fourier Fourier transforms and transform properties. These are seen to be extensions of the Fourier transforms and properties in 1-D, with the addition of the issue of separability.

The analytic form of the 2-D Fourier transform

Following Bracewell's The Fourier Transform and Its Applications[13], the relationships between a 2-D function $f(x,y)$ and its 2-D Fourier transform are:

$$\begin{split}F(u,v) &= \iint^{\infty}_{-\infty}f(x,y) exp\big... ...infty}F(u,v) exp\big{(}j 2\pi(ux+vy)\big{)} du dv\ \end{split}$$ (2.12)

The convolution integral is another extension to 2-D of a familiar relationship:

$$f(x,y) \ast g(x,y) = \iint^{\infty}_{-\infty}f(\eta,\nu) g(x-\eta,y-\nu) d\eta d\nu$$ (2.13)

A 2-D function is separable if it can be described as the product of two 1-D functions along orthogonal dimensions:

$$\begin{split}f(x,y) &= f(x)f(y)\ F(u,v) &= F(u)F(v)\ \end{split}$$ (2.14)

Thus if a function is separable in the space domain, it is also separable in the frequency domain, and vice versa. Identifying a function as separable can simplify its analysis by reducing its 2-D transform to a product of 1-D transforms. Identifying a function as separable can simplify its analysis by reducing its 2-D transform to a product of 1-D transforms.

One class of functions that are not separable are those having circular symmetry. The 2-D Fourier transforms of these functions are also circularly symmetric:

$$\begin{split}\text{If  } f(x,y) &= f(r)\text{, where } r^2=x^... ...then  } F(u,v) &= F(q)\text{, where } q^2=u^2+v^2\ \end{split}$$ (2.15)

A series of Fourier transform pairs are shown in Figures 2.2 through 2.10 using Bracewell's notation for the delta ($\delta$), comb (III), rectangle ($\Pi$), and triangle ($\Lambda$) functions.

Figure 2.2:

Figure 2.3:

Figure 2.4:

Figure 2.5:

Figure 2.6:

Figure 2.7:

Figure 2.8:

Figure 2.9:

Figure 2.10:

Important properties of the 2-D Fourier transform

$$\begin{split}\text{Scaling:}&  f(ax,by)  \Longleftrightarrow ... ...  \int^{\infty}_{-\infty}F(u,v) G^{\ast}(u,v) du dv \end{split}$$ (2.16)

Transforms and properties in the context of a typical aperture

The basis set of the Fourier transform is the complex exponential of general form $A\exp(j2\pi f x + \theta)$. In systems with circular symmetry, such as most optical systems constructed with circular lenses, a different basis set and transform are usually applied, i.e. Bessel functions and the Hankel transform. These allow a circularly symmetric 2-D problem to be solved as a 1-D function of radius $r$. In ultrasound, functions with this symmetry are encountered in work with circular pistons. In this text, we restrict our application of the Fraunhofer approximation to the Fourier transform.

The aperture of a typical 1-D linear array ultrasound transducer is rectangular with electronic focusing applied along the longer dimension (lateral) and an acoustic lens providing a fixed focus in the second dimension (elevation). With no apodization (amplitude weighting of the elements) this aperture function can be represented as a separable product of two rectangle functions, $\Pi(x)\Pi(y)$. As seen in Figure 2.9, the far-field amplitude pattern of this aperture on transmit or receive is approximated by sinc$(u)$sinc$(v)$.

When this aperture is used on both transmit and receive, as is almost always the case, the corresponding function is sinc$^2(u)$sinc$^2(v)$, seen in Figure 2.10. In both cases, the far-field pattern is approximated by these Fourier transforms evaluated at spatial frequencies $f_{u}=\lambda z / x_{o}$ and $f_{v}=\lambda z / y_{o}$. Intuitively, the transmit beam pattern shows where sound is sent into the field, while the receive beam pattern shows the spatial sensitivity of the transducer to returning echoes.

To achieve electronic focusing, one or more dimensions of an aperture is typically divided into discrete elements that are themselves small rectangles. The subdivided dimension is more accurately modeled as a small rectangle convolved with a comb function the length of the array. Thus an aperture of size $D$ composed of elements of size $d$ is modeled as $\Pi(D) [\Pi(d)\ast$III$(x/d)]$, where III$(x)=\Sigma^{\infty}_{-\infty}\delta(x-n)$. The transform of this function is a narrow sinc function convolved with a comb, all multiplied by a broad sinc function:

$$\Pi(D) [\Pi(d)\ast (x/d)]  \Longleftrightarrow (u/d) [ (u/D)\ast (d u)]$$ (2.17)

This is shown schematically in Figure 2.11.

Figure: (Left) An aperture of size $D$ composed of elements of size $d$ is modeled as $\Pi(D) [\Pi(d)\ast$III$(x/d)]$. The Fourier transform of this pattern (right) is sinc$(u/d) [$sinc$(u/D)\ast$III$(d u)]$. The subdivision of the array into elements produces grating lobes.

The transforms of the individual element and of the entire array have a number of implications. First, each individual element has its own far-field amplitude pattern. This pattern in part determines the angular sensitivity of the element, and by extension, of the entire array. This shows up as the broad sinc pattern in Figure 2.11. Also, the periodicity created by subdividing the array introduces a periodicity in the beam pattern. The desired main lobe of the beam pattern is now duplicated off axis as a series of grating lobes. These can be a source of image artifacts, as they literally correspond to regions of unintended off-axis emission and sensitivity. The magnitude of these grating lobes is an important design parameter. When the element spacing, or pitch, is reduced to $\lambda$, the side lobes are located at $90^\circ$ off the main lobe. With a pitch of $\lambda/2$ , the grating lobes are practically eliminated. This approach is necessary when the beam is to be aggressively steered off axis, such as in a cardiac probe. The array's angular sensitivity and the distance of the grating lobes off-axis are both dependent on the element size, so it is common to increase the element size beyond $\lambda/2$ and rely on the reduced angular sensitivity of the elements to reduce the impact of the grating lobes. This design approach is appropriate for linear arrays that will not be used with aggressive beam steering.

Propagation as a spatial transfer function

Under most circumstances in medical imaging, propagation is a linear operation. It has a transfer function associated with it, which is in turn the Fourier transform of the impulse response. (Recently non-linear imaging has also been introduced in clinical scanners, but this modality is beyond the scope of this discussion.)

The free-space Green's function is an impulse function. We can use it to represent a point source, and therefore its Fourier transform is the transfer function of propagation, $H(f_x,F_y) = \exp(jkz)\exp(-j\pi \lambda z(f_{x}^2+f_{y}^2))$. Note that only a change in phase occurs in this operation. There is no change in amplitude, i.e. a pattern at the Fraunhofer distance is the Fourier transform of the aperture, and its Fourier transform gives the aperture distribution back again.

Remember that for an unfocused transducer the FA is only valid beyond the Fraunhofer distance, $z» k(x_1^2+y_1^2)_{max}$, such that the greatest off-axis point in the aperture function satisfies the small angle approximation. For a focused transducer, the FA is valid in the focal plane. This is because at the focal point, the spherically converging waves emanating from the focused aperture are indistinguishable from plane waves modulated by a much more distant, unfocused aperture. Focusing brings the far-field pattern to the focal plane.

To summarize, k-space analysis is based on representing an optical or ultrasound system as a linear system in complex amplitude. Propagation itself can be represented as a transfer function. Under the conditions of linear propagation, only the phase and amplitude of frequency components are modified.