CSL 5.2

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

Ambisonic.cpp
Go to the documentation of this file.
1 //
2 // Ambisonic.cpp -- Higher Order Ambisonic abstract base class for all Ambisonic effects and panners.
3 // See the copyright notice and acknowledgment of authors in the file COPYRIGHT
4 // Higher Order Ambisonic classes written by Jorge Castellanos, Graham Wakefield, Florian Hollerweger, 2005
5 //
6 // AmbisonicUnitGenerator.cpp -- Higher Order Ambisonic superclass for all Ambisonic encoded framestreams
7 // See the copyright notice and acknowledgment of authors in the file COPYRIGHT
8 // Higher Order Ambisonic classes written by Jorge Castellanos, Graham Wakefield, Florian Hollerweger, 2005
9 //
10 // Mixin superclass for the ambisonic framestream classes (e.g. HOA_Encoder, HOA_Rotator, AmbisonicDecoder).
11 // Encapsulates parameters regarding the ambisonic order of the class (including hybrid orders)
12 //
13 // See HOA_AmbisonicFramestream.h for descriptions of the the class members.
14
15 #include "Ambisonic.h"
16
17 #include <stdlib.h>
18 #include <math.h>
19 #include <string.h>
20 #include <iostream>
21
22 #include "AmbisonicUtilities.h"
23
24 #define AMBI_INVSQRT2 (1/1.414213562)
25 #define INV_SQRT2 (1/1.414213562)
26 #define AMBI_INVSQRT2 (1/1.414213562)
27
28 using std::cout;
29 using std::endl;
30
31 using namespace csl;
32
33
34 /*******************CONSTRUCTORS & DESTRUCTOR*******************/
35
36 AmbisonicUnitGenerator::AmbisonicUnitGenerator(unsigned torder) : mOrder(torder, torder) {
37 initOrder(); // initialize
38 }
39
40 AmbisonicUnitGenerator::AmbisonicUnitGenerator(unsigned hOrder, unsigned vOrder) : mOrder(hOrder, vOrder) {
41 mOrder.verticalOrder = vOrder;
42 initOrder();
43 }
44
45 AmbisonicUnitGenerator::AmbisonicUnitGenerator(AmbisonicOrder torder) : mOrder(torder) {
46 initOrder(); // initialize
47 }
48
49 AmbisonicUnitGenerator::~AmbisonicUnitGenerator() { }
50
51
52 void AmbisonicUnitGenerator::initOrder() {
53
54 setCopyPolicy(kIgnore); // This is needed so the default kCopy doesn't overide the multiple channel panning done here.
55
56 mOrder.isUniform = (mOrder.horizontalOrder == mOrder.verticalOrder); // false for hybrid orders
57 mNumChannels = orderToChannels(mOrder); // derive necessary output channels from ambisonic order
58
59 }
60
61 void AmbisonicUnitGenerator::setOrder(AmbisonicOrder torder) {
62 mOrder = torder;
63 initOrder();
64 }
65
66 /*******************UTILITY FUNCTIONS*******************/
67
68 unsigned AmbisonicUnitGenerator::channelsToUniformOrder(const unsigned channels) {
69 // M = floor(sqrt(N) - 1)
70 return (unsigned)floor(sqrt((double)channels)-1);
71 }
72
73 unsigned AmbisonicUnitGenerator::greaterOrder(const AmbisonicOrder torder) {
74 return torder.horizontalOrder > torder.verticalOrder ? torder.horizontalOrder : torder.verticalOrder;
75 }
76
77 unsigned AmbisonicUnitGenerator::orderToChannels(const AmbisonicOrder torder) {
78 // hybrid system version: N = 2*M_h + 1 + (M_v + 1)^2 - (2*M_v + 1)
79 return (2 * torder.horizontalOrder + 1) + ((unsigned)pow(torder.verticalOrder + 1.f, 2) - (2 * torder.verticalOrder + 1));
80 }
81
82 unsigned AmbisonicUnitGenerator::orderToChannels(unsigned torder) {
83 // uniform system version: N = (M+1)^2
84 return (unsigned)pow(torder + 1.f, 2);
85 }
86
87 unsigned AmbisonicUnitGenerator::orderToHorizontalChannels(const AmbisonicOrder torder) {
88 // N_h = 2*M_h + 1
89 return (2 * torder.horizontalOrder + 1);
90 }
91
92 unsigned AmbisonicUnitGenerator::orderToVerticalChannels(const AmbisonicOrder torder) {
93 // N_v = (M_v + 1)^2 - (2*M_v + 1)
94 return ((unsigned)pow(torder.verticalOrder + 1.f, 2) - (2 * torder.verticalOrder + 1));
95 }
96
97
98 /*******************CHANNEL INDEXER FUNCTIONS*******************/
99
100 void AmbisonicUnitGenerator::channelIndexer(unsigned *indexArray) {
101 AmbisonicOrder torder = mOrder;
102 unsigned *channelIndex = indexArray;
103 // how many channels it would need *if* it was a uniform order (hOrder == vOrder)
104 unsigned numChannelsGreaterOrder = orderToChannels(greaterOrder(torder));
105
106 // how many channels it actually needs for the hybrid order
107 unsigned nnumChannels = orderToChannels(torder);
108
109 unsigned n = 1;
110 unsigned currentOrder = 1;
111 unsigned hoaChannel = 1;
112 unsigned numChannelsCurrentOrder;
113
114 // fill the channel indexer with zeroes
115 for (unsigned i = 0; i < numChannelsGreaterOrder; i++)
116 channelIndex[i] = 0;
117
118 // work our way up the orders until we have exhausted all the required channels
119 while (n < nnumChannels) {
120 // horizontal indexing required?
121 if (currentOrder <= torder.horizontalOrder) {
122 // there are always 2 channels for horizontal encoding at each order
123 channelIndex[hoaChannel++] = n++;
124 channelIndex[hoaChannel++] = n++;
125 } else {
126 hoaChannel += 2; // skip them
127 }
128
129 // how many channels are required up to the current order?
130 numChannelsCurrentOrder = orderToChannels(currentOrder);
131
132 // vertical indexing required?
133 if (currentOrder <= torder.verticalOrder) {
134 // do until all required vertical channels up to the current order are exhausted
135 while(hoaChannel < numChannelsCurrentOrder) {
136 channelIndex[hoaChannel++] = n++; // set and increment
137 }
138 } else { // skip them
139 hoaChannel = numChannelsCurrentOrder;
140 }
141 // finished with the current order, increment and do again
142 currentOrder++;
143 }
144 return;
145 }
146
147
148 void AmbisonicUnitGenerator::invChannelIndexer(unsigned *indexArray) {
149
150 AmbisonicOrder torder = mOrder;
151 unsigned *channelIndex = indexArray;
152
153 // how many channels it actually needs for the hybrid order
154 unsigned nnumChannels = orderToChannels(torder);
155 channelIndex[0] = 0; // w channel
156
157 unsigned n = 1;
158 unsigned currentOrder = 1;
159 unsigned numChannelsCurrentOrder;
160
161 // work our way up the orders until we have exhausted all the required channels
162 while (n < nnumChannels) {
163 // horizontal indexing required?
164 if (currentOrder <= torder.horizontalOrder) {
165 // there are always 2 channels for horizontal encoding at each order
166 channelIndex[n] = n++;
167 channelIndex[n] = n++;
168 } else {
169 n += 2; // skip them
170 }
171 // how many channels are required up to the current order?
172 numChannelsCurrentOrder = orderToChannels(currentOrder);
173 // vertical indexing required?
174 if (currentOrder <= torder.verticalOrder) {
175 // do until all required vertical channels up to the current order are exhausted
176 while(n < numChannelsCurrentOrder) {
177 channelIndex[n] = n++; // set and increment
178 }
179 } else {
180 n = numChannelsCurrentOrder; // skip them
181 } // finished with the current order, increment and do again
182 currentOrder++;
183 }
184 return;
185 }
186
187
188 //////////////////////////////////////////////
189 /* AMBISONIC ENCODER IMPLEMENTATION */
190 //////////////////////////////////////////////
191
192 /*******************CONSTRUCTORS & DESTRUCTOR*******************/
193 AmbisonicEncoder::AmbisonicEncoder() : AmbisonicUnitGenerator() { }
194
195 AmbisonicEncoder::AmbisonicEncoder(SpatialSource &input, unsigned order) : AmbisonicUnitGenerator(order), mInputPort(NULL) {
196 setInput(input);
197 initialize();
198 }
199
200 AmbisonicEncoder::AmbisonicEncoder(SpatialSource &input, unsigned horder, unsigned vorder) : AmbisonicUnitGenerator(horder, vorder), mInputPort(NULL) {
201 setInput(input);
202 initialize();
203 }
204
205 AmbisonicEncoder::~AmbisonicEncoder() {
206 delete [] mWeights; // de-allocate memory
207
208 if (mInputPort != NULL) {
209 delete mInputPort;
210 }
211 }
212
213
214 /*******************INITIALIZE METHOD*******************/
215
216 void AmbisonicEncoder::initialize() {
217
218 // allocate memory for encoding weights
219 mWeights = new sample[mNumChannels];
220 // initialize them
221 mWeights[0] = AMBI_INVSQRT2; // W channel weight is constant
222 for (unsigned i = 1; i < mNumChannels; i++) {
223 mWeights[i] = 0.f;
224 }
225 }
226
227 /// Use to set the input to be encoded.
228 void AmbisonicEncoder::setInput(SpatialSource &input) {
229
230 if (mInputPort == NULL)
231 mInputPort = new Port(&input);
232
233 if (mInputPort->mUGen != &input) { // Make sure no setting the same input.
234 mInputPort->mUGen->removeOutput(this); // Remove the encoder (slef) from the previous input ugen.
235 mInputPort->mUGen = &input; // mInputPort->setInput(input);
236
237 input.addOutput(this); // be sure to add me as an output of the other guy
238 }
239 #ifdef CSL_DEBUG
240 logMsg("AmbisonicEncoder::setInput");
241 #endif
242 }
243
244 /*******************NEXT_BUFFER METHOD*******************/
245
246 // To get my next buffer, get a buffer from the input, and then "process" it...
247 void AmbisonicEncoder::nextBuffer(Buffer &outputBuffer, unsigned outBufNum) throw(CException) {
248
249 unsigned numFrames = outputBuffer.mNumFrames; // Block size
250 Buffer *inputBuffer;
251 outputBuffer.zeroBuffers(); // Initialize output buffer with zeros
252 #ifdef CSL_DEBUG
253 logMsg("AmbisonicEncoder::nextBuffer");
254 #endif
255 if (outputBuffer.mNumChannels < mNumChannels) // Make sure the buffer passed has enough channels
256 logMsg(kLogError, "The AmbisonicEncoder requires %d output channels, found only %d ???\n", mNumChannels, outputBuffer.mNumChannels);
257
258 // Get a pointer to the buffer of the sound source to be encoded
259 // inputBuffer = mInputPort->pullInput(numFrames);
260 // THE BLOCK BELOW WAS PULL_INPUT OF THE UGEN PORT...
261 inputBuffer = mInputPort->mBuffer; // else get the buffer
262 inputBuffer->mNumFrames = numFrames;
263 inputBuffer->mType = kSamples;
264 mInputPort->mUGen->nextBuffer(*inputBuffer); // and ask the UGen for nextBuffer()
265 inputBuffer->mIsPopulated = true;
266
267 SpatialSource *input = (SpatialSource *)mInputPort->mUGen;
268
269 // if (input->positionChanged()) {
270 // Get the encoding weights for this encoder
271 fumaEncodingWeights(mWeights, mOrder, input->azimuth(), input->elevation());
272 // }
273 // Start encoding audio from zeroth order; higher orders cascade from this function
274 for (unsigned i = 0; i < mNumChannels; i++) {
275 SampleBuffer inPtr = inputBuffer->buffer(0); // Reset the input pointer
276 SampleBuffer outPtr = outputBuffer.buffer(i); // Set a pointer to the output buffer
277
278 for (unsigned j = 0; j < numFrames; j++) {
279 // Actual audio processing: encode sound sources following the Furse-Malham encoding convention
280 *outPtr++ += (*inPtr++) * mWeights[i]; // encoding W channel
281 }
282 }
283 }
284
285
286 ////////////////////////////////////////////
287 /* AMBISONIC DECODER IMPLEMENTATION */
288 ////////////////////////////////////////////
289
290
291 /*******************CONSTRUCTORS & DESTRUCTOR*******************/
292
293 AmbisonicDecoder::AmbisonicDecoder(AmbisonicUnitGenerator &input, SpeakerLayout *layout,
294 AmbisonicDecoderMethod method, AmbisonicDecoderFlavour flavour)
295 : AmbisonicUnitGenerator(input.order()), mInputPort(NULL), mSpeakerLayout(layout),
296 mDecodingMethod(method), mDecoderFlavour(flavour) {
297 initialize(input, method, flavour);
298 }
299
300 AmbisonicDecoder::AmbisonicDecoder(UnitGenerator &input, unsigned order, SpeakerLayout *layout,
301 AmbisonicDecoderMethod method, AmbisonicDecoderFlavour flavour)
302 : AmbisonicUnitGenerator(order), mInputPort(NULL),
303 mDecodingMethod(method), mDecoderFlavour(flavour) {
304 initialize(input, method, flavour);
305 }
306
307 AmbisonicDecoder::AmbisonicDecoder(UnitGenerator &input, unsigned hOrder, unsigned vOrder,
308 SpeakerLayout *layout, AmbisonicDecoderMethod method, AmbisonicDecoderFlavour flavour)
309 : AmbisonicUnitGenerator(hOrder, vOrder), mInputPort(NULL),
310 mDecodingMethod(method), mDecoderFlavour(flavour) {
311 initialize(input, method, flavour);
312 }
313
314 AmbisonicDecoder::~AmbisonicDecoder() {
315 mInputPort->mUGen->removeOutput(this);
316 delete mInputPort;
317 delete [] mIOChannelMap;
318
319 for (unsigned s = 0; s < mSpeakerLayout->numSpeakers(); s++)
320 delete [] mDecodingMatrix[s];
321 free(mDecodingMatrix);
322 }
323
324
325 /*******************INITIALIZE METHOD*******************/
326
327 void AmbisonicDecoder::initialize(UnitGenerator &input, AmbisonicDecoderMethod method, AmbisonicDecoderFlavour flavour) {
328
329 // if (mSpeakerLayout == NULL) // if the user didn't set a layout then ask for the default layout.
330 // mSpeakerLayout = SpeakerLayout::defaultSpeakerLayout();
331
332 mInputPort = new Port(&input);
333 input.addOutput(this); // be sure to add me as an output of the other guy
334
335 unsigned numSpeakers = mSpeakerLayout->numSpeakers(); // get a local copy of the number of speakers.
336 AmbisonicOrder decodingOrder = mOrder;
337
338 // are there enough speakers to decode this order?
339 // if (orderToChannels(decodingOrder) > numSpeakers) {
340 // // less speakers than necessary, so reduce decoding order by successively decrementing the vertical,
341 // // then horizontal order, until L >= N criterium applies
342 // while (orderToChannels(decodingOrder) > numSpeakers) {
343 // if (decodingOrder.verticalOrder >= decodingOrder.horizontalOrder) {
344 // --decodingOrder.verticalOrder;
345 // if (orderToChannels(decodingOrder) <= numSpeakers)
346 // break;
347 // }
348 // --decodingOrder.horizontalOrder;
349 // }
350 //
351 // logMsg(kLogWarning,
352 // "Reducing decoding order because of available number of speakers (%d): horizontal order = %i, vertical order = %i\n",
353 // numSpeakers, decodingOrder.horizontalOrder, decodingOrder.verticalOrder);
354 // }
355
356 // are there enough encoded channels in the input to decode the specified order?
357 if (orderToChannels(decodingOrder) > input.numChannels()) {
358 // reduce order if not enough channels in the input ugen
359 while (orderToChannels(decodingOrder) > input.numChannels()) {
360 if (decodingOrder.verticalOrder >= decodingOrder.horizontalOrder) {
361 --decodingOrder.verticalOrder;
362 if (orderToChannels(decodingOrder) <= input.numChannels())
363 break;
364 }
365 --decodingOrder.horizontalOrder;
366 }
367 logMsg(kLogWarning,
368 "Reducing decoding order because encoded input does not contain enough channels for specified order.\n Input channels:(%d): horizontal order = %i, vertical order = %i\n",
369 input.numChannels(), decodingOrder.horizontalOrder, decodingOrder.verticalOrder);
370 }
371
372
373 // recalculate the greater order & the uniformity
374 decodingOrder.isUniform = (decodingOrder.horizontalOrder == decodingOrder.verticalOrder);
375 mNumChannels = orderToChannels(decodingOrder); // how many input ambisonic channels needed at this order
376 // how many channels it would need *if* it was a uniform order (hOrder == vOrder)
377 unsigned numChannelsGreaterOrder = orderToChannels(greaterOrder(decodingOrder));
378 unsigned channelIndex[numChannelsGreaterOrder];
379 unsigned invChannelIndex[numChannelsGreaterOrder];
380
381 // The invChannelIndexer doesn't go thru all channels, so make sure to set others to zero.
382 memset(invChannelIndex , 0, numChannelsGreaterOrder * sizeof(unsigned));
383
384 channelIndexer(channelIndex); // get a channel indexer based on the input order
385 invChannelIndexer(invChannelIndex); // get an indexer from output framestream channel to 'ideal' ambisonic channel (after fixing the mOrder)
386
387 mIOChannelMap = new int[mNumChannels];
388 for (unsigned n = 0; n < mNumChannels; n++) {
389 // cout << n << ": " << invChannelIndex[n] << endl;
390 mIOChannelMap[n] = channelIndex[invChannelIndex[n]];
391 }
392
393 mDecodingMatrix = (sample**) malloc(numSpeakers * sizeof(sample*));
394 for (unsigned s = 0; s < numSpeakers; s++)
395 mDecodingMatrix[s] = new sample[mNumChannels];
396
397 if (mDecodingMethod == kPSEUDOINVERSE)
398 asPseudoInverse(); // build the Decoding matrix D using the pseudo inverse method
399 else
400 asProjection(); // build the Decoding matrix using the projection method
401
402
403 // process the decoding matrix according to the desired decoder flavour
404 if (mDecoderFlavour == kINPHASE) {
405 if ( ! decodingOrder.isUniform)
406 logMsg(kLogWarning, "AmbisonicDecoder: the in-phase decoding method is only well defined for non-hybrid systems - you have asked for an in phase decoder for a hybrid ambisonic decoding.\n");
407 makeInPhase(greaterOrder(decodingOrder));
408 }
409 else if (mDecoderFlavour == kMAXRE) {
410 if ( ! decodingOrder.isUniform)
411 logMsg(kLogWarning, "AmbisonicDecoder: the max-rE decoding method is only well defined for non-hybrid systems - you have asked for a max rE decoder for a hybrid ambisonic decoding.\n");
412 makeMaxRE(greaterOrder(decodingOrder));
413 }
414
415 #ifdef CSL_DEBUG
416
417 printf("Decoding from %i ambisonic channels to %i speakers", mNumChannels, numSpeakers);
418 printf(mDecodingMethod == kPSEUDOINVERSE ? ", using pseudo inverse method" : ", using projection method");
419 switch(mDecoderFlavour) {
420 case 0: printf(" of basic flavour\n"); break;
421 case 1: printf(" of in-phase flavour\n"); break;
422 case 2: printf(" of max rE flavour\n"); break;
423 }
424
425 printf("Decoding Matrix:\n");
426 for (unsigned s=0; s< numSpeakers; s++ ) {
427 printf("speaker %d\t",s);
428 for (unsigned n=0; n < mNumChannels; n++) {
429 printf("%f\t\t", mDecodingMatrix[s][n]);
430 }
431 printf("\n");
432 }
433 printf("\n");
434
435 printf("Created AmbisonicDecoder with horizontal order = %d, vertical order = %d\n", mOrder.horizontalOrder, mOrder.verticalOrder);
436
437 #endif
438
439 }
440
441
442 /*******************DECODING MATRIX METHODS*******************/
443 /*
444 B ... column vector of encoded Ambisonic channels, dimension N x 1
445 C ... re-encoding matrix, dimension N x L
446 C' ... transposed re-encoding matrix, dimension L x N
447 pinv(C) ... pseudoinverse of the re-encoding matrix C, dimension L x N
448 C * C' ... matrix product of the matrices C and C', dimension N x N
449 inv (C * C') ... inverse of C * C', dimension N x N
450 D ... decoding matrix, dimension L x N
451 p ... column vector of decoded loudspeaker driving signals, dimension L x 1
452
453 The matching condition (reference soundfield = synthesized soundfield): B = C * p
454 gives the decoding equation: p = D * B
455 Where D can be calculated by the projection method, or by the pseudoinverse method
456 */
457
458 /// build the Decoding matrix D using the projection method D = (1/L) * C'
459 void AmbisonicDecoder::asProjection() {
460
461 makeTransposedReEncodingMatrix(mDecodingMatrix); // get C' before scaling it down by 1/L
462
463 unsigned numSpeakers = mSpeakerLayout->numSpeakers();
464 sample invNumSpeakers = 1.f/((float)numSpeakers); // inverse of number of speakers (1/L)
465
466 // for each element in the matrix...
467 for (unsigned s = 0; s < numSpeakers; s++) {
468 for (unsigned i = 0; i < mNumChannels; i++) {
469 mDecodingMatrix[s][i] *= invNumSpeakers; // ...scale down matrix element by 1/L:
470 }
471 }
472 }
473
474 /// Build the Decoding matrix D using the pseudo inverse method \n D = pinv(C) = C' * inv(C * C'). \n
475 /// Pseudo inverse code based on the matrix library found at: http://home1.gte.net/edwin2/Matrix/
476 void AmbisonicDecoder::asPseudoInverse() {
477 makeTransposedReEncodingMatrix(mDecodingMatrix); // get C' before scaling it down by 1/L
478
479 unsigned nnumChannels = mNumChannels;
480 unsigned numSpeakers = mSpeakerLayout->numSpeakers();
481
482
483 // build temporary 2D array of dimension N x N to be used in calculating the psuedo-inverse: (C * C') and inv(C * C')
484 SampleBufferVector squareMatrix = (sample**) malloc(nnumChannels * sizeof(sample*));
485 for (unsigned n = 0; n < nnumChannels; n++) {
486 squareMatrix[n] = new sample[nnumChannels];
487 }
488
489 SampleBufferVector invSquareMatrix = (sample**) malloc(nnumChannels * sizeof(sample*));
490 for (unsigned n = 0; n < nnumChannels; n++) {
491 invSquareMatrix[n] = new sample[nnumChannels];
492 }
493
494 // create temporary matrices
495 sample** transposedReEncodingMatrix = (sample**) malloc(numSpeakers * sizeof(sample*));
496 for (unsigned s = 0; s < numSpeakers; s++)
497 transposedReEncodingMatrix[s] = new sample[nnumChannels];
498
499 makeTransposedReEncodingMatrix(transposedReEncodingMatrix);
500
501 // calculate (C * C'):
502 for (unsigned n = 0; n < nnumChannels; n++) { // for each row of (C * C')...
503
504 for (unsigned j = 0; j < nnumChannels; j++ ) { // for each column of (C * C')...
505 squareMatrix[n][j] = 0.f; // ... zero the new cell
506 // matrix multiplication: stepping along the row of C and down the column of C', summing as it goes
507 for (unsigned s = 0; s < numSpeakers; s++ ) {
508 // (C * C') - we read C' as C by swapping the rows and columns of C' in the first term
509 squareMatrix[n][j] += transposedReEncodingMatrix[s][n] * transposedReEncodingMatrix[s][j];
510 }
511 }
512 }
513
514 // ***start of matrix inversion: inv(C * C')***
515
516 SampleBuffer uu = new sample [nnumChannels * nnumChannels]; // more temporary variables
517 SampleBuffer vv = new sample [nnumChannels * nnumChannels];
518 SampleBuffer w = new sample [nnumChannels];
519 SampleBufferVector u = new sample*[nnumChannels];
520 SampleBufferVector v = new sample*[nnumChannels];
521
522 for (unsigned i = 0; i < nnumChannels; i++) {
523 u[i] = &(uu[nnumChannels*i]);
524 v[i] = &(vv[nnumChannels*i]);
525 for (unsigned j = 0; j < nnumChannels; j++)
526 u[i][j] = squareMatrix[i][j];
527 }
528
529 singularValueDecomposition(u, nnumChannels, nnumChannels, w, v); // singular value decomposition
530 sample wmax = 0.0; // maximum singular value.
531
532 for (unsigned j = 0; j < nnumChannels; j++) {
533 if (w[j] > wmax)
534 wmax = w[j];
535 }
536
537 sample wmin = wmax*1e-12; // epsilon
538 for (unsigned k = 0; k < nnumChannels; k++)
539 if (w[k] < wmin)
540 w[k] = 0.0;
541 else
542 w[k] = 1.0/w[k];
543
544 for (unsigned i = 0; i < nnumChannels; i++) {
545 for (unsigned j = 0; j < nnumChannels; j++) {
546 invSquareMatrix[i][j] = 0.0;
547 for (unsigned k = 0; k < nnumChannels; k++)
548 invSquareMatrix[i][j] += v[i][k]*w[k]*u[j][k];
549 }
550 }
551 // clean up
552 delete [] w;
553 delete [] u;
554 delete [] v;
555 delete [] uu;
556 delete [] vv;
557
558 // ***end of matrix inversion***
559
560 // get decoding matrix C by pseudoinverting it, i.e. multiply transposed re-encoding matrix by matrix
561 // we just inverted: D = C' * inv(C * C')
562
563 for (unsigned s = 0; s < numSpeakers; s++ ) { // for each row of mDecodingMatrix...
564 for (unsigned n = 0; n < nnumChannels; n++) { // for each column of mDecodingMatrix...
565 mDecodingMatrix[s][n] = 0.f; // ...zero the new cell
566 for (unsigned j = 0; j < nnumChannels; j++) { // for each column of C'...
567 // ... do operations to obtainD = C' * inv(C * C')
568 mDecodingMatrix[s][n] += transposedReEncodingMatrix[s][j] * invSquareMatrix[j][n];
569 }
570 }
571 }
572
573 // Cleanup
574 for (unsigned n = 0; n < nnumChannels; n++) {
575 delete [] squareMatrix[n];
576 delete [] invSquareMatrix[n];
577 }
578 for (unsigned s = 0; s < numSpeakers; s++)
579 delete [] transposedReEncodingMatrix[s];
580
581 free(squareMatrix);
582 free(invSquareMatrix);
583 free(transposedReEncodingMatrix);
584 }
585
586
587 /*******************TRANSPOSED RE-ENCODING MATRIX GENERATION METHOD*******************/
588
589 void AmbisonicDecoder::makeTransposedReEncodingMatrix(float **transposeMatrix)
590 {
591 unsigned numSpeakers = mSpeakerLayout->numSpeakers();
592 Speaker *speaker;
593 AmbisonicOrder order = mOrder;
594 float **transposedReEncodingMatrix = transposeMatrix;
595
596 // for each speaker
597 for (unsigned s = 0; s < numSpeakers; s++) {
598 speaker = mSpeakerLayout->speakerAtIndex(s); // get the details of this speaker
599
600 // create a row in the LxN matrix containing the encoding weights (FuMa) for this speaker
601 transposedReEncodingMatrix[s][0] = AMBI_INVSQRT2; // W channel is always the same
602 fumaEncodingWeights(transposedReEncodingMatrix[s], order, speaker->azimuth(), speaker->elevation());
603 }
604
605 return;
606 }
607
608
609 /*******************DECODING MATRIX FLAVOUR ADOPTION METHODS*******************/
610
611 /// Scales the decoding matrix according to the factors for max rE decoding
612 void AmbisonicDecoder::makeMaxRE(unsigned greaterOrder) {
613 unsigned numSpeakers = mSpeakerLayout->numSpeakers();
614
615 // hard-coded gain coefficients for in max-rE decoding (parameters specified up to 4th order):
616 sample maxREParameters[][5] = { {1.0, 1.0, 1.0, 1.0, 1.0},
617 // n = 0, M = 0, 1, 2, 3, 4
618 {0.f, 0.577, 0.775, 0.861, 0.906},
619 // n = 1, M = 0, 1, 2, 3, 4
620 {0.f, 0.f, 0.4, 0.612, 0.732},
621 // n = 2, M = 0, 1, 2, 3, 4
622 {0.f, 0.f, 0.f, 0.305, 0.501},
623 // n = 3, M = 0, 1, 2, 3, 4
624 {0.f, 0.f, 0.f, 0.f, 0.246},
625 // n = 4, M = 0, 1, 2, 3, 4
626 };
627
628 // calculate gain factors for each order
629 unsigned M = greaterOrder; // the order we are assuming as the context
630
631 for (unsigned i = 1; i < mNumChannels; i++) { // skip W channel as it's factor is always 1
632 for (unsigned l = 0; l < numSpeakers; l++) {
633 // channelsToUniformOrder(i)+1 is a convenient way to get n (the order of each ambisonic channel), but it is not defined for W
634 mDecodingMatrix[l][i] *= maxREParameters[channelsToUniformOrder(i)+1][M];
635 }
636 }
637 }
638
639 /// Scales the decoding matrix according to the factors for in-phase decoding
640 void AmbisonicDecoder::makeInPhase(unsigned tgreaterOrder)
641 {
642 unsigned numSpeakers = mSpeakerLayout->numSpeakers();
643
644 // hard-coded gain coefficients for in phase decoding (parameters specified up to 4th order)
645
646 sample inPhaseParameters[][5] = { // n = 0, 1, 2, 3, 4, M = 0, 1, 2, 3, 4
647 {1.0, 1.0, 1.0, 1.0, 1.0},
648 {0.f, 0.333, 0.5, 0.6, 0.667},
649 {0.f, 0.f, 0.1, 0.2, 0.286},
650 {0.f, 0.f, 0.f, 0.029, 0.071},
651 {0.f, 0.f, 0.f, 0.f, 0.008},
652 };
653
654 // calculate gain factors for each order
655 unsigned M = tgreaterOrder; // the order we are assuming as the context
656
657 for (unsigned i = 1; i < mNumChannels; i++) { // skip W channel as its factor is always 1
658 for (unsigned l = 0; l < numSpeakers; l++) {
659 // channelsToUniformOrder(i)+1 is a convenient way to get n (the order of each ambisonic channel), but it is not defined for W
660 mDecodingMatrix[l][i] *= inPhaseParameters[channelsToUniformOrder(i)+1][M];
661 }
662 }
663 }
664
665
666 /*******************NEXT_BUFFER METHOD*******************/
667
668 // To get my next buffer, get a buffer from the input, and then "process" it...
669 void AmbisonicDecoder::nextBuffer(Buffer &outputBuffer, unsigned outBufNum) throw(CException) {
670
671 unsigned numChannels = mNumChannels;
672 unsigned numSpeakers = mSpeakerLayout->numSpeakers();
673 unsigned numFrames = outputBuffer.mNumFrames; // block size
674 Buffer *inputBuffer;
675
676 if (outputBuffer.mNumChannels < numSpeakers) { // USED TO BE numChannels. I think it makes more sense to be numSpeakers.
677 logMsg(kLogError, "AmbisonicDecoder needs a buffer with %d channels, found only %d\n", numChannels, outputBuffer.mNumChannels);
678 throw RunTimeError("Insufficient number of channels in buffer passed.");
679 }
680
681 outputBuffer.zeroBuffers(); // fill the output buffers with silence
682
683 // inputBuffer = mInputPort->pullInput(numFrames);
684 // THE BLOCK BELOW WAS PULL_INPUT OF THE UGEN PORT...
685 inputBuffer = mInputPort->mBuffer; // else get the buffer
686 inputBuffer->mNumFrames = numFrames;
687 inputBuffer->mType = kSamples;
688 mInputPort->mUGen->nextBuffer(*inputBuffer); // and ask the UGen for nextBuffer()
689 inputBuffer->mIsPopulated = true;
690
691
692 // since all equations for all speakers & encoded channels are the same, they can be merged into a single loop:
693 for (unsigned n = 0; n < numChannels; n++) { // for each input ambisonic channel to be used
694 // get pointers to the encoded input channels (according to the channel indexer to accomodate hybrid orders)
695 SampleBuffer inPtr = inputBuffer->buffer(mIOChannelMap[n]);
696
697 // for each speaker, for each order, for each sample & encoded channel
698 for (unsigned l = 0; l < numSpeakers; l++) {
699 SampleBuffer outPtr = outputBuffer.buffer(l); // get a pointer to the output channel
700
701 // for each frame of the buffer
702 for (unsigned i = 0; i < numFrames; i++) {
703
704 *outPtr++ += mDecodingMatrix[l][n] * inPtr[i];
705 }
706 }
707 }
708
709 return;
710 }