# VECTOR-BASE AMPLITUDE PANNING LIBRARY

## Contents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % 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

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.