beamingExample.scd

/*

Example code illustrating matrices suitable for decomposing a soundfield into two parts:

a) beam
b) residual, aka null

Three classic Ambisonic beam patterns are available:

1) planewave
2) maximum energy
3) minimum energy

The beam is re-encoded, as is the residual (null). In other words, both the resulting beaming and nulling matrices return B-format signals.
Summing the beam and null will return the original soundfield!

Joseph Anderson, 2016

*/


// Version I: Using the ATK's A-format encoder & decoder
(
var weightDict, k, patternDict;
var orientation;
var totalA, beamA;
var totalB, beamB, nullB;
var beamDecoding;

weightDict = Dictionary.with(*[
	'velocity' -> 'dec',    // planewave
	'energy' -> 'uns',     // maximum energy
	'controlled' -> 'car'  // minimum energy
]);

patternDict = Dictionary.with(*[
	'velocity' -> 0.75,                // planewave
	'energy' -> ((3-sqrt(3))/2),  // maximum energy
	'controlled' -> 0.5              // minimum energy
]);

// choose a beam shape, named in terms of decoder k
k = 'velocity';
// k = 'energy';
// k = 'controlled';

// choose an orientation: pick 'fbd' or 'fbu' for beam at [0, 0]
orientation = 'fbu';

// design the A-format decoding matrix
totalA = FoaDecoderMatrix.newBtoA(orientation, weightDict.at(k)).matrix;
beamA = Matrix.newClear(totalA.rows, totalA.cols);
beamA.putRow(0, totalA.getRow(0));

// design the A-format encoding matrix
totalB = FoaEncoderMatrix.newAtoB(orientation, weightDict.at(k)).matrix;

// design the beam matrix
beamB = totalB.mulMatrix(beamA);

// design the residual matrix
nullB = Matrix.newIdentity(beamB.size) - beamB;

// design the beam decoding matrix
beamDecoding = FoaDecoderMatrix.newMono(0, 0, patternDict.at(k)).matrix;


// post the resulting matrices
"-----------------".postln;
"Beam shape : ".post;
k.postln;

"-----------------".postln;
"Beaming ".post;
beamB.postln;

"-----------------".postln;
"Nulling ".post;
nullB.postln;
"-----------------".postln;
"".postln;

// FYI compare to weights for just decoding
"-----------------".postln;
"Beam Decoding ".post;
beamDecoding.postln;
"".postln;
)


// Version II: Using the ATK's diametric decoder
(
var k, patternDict;
var directions;
var totalA, beamA;
var totalB, beamB, nullB;
var beamDecoding;

patternDict = Dictionary.with(*[
	'velocity' -> 0.75,                // planewave
	'energy' -> ((3-sqrt(3))/2),  // maximum energy
	'controlled' -> 0.5              // minimum energy
]);

// choose a beam shape, named in terms of decoder k
k = 'velocity';      // planewave
// k = 'energy';       // maximum energy
// k = 'controlled';  // minimum energy

// specify 1/2 of directions, on +X, +Y, +Z
// the first beam is on +X,  at [0, 0]
directions = [[0.0, 0.0], [pi/2, 0.0], [0.0, pi/2]];

// design the A-format decoding matrix
totalA = FoaDecoderMatrix.newDiametric(directions, k).matrix;
beamA = Matrix.newClear(totalA.rows, totalA.cols);
beamA.putRow(0, totalA.getRow(0));

// design the A-format encoding matrix
totalB = totalA.pseudoInverse;

// design the beam matrix, and rescale
beamB = totalB.mulMatrix(beamA);
beamB = beamA.rows / beamA.cols * beamB;

// design the residual matrix
nullB = Matrix.newIdentity(beamB.size) - beamB;

// design the beam decoding matrix
beamDecoding = FoaDecoderMatrix.newMono(0, 0, patternDict.at(k)).matrix;


// post the resulting matrices
"-----------------".postln;
"Beam shape : ".post;
k.postln;

"-----------------".postln;
"Beaming ".post;
beamB.postln;

"-----------------".postln;
"Nulling ".post;
nullB.postln;
"-----------------".postln;
"".postln;

// FYI compare to weights for just decoding
"-----------------".postln;
"Beam Decoding ".post;
beamDecoding.postln;
"".postln;
)


// Version III: Direct beam encoder from planewave prototype - decoder from virtual mic decoder
(
var k, patternDict, orderGains, normOrderGains;
var order;
var planewaveEncoding;
var beamEncoding;
var beamDecoding;
var beamB, nullB;

// FOA
order = 1;

patternDict = Dictionary.with(*[
	'velocity' -> 0.75,                // planewave
	'energy' -> ((3-sqrt(3))/2),  // maximum energy
	'controlled' -> 0.5              // minimum energy
]);

// choose a beam shape, named in terms of decoder k
k = 'velocity';
// k = 'energy';
// k = 'controlled';

// calculate order gains
orderGains = [1-patternDict.at(k), patternDict.at(k)];
orderGains = orderGains / orderGains.at(0);
orderGains = orderGains / Array.fill(order+1, { arg i; 2*i+1 });

// calculate normalized order gains
normOrderGains = orderGains.reciprocal * ((Array.fill(order+1, { arg i; 2*i+1 }) * orderGains).sum / 2.pow(order+1));


// form a planewave encoding matrix,  at [0, 0]
planewaveEncoding = FoaEncoderMatrix.newDirection.matrix;
planewaveEncoding = planewaveEncoding.addRow([0]);   // make 3D

// apply normalized order gains to form beam encoding matrix,  at [0, 0]
beamEncoding = Array.fill(4, {
	arg i;
	normOrderGains.at(i.sqrt.floor) * planewaveEncoding.at(i, 0);
}).reshapeLike(planewaveEncoding);


// design the beam decoding matrix (normalized)
beamDecoding = FoaDecoderMatrix.newMono(0, 0, patternDict.at(k)).matrix;


// design the beam matrix
beamB = beamEncoding.mulMatrix(beamDecoding);

// design the residual matrix
nullB = Matrix.newIdentity(beamB.size) - beamB;

// post the resulting matrices
"-----------------".postln;
"Beam shape : ".post;
k.postln;

"-----------------".postln;
"Beaming ".post;
beamB.postln;

"-----------------".postln;
"Nulling ".post;
nullB.postln;
"-----------------".postln;
"".postln;

// FYI compare to weights for just decoding
"-----------------".postln;
"Beam Decoding ".post;
beamDecoding.postln;
"".postln;

// FYI compare to weights for just encoding
"-----------------".postln;
"Beam Encoding ".post;
beamEncoding.postln;
"".postln;
)


// Version IV: Direct beam encoding and decoding from planewave prototype
(
var k, patternDict, orderGains, normEncOrderGains, normDecOrderGains;
var order;
var planewaveEncoding;
var beamEncoding;
var beamDecoding;
var beamB, nullB;

// FOA
order = 1;

patternDict = Dictionary.with(*[
	'velocity' -> 0.75,                // planewave
	'energy' -> ((3-sqrt(3))/2),  // maximum energy
	'controlled' -> 0.5              // minimum energy
]);

// choose a beam shape, named in terms of decoder k
k = 'velocity';
// k = 'energy';
// k = 'controlled';

// calculate order gains
orderGains = [1-patternDict.at(k), patternDict.at(k)];
orderGains = orderGains / orderGains.at(0);
orderGains = orderGains / Array.fill(order+1, { arg i; 2*i+1 });

// calculate normalized order gains
normEncOrderGains = orderGains.reciprocal * ((Array.fill(order+1, { arg i; 2*i+1 }) * orderGains).sum / 2.pow(order+1));

normDecOrderGains = orderGains / ((Array.fill(order+1, { arg i; 2*i+1 }) * orderGains).sum);
normDecOrderGains = normDecOrderGains * [2, 3];  // (MaxN to N3D).pow(2)


// form a planewave encoding matrix,  at [0, 0]
planewaveEncoding = FoaEncoderMatrix.newDirection.matrix;
planewaveEncoding = planewaveEncoding.addRow([0]);   // make 3D

// apply normalized order gains to form beam encoding matrix,  at [0, 0]
beamEncoding = Array.fill(4, {
	arg i;
	normEncOrderGains.at(i.sqrt.floor) * planewaveEncoding.at(i, 0);
}).reshapeLike(planewaveEncoding);

// apply normalized order gains to form beam decoding matrix,  at [0, 0]
beamDecoding = Array.fill(4, {
	arg i;
	normDecOrderGains.at(i.sqrt.floor) * planewaveEncoding.at(i, 0);
}).reshapeLike(planewaveEncoding).flop;


// design the beam matrix
beamB = beamEncoding.mulMatrix(beamDecoding);

// design the residual matrix
nullB = Matrix.newIdentity(beamB.size) - beamB;

// post the resulting matrices
"-----------------".postln;
"Beam shape : ".post;
k.postln;

"-----------------".postln;
"Beaming ".post;
beamB.postln;

"-----------------".postln;
"Nulling ".post;
nullB.postln;
"-----------------".postln;
"".postln;

// FYI compare to weights for just decoding
"-----------------".postln;
"Beam Decoding ".post;
beamDecoding.postln;
"".postln;

// FYI compare to weights for just encoding
"-----------------".postln;
"Beam Encoding ".post;
beamEncoding.postln;
"".postln;
)