Biomedical Engineering Reference
In-Depth Information
where
k
= ω/
c
is the wavenumber for a monochromatic excita-
tion with radian frequency ω, the spatial frequency components
corresponding to propagation in the
x
,
y
, and
z
directions are
k
x
,
k
y
, and
k
z
, respectively, and
kkkk
evaluated and multiplied by
H
p
(
k
x
,
k
y
,Δ
z
) from Equation 6.28,
and then an inverse 2D FFT is applied to the result to obtain the
output pressure plane
p
(
x ,y,z
) at
z
=
z
0
+ Δ
z
.
In pressure calculations with the angular spectrum approach,
each transfer function is associated with a specific error distri-
bution. For example, angular spectrum calculations that eval-
uate the Green's function
h
u
(
x ,y,z
) in the spatial domain and
then numerically calculate the 2D FFT of the result to obtain
H
u
(
k
x
,
k
y
,
z
) are accurate in the central portion of the computa-
tional grid; however, the errors are large near the edges of the
computed field due to wraparound errors that occur in the spa-
tial domain (Zeng and McGough 2008). Computing
H
u
(
k
x
,
k
y
,
z
)
directly in the spatial frequency domain with Equation 6.26
avoids the problems with wraparound errors near the edge, but
when the analytical expression
H
u
(
k
x
,
k
y
,
z
) is sampled in the
spatial frequency domain, the spectrum is inherently aliased
due to the rapid oscillations that occur near
kkk
x
2 2 2 2
++=
. In Equation 6.26,
the term in the numerator is purely real and exponentially
decaying for combinations of spatial frequencies outside of the
circle with radius
k
, and the numerator is complex and oscilla-
tory for combinations of spatial frequencies inside of the circle
with radius
k
. Thus, spatial frequencies inside this circle corre-
spond to propagating waves, and spatial frequencies outside this
circle correspond to evanescent waves.
For certain transducer geometries, the Fourier transform
U
(
k
x
,
k
y
,0) of the normal particle velocity
u
(
x ,y
,0) can be evalu-
ated analytically, but in general,
U
(
k
x
,
k
y
,0) is obtained from
the 2D FFT of
u
(
x ,y
,0). The expression for
H
u
(
k
x
,
k
y
,
z
) can also
be obtained from the 2D FFT of the Green's function
h
u
(
x ,y,z
).
Ater
U
(
k
x
,
k
y
,0) and
H
u
(
k
x
,
k
y
,
z
) are computed, the product is
evaluated, and the pressure is then calculated with an inverse 2D
FFT. The result describes the diffracted pressure within a single
plane. If additional planes are required (for example, to model
ultrasound propagation in a large 3D volume), then
H
u
(
k
x
,
k
y
,
z
)
is calculated for each additional plane, and then the inverse 2D
FFT is evaluated for the product of
U
(
k
x
,
k
y
,0) and
H
u
(
k
x
,
k
y
,
z
) in
as many planes as needed.
The angular spectrum approach also defines a transfer func-
tion in the spatial frequency domain that relates the Fourier
transform of the pressure in one plane to the Fourier transform
of the pressure in another parallel plane. This transfer function
is designated by
H
p
(
k
x
,
k
y
,Δ
z
), where the subscript
p
indicates that
the transfer function input is the Fourier transform of the pres-
sure (as opposed to the Fourier transform of the normal particle
velocity for the transfer function
H
u
),
k
x
and
k
y
are the spatial fre-
quency components defined previously, and Δ
z
is the orthogonal
distance between the two parallel planes. The transfer function
H
p
(
k
x
,
k
y
,Δ
z
) describes forward propagation between an input
pressure plane located at
z
=
z
0
and an output pressure plane
located at
z
=
z
0
+ Δ
z
according to
x
y
z
2 2 2
=+
. This
frequency-domain aliasing causes errors in the spatial domain
that appear as unwanted high-frequency signals in the simu-
lated pressure field. Thus, instead of errors that are confined
to the edge for angular spectrum calculations with the Green's
function
h
u
(
x ,y,z
), the errors are distributed across the entire
computed pressure field for angular spectrum calculations with
the analytical expression
H
u
(
k
x
,
k
y
,
z
). These same errors are also
caused by the rapid oscillations near
kkk
x
y
2 2 2
=+
that occur in
the analytical transfer function
H
p
(
k
x
,
k
y
,Δ
z
) in Equation 6.28.
In calculations with the analytical expressions
H
u
(
k
x
,
k
y
,
z
) and
H
p
(
k
x
,
k
y
,Δ
z
), these errors are reduced with a low-pass filter in
k
-space (Christopher and Parker 1991a; Wu et al. 1997). The
errors are also reduced when the medium is sufficiently attenu-
ative or if the aperture is apodized (Zeng and McGough 2008).
Comparisons between angular spectrum simulations evalu-
ated with
h
u
(
x ,y,z
) in Equation 6.23,
H
u
(
k
x
,
k
y
,
z
) in Equation 6.26,
and
H
p
(
k
x
,
k
y
,Δ
z
) in Equation 6.28 indicate that each approach has
distinct advantages and disadvantages. Of these, angular spec-
trum calculations using the analytical expression for
H
u
(
k
x
,
k
y
,
z
)
in Equation 6.26 are the fastest, requiring minimal computa-
tional effort during initialization and the fewest 2D FFT cal-
culations. The other two approaches are more time consuming,
where simulations with
H
p
(
k
x
,
k
y
,Δ
z
) in Equation 6.28 require
additional time to compute
p
0
(
x ,y,z
0
) in the input pressure plane,
and angular spectrum calculations with
h
u
(
x ,y,z
) in Equation
6.23 require an additional 2D FFT in each plane where the pres-
sure is computed. In contrast, the smallest errors are achieved
by the analytical expression for
H
p
(
k
x
,
k
y
,Δ
z
) in Equation 6.28
when the errors due to frequency domain aliasing are eliminated
with a low-pass filter applied in
k
-space (or through attenuation
or apodization) and the input pressures are computed with the
fast near-field method in a sufficiently large, well-sampled plane
located near the transducer or array face (Zeng and McGough
2009). The errors produced by angular spectrum calculations
with normal particle velocity inputs and the Green's function
h
u
(
x ,y,z
) or the transfer function
H
u
(
k
x
,
k
y
,
z
) are reduced with
increased spatial sampling rates; however, this significantly
y
Pk kz Pk kzHkkz
(,
,)
=
(,
,
) (,
,
,
(6.27)
xy
0
xy
0
px y
where
P
0
(
k
x
,
k
y
,
z
0
) is the 2D Fourier transform of the input pres-
sure plane
p
0
(
x ,y,z
0
) located at
z
=
z
0
, and
P
(
k
x
,
k
y
,
z
0
) is the 2D
Fourier transform of the output pressure plane
p
(
x ,y,z
) located at
z
=
z
0
+ Δ
z
. In Equation 6.27, the transfer function
H
p
(
k
x
,
k
y
,Δ
z
) is
defined by the analytical expression
222
−
jz kkk
−
−
2
2
2
xy
e
for
k
+≤
kk
x
y
Hkkz
(,
,
)
=
,
(6.28)
px y
222
−
zk kk
+
−
2
2
2
e
xy
for
k
+>
kk
x
y
and the input pressure plane
p
0
(
x ,y,z
0
) is typically calculated
with the Rayleigh-Sommerfeld integral or the fast near-field
method. Therefore, to compute the forward propagating ultra-
sound field, the 2D FFT of the input pressure plane
p
0
(
x ,y,z
0
) is
