1 Introduction
There is an extensive literature on microphone and speakerarrays for
audio measurement and reconstruction [Pulkki, 2017, Ahrens, 2012]. Let denote the
number of microphones, and the number of speakers. When ,
we have monaural recording and playback, while describes
stereo, etc.
There are many approaches to making larger microphone and/or
speakerarrays ( and/or greater than 2). Since only a small
number of speakers is affordable in typical practice, we are normally
very concerned with human perception of spatial sound
[Blauert, 1997], informing stereophonic, quadraphonic, and more generally
ambisonic sound systems [Cooper, 1972, Gerzon, 1974, 1985].
Ambisonics extends stereo and quad with an expansion of the soundfield
in terms of spherical harmonic basis functions centered on
one listening point.^{1}^{1}1A nice overview of ambisonics
references appears on the
Web:
http://www.york.ac.uk/inst/mustech/3d_audio/gerzonrf.htm
Such systems must deal with the psychoacoustics of direction and timbre
perception in various frequency ranges and for various geometries.
Given a very large number of microphones and speakers, it is possible to approximate complete reconstruction of the soundfield in a given space. The best known approach to this problem is Wave Field Synthesis (WFS) [Berkhout et al., 1993],^{2}^{2}2https://en.wikipedia.org/wiki/Wave_field_synthesis also called “acoustic holography,” or “holophony” [Berkhout, 1988]. WFS reproduces (or synthesizes) a recorded soundfield physically, so that psychoacoustic questions can in principle be avoided.^{3}^{3}3Ambisonics can also approach physical completeness of the soundfield, as the order (number of spherical harmonics) increases, but most practical systems are fairly low order. (The highest order used at CCRMA is presently seven.) Firstorder ambisonics is essentially stereo (representable as a monopole plus one leftright dipole) augmented to include frontback and topdown dipoles. However, for best results at minimum expense, psychoacoustic considerations remain important.
The basic idea of an “acoustic curtain” for reconstructing soundfields in principle was described by Harvey Fletcher [1934], and at that time, two or three speakers was considered an adequate psychoacoustic approximation [Steinberg and Snow, 1934]. Generating wave propagation from spherical waves (“secondary sources”) emitted along the wavefront is the essence of Huygens’ Principle (1690).^{4}^{4}4https://www.britannica.com/science/Huygensprinciple Both Huygens and Fletcher called for a continuum of wavefront samples. The theory of bandlimited sampling was introduced by Nyquist [1928], which, together with basic wave theory, can be considered the basis of this paper.
Deriving WFS begins with the KirchhoffHelmholtz integral (or in simplified form from the Rayleigh integral), which expresses any sourcefree acoustic field as a sum of contributions—called secondary sources—from the boundary of any enclosing surface [Firtha, 2018, Pierce, 1989, Berkhout et al., 1993, Ahrens, 2012]. The same basic approach is used by the well known Boundary Element Method (BEM) for numerically computing a wavefield from its values along a boundary surface [Kirkup, 2007]. The secondary sources in WFS aim to reconstruct (in the listening zone) the same soundfield produced by the original (primary) sources “on stage.” In practice, the secondary sources are simplified from an enclosing sphere down to (typically) a ring of loudspeakers in a line array around the listening space (which should not be reverberant and ideally anechoic—a major goal of Berkhout’s WFS formulation was to include the reverberant as well as the direct soundfield). There are many variations on the details of deriving a practical WFS system, and some of them get close to the samplingbased pointofview taken here. However, there does not appear to be a WFS paper which formulates the problem as basic soundfield sampling and reconstructionfromsamples problem (spatial analogtodigital and digitaltoanalog conversion). As a result, differences in final implementation do emerge, as will be brought out below.
FarField WFS (FFWFS) is the limiting form of WFS in which the sources are many wavelengths away from the recording mics and listening audience. By adopting this simplifying assumption, which is not restrictive in many applications, we can derive FFWFS very easily and clearly from sampling theory, and this is where we begin below.
Overview
This paper defines several speakerarray systems based on a samplingtheory approach:

PlanewaveBased Angle Panning (PBAP) consists of at least one linear speaker array implementing an approximation to FarField Wave Field Synthesis (FFWFS).

Polygonal PBAP (PPBAP) uses multiple PBAP arrays to define an interior acoustic space which is preferably polygonal in the 2D case. Thus, multiple line arrays are combined to form a polygonal listening area, such as a square, rectangle, octagon, etc.

Truncated Polygonal PBAP is a variation on PPBAP which trims each line array down to the polygon side it includes. In the limit of many polygon sides, a circular array is obtained, giving the simplest form of VectorBased AmplitudePanning (VBAP) [Pulkki, 1997, 2001] in which only one speaker (the one closest to the desired angle) emits each source.

Huygens Array (HA) is a generalized approach allowing sources and speakers/microphones to be placed anywhere on the listenerside of the source(s). In this case, samplingbased insights provide guidance for what to expect, and how to avoid artifacts by maintaining a valid (nonuniform) sampling grid.

Various straightforward extensions to 3D are discussed.

Multiband arrays and integrations with VBAP, stereo, and subwoofers are discussed.
When conditions for samplebased reconstruction are not met, we transition as gracefully as possible to some form of Vector Based Amplitude Panning (VBAP) [Pulkki, 2001].
A benefit of extending VBAP to PBAP is providing a larger “sweet spot,” since the quasi spherical waves emanating from VBAP speakers become upgraded to approximations of sampled plane waves in PBAP, and plane waves are the same for all listeners to within a delay. In other words, it is well known that point sources, such as ideal speakers used for ordinary stereo, are sensitive to listener proximity, while the source mix remains constant in soundfields composed of a plane wave from each farfield source.
Background
This paper evolved from 2011 to now as a backburner interest by the author who has never published previously in the field of spatial audio, and who has limited experience in the field, but much experience with bandlimited sampling theory. This paper contains the author’s accumulated reactions to drilling down on currently used spatial audio methods such as Wave Field Synthesis (WFS) and VectorBased Amplitude Panning (VBAP), without the level of background research (i.e., reading all known relevant papers) expected of a typical graduate student in the field. As a result, it is possible and even likely that some or many of the speaker arrays “introduced” here exist already in the literature and/or patent record in some form. The author welcomes citations that can be incorporated into an updated version of this paper.
Another result of the evolving nature of this paper is that it unfolds progressively from basic soundfield sampling to PBAP and its variations and finally to Huygens Arrays (HA). While it could make sense to write a fresh, shorter paper on one or more of these subtopics, the full story has greater tutorial value, and is less work to write, so here it all is in one relatively long paper.
Overview of PlanewaveBased Angle Panning (PBAP)
Studying Wave Field Synthesis (WFS) led to the idea of the considerably simpler PlanewaveBased Angle Panning (PBAP), which could just as well be called FarField Wave Field Synthesis (FFWFS). The simplifying assumptions of PBAP are

the microphonearray is positioned between the primary sources and the listeners, e.g., between a “stage” and all the “seats” in an auditorium, for example, as shown in Fig. 1;

the microphonearray is at least several wavelengths away from the nearest source, so that the micarray sees a superposition of plane waves to a good approximation;^{5}^{5}5Since every sourcefree soundfield can be constructed as a superposition of plane waves, it follows that solving the planewave synthesis problem is quite general. and

the speakerarray lies along a line parallel to or coincident with the mic array in the direction of the audience.
The simplest case of assumption 3 is when the speakers and microphones are coincident, i.e., and the speaker array lies along the same line as the mic array, so that each speaker reproduces the signal recorded at its microphone. The mic array can be viewed as a spatial A/D converter, sampling the incident wavefront, and the coincident speaker array gives the corresponding spatial D/A converter for these microphone samples. When or when not coincident and not spaced in a compensating manner, we must compute the speaker signals from spatialresampling of the mic signals.
The main point is that PBAP is obtained from WFS by moving all acoustic sources to “infinity”—i.e., onto the “celestial sphere”.^{6}^{6}6This analogy for PBAP sources as pointsources (stars) on the celestial sphere motivates the alternate name “Star Field Synthesis (SFS).” Viewing a starfield through a window is a vivid analogy for PBAP using a rectangular speaker array. However, diffraction effects are much more significant in SFS (especially at low frequencies) than when viewing light through a rectangular aperture, so the analogy can be misleading. In practical acoustics, “infinity” means “many wavelengths,” and thus depends on frequency. Huygens Arrays, developed next, generalize the linear/planar sampling arrays to more geometries while keeping a separation plane between primary sources and listeners.
Figure 1 illustrates the basic geometry assumed for singlelinearray PBAP. All sources are confined to a “stage area,” which can be a 3D distribution of sources when the microphonearray is planar, but such a distribution is effectively 2D due to assumption 2 (the microphones are in the “far field” of each source). All listeners are confined to the “audience area” which can also be a 3D distribution of “seats,” when the speakerarray is planar, but again the farfield assumption implies that each seat hears substantially the same planar distribution of sources far away. The microphone and speakerarrays thus form a separation plane between the stage and audience areas. For the linearray case, the sources are considered as being arranged on a flat stage (same height—only azimuth varying—the same case addressed by ordinary stereo).
In summary, our assumptions allow us to reconstruct a superposition of plane waves given pressure samples along any plane separating the source and audience area. Since the microphonearray is our pressuresampling array, it is natural to choose uniform spacing of the microphones along each coordinate dimension. (It is reasonable to choose a smaller spacing for the horizontal direction since people are more sensitive to azimuth than to elevation of source angleofarrival.) Furthermore, we need to ensure that the incident wavefield is properly bandlimited^{7}^{7}7https://ccrma.stanford.edu/~jos/resample/ [Smith, 2010]. In normal signal sampling, an antialiasing lowpassfilter removes highfrequency components that would alias. In the spatial sampling array, we can simply require that the incoming plane waves have a known maximum angle that does not “alias” spatially. This limits the angular “stage width.”
2 Plane Wave Sampling Theory
This section develops a special case of Wave Field Synthesis (WFS) by spatially sampling simple plane waves. Sampling plane waves is much simpler than the traditional WFS formulation which begins with the classical KirchhoffHelmholtz integral [Pierce, 1989, Firtha, 2018]. In return for this simplicity, we are restricted to virtual primary sources that are many wavelengths away from the speaker array, and on the other side of it from the listener. As we shall see, we can relax these restrictions in various ways, and the remaining sampling conditions are generally equally binding for WFS systems. In other words, spatial sampling theory is fundamental to all spatial audio systems using discrete drivers arranged in linear, planar, or even more general array geometries. What does not seem to be generally known, however, is that a samplingbased approach is sufficient (and much more to the point) for deriving and optimizing the system, as pursued in this paper.
Figure 2 shows a monochromatic plane wave traveling down and to the right at a 45 degree angle. The solid black line across the middle represents the microphonearray, ideally a uniformly spaced grid of tiny omnidirectional pressure microphones having no “acoustic shadow” at all; these microphones serve to sample the plane wave along the line. In the 3D case, the line represents one cut along a planar microphonearray. The sinusoid drawn along the microphone line indicates the pressure seen by each microphone. By the sampling theorem (applied now to spatial sampling using a microphonearray), we must have more than two microphones per wavelength along the line array. Thus, the required microphone density is determined by the minimum incident wavelength and the maximum angle of incidence , as derived below.
Figure 3 illustrates the geometry of the wavelengths involved. The wavelength of the incident sinusoidal plane wave is denoted , and denotes the wavelength of the sinusoidal pressure fluctuation seen by the microphone line array. As Fig. 3 makes clear, from the angle of incidence and incident wavelength , we have
(1) 
Let denote the microphone spacing along the axis. Then the sampling theorem requires
(2)  
(3) 
where is the speed of sound (m/s), is the maximum temporal frequency in Hz (typically 20 kHz for audio), and (radians) is the maximum planewave angle allowed.
For example, choosing kHz and (stage angle 90 degrees), and using m/s for sound speed, we obtain mm, or about halfinch spacing for the microphones. (The coincident speakerarray has the same samplingdensity requirement as the microphonearray.)
Reducing either or relaxes the spatial sampling density requirement. For example, if the “stage width” is reduced from 90 degrees () down to 40 degrees (), then oneinch spacing of the microphones (and speakers) is allowed. If we bandlimit our reconstruction bandwidth to 5 kHz, then we get by with fourinch spacing, as pursued below in a practical PBAP design (§2.15).
If we don’t bandlimit to below the spatial Nyquist limit, then we obtain “spatial angle aliasing” at very high frequencies for sources near the left or right edge of the “stage viewing window”. That is, for sources near the left or right edges of the “stage”, the highestfrequency components may not appear to come from the same direction as components below the cutoff frequency of 5 kHz. On the other hand, perception is such that the apparent angleofarrival typically may not alias at high frequencies because the desired angle remains a psychoacoustic choice and keeps the whole source spectrum in one place. Sources near the center of the stage are spatially oversampled by the microphonearray at all frequencies, so they are never a problem. In lowpassedwidebandnoise tests (see Appendix A), highfrequency spatial aliasing has been observed to break up a formerly coherent wideband virtual noise source, but not in a manner indicating “folding” as in aliasing of tones due to temporal sampling, but instead as the sound of a spurious new noise source somewhere along the array. The psychoacoustics of spatial aliasing perception is a fascinating topic for ongoing research.
2.1 Only Planar and Spherical Arrays Can Work Perfectly
It is well known that only spherical waves propagate in 3D without a wake, which includes plane waves as a limiting case. A line array produces a cylindrical wave, which has a wake [Morse and Ingard, 1968, p. 243]. This takes the form of temporal dispersion behind the wavefront in the soundfield reconstruction from linear arrays. Interestingly, nobody seems to complain about it, so it can’t be a particularly audible distortion. However, for best results, planar (or spherical) sample distributions are preferred.
Another disadvantage of cylindrical waves relative to plane waves is amplitude decay by (consider energy conservation), where is distance from the cylinder axis. This is less of a problem, since it preserves audio fidelity and spatial locations, and is potentially even desirable since it gives the listener a way to alter listening level by merely moving closer to or farther from the array (3dB per distance doubling).
We will continue to consider primarily linear arrays, in which soundfield reconstruction from samples happens only along the leftright dimension, and VBAP or stereo is used for the vertical dimension, if anything. There will be other approximation errors as well, and the key question as always is the audibility of those errors.
2.2 Fractional Delay
A nice property of ideal plane waves is that they do not decay, i.e., the amplitude of a given source is the same at every spatial sampling point (microphone or speaker). The only change from one spatial sample to the next is a relative delay corresponding to the different propagation distance. Thus, the only processing needed for each source for each speaker is the fractional delay [Franck, 2008, 2011, Välimäki, 1995] corresponding to the sourcespeaker distance. Furthermore, in a line or plane array, the change in delay from one spatialsample to the next is constant along any straight line, and this can be used to save overall computation [Franck, 2011].
2.3 A Synthesis Scenario
Let denote the number of virtual acoustic sources and let be the number of speakers. Then in PBAP, the th speaker signal is given by a sum of delayed source signals:
(4) 
where is the th source signal, and is the timedelay, in samples, from virtual source to speaker . The signal is obtained as the output of a fractionaldelay filter driven by . The same delay line can be used for source to all speakers, using “taps” to interpolate [Smith et al., 2002].
2.4 Quantized Angles of Arrival
Since highquality fractionaldelay filtering is expensive, it is worth considering restriction to anglesofarrival corresponding to integer delays (in samples). If the speakertospeaker spacing along a line array is , then the speakertospeaker delay for a plane wave at angleofincidence is , where denotes sound speed. Thus, an angleofarrival corresponds to an integer speakertospeaker delay (in samples) when
(5) 
where denotes the temporal sampling interval in seconds, and
Note that increasing the speaker spacing for a given temporal sampling rate gives more integerdelay angles . However, doing this also decreases the stagewidth (or supported bandwidth) by the same factor.
It is clearly inaudible to shift the location of each virtual source so that the time delay to the nearest speaker is an integer number of samples. Then having an integer number of samples for each interspeaker delay makes all the delays integer. Finally, this can all be implemented as a single delay line with a tap (noninterpolating) for each speaker signal. For moving sources, to avoid clicks, moving taps should be crossfaded from one integer delay to the next in the usual way [Smith, 2010].^{8}^{8}8https://ccrma.stanford.edu/~jos/pasp/Delay_Line_Signal_Interpolation.html
Solving Eq. (5), the collection of angles corresponding to integer interspeaker delays (in samples) is
(6) 
For example, with m (), m/s, and kHz, the first 11 available angles are
(7) 
degrees, to two digits precision. Azimuth perception is accurate to approximately 1 degree at centerfront.^{9}^{9}9http://acousticslab.org/psychoacoustics/PMFiles/Module07a.htm
Figure 4 depicts the available geometric rays of planewave propagation for this example. Thicker rays are drawn for 0 degrees (directly in front) and degrees (full left and right).
If the 21 anglesofarrival across a stage listed in Eq. (7) are deemed sufficient, then PBAP is essentially free: just provide the appropriate integer adjacentspeaker delays for each source in the sum for each speaker. As is well known, an integer delay is an computation, requiring only a single read, write, and circularbuffer pointerincrement each sampling instant [Smith, 2010].^{10}^{10}10https://ccrma.stanford.edu/~jos/pasp/Software_Delay_Line.html
2.5 FourQuadrant PBAP
If four line arrays are arranged in a “tic tac toe board” configuration (or perhaps just a square) enclosing the listener in its central square, as shown in Fig. 5, then each line array need only cover a degree range, which is more uniformly sampled in angle.
FourQuadrant PBAP uses only a fourth of the speakers for each plane wave. Infinitely long line arrays emit cylindrical waves, which are equivalent (ignoring the wake) to plane waves in one listening plane passing through the cylindrical axis. However, more practical truncated line arrays can benefit from using more than a fourth of the speakers. It is intuitively obvious that at least half of the speakers could help construct a desired plane wave at the listener—all speakers having a radial component in the desired direction.^{11}^{11}11The analogous condition for Wave Field Synthesis is to use only secondary sources that are “illuminated” by the virtual source being rendered [Ahrens, 2012]. At a single listening point, as in ambisonics, all of the speakers can be put to work toward approximating the desired soundfield pressure and velocity versus frequency at that point.
2.6 Polygonal PBAP (PPBAP)
The FourQuadrant PBAP of the previous section readily extends to sided polygons. Figure 6 shows a progression from (left), to (middle), and (right). Along the top of the figure is the result of using line arrays (as long as available, yielding the best planewave slices in the listening plane). Along the bottom is the compromise obtained by truncating each line array to the portion that borders the interior polygonal listening area (Truncated PPBAP).
2.7 Circular PBAP Truncates to VBAP
In the limit as the number of polygon sides becomes large, we obtain a circular array, having only one speaker representing each line array in the truncated PPBAP case. Furthermore, the available angles (when avoiding interpolation) are simply the speaker angles. This coincides with zerothorder VBAP (Vector Based Amplitude Panning) [Pulkki, 2001, 1997].
In practical VBAP, sources are typically enlarged to more than a single speaker, both to make a sonically larger source, and so that they pan more smoothly from one location to the next around the ring (or dome, etc.). This of course gets us into interpolation, and can be viewed as such. The multispeaker interpolation strategies of VBAP can be applied to PBAP, both to the interior polygon/circle as in VBAP, and more generally to the line arrays creating an sided polygon. In the latter case, PBAP can be viewed as a sweetspot enlarging strategy for VBAP. The longer the tangential line arrays, the straighter the plane waves emitted, and the less sensitivity to listener position downstream from the plane waves.
2.8 Wrapped Polygonal PBAP (WPPBAP)
A problem with VBAP’s sweetspot size is that each speaker is approximately a spherical wave source. Therefore, Truncated Polygonal PBAP produces quasi spherical waves from each polygon side (when the wavelength is large compared with polygon side length). To address this, we can use nearly half of all speakers surviving the truncation to participate in the planewave generation from each original line array.
In particular, starting with Polygonal PBAP, the process of truncating to the interior polygon can include summing the signal “seen” from each truncated speaker arriving at each surviving speaker. This is essentially just applying the sampling principles used to derive PBAP in the first place. Thus, each speaker signal will include the signal from its original line array, plus a delayed and attenuated signal from each truncated speaker that is “behind it” relative to the listener, and reasonably close. The truncated speakers are now treated as point sources, so the attenuation is proportional to inversedistance as usual for spherical waves. A maximum distance is set beyond which speakers from the (imagined) extended line array are not heard. This leaves less than half of the surviving speakers receiving contributions from the truncated speakers on any given line array.
2.9 Combining Line Arrays to make Polygons
When avoiding delayline interpolation and accepting the angles given by integer interspeaker delays, we should choose the sampling rate and speaker linearray spacing so that the angles available from each line array include the angles to the polygon vertices.
For an sided polygon, the two needed angles are . The set of all vertex angles is , , where is sound speed and is the sampling interval. Thus, we need some integer to give , or
for some .
In the fourquadrant case (four line arrays defining a square) with speaker spacing in (example from §2.15 below), the sampling rate wants to be a multiple of , and kHz happens to be , so the 10th angle is very close to 45 degrees. For polygon sides and fourinch speaker spacing, we need to be a multiple of , and it so happens that 44.1 kHz is close to five times that ().
2.10 Increasing the Number of Source Angles
The number of available ray angles in PBAP can be increased by various means.
Oversampling in time to kHz doubles the number of integerdelay angles over the same range. The example of Eq. (7) on page 7 goes from 21 to 41 anglesofarrival across a frontal stage.
Another avenue for doubling the number of source angles is to implement a halfsample fractionaldelay filter. In that case, half of the available angles require use of the filter while the other half require no interpolation filter.
Generalizing, given fractionaldelay filters, with the th filter providing fractional delay samples, , the number of available arrival angles is increased by the factor . A benefit of a fixed grid of available angles is that each fractionaldelay filter can be individually designed and optimized for the delay it provides. Three such filters () suffice to provide approximately complete perceptual resolution in the example of Eq. (7).
2.11 Continuous Angles of Arrival
Quantized anglesofarrival do not suffice when sources move over time. Also, continuous anglesofarrival allow a less quantized “source width” parameter for each source—as if the source were coming from a distant solidangle region (like a nebula) instead of one point (star).
There are various methods for continuously variable fractionaldelay filtering [Smith, 2010]. Perhaps the simplest is by means of Lagrange interpolation [Välimäki, 1995, Franck, 2008, 2011, Smith, 2010].
The case of firstorder Lagrange interpolation is especially simple, being simple linear interpolation. Thus, one can linearly “crossfade” in amplitude from one source angle to the next to implement a moving source. Distributed sources can be formed as a linear combination of adjacent source angles.
2.12 Covering the Spatial Hearing Frequency Range
Spatial hearing [Blauert, 1997] is accomplished by two ears sampling
the acoustic field through small apertures (ear canals) having
diameters smaller than the wavelength across almost the entire audio
band. As a result, the directionality of a sound is inferred primarily
from the relative intensity () and timeofarrival () at the two
ears. There is also directionality information impressed on the
signal by pinnae filtering and shoulder reflections, etc., that are
especially important for elevation perception.^{12}^{12}12http://en.wikipedia.org/wiki/Headrelated_transfer_function
For azimuth perception, is the dominant cue below about 800 Hz, and
dominates above 1600 Hz or so.^{13}^{13}13http://en.wikipedia.org/wiki/Sound_localization#ITD_and_ILD.
In [Gerzon, 1974], citing
Rayleigh from 1907, the lowfrequency crossover is given as 700
Hz. It is also noted in [Gerzon, 1974] that pinnae filtering is
thought to be important above 5 kHz.
In the octave between these
limits, both and are used. Also, is picked up as
phase delay for low frequencies, and group delay at high
frequencies (ibid.).
Perceptual accuracy is on the order of 1 for azimuth in
front of the listener. The lower limit of azimuth perception based on
is approximately 80 Hz, below which phase differences become
imperceptible. Thus, our spec is to synthesize correct s down to
80 Hz. Note, however, that since the acoustic wavelength at 80 Hz is
over 4 meters long, we could get by with reduced spatial resolution in
this frequency range, such as simple stereo. Multiresolution speaker
arrays are discussed starting in §4 below.
2.13 Delaying High Frequencies to Suppress Aliasing via Precedence Effect
As mentioned in §2, spatial aliasing limits the highest frequency and/or the widest source angle supported by a uniformly sampled line array. As also mentioned there, spatial aliasing may not be perceived because the spectrum as a whole may lock in perceptually at the correct angle. In other words, the ambiguity of the spatial angle of the highestfrequency components may be resolved perceptually by the brain’s natural desire to “make sense” of an auditory scene. This effect can be enhanced by slightly delaying the highfrequency components relative to the lowfrequency components that have no aliased interpretation. The idea is to force perception to hear the desired angle before the ambiguous spectral components are heard, so that they will all fuse at correct angle of arrival. This is of course altering the timbre of the sound, and may be considered off limits for that reason.
2.14 FiniteArray Correction
In practice, it is necessary to truncate an array to finite bounds. This causes reconstruction error analogous to the error obtained when restoring a continuoustime signal from a finite segment of its samples. Thus, most of the error occurs for sources near the edges of the array (i.e., arriving from the maximum anglesofarrival supported). This error can be reduced by compensating for the missing contributions from the truncated “sampling kernels”. From this point of view, the error is equivalent mathematically to the “Gibbs phenomenon,” and many forms of “windowing” and “apodization” have been advanced to address the issue [Smith, 2011].^{14}^{14}14https://ccrma.stanford.edu/~jos/sasp/Spectrum_Analysis_Windows.html
One can also formulate a customized optimization that maximizes perceptual criteria; for this problem, linear programming formulations may suffice for correcting amplitude error.
^{15}^{15}15https://ccrma.stanford.edu/~jos/sasp/Window_Design_Linear_Programming.html A Hann window is used for array windowing in the Sound Field Synthesis Toolbox [Wierstorf and Spors, 2012].^{16}^{16}16https://github.com/sfstoolbox/sfs2.15 A FourInch Grid Implementation
The Meyer Sound MM4XP Miniature Loudspeaker (Fig. 7) is a selfpowered speaker that provides an approximate fourinch by fourinch cell. Thus, with this speaker we can make either a line array or a 2D grid with fourinch spacing. While this speaker is relatively expensive, it exhibits excellent power and linearity for its size.
The MM4XP power output is 113 dB at 120 Hz. Extending the low end beyond the spatialhearing range down one octave requires a 6 to 12 dB sacrifice in power for the same peak diaphragm excursion.^{17}^{17}17To see this, consider that the speaker is much smaller than the wavelength produced at its low end, so it can be regarded as an approximate “point source” in that frequency range. From the theory of a point source [Morse and Ingard, 1968, p. 310], the peak pressureamplitude from a sinusoidally oscillating point source is proportional to the peak volume acceleration from the source, which is in turn proportional to the radial acceleration (second timederivative of spherical radius) for any small sphere used to model the simple source. (Any sphere much smaller than a wavelength in diameter will do.) Let the spherical radius acceleration be denoted , which is proportional to the farfield pressure a fixed distance away. Then the peak radial excursion of the spherical surface is given by , and keeping the excursion fixed while decreasing by one octave reduces the farfield pressure by a factor of four, or dB. When operating as a cell of a planar array, on the other hand, farfield pressure is proportional to the driver surface velocity instead of acceleration. In that case, only 6 dB per octave is lost integrating velocity to get displacement. Thus, the extendedlowend array is plenty loud—and we may cross over to a subwoofer below the spatial hearing range.
2.16 Extension from a 2D Listening Plane to 3D
Since PBAP converges to VBAP when using many line arrays truncated to the enclosed polygon, a simple VBAPstyle extension to 3D is to place a new speaker directly overhead and a second new speaker directly below. Then, elevation cues can be imparted by mixing in the above or below speaker according to a psychoacoustically measured panning law between the array and an outofplane speaker above or below.
The proper extension of PBAP to 3D is of course obtained using sampled plane waves arriving at the correct 3D angles, instead of cylindrical waves sampled in the listening plane for 2D PBAP from a line array. That means our line arrays must be replaced by planar speaker arrays, and the polygonal listening space becomes a polyhedron, or sphere in the limit. This of course also reduces to noninterpolating spherical VBAP when the planar arrays are truncated down to one point in each plane, and normal VBAP interpolation can be used as described above here as well. As in the 2D case (§2.7 on page 2.7), 3D PBAP can be used as a sweetspot enlarger for 3D VBAP.
3 Huygens Arrays (HA)
We have so far considered only PBAP and its variants, which can be considered FarField Wave Field Synthesis (FFWFS), in which each source contributes a plane wave to each listening point. This is pretty general, in principle, because, as is well known, every sourcefree soundfield can be expressed as a sum of plane waves at various amplitudes, phases, and directions of arrival. In fact, Fourier transform methods can be used for this purpose.
^{18}^{18}18https://ccrma.stanford.edu/~jos/pasp/Vector_Wavenumber.html Therefore, a straightforward path from FFWFS to fullfledged WFS is to decompose any desired soundfield into a sum of plane waves, and then generate those plane waves using FFWFS/PBAP. There are many known methods for socalled Planewave Decomposition (PWD) [Pulkki, 2017].A more direct extension of PBAP toward WFS is based again on simple sampling of the acoustic source wave, but now allowing spherical waves instead of only plane waves in the soundfield reconstruction. We could call this SphereBased Range and Angle Panning (SBRAP). However, reconstructing a wavefront as a superposition of spherical waves is essentially the idea of Huygens’ Principle. We therefore choose the name Huygens Array (HA) for the extension of PBAP to include spherical as well as planar wavefronts.
For constructing a Huygens Array, each virtual source is at a known location in 3D space:
According to the basic sampling assumption in PBAP, each speaker location is also a microphone location, so we can denote the th speaker/mic location by , . Different mic distributions are obtainable via spatial resampling as before.
In PBAP, each virtual source was characterized by an angle of arrival , which determined the interspeaker delay
in seconds along the line array, where denotes the centertocenter speaker spacing. To generalize from plane waves to spherical waves, we need both a delay and a gain describing the acoustic ray from source to speaker .
Let denote the amplitude of source at a distance one meter from its center. (Each source is assumed to be a point source for now; distributed sources can be modeled as weighted sums of point sources.) Then for the delays we have
where denotes the Euclidean norm of , denotes sound speed, and is the digital audio sampling interval, as before. For the gains we have
and, if desired, lowpassfiltering due to air absorption can be included:
where is a standard air absorption filter corresponding to propagating meters through air at some assumed standard conditions (humidity level being the most important) [Smith, 2010].
3.1 Linear Huygens Arrays
We have so far not used any assumptions regarding the microphone/speaker array to be used.
The sampling analysis of §2 on page 2 made use of the farfield assumption in obtaining a spatial aliasing limit that depended only on the source angle and spatial frequency . Generalizing to the near field (arbitrary source distances) means that the sampling analysis is applicable only locally along the array. That is, the wavelength seen by the line array depends on both the source angle and the distance of the source to the array (or equivalently, the relative distance of the source to the array and to the listening point). For example, the line from the source normal to a horizontal array (see Fig. 3 on page 3) is at angle , which is always oversampled by the array. Points on the array far away from the normal line, however, see an angle approaching degrees to the right and degrees to the left. If a source touches the array at , then all of array points other than the point at see a right angle ( degrees). This behavior means we cannot set a limited stageangle to avoid spatial aliasing like we did in the farfield case (§2). We can now accept a 180degree stage, or limit the closeness and layoutwidth of the sources to obtain a worstcase angle limit (maximally close to the array at the edge of the allowed stage), and treat that as before.
3.2 Sampling Spreading Loss
In addition to spatial audio oscillations from a point source, there is amplitude change due to “spreading loss” away from the source. To look at this, consider that a unit pressure pointsource at the origin can be expressed as the real part of [Morse and Ingard, 1968]
where is the radial coordinate axis. There is no change in the complex amplitude along directions with constant radius . Along we observe the maximum amplitude changerate. This rate of change is approached asymptotically along the array in both directions. The pressure gradient is given by
Intuitively, a sampling grid that is adequate for sampling spatial frequencies should be adequate for sampling spherical spreading (decay by ) when
or
(8) 
That is, to
keep the rate of amplitudechange due to spreading loss much less than
that due to acoustic vibration, we can keep all sources a few
minimumwavelengths or more away from the line array. Note that this
strategy only provides approximately valid sampling of the
spreadingloss decay, because is not a bandlimited function and
therefore cannot be sampled without some error in the
reconstruction.^{19}^{19}19The 2D Fourier transform of can be
shown to be :
http://sepwww.stanford.edu/public/docs/sep103/jon3/paper_html/node3.html
Fortunately, the error can be made zero psychoacoustically at a
reasonable sampling density. Amplitude error perception generally
requires at least a quarterdB difference, and that’s in the most
demanding case of comparing alternating amplitude levels.
Since the speaker spacing must be smaller than half the minimum wavelength , we can stipulate that all sources should stay at least several speakerspacings away from the array.
An alternative strategy to Eq. (8) is to double the linear sampling density of the array and allow amplitudechange due to spreading loss become comparable to that due to vibration. In this case, the minimum approach distance becomes (), allowing sources get to within of the minimum wavelength from the array, or about a third of the centertocenter speaker spacing. Intuitively, thinking of the speaker array as a sampling grid, it makes sense to keep virtual pointsources on the order of a sample away or more.
We learn in a first course on digital signal processing that a signal must be bandlimited to less than half the sampling rate in order to avoid aliasing [Smith, 2007b]. Setting a minimum on how close a virtual source may approach the sampling array effectively spatially bandlimits the wavefront geometry. This enables soundfield sampling to work as intended. Analogous bandlimiting happens along the time dimension when we apply an A/D lowpass filter prior to sampling in time.
3.3 Virtual Sources in Front of the Array
In all cases considered so far, the virtual sources (primary sources) have been restricted to be behind the speaker array by some minimum distance for valid sampling. We can now extend as in WFS to allow sources between the array and the listener, but we still must maintain the same minimum distance, but now from the other side of the array. There are differences, however, to keep in mind relative to the behindthearray case. For simplicity, consider a line or plane array as a starting point.

Unlike the primary sources behind the array, those in front of the array have to first propagate the “wrong way” to the array to provided recorded signal components (spatial samples) that can be played back from the array to create a converging wavefront back to the virtual source and then on to the listener. This means there is an added delay between the source and the listener, as if the listener can only hear the firstorder reflection of the source bouncing off the array, with no direct signal from the virtual source. The listener hears the reflection from the array after it passes through the point of convergence at the virtual source and then propagates to the listener.

The listener receives a mirrored reflection of the rear of the virtual source. This is no problem for an isotropic source, like any monopole, but it can be an error for distributed sources trying to achieve a specific natural radiation pattern. As a result, primary sources should face the array instead of the listener, and be flipped as needed.

In offline applications the extra delay for interior sources is easily removed in postprocessing.
3.4 Undersampled Huygens Arrays become VBAP
A practical issue that arises when an array is allowed to undersample the soundfield is level normalization. An extreme example is the amplitude () of a point source when it touches a microphone/speaker point on the array, which is well out of bounds for any reasonable practical system.
Level regularization is another reason to keep sources a few wavelengths or more away from the microphone/speaker array.
When the array becomes undersampled, say because a wideband source is approaching an extreme stage angle, the array can be regarded as transitioning from soundfield reconstruction by summing interpolation kernels to simple panning between/among available speakers (i.e., what we’re normally always doing). Note that the highfrequencies must go into “panning mode” first, while lower frequencies may remain adequately sampled. We can implement an adaptive spectral partition between samplebased reconstruction below and panning above.^{20}^{20}20This spectral partition issue is related to the classic “panning problem” in which low frequencies see a 3dB boost relative to high frequencies, due to coherent versus noncoherent summation from stereo speakers for an offaxis listener: https://ccrma.stanford.edu/~jos/sasp/Panning_Problem.html
3.5 More General Huygens Arrays
We can also drop the restriction to microphone/speaker line/plane arrays and allow the speakers to be more generally laid out, such as on a sphere (a typically available layout for ambisonics systems). The basic Huygens principle remains the same: when a wavefront reaches a speaker, it fires out a spherical secondary wave (or hemispherical wave aimed toward the sphere’s interior). This configuration was simulated with a plane wave excitation, and a decent planewave reconstruction was obtained inside the spherical space. For a spherical speaker array, only half of the speakers become activated by a source outside the sphere.
More generally, when virtual sources are kept away from uniformly laid out speakers, Huygens’ Principle tends to hold up pretty well, according to simulation results to date.
It remains preferable to have a reasonable sampling grid for most frequencies, and we continue to prefer a separation plane between the sources and listeners (Fig. 1 on page 1).
One argument in favor of a separation plane has to do with the fundamental limitation of Huygens’ Principle, which only considers pressure (a scalar) and not velocity (a 3D vector). Without velocity matching, a line array of point sources creates a cylindrically symmetric output. A plane array of point sources similarly emits identical wavefronts in both directions away from the plane. If all listeners are on one side of the speaker array, then we don’t care what happens on the other side. However, this argument can be overcome.
Since audio loudspeakers are normally baffled, we get an approximate hemispherical source from them, which is ideal in the limit of a continuous distribution of point sources. Also, matching pressure across both time and space implies velocity matching, since ultimately the two are tied together (in a sourcefree region) by the wave equation, as discussed further in the next subsection.
From our sampling point of view, what we want are speakers having a radiation pattern that serves well as an interpolation kernel (the spatial shape of one sample) for reconstructing a soundfield from its samples in one direction leaving the speaker array. Like Huygens, we want to neglect velocity and work only with pressure samples, but generate the correct velocities indirectly using pressure differentials along the array and across time, when the speakers are close enough for this to work, as they are in a valid sampling grid.
Another argument in favor of a sourcelistener separation plane is that there can be no standing waves when sound generally propagates from a set of sources to a set of listeners. Standing waves could pose problems if/when our microphone/speaker array happens to line up with a node line. A soundfield is uncontrollable and unobservable at nodes of vibration, leading to possible degeneracies in implementation. A progressive wavefield cannot have these problems.
It can also be taken as a simple design decision that we want our wavefronts to cross a separation plane from sources to listeners. This implies we are not trying to synthesize the reverberant field like WFS, and we set up our arrays in normal acoustic environments, as opposed to the anechoic environments called for by WFS. WFS solves for signals at all speakers enclosing the listening space to produce the desired interior field, even when it contains sources, and including any reverberation. We are less ambitious with Huygens Arrays, and we can expect more robustness and intuitionguided design freedom as a result.
3.6 Interpolation Accuracy
In a uniform speaker array (line or plane) that can be considered a “spatial digital to analog converter (D/A)”, the speaker’s radiation pattern plays the role of “sampling kernel,” or “reconstruction lowpassfilter impulseresponse” used in D/A conversion.
A speaker’s radiation pattern is naturally given as a function of angle relative to the speaker’s central axis. It may therefore appear to be a problem that our acoustic “samples” are changing (diffracting) as they propagate away from the speaker. However, this is ok as long as there is any distance from the speaker array at which the soundfield has been reconstructed (both pressure and velocity versus position). Beyond that, the wavefront “takes care of itself” (consider the more complete HuygensFresnel principle [Miller, 1991]).
In the farfield case (PBAP), we observe only traveling plane waves, where pressure and velocity are proportional to each other. In the case of a plane wave impinging on the array at angle , the velocity is obviously in the correct direction by symmetry. For an angled plane wave, it is easy to show that the velocity angle is correct far downstream from the array, since the contribution of each pointsource along the array becomes a plane wave, and the lines of constant phase are along the desired angle due to the relative timing of the point sources. Having the pressure of a plane wave propagating in the desired angle means that its velocity points in that direction as well.
It is more difficult to show velocity reconstruction in the nearfield case (general HA), where diverging wavefronts are sampled and reemitted by the array. In this case, we need to show that the pressure samples add up to give the both the correct pressure and velocity for continuing the spherical wavefront expansion from the source. We have found simulation results to be helpful in the absence of analytical results (see Appendix A). There is no error in the timing of the secondary wavefronts, making the simulation results look great (a reconstructed plane wave looks very planar), but there is in principle amplitude error along the synthesized wavefront caused by the individualized factors. This error declines in relative terms downstream from the array.
While the instantaneous pressure along a line or plane does not determine the corresponding velocity, unless we can assume a traveling plane wave, etc., the wavefront pressure over a nonzero time interval does determine the particle velocity . In fact, the wave equation itself can be integrated to compute it. The time interval creates an interval of pressure history which allows the pressure gradient to be calculated (in a progressive wave), and the pressure gradient drives the velocity in the absence of a coincident source. Specifically, the sound velocity is given by the timeintegral of the pressuregradient divided by the air’s massdensity (Newton’s second law of motion , which, in the wave equation, appears as , neglecting highorder terms [Morse and Ingard, 1968, p. 243]).
It thus suffices to reconstruct sound pressure as a function of time and position along any parallel line (or plane) in front of the array. Since we assume no sources along that line or plane, the velocity is determined by the pressure in any spatial neighborhood.
In the case of the line array, we are only concerned about the velocity vector lying in the plane determined by the line array and any vector from the line array to the sources or listeners, all of which are assumed to lie in one plane.
For a planar microphone/speaker array, the velocity vector can point anywhere in the halfspace from the array toward the listeners.
3.7 Evolving Radiation Pattern (Sampling Kernel) due to Diffraction
It is well known that the farfield radiation pattern of an acoustic source is proportional to the spatial Fourier transform of the source’s radiating amplitude distribution. In optics, this is called Fraunhofer diffraction theory [Goodman, 2005]. By linearity of the Fourier transform, it follows that if the speaker array emits valid sampling kernels near the array (i.e., the radiation patterns overlap and add to a constant overall gain for a plane wave), then the farfield patterns will similarly overlap and add to a constant gain. All points in between must be valid as well, being progressive diffractions of the source distribution, but a rigorous proof with quantified approximation error would be nice to see (there are always terms to neglect, and it is good to be mindful of them).
3.8 Specific Sample Shapes (Speaker Radiation Patterns)
Perhaps the most obviously ideal sampling kernel for a speaker centered at within a linearray with cellwidth is the rectangular pulse:
That is, each speaker in the array pushes a unit volume of air in one second (unit area). Such a speaker array can be regarded as a series of contiguous pistons, each of width .
The unitarea pulse function is not bandlimited. A more graceful choice is the “spatial sinc function,” which is the ideal sampling kernel used in bandlimited audio sampling:
where
It is well known that the sinc function is the Fourier transform of a rectangular pulse, and vice versa. Therefore, rectangularpulse samples at the line array propagate and diffract to become overlapping sinc functions in the far field, while sincshaped radiation patterns diffract to become adjacent or overlapping rectangular spatial bands in the far field. Of course, a true sinc function can never be implemented precisely because it is infinite in spatial extent. Still, it is interesting to imagine disjoint spatial samples in the far field, and consider whether approaching that might be desirable in some situations, such as delivering different program material to different listening positions (a current problem in automotive sound).
The most typical example in practice is a circular driver—an ordinary circular speaker cone. In this case, the far field shape is the socalled “Airy function” involving Bessel function .^{21}^{21}21 https://adriftjustoffthecoast.wordpress.com/2013/06/06/2dfouriertransformoftheunitdisk/, http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf
In addition to the rectangular piston and sinc function, or circular driver and Airy pattern, any spectrumanalysis window , or its Fourier transform , having the constantoverlapadd (COLA) property
can be used as a spatial sampling kernel, with the added stipulation that should have an effective width on the order of (one spatial sample) in order that spatial resolution be maximized. In other words, it suffices for the speaker radiation gains at a particular distance away from the array to overlapandadd to a constant gainversusposition at the speaker spacing used. In the limit of infinite sampling density, all windows have the COLA property, including hemispherical, which is in part why Huygens’ Principle works as a spatial D/A converter with its spherical/circular wavefronts. A review of COLA windows and their design is given in [Smith, 2011].
Spatial oversampling gives more flexibility in the choice of radiation pattern. For example, doubling the spatial sampling rate allows the radiation patterns to overlap by an additional factor of two with no loss of spatial resolution along the line array. On the other hand, that factor of two could be given to the aliasing cutoff frequency, which is probably a better use of it.
As in audio interpolation, windowed sinc interpolation can be used in practice [Smith and Gossett, 1984], if a speaker driver can be devised to generate that radiation pattern at some distance from the speaker.
The spatial sampling kernel should ideally be frequency independent, but this is never the case for typical speaker systems. Instead, typical speakers look like point sources at low frequencies, radiate efficiently and widely at wavelengths comparable to the speaker diameter, and begin to “spotlight” increasingly at higher frequencies (where the diameter is multiple wavelengths).
Due to the naturally narrowing spatial beamwidth with increasing frequency for typical speaker drivers, a sufficient sampling density for high frequencies corresponds to heavy oversampling (spatially) at low frequencies. PBAP and its extensions therefore should either be implemented in separate frequency bands (multiband PBAP is discussed below), or using a new kind of speaker having a frequencyindependent radiation pattern that sums to a constant when the speaker outputs are all added together at any point of the listening region. One solution is to approximate a point source (see Appendix B), for which the sampling kernel is a substantially identical sphere for all speaker drivers much smaller than a wavelength. Such speakers radiate inefficiently, but they are already widely used at the lowfrequency end in practice (subwoofers are smaller than most of the wavelengths they must produce). The main problem with a single set of pointsourceapproximation drivers is that they must be packed very densely for the high frequencies and also have longthrow excursion for low frequencies—expensive. Furthermore, there is always intermodulation distortion in any wideband driver (e.g., Doppler shift of highfrequency components by lowfrequency excursion, which nobody apparently compensates).
4 Multiresolution Spatial Sampling Arrays
Since typical speakers have a frequencydependent radiation pattern, it makes sense to take a multiband approach and combine an array of “woofers” with a denser grid of midrange speakers and a yet denser array of “tweeters,” for example. In principle, each speaker is assigned a frequency band containing wavelengths up to its diameter or so, and larger wavelengths can be pushed using additional volumevelocity drive. Of course, one large 17meter (56’) diameter diaphragm can handle the entire audio band, but then we have poor spatial resolution at high frequencies.
4.1 Multiresolution Coaxial Drivers Array
One approach suitable for 2D arrays is coaxial heterogeneous drivers, as depicted in Fig. 9. The Kenwood KFC1695PS, for example, provides three concentric drivers with diameters 6 1/2, 1 9/16, and 1/2 inches (75 Watts RMS). This could be packed into a rectangular panel and treated as a multiresolution Huygens Array or PBAP/VBAP system.
One can also imagine this type of array implemented using pipes of various diameters, with the pistons driving the larger pipes operating behind the pistons of the smaller pipes, where each piston is driven at the end of a narrow rod that passes through a small hole in the piston(s) behind.
4.2 Multiresolution Line Array
A simpler geometry for the linearray case is parallel strips covering different bands (Fig. 10).
Of course, the drivers do not have to be circular, and some tweeter geometries are noncircular. Also, circular drivers can drive the base of horns having rectangular exit apertures that tessellate a linear or planar region.
4.3 Multiresolution Sampling Considerations
Intuitively, spatial sample reconstruction requires that each driver be capable of pressurizing a fraction of one wavelength (in its band) to the desired pressure level. Thus, while the speakers are ideally very small, on the order of drinkingstraw diameters at high frequencies, the drivers need a long excursion in order to push enough air down the straw to achieve the desired pressure within the subwavelength zone being served. It is straightforward to calculate the maximum piston excursion needed for a given sound pressure level and lowest sinusoidal frequency. Dividing up the spectrum into frequency bands makes this easier.
The exact shape of the drivers is not important when they are smaller than a wavelength, only that they can pressurize their subwavelength zone as needed. Perhaps the easiest solution conceptually is a grid of contiguous square pistons. In that case it is easy to see that it must work very well, because the pistons can generate the wave propagation leaving the surface in great detail.
In classical “critical sampling,” there would two pistons per wavelength, one to push while the adjacent piston pulls. In practice, critical sampling may cause undesired noise due to turbulence, since there is no guarantee that laminar flow is maintained. Obtaining silent pressurization of half a wavelength may prove difficult at high sound pressure levels, so spatial oversampling helps.
4.4 Multiresolution Speaker Systems
We are all familiar with speaker cabinets containing woofers, midrange, and tweeters, etc. Each cabinet can be considered a monaural multiresolution speaker system, typically threeway. Additionally there is often a subwoofer somewhere putting out the deep bass.
The Kenwood JL840W speaker systems are fourway:^{22}^{22}22https://www.acs.psu.edu/drussell/Demos/BaffledPiston/BaffledPiston.html They use four circular drivers having diameters 30, 12, 6, and 3 cm.^{23}^{23}23These are “octave spaced” from midrange to supertweeter, with an extra large woofer. To adhere to octave spacing, the woofer diameter could be changed to 24 cm and the lower frequencies could be taken over by a subwoofer where the woofer leaves off. The crossover frequencies are at 2, 5, and 10 kHz, which is at for the speaker driving below crossover, where is wave number in radians per meter, as usual, and is speaker radius.^{24}^{24}24
This implies the lower speaker diameter is 1.75 wavelengths at crossover while the upper speaker diameter is 0.7 (midrange) or 0.88 (tweeter and supertweeter) wavelengths at crossover. The geometric means of these two diameters in wavelengths are 1.1 (woofermidrange) and 1.24 (other two crossovers) wavelengths. The general tradeoff is that driving with diameter less than a wavelength is inefficient (below cutoff), but yields nicely omnidirectional radiation, while driving with diameter much larger than a wavelength becomes highly directional.
The nominal total frequency range of the system is 20–20 kHz, but amplitude drops off in the woofer (30 cm) for wavelengths much greater than the diameter, which must be compensated by extra drive. The supertweeter, tweeter, and midrange drivers have diameters on the order of one to two wavelengths. The woofer highend is near that range, but must handle all lower frequencies as well.We can extend existing way speaker systems to multiresolution line arrays in a straightforward manner. We use the term Huygens Array (HA) to refer to any collection of drivers used as a spatial imaging array, and in particular, Huygens Octave Panels (HOP) will refer to multiresolution line arrays having octave divisions.^{25}^{25}25We avoid the term “Huygens Octave Array” (HOA), which would be nice to use for planar arrays as in Fig. 9, due to the common use of HOA as “Higher Order Ambisonics”. It is convenient to organize multiresolution line arrays into panels, so not much is lost.
4.5 FourWay Huygens Arrays
A fourway line array inspired by the Kenwood JL840W fourway speaker system can be made using four parallel rows of contiguous speaker drivers, packed together as closely as possible, as indicated in Fig. 11.
The extra large woofer breaks the pattern, suggesting using two of them as a stereo pair. Doing this makes room for a center channel, as shown in Fig. 12.
4.6 Huygens Octave Panels (HOP)
In audio, octavebased resolution is very common. If the upper limit of the audio bandwidth is taken to be 20 kHz, then the top octave covers 10–20 kHz, the next octave down is 5–10 kHz, and so on, giving crossover frequencies at 10, 5, 2.5, and 1.25 kHz, and 612, 312, 156, 78, and 39 Hz, which can be taken as the crossover to the lowest octave including 20 Hz. Thus, the complete audio spectrum spans 10 octaves, and the 9 crossover frequencies are given in kHz by , where is the octave number counting down from the top.
It is common to use a subwoofer to take the low end (spanning 20–80 Hz for THX, 20–100 Hz for “pro”, or 40–200 for “consumer” quality level, etc.^{26}^{26}2620190930: https://en.wikipedia.org/wiki/Subwoofer), so that only 8 (starting from 78 Hz) or 7 (from 156 Hz) octave bands are needed from a multiresolution array.
It is instructive to look at the wavelengths in each band, since each speakerdriver diameter needs to be on the order of a wavelength. Let’s define the centerfrequency of each octave as the geometric mean of its limits, so , , which gives
for the center frequencies in Hz. Then for a speed of sound m/s, using , we obtain the centerfrequency wavelengths to be
or
or
We see that even consumer quality level wants a lowestoctave driver diameter on the order of 10 feet, and THX and pro quality want a 20foot cone! Practical systems rarely go for such large drivers. Instead, we settle for the top five or six octaves, and drive the lowest octave with additional gain to get the desired power level. In other words, the lowend speaker(s) operate in a “rolling off” region, driving only a fraction of a wavelength, and so they require a 6 dB boost for each halving of frequency in that zone. We can make up for driving less than a wavelength by using extra power, to the extent no audible turbulence is generated.
Another point to consider is that sound localization is only meaningful down to a fraction of a wavelength, and wavelengths much larger than our heads cannot be localized based on the steadystate soundfield, because our ears get nearly identical signals from the field.^{27}^{27}27http://acousticslab.org/psychoacoustics/PMFiles/Module07a.htm We generally localize a sound based on its higher frequency components, such as above 500 Hz where azimuth changes cause noticeable Interaural Intensity Differences (IID), i.e., audible “head shadowing”. The localization determined during a sound’s onset, which normally has the most highfrequency content, tends to persist psychologically even after the highfrequencies that localized it have faded away.^{28}^{28}28A great demonstration of this effect in stereo, is to turn on a gated sinusoid in, say, the right speaker, crossfade the signal over to the left speaker, then pull out the right speaker cable and hand it to the befuddled listener who still hears the tone coming out of the right speaker. Thanks to Bill Putnam for showing me this one.
4.7 FiveBand HOPs
Extending the fourway Huygens Array of §4.5 to five bands, and enforcing octaveband design all the way down yields a fiveband Huygens Octave Panel for our consideration.
Figure 13 shows the driver outlines for two fiveoctave HOPs side by side. Each panel is about half a meter wide, so a 4m wall needs 8 HOPs.
The top three rows use drivers of the same sizes used in the Kenwood JL840W fourway speakers (3, 6, and 12 cm), while the fourth row is reduced to 24 cm (down from 30 cm) in order to continue the octave divisions, and a fifth row is added employing 48 cm (19 in) drivers. Fifteen inch speakers are very common, and could be used as a compromise for the bottom row, and subwoofer drivers of diameter 18 and 20 inches appear to be available, but probably quite expensive to use in an array system.^{29}^{29}29https://www.svsound.com/blogs/svs/strengthsandpitfallsofbigsubwooferdrivers
4.7.1 HOP Stage Angle
An immediate question that arises is how large is the maximum angleofarrival for a PBAP array made using such a HOP? From §2 on page 2, we have
Assuming it is possible to pack the speakers contiguously, and setting m/s, we obtain for a top row using driver diameters 0.03 m as in Fig. 13
or 16.6 degrees, giving a 33.2degree stage at . The same result is obtained for the highest frequency in each band of this example HOP, due to its strict octave scaling. The aliasfree angleofarrival range is unfortunately narrow, but it is calculated at the high edge of each band. At the bottom of each band, we obtain a much nicer 70 degree soundstage. At the centerfrequency of each band, defined as the geometric mean of the band edges, we find that a 47.7degree soundstage is supported without spatial aliasing, which is not too bad. Again, when pushing the angleofarrival limits, only certain frequency bands near the top of each octave band start to spatially alias, and the spectrum as a whole is likely to keep the brain fusing everything together at the intended angle of arrival.
To reduce spatial aliasing or increase the stage width of this HOP, it is tempting to shift the operating range of each speaker down by some fraction of an octave, thereby increasing the stagewidth while making each speaker a less efficient radiator, a better approximation to a point source, and a finer spatial sample within the array. Decreasing driver diameter is of course equivalent to decreasing the frequency range over which it operates. However, since most people cannot hear frequencies near 20 kHz, it makes sense to push first on downward frequency scaling. This point is pursued further in §4.9 below.
The need to downscale the frequency band served by each speaker driver is a fundamental problem with circular drivers laid out in a row. Alternatives driver geometries include (1) overlapping driver cones (perhaps using a shared membrane across multiple drivers instead of individual cones), (2) staggered packing of two or more identical rows, as in the hexagonal mesh, or (3) square or rectangular horn drivers with contiguous exitapertures driven by longthrow pistons.
To summarize the situation from an elementary spatialsampling perspective, the maximum stage width determines the minimum wavelength seen by the array, and we need two or more drivers per wavelength for proper spatial sampling. The maximum width (diameter) of each driver must be less than half the minimum spatial wavelength seen by the array (which can be made arbitrarily large by restricting the stage angle). If wide stage angles are to be supported, then every driver must be smaller than spatial wavelengths it is emitting. They’re all approaching “simple sources” emitting omnidirectional radiation patterns.
The crossover frequency 612 Hz for the low end of the 48cm woofer in Fig. 13 is based on the preferred driving frequencies and dispersion pattern for a conventional loudspeaker of that size, as chosen in the Kenwood JL840W fourway speaker system. We see that stagewidth and spatialaliasing considerations argue for a smaller crossover frequency, perhaps even an octave lower, if adequate undistorted power can be obtained covering that range.
The nominal lowest frequency in the fifth row (612 Hz) is well above the typical subwoofer crossover range of 80–200 Hz. This means we have almost three missing octaves between this particular fiveband HOP and the ideal subwoofer taking over at 80 Hz. Even downscaling the frequency bands by an octave (discussed above), two missing octaves remain. Raising the subwoofer cutoff from 80 to 160 gets one more, leaving almost a oneoctave gap, which can be reduced further from either side. A sixth row in the HOP would call for a 96 cm speaker diameter (38 inches or 3.15 feet), which does not appear to be practical using conventional driver technology. We therefore presently allow the bottom row to go undersampled, and crossfade over to VBAP signals, or even simple stereo using the left speakers as a left channel and the right speakers as a right channel (panning each virtual source accordingly in the stereo field). Fewer than speakers can be used for wider separation but less power. Experiments are needed to determine best practices in this (undersampled) frequency range. At frequencies immune to spatial aliasing, we expect to be limited only by driver quality and noise.
4.8 Upper, Middle, and Lower HOPs
Figure 14 illustrates a simple extension from one to three multiband linearrays for the purpose of providing three elevation levels. In the vertical dimension, it can be driven as threechannel VBAP, or stereo with a center channel. While azimuth perception is accurate to degree straight ahead, elevation perception is an order of magnitude less accurate, so a cruder sampling of the vertical axis is appropriate.
4.9 Nominal Design Guidelines
In general, we want

Subwoofer(s) up to Hz

VBAP from to Hz, or some simple stereo reduction

HOPs from to , where depends on listener age, etc.
Crossover frequency is normally set near the lower limit of stereo perception, which we are taking to be around 80 Hz (§2.12 on page 2.12).
Crossover frequency is chosen to be the lowest frequency at which the lowest octave of the HOP is adequately sampled, following the analysis of §2 on page 2. For the fiveband HOP example, it is Hz (§4.7 on page 4.7). From a spatial sampling point of view, depends on listening geometry, particularly the stagewidth, minimumsourcedistance, and listenerdistance from the array.
The setting for the upper limit can be based on how much highfrequency spatialization is desired. For older listeners,^{30}^{30}30See International Standard ISO7029:2017 (3rd edition): https://www.iso.org/standard/42916.html the top row is normally inaudible in arrays such as the previous example and in the fiveband HOP examples above, such as Fig. 13 on page 13. Therefore, the top row can either be simply omitted for older listeners, or, as a compromise preserving at least the presence of the high end for the occasional younger listener, replaced by conventional stereo tweeters (twochannel VBAP), or any number of tweeters in an undersampled VBAP row, or coaxial mounts along any lower row.
4.10 An Equilateral Triangle Design Example
Suppose we build a linear array forming an equilateral triangle with the centered listener, as in typical stereo speaker placements. Denote the sidelength by . The listener distance is then from the array. Suppose we want to sit back three meters from the array for normal listening. Then , which is about 7 HOPs across, following the example of Fig. 13. To reduce edge effects, we can add one more HOP on each side to obtain nine HOPs, as shown in Fig. 14. As discussed in §4.6 on page 4.6, the fiveband HOP places the crossover to VBAP/stereo at Hz, or lower as discussed in §4.7 above, and the subwoofer nominally crosses over at 80 Hz, or wherever the VBAP/stereo band lower limit may be.
4.10.1 Stereo As Two Spatial Samples
It is interesting to consider the woofer upper frequency in this example from a spatial sampling point of view. Suppose for simplicity that stereo will be used immediately above . The array width is , and for simplicity, consider the effective stereo speaker separation to be when half the speakers comprise the left channel and the other half comprise the right, or when using only the extreme left and right drivers, or various choices in between. Considering the stereo pair as two soundfield samples, stereo should take over for the subwoofer before the wavelength shrinks to (critical spatial sampling). The formula for is then
In the above example with , we obtain Hz, which includes the desirable 80 Hz setting. In the much smaller fourway Huygens Array of Fig. 11 on page 11, the stereo speaker separation (center to center) is close to meter, where we obtain the range Hz, which is quite high compared to normal subwoofer crossovers. It is apparent that in typical listening geometries the subwoofer cutoff due to perceptual limits is comparable to that obtained by considering the stereo speakers as a pair of spatial samples needing to be outside each other’s “highcorrelation zone.”^{31}^{31}31A singlefrequency soundfield is always highly correlated at distances less than a quarter wavelength or so, even when it is randomly constructed as a sum of plane waves from all directions with random phases (a “diffuse field” [Pierce, 1989, Beranek, 1986, Smith, 2010]). Even in a richly reverberant environment, one can imagine “correlation bubbles” on the order of the wavelength at each frequency.
5 Conclusions and Future Work
We have approached the problem of soundfield synthesis from the point of view of spatial sampling theory. Starting with the simplest case of farfield sources, we derived PlanewaveBased Angle Panning (PBAP) which can be considered both a specialcase of wave field synthesis (WFS) for distant sources without reverberation, and a path for enlarging the “sweetspot” in VectorBased Amplitude Panning (VBAP). Linear, planar, and more generally distributed arrays were considered, with the most general case termed Huygens Arrays (HA). Considerations of basic sampling theory led to straightforward practical guidelines on how to use the arrays, such as determining speaker spacing and size, speaker radiation pattern (interpolation kernel), maximum source angle, and how closely a source can approach the array. Finally, various multiband systems were considered, allowing each driver to focus on a particular band, such as an octave band. Multiband linear arrays, for use with conventional stereo and/or subwoofer(s) at low frequencies, were found to be particularly attractive and relatively practical.
Interesting next steps include studying the overlapadd properties of the near and farfield radiation patterns for circular and other commonly used drivers. The SFS software (linked in Appendix A) could be adapted for this purpose. The simple case of pointsource drivers is given a start in Appendix B.
The shape of a spatial sample (speaker radiation pattern) can be optimized for various purposes. One is to make it easy to produce using conventional drivers, such as the simple matrix of square pistons considered in §3.8. Another is to optimize farfield considerations, such as minimizing crosstalk between angular directions (also considered in §3.8). Other optimization criteria can be imagined. For example, if the center of the array at a particular listening distance is given special status, then the array can be optimized for that region in various ways, including minimizing the spatial width of the radiation pattern and thence array truncation effects.
6 Acknowledgments
Thanks to HyungSuk Kim for fruitful discussions on the topic of this paper, and to Madeline Huberth for reporting errata and providing helpful suggestions.
Appendix A Relevant Software
a.1 Soundfield Synthesis (SFS) in Matlab
The Soundfield Synthesis (SFS) Toolbox for Matlab [Wierstorf and Spors, 2012, Ahrens, 2012]:
https://github.com/sfstoolbox/sfs
Note that there are many useful pointers to references in the software documentation.^{32}^{32}32https://sfs.readthedocs.io/en/3.2/problem/
a.2 PBAP in Faust, Integer Delays
Below is the initial interactive test program in the Faust language for integerdelay PBAP on the eightchannel array described in §2.15 on page 2.15. Integer delays are extremely efficient computationally, but assume fixed source positions. The lowpass filter is provided to facilitate experiments with the perception of spatial aliasing at various cutoffs.
N = 8; // number of channels (speakers) MAXDELAY = 1024; // maximum delayline length for each channel lpf = lowpass(3); // Butterworth lowpass filter to use (filter.lib) declare name "Far Field Wave Field Synthesis (FFWFS) Simple Tests"; declare author "Julius O. Smith (jos at ccrma.stanford.edu)"; import("oscillator.lib"); // saw2, pink_noise import("filter.lib"); // lowpass, smooth, ... // GUI: level = hslider("v:FFWFS/[0]Level (dB)", 10, 70, 10, 0.1); del = hslider("v:FFWFS/[1]Interspeaker Delay (samples)",0,misd,misd,1) with { misd = int(0.5*MAXDELAY/N);}; // max interspeaker delay freq = hslider("v:FFWFS/[2]Frequency (Hz)",440,20,10000,1); cutoff = hslider("v:FFWFS/[3]Lowpass Cutoff (Hz)",1000,20,20000,1); nsw = checkbox("v:FFWFS/[4]Sawtooth (instead of Pink Noise)"); // Signal Processing: amp = level : db2linear : smooth(0.999); signal = _ + select2(nsw,pink_noise,saw2(freq)) : lpf(cutoff); ffwfs(i) = delay(MAXDELAY,int(MAXDELAY/2)i*del); process = signal * amp <: par(i,N,ffwfs(i));
Appendix B Speakers as Spatial Samples
The polar pattern for a microphone or loudspeaker is its gain along a circle of constant radius away from the diaphragm/driver. For a sphericalwave “pointsource”, the polar pattern is simply a constant at each radius, e.g., , where denotes the pressurescaling at and denotes the distance from the center of the source.
Since our speaker arrays are typically flat, we need to calculate a slice through the polar pattern along a listening line (or plane) which we will take to be parallel to the array and to the axis, as shown in Fig. 15. The polarpattern slice is then be considered as one sample (interpolation kernel) used to reconstruct the soundfield at a distance from the array.
b.1 PointSource Speakers
Consider a spherical secondary source (“point source”) located at as shown in Fig. 15. Then the pressure amplitude along the listening line for all is given by
(9) 
where denotes the pressure amplitude observed from the point source at , and denotes the angle of the line from the source center to the linearray point at , i.e., . This polarpattern slice is plotted in Fig. 16 for the triangular array+listener geometry such as described in §4.10 on page 4.10, with . The top plot shows over a 12 m width centered about a 6 m wide array, and the bottom plot shows for , thereby covering the entire axis containing array.
In addition to the gain variation in Eq. (9), we also have phase effects that we don’t have when keeping radius constant and only varying . The envelope in Fig. 16 is due only to spherical spreading loss according to . The difference in propagation distance between and is
(10) 
The resulting frequencyresponse is plotted in Fig. 17 for 1 kHz (wavelength about a foot (1.126 ft)) and , where is the wavenumber (spatial radian frequency), which has made its appearance for the first time in these formulas.^{33}^{33}33For a definition of “frequency response” and related terms, see, e.g., [Smith, 2007a].
For small , we can use the approximation
b.2 OverlapAdd of PointSource Frequency Responses
The listening point receives a sum of the radiating pointsources along the array. For an incident planewave traveling toward the listener, the sources are all “in phase”, emitting identical signals. Let’s focus on two adjacent sources for this case.
Figure 18 shows the magnitude of the sum of two identical pointsources separated by half a wavelength, and observed along a listening line one wavelength away from the twosource line array. Also shown is an overlay of the real parts of the component pointsources being summed, as in Fig. 17. Notice how the pressure sums coherently near the center, but largely cancels far away from the center, especially to the far left and right where the halfwavelength spacing of the sources leads to almost complete cancellation.
Consider now a larger array spanning 6 meters centered within a 12 m
listening line, again one wavelength
away. Figure 19(a) shows the OverLapAdd
(OLA) result for the pointsource frequencyresponse
of Fig. 17, again using sourceseparation
,
corresponding to “critical sampling” for extreme stage angles
at 1 kHz. The Gibbs phenomenon is highly visible, strongly
suggesting a tapering window to reduce array truncation effects
(§2.14). Tapering sinusoidally to zero (Hann halfwindows) over
3 wavelengths (6 source points) produces the result shown
in Fig. 19(b). The normalized driving
amplitude from each source is indicated by a small black
circle.^{34}^{34}34Matlab for these figures is available
at
https://ccrma.stanford.edu/~jos/huygens/matlab.tgz
While spherical waves do not give perfectreconstruction at halfwavelength overlapadd, one wavelength out, we see from the dB plots in Fig. 19 that the error is quite small in audio terms. The “ripple” stays small down to well under wavelength before source proximity effects increase the ripple significantly. Note, incidentally, that we are ignoring the reactive “mass” component of the pointsource field near its center, which is nonpropagating [Morse and Ingard, 1968].
The ripple reduces as we get farther out from the array () because the interpolation kernels expand, giving more overlap. (Increase the number of sources by the same factor to see this more clearly without additional arraytruncation error.) Since we normally listen to arrays at some distance away, we see that a larger issue than sampling density is array extent; looking in the direction a wave is coming from, we should see plenty of samples surrounding the point that our “look direction” intersects along the array. The commonly used surrounding ring/sphere architectures (often set up for ambisonics or VBAP) are excellent in this respect: Every look direction has samples uniformly about/around it. Linear/planar arrays are at a disadvantage because they must be truncated or windowed, limiting the virtual stage view.
References
 Ahrens [2012] J. Ahrens. Analytic Methods of Sound Field Synthesis. Springer, New York, 2012. ISBN 9783642257438. doi: 10.1007/9783642257438. URL https://www.springer.com/gp/book/9783642257421.
 Beranek [1986] L. L. Beranek. Acoustics. American Institute of Physics, for the Acoustical Society of America, http:http://asa.aip.org/publications.html//asa.aip.org/publications.html, 1986. 1st ed. 1954.
 Berkhout [1988] A. J. Berkhout. A holographic approach to acoustic control. Journal of the Audio Engineering Society, 36(12):977–995, 1988.
 Berkhout et al. [1993] A. J. Berkhout, D. de Vries, and P. Vogel. Acoustic control by wave field synthesis. Journal of the Acoustical Society of America, 93:2764–2778, 1993.
 Blauert [1997] J. Blauert. Spatial Hearing. The Psychophysics of Human Sound Localization. MIT Press, Cambridge, MA, USA, 1997.
 Cooper [1972] T. Cooper, Duane H.; Shiga. Discretematrix multichannel stereo. J. Audio Eng. Soc, 20(5):346–360, 1972. URL http://www.aes.org/elib/browse.cfm?elib=2070.
 Firtha [2018] G. Firtha. Unified Wave Field Synthesis Framework with Application for Moving Virtual Sources  DRAFT. Budapest University of Technology and Economics, Hungary, Febuary 7, 2018. Accessed October 19, 2019 from http:http://last.hit.bme.hu/download/firtha/WFS_system/Firtha_Dissertation_DRAFT.pdf//last.hit.bme.hu/download/firtha/WFS_system/Firtha_Dissertation_DRAFT.pdf.
 Fletcher [1934] H. Fletcher. Auditory perspective — basic requirements. Electrical Engineering, 53(1):9–11, Jan. 1934. doi: 10.1109/EE.1934.6540356.
 Franck [2008] A. Franck. Efficient algorithms and structures for fractional delay filtering based on lagrange interpolation. Journal of the Audio Engineering Society, 56(12):1036–1056, Dec. 2008. URL http://www.aes.org/elib/browse.cfm?elib=14647.
 Franck [2011] A. Franck. Efficient Algorithms for Arbitrary Sample Rate Conversion with Application to Wave Field Synthesis. PhD thesis, Technischen Universität Ilmeneau, Nov. 2011.
 Gerzon [1974] M. A. Gerzon. Surroundsound psychoacoustics. Wireless World, 80:483–486, Dec. 1974. available online from The Gerzon Archive at http:http://www.audiosignal.co.uk//www.audiosignal.co.uk.
 Gerzon [1985] M. A. Gerzon. Ambisonics in multichannel broadcasting and video. Journal of the Audio Engineering Society, 33(11):859–871, 1985.
 Goodman [2005] J. W. Goodman. Introduction to Fourier Optics, 3rd Edition. Roberts & Company, Englewood, CO, 2005. URL https://books.google.com/books?id=ow5xs_Rtt9AC.
 Kirkup [2007] S. Kirkup. The Boundary Element Method in Acoustics, volume 8. World Scientific Publishing, 01 2007. ISBN 0953403106. doi: 10.1016/S0218396X(00)000169. URL http://clok.uclan.ac.uk/7333/1/tbemia07.pdf.
 Miller [1991] D. A. B. Miller. Huygens’s wave propagation principle corrected. Optics Letters, 16:1370–2, 1991. doi:10.1364/OL.16.001370https://www.osapublishing.org/ol/abstract.cfm?uri=ol16181370.
 Morse and Ingard [1968] P. M. Morse and K. U. Ingard. Theoretical Acoustics. McGrawHill, New York, 1968.
 Nyquist [1928] H. Nyquist. Certain topics in telegraph transmission theory. Transactions of the American Institute of Electrical Engineers, 47(2):617–644, April 1928. doi: 10.1109/TAIEE.1928.5055024.
 Pierce [1989] A. D. Pierce. Acoustics. American Institute of Physics, for the Acoustical Society of America, 1989. http:http://asa.aip.org/publications.html//asa.aip.org/publications.html.
 Pulkki [1997] V. Pulkki. Virtual sound source positioning using vector base amplitude panning. Journal of the Audio Engineering Society, 45(6):456–, June 1997.
 Pulkki [2001] V. Pulkki. Tech. rep. 62. spatial sound generation and perception by amplitude panning techniques. Technical report, Lab. of Acoustics and Audio Signal Processing, Helsinki University of Technology, 2001. http:http://lib.tkk.fi/Diss/2001/isbn9512255324/isbn9512255324.pdf//lib.tkk.fi/Diss/2001/isbn9512255324/isbn9512255324.pdf.

Pulkki [2017]
V. Pulkki, editor.
Parametric Time‐Frequency Domain Spatial Audio
. John Wiley and Sons, Inc., New York, 2017. URL https://onlinelibrary.wiley.com/doi/book/10.1002/9781119252634.  Smith [2007a] J. O. Smith. Introduction to Digital Filters with Audio Applications. https:https://ccrma.stanford.edu/ jos/filters///ccrma.stanford.edu/~jos/filters/, Sept. 2007a. online book.
 Smith [2007b] J. O. Smith. Mathematics of the Discrete Fourier Transform (DFT), with Audio Applications, Second Edition. https:https://ccrma.stanford.edu/ jos/mdft///ccrma.stanford.edu/~jos/mdft/, Apr. 2007b. online book.
 Smith [2010] J. O. Smith. Physical Audio Signal Processing. https:https://ccrma.stanford.edu/ jos/pasp///ccrma.stanford.edu/~jos/pasp/, Dec. 2010. ISBN 9780974560724. online book.
 Smith [2011] J. O. Smith. Spectral Audio Signal Processing. https:https://ccrma.stanford.edu/ jos/sasp///ccrma.stanford.edu/~jos/sasp/, Dec. 2011. ISBN 9780974560731. online book.
 Smith and Gossett [1984] J. O. Smith and P. Gossett. A flexible samplingrate conversion method. In Proc. 1984 Int. Conf. Acoustics, Speech, and Signal Processing (ICASSP84), San Diego, volume 2, pages 19.4.1–19.4.2, New York, Mar. 1984. IEEE Press. expanded tutorial and associated free software available at the Digital Audio Resampling Home Page: https:https://ccrma.stanford.edu/ jos/resample///ccrma.stanford.edu/~jos/resample/.
 Smith et al. [2002] J. O. Smith, S. Serafin, J. Abel, and D. Berners. Doppler simulation and the leslie. In Proceedings of the COSTG6 Conference on Digital Audio Effects (DAFx02), Hamburg, Germany, pages 13–20, September 26 2002. https://ccrma.stanford.edu/~jos/doppler/.
 Steinberg and Snow [1934] J. C. Steinberg and W. B. Snow. Auditory perspective — physical factors. Electrical Engineering, 13(2):12–17, Jan. 1934. doi: 10.1002/j.15387305.1934.tb00661.x.
 Välimäki [1995] V. Välimäki. DiscreteTime Modeling of Acoustic Tubes Using Fractional Delay Filters. PhD thesis, Report no. 37, Helsinki University of Technology, Faculty of Electrical Engineering, Laboratory of Acoustic and Audio Signal Processing, Espoo, Finland, Dec. 1995. http:http://www.acoustics.hut.fi/ vpv/publications/vesa_phd.html//www.acoustics.hut.fi/~vpv/publications/vesa_phd.html.
 Wierstorf and Spors [2012] H. Wierstorf and S. Spors. Sound Field Synthesis Toolbox. In 132nd Convention of the Audio Engineering Society, 2012.