VECTOR-BASE AMPLITUDE PANNING LIBRARY

EXAMPLE 1: Triangulation
EXAMPLE 2: VBAP Gains
EXAMPLE 3: Spreading by amplitude panning
EXAMPLE 4: Gain Tables
EXAMPLE 5: Moving source
EXAMPLE 6: VBAP filters with frequency-dependent P-value
REFERENCES

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Archontis Politis, 2015
%   Department of Signal Processing and Acoustics, Aalto University, Finland
%   archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This is a compact Matlab/Octave library implementing vector-base amplitude panning (VBAP) [ref.1], VBAP-base spreading of a panned source [ref.2&3], and Multiple-direction amplitude panning (MDAP) [ref.3]. A function implementing the variant vector-base intensity panning [ref.4] is also included. Recently, both VBAP and VBIP have been additionaly used for the design of robust ambisonic decoding matrices, see [ref.5&6].

The code is written by Archontis Politis, except the core vbap() function contributed by Ville Pulkki, with small modifications by Archontis Politis. The following code examples are meant to give a quick idea how to use the library for common operations in amplitude panning, such as triangulation of a 3D loudspeaker setup into loudspeaker triplets, spreading of a panned source, construction of panning gain tables, and panning a moving source in a real-time block processing context.

The library contains the following main functions:

findLsPairs:    find sorted loudspeaker pairs from loudspeaker
                directions (for 2D layouts)
findLSTriplets: find valid loudspeaker triangles from loudspeaker
                directions (for 3D layouts)
invertLsMtx:    precompute inversion of matrix of loudspeaker triplets
                or pairs, for use in VBAP
getSpreadSrcDirs:   get auxiliary source directions around panning
                    direction, for source spreading and MDAP
vbap:   Return VBAP panning gains for multiple panning directions, with
        spread control if needed

Additionaly:

plotTriangulation:  Plots the loudspeaker triangulated mesh
getGainTable:   Construct a look-up VBAP gain table of VBAP for a
                specified regular grid
vbip:   Similar to VBAP, but implementing its energy-based variant
        (see [ref.4])
getPValueResponse:  Returns VBAP frequency-dependent normalization
                    values, for approximate flat perceived response of
                    a panned source in dry playback environments
                    (see [ref.7])

For any questions, comments, corrections, or general feedback, please contact archontis.politis@aalto.fi

EXAMPLE 1: Triangulation

Triangulation of 3D setups is done using the convex hull of the loudspeaker points on the unit sphere (equivalent to Delaunay triangulation)

% FULL SPHERE EXAMPLE OF 29 SPEAKERS
ls_dirs_full = [-18 -54 -90 -126 -162 -198 -234 -270 -306 -342  0    -72  -144 -216 -288 -45  -135 -225 -315 0;
                 0   0   0   0    0    0    0    0    0    0   -10   -10  -10  -10  -10   45   45   45   45  90]';
% Delaunay triangulation
[~, ls_full] = findLsTriplets(ls_dirs_full);

OMITTING LARGE TRIANGLES:

In many cases (like the setup above), it is advantageous to remove very large triangles (with aperture>~100deg) that are going to result in unstable localization, such as the ones below -10deg elevation in the example configuration above (see next plot). If you use this, always check visually the triangulation to be sure that you haven't omitted by mistake some triangle that you would like to keep. An alternative strategy is to introduce virtual loudspeakers that their channels are then either discarded or mixed with the actual ones.

% omit large triangles
OMIT_LARGE_TRI = 1;
large_tri_aperture = 100; % largest allowed triangle side in degrees
[~, ls_full_omit] = findLsTriplets(ls_dirs_full, OMIT_LARGE_TRI, large_tri_aperture);

% plot triangulations
figure
subplot(121)
plotTriangulation(ls_full)
title('full'), set(findall(gca, 'type', 'text'), 'visible', 'on', 'fontsize',16) % make title visible
view(60,-20), zoom(2)
subplot(122)
plotTriangulation(ls_full_omit)
title(['full' char(10) 'with large triangles discarded']), set(findall(gca, 'type', 'text'), 'visible', 'on', 'fontsize',16)
view(60,-20), zoom(2)
h = gcf; h.Position(3:4) = 2*h.Position(3:4);

PARTIAL SETUP OF LOUDSPEAKERS

Partial setups can be domes, or setups around a screen, etc. In this case loudspeaker triplets should be discarded when their normals do not point outwards. This is handled automatically by the findLsTriplets() function.

% define a partial frontal setup
ls_dirs_partial = [-80 -45 0 45 80 -60 -30 30 60;
                    0   0  0 0  0   60  60 60 60]';
% triangulate directly, without considering invalid triangles
ls_dirs_rad = ls_dirs_partial*pi/180;
[tempx, tempy, tempz] = sph2cart(ls_dirs_rad(:,1), ls_dirs_rad(:,2), 1);
ls_part_invalid.vert = [tempx, tempy, tempz];
ls_part_invalid.faces = convhulln(ls_part_invalid.vert);
% trinagulate and discard invalid faces
[~, ls_part_valid] = findLsTriplets(ls_dirs_partial);

% plot triangulation
figure
subplot(121)
plotTriangulation(ls_part_invalid)
title('partial - including invalid triangles'), set(findall(gca, 'type', 'text'), 'visible', 'on', 'fontsize',16)
view(3), zoom(1.5)
subplot(122)
plotTriangulation(ls_part_valid)
title('partial - after omitting invalid triangles'), set(findall(gca, 'type', 'text'), 'visible', 'on', 'fontsize',16)
view(3), zoom(1.5)
h = gcf; h.Position(3:4) = 1.5*h.Position(3:4);

EXAMPLE 2: VBAP Gains

The steps to obtaining the VBAP gains for a set of directions are:

% a) define loudspeaker setup ls_dirs = [azi1 azi2 ... aziK] for 2D
%    or  ls_dirs = [azi1 elev1; azi2 elev2; ...; aziK elevK] for 3D
ls_dirs = [30 -30 0 110 -110]; % define a 2D 5.0 setup in degrees
%
% b) find valid loudspeaker pairs or triplets:
%    findLsPairs(ls_dirs) for 2D, or
%    findLsTriplets(ls_dirs) for 3D
ls_groups = findLsPairs(ls_dirs);
%
% c) compute inverse matrices for loudspeaker pairs or triplets, needs to
%    be done once for any panning direction
layoutInvMtx = invertLsMtx(ls_dirs, ls_groups);
%
% d) compute vbap gains for the required source directions, in degrees
src_dirs2D = (0:359)'; % 2D panning directions at every 1deg
gains2D = vbap(src_dirs2D, ls_groups, layoutInvMtx); % compute vbap gains

% Plot panning gains
figure
polar(src_dirs2D*ones(1,5)*pi/180, gains2D)
title('5.0 VBAP gains')

% Repeat process for an 11.0 3D setup
ls_dirs = [30 -30 0 120 -120 90 -90 45 -45 135 -135; 0 0 0 0 0 0 0 45 45 45 45]';
[ls_groups, layout] = findLsTriplets(ls_dirs); % return also triangulation mesh for plotting
figure, subplot(121), plotTriangulation(layout); view(50,30), zoom(2) % plot triangulation
layoutInvMtx = invertLsMtx(ls_dirs, ls_groups);
% Generate a regular 2D grid of panning directions covering the sphere
aziRes = 5;
elevRes = 5;
[Elev, Azi] = meshgrid(-90:elevRes:90, 0:aziRes:360);
src_dirs3D = [Azi(:) Elev(:)];
% Get VBAP gains
gains3D = vbap(src_dirs3D, ls_groups, layoutInvMtx);

% Plot panning gains
ls_num = size(ls_dirs,1);
[nAzi, nElev] = size(Azi);
[X,Y,Z] = sph2cart(Azi*pi/180,Elev*pi/180,1);
subplot(122)
hold on
for nl = 1:ls_num
    gains_grid_nl = reshape(gains3D(:,nl), nAzi, nElev);
    surf(gains_grid_nl.*X,gains_grid_nl.*Y,gains_grid_nl.*Z,gains_grid_nl);
end
axis([-1 1 -1 1 -0.5 1]), axis equal
colorbar, view(50,30), zoom(2), grid
title('11.0 VBAP gains')
h = gcf; h.Position(3:4) = 1.5*h.Position(3:4);

EXAMPLE 3: Spreading by amplitude panning

Spreading of sources by means of amplitude panning has two useful applications: a)creating synthetic sounds at arbitrary directions that have perceptually some spatial extent (see [ref.2&3]), and b)creating panning gains for multichannel systems that avoid the "loudspeaker detent" effect of sources collapsing to the loudspeakers or appearing spread at certain directions and very point-like at others (see [ref.3]). This second application is useful also on the design of stable hybrid ambisonic decoders (see [ref.6]).

Spread can be controlled if an extra spread parameter is added to the vbap() function. This parameter determines the extent of the source in degrees. The spread effect is created by using auxiliary spread sources around the main panning direction. If no additional arguments are passed along the desired spread, a default of 8 auxiliary sources are used. Otherwise, an arbitrary number of auxiliary sources can be generated by an additional argument num_spread_src. For 3D cases, the spread sources are arranged on a ring around the panning direction. By default, a single ring is used. More rings can be generated by an extra argument num_spread_rings.

% plot spread source directions for a 2D setup
spread = 60;
U_spread1 = getSpreadSrcDirs(-90, spread);
spread = 100;
num_spread_src = 32;
U_spread2 = getSpreadSrcDirs(45, spread, num_spread_src);
U_circle = [cos(0:pi/100:2*pi)' sin(0:pi/100:2*pi)'];
figure
plot(U_circle(:,1), U_circle(:,2), 'k')
hold on
plot(U_spread1(:,1), U_spread1(:,2), 'ro') % spread source 1
plot(U_spread2(:,1), U_spread2(:,2), 'bo') % spread source 2
axis equal, title('2D spread sources')

% plot spread source directions for a 3D setup
spread = 60;
U_spread1 = getSpreadSrcDirs([0 0], spread);
spread = 90;
num_spread_src = 64;
U_spread2 = getSpreadSrcDirs([0 90], spread, num_spread_src);
spread = 140;
num_spread_src = 32;
num_spread_rings = 4;
U_spread3 = getSpreadSrcDirs([-135 -30], spread, num_spread_src, num_spread_rings);

[Xs, Ys, Zs] = sphere(30);
figure
surf(0.95*Xs, 0.95*Ys, 0.95*Zs, 1, 'facecolor','c')
hold on
plot3(U_spread1(:,1), U_spread1(:,2), U_spread1(:,3), 'ro') % spread source 1
plot3(U_spread2(:,1), U_spread2(:,2), U_spread2(:,3), 'bo') % spread source 2
plot3(U_spread3(:,1), U_spread3(:,2), U_spread3(:,3), 'mo') % spread source 3
view(36,15), axis equal, title('3D spread sources')
h = gcf; h.Position(3:4) = 1.5*h.Position(3:4);

% compute 2D vbap gains for a 5.0 layout, with and without spreading
ls_dirs = [30 -30 0 110 -110]; % define a 2D 5.0 setup in degrees
ls_groups = findLsPairs(ls_dirs);
layoutInvMtx = invertLsMtx(ls_dirs, ls_groups);
src_dirs2D = (0:359)';
gains_nospread = vbap(src_dirs2D, ls_groups, layoutInvMtx);
spread = 45;
gains_spread = vbap(src_dirs2D, ls_groups, layoutInvMtx, spread);

% Plot panning gains
figure
subplot(121)
polar(src_dirs2D*ones(1,5)*pi/180, gains_nospread)
title('VBAP - no spreading')
subplot(122)
polar(src_dirs2D*ones(1,5)*pi/180, gains_spread)
title('MDAP - 45deg spreading')
h = gcf; h.Position(3) = 2*h.Position(3);

EXAMPLE 4: Gain Tables

For some applications it may be advantageous to precalculate a gain table of panning gains with some required spatial resolution, and then just look-up the gains for the panning direction. This is faster than computing the gains online but it requires more memory. Gain tables can be computed with the getGainTable() function.

% 2D gain table example and look-up
azi_res = 3; % compute gains at every 3 degrees (default 1deg if not defined)
ls_dirs = [30 -30 0 110 -110];
gtable2D = getGainTable(ls_dirs, azi_res);
% look-up gains for a certain direction
azi = 97.65;
idx2D = round(mod(azi+180,360)/azi_res)+1;
gains2D = gtable2D(idx2D,:);

% 3D gain table example and look-up
azi_res = 5;
elev_res = 5; % compute gains at every 5deg azimuth and elevation (default 2deg and 5deg respectively)
N_azi = round(360/azi_res) + 1;
ls_dirs = [30 -30 0 120 -120 90 -90 45 -45 135 -135; 0 0 0 0 0 0 0 45 45 45 45]';
gtable3D = getGainTable(ls_dirs, [azi_res elev_res]);
% look-up gains for a certain direction
azi = -123.6;
elev = 37.65;
aziIndex = round(mod(azi+180,360)/azi_res);
elevIndex = round((elev+90)/elev_res);
idx3D = elevIndex*N_azi+aziIndex+1;
gains3D = gtable3D(idx3D,:);

EXAMPLE 5: Moving source

A simple example of panning a moving source with block processing (one-and-a-half horizontal rotation of a 500Hz sinusoid, at 3 seconds).

% initial parameters
fs = 48000; % samplerate
blocksize = fs/20; % (~50msec)
hopsize = blocksize/2; % panning hopsize (update the panning values twice per bocksize)
ls_dirs = [0  45 90  135 180 -135 -90 -45 0;    % define an octagon with a top
           0  0  0   0   0    0    0   0  90]'; % loudspeaker
ls_num = size(ls_dirs,1);

% define signal
sig = sin(2*pi*250*(1:3*fs)/fs)'; % 3sec of 250Hz sinewave as example
Lsig = length(sig);
Nhop = ceil(Lsig/hopsize) + 2;
padsig = [zeros(hopsize,1); sig; zeros(Nhop*hopsize - Lsig - hopsize,1)]; % zero padding
pansig = zeros(size(padsig,1), ls_num);

% define the trajectory for panning, one and a half horizontal rotation
azis = (0:(Nhop-1)-1)'*(1.5*360)/(Nhop-1);
eles = 0*azis;

% precompute VBAP triplet inversion
ls_groups = findLsTriplets(ls_dirs); % return also triangulation mesh for plotting
layoutInvMtx = invertLsMtx(ls_dirs, ls_groups);

% do panning of signal to defined trajectory with overlap-add method
counter = 1;
window = hanning(blocksize);
spread = 30;
for idx = 0:hopsize:(Nhop-2)*hopsize
    winsig = padsig(idx+(1:blocksize),1) .* window;
    azi = azis(counter);
    ele = eles(counter);
    gains = vbap([azi ele], ls_groups, layoutInvMtx, spread);
    panwinsig = winsig*gains;
    pansig(idx+(1:blocksize),:) = pansig(idx+(1:blocksize),:) + panwinsig;
    counter = counter+1;
end
% truncate loudspeaker signals to original length (omit zeropadding)
pansig = pansig(hopsize+(1:Lsig),:);

% plot original signal and panned front and front-left channel, half rotation only
figure, plot(1:Lsig/3, sig(1:Lsig/3), 1:Lsig/3, pansig(1:Lsig/3,1), '--m', 1:Lsig/3, pansig(1:Lsig/3,2), '--y');
legend('original','panned - front', 'panned - left-front'), xlabel('samples'), ylabel('amplitude')
h = gcf; h.Position(3) = 2*h.Position(3);

% write audio output to file
audiowrite('panning_example.wav', pansig, fs)

EXAMPLE 6: VBAP filters with frequency-dependent P-value

Commonly VBAP is used with a frequency-independent power normalization of the amplitude panning gains, assuming incoherent summation of the speaker channels at the listener's ears, in order to result in a constant perceived loudness with panning direction. This assumption is on average valid for most domestic listening environments and for mid-to high frequencies.

More generally, the gain normalization factor for L loudspeakers can be defined as

$N = \sum_{l=1}^{L} \sqrt[p]{ g_l^p}, \quad p\in[0,1]$

where p=1 corresponds to amplitude normalization and p=2 to standard power normalization. Amplitude normalization is more appropriate at low frequencies and in dry playback environments due to coherent summation of the loudspeaker channels. Power normalization can cause a clearly perceived bass-boosting effect in these cases.

A solution is proposed in [ref.7], where the p-value becomes frequency-dependent, with respect to a room-related parameter corresponding roughly to the direct-to-total energy ratio. The function getPValueResponse() returns these normalization values at specified frequencies, and the code below demonstrates how it can be used to construct VBAP gain filters.

fs = 40000; % sample rate
Lfilt = 512; % filter length (unnecessary long - used for smooth plots below)

f = (0:Lfilt/2)*fs/Lfilt; % frequency vector
DTT = 1; % 1 for anechoic conditions, ~0.5 for listening rooms, 0 for standard power normalization
pValue = getPValueResponse(f, DTT);
ls_dirs5 = [30 -30 0 110 -110];
gains = vbap(15, findLsPairs(ls_dirs5),invertLsMtx(ls_dirs5, findLsPairs(ls_dirs5))); % VBAP gains for source at 15deg
ls_num = length(ls_dirs5);
for nf=1:Lfilt/2+1
    pv_f = pValue(nf); % p-value for this frequency
    H_filt(nf,:) = gains./((sum(gains.^pv_f)).^(1/pv_f) * ones(1,ls_num));
end
% get impulse response
h_filt = fftshift(ifft([H_filt; H_filt(end-1:-1:2,:)]),1);

% plot the p-value curve and the magnitude response of the filter for the
% left speaker
figure
subplot(121), semilogx(f, pValue), grid, axis([100 2e4 1 2]), title(['p-value for DTT=' num2str(DTT)])
subplot(122), semilogx(f, H_filt(:,1)), grid, axis([100 2e4 0 1]), title('|H_{vbap}| for left speaker')
h = gcf; h.Position(3) = 2*h.Position(3);

Note that the experimental curves of [ref.7] are derived for listening with pairs of loudspeakers, however the authors have used them with success also in 3D loudspeaker setups.

REFERENCES

[1] Pulkki, V. (1997). Virtual Sound Source Positioning Using Vector Base Amplitude Panning. Journal of the Audio Engineering Society, 45(6), 456-466.

[2] Pulkki, V. (2000). Generic panning tools for MAX/MSP. International Computer Music Conference (ICMC), Berlin, Germany

[3] Pulkki, V. (1999). Uniform Spreading of Amplitude Panned Sources. IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA), New Paltz, NY, USA

[4] Jot, J.-M., Larcher V., Pernaux, J.-M. (1999). A comparative study of 3-D audio encoding and rendering techniques. 16th International Conference of the AES, Rovaniemi, Finland

[5] Zotter, F., Frank, M. (2012). All-Round Ambisonic Panning and Decoding. Journal of the Audio Engineering Society, 60(10), 807-820.

[6] Epain, N., Jin, C.T., Zotter, F. (2014). Ambisonic Decoding With Constant Angular Spread. Acta Acustica united with Acustica, 100(May), 928-936.

[7] Laitinen, M., Vilkamo, J., Jussila, K., Politis, A., Pulkki, V. (2014). Gain normalization in amplitude panning as a function of frequency and room reverberance. 55th International Conference of the AES. Helsinki, Finland.