CSL 5.2

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

- File List
- File Members

AmbisonicUtilities.cpp
Go to the documentation of this file.
1 //
2 // HOA_Utilities.cpp -- Higher Order Ambisonic utility classes
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 // Utility classes used by HOA_Encoder, MOA_Encoder, AmbisonicRotator, MOA_Rotator, HOA_Decoder
7
8 #include <stdlib.h>
9 #include <stdio.h>
10 #include <math.h>
11 #include "AmbisonicUtilities.h"
12
13 using namespace csl;
14
15 ////////////////////////////////////////////
16 /* AMBISONIC MIXER IMPLEMENTATION */
17 ////////////////////////////////////////////
18
19 /*******************CONSTRUCTORS & DESTRUCTOR*******************/
20
21 AmbisonicMixer::AmbisonicMixer(unsigned norder) : AmbisonicUnitGenerator(norder) {
22 initialize();
23 }
24
25 AmbisonicMixer::AmbisonicMixer(unsigned hOrder, unsigned vOrder) : AmbisonicUnitGenerator(hOrder, vOrder) {
26 initialize();
27 }
28
29 AmbisonicMixer::~AmbisonicMixer() { }
30
31
32 /*******************INITIALIZE METHOD*******************/
33
34 void AmbisonicMixer::initialize()
35 {
36 // allocate memory for source signal buffers according to number of channels required at this order
37 mInBuffer = new Buffer(mNumChannels, CGestalt::blockSize());
38 mInBuffer->allocateBuffers();
39
40 #ifdef CSL_DEBUG
41 logMsg("Created an AmbisonicMixer with horizontal order = %d, vertical order = %d\n", mOrder.horizontalOrder, mOrder.verticalOrder);
42 #endif
43
44 }
45
46 void AmbisonicMixer::addInput(UnitGenerator &input) {
47
48 if (mNumChannels <= input.numChannels()) { // Make sure the input has at least the number of channels needed for the order specified
49
50 mInputs.push_back(&input); // add the new input
51
52 // calculate inverse of the number of input sound sources (used for normalization)
53 mInvNumInputs = 1.f/(float)mInputs.size();
54
55 #ifdef CSL_DEBUG
56 logMsg("AmbisonicPanner: added %ith input\n", mInputs.size());
57 #endif
58 } else {
59 logMsg(kLogError, "AmbisonicMixer: cannot add input to a mixer of %ix%i order.\n", mOrder.horizontalOrder, mOrder.verticalOrder);
60 }
61 return;
62 }
63
64 void AmbisonicMixer::addInput(AmbisonicUnitGenerator &input) {
65
66 AmbisonicOrder inputOrder = input.order();
67 if (inputOrder.horizontalOrder >= mOrder.horizontalOrder && inputOrder.verticalOrder >= mOrder.verticalOrder) {
68 mInputs.push_back(&input); // add the new input
69
70 // calculate inverse of the number of input sound sources (used for normalization)
71 mInvNumInputs = 1.f/(float)mInputs.size();
72
73 #ifdef CSL_DEBUG
74 logMsg("AmbisonicPanner: added %ith input\n", mInputs.size());
75 #endif
76
77 } else {
78 logMsg(kLogError, "AmbisonicMixer: cannot add input of %ix%i order to a mixer of %ix%i order.\n", inputOrder.horizontalOrder, inputOrder.verticalOrder, mOrder.horizontalOrder, mOrder.verticalOrder);
79 }
80 return;
81
82 }
83
84 /*******************NEXT_BUFFER METHOD*******************/
85
86 // To get my next buffer, get a buffer from the input, and then "process" it...
87 void AmbisonicMixer::nextBuffer(Buffer &outputBuffer, unsigned outBufNum) throw(CException) {
88 SampleBuffer outPtr, inPtr;
89 unsigned i, c, k, numFrames;
90 unsigned numInputs = mInputs.size();
91
92 mInBuffer->mNumFrames = numFrames = outputBuffer.mNumFrames; // block size
93
94 outputBuffer.zeroBuffers(); // initialize output buffer with zeros
95
96 #ifdef CSL_DEBUG
97 logMsg("AmbisonicMixer::nextBuffer");
98 #endif
99
100 if (numInputs) { // if there are no inputs, then the outputbuffer can remain as zeroes.
101
102 if (outputBuffer.mNumChannels < mNumChannels)
103 logMsg(kLogError, "AmbisonicMixer requires %d output channels, found only %d ???\n", mNumChannels, outputBuffer.mNumChannels);
104
105 // loop through each input in turn, adding it's next_buffer output onto the outputbuffer channel by channel
106 for (i = 0; i < numInputs; i++) {
107 // pointer to the current input framestream
108 UnitGenerator * input = mInputs[i];
109
110 if (input->isActive()) { // if active input found, get a buffer and sum it into the output
111
112 // fill up mInBuffer with the ambisonic encoded audio from the current input
113 input->nextBuffer(*mInBuffer);
114
115 for (c = 0; c < mNumChannels; c++) { // loops through channels
116 outPtr = outputBuffer.buffer(c);
117 inPtr = mInBuffer->buffer(c);
118 for (k = 0; k < numFrames; k++) // loops through sample buffers
119 *outPtr++ += *inPtr++;
120 }
121 }
122 } // end of input loop
123
124 // after mixing, scale down our encoded output by the number of encoded sources to avoid clipping:
125 // for each actual audio output channel (not each ambisonic channel of the greater order)...
126 for (c = 0; c < mNumChannels; c++) {
127 outPtr = outputBuffer.buffer(c); // reset the outPtr to the beginning of this channel in the outputbuffer
128
129 for (i = 0; i < numFrames; i++) // ...scale this channels' output buffer down
130 *outPtr++ *= mInvNumInputs;
131 }
132 } // end of if (numInputs)
133
134 return;
135 }
136
137
138 ////////////////////////////////////////////
139 /* AMBISONIC ROTATOR IMPLEMENTATION */
140 ////////////////////////////////////////////
141
142
143 /*******************CONSTRUCTORS & DESTRUCTOR*******************/
144
145 AmbisonicRotator::AmbisonicRotator(AmbisonicUnitGenerator &input) : AmbisonicUnitGenerator(input.order()) {
146
147 initialize(input);
148 }
149
150 AmbisonicRotator::AmbisonicRotator(UnitGenerator &input, unsigned order) : AmbisonicUnitGenerator(order) {
151
152 initialize(input);
153 }
154
155 AmbisonicRotator::AmbisonicRotator(UnitGenerator &input, unsigned horder, unsigned vorder) : AmbisonicUnitGenerator(horder, vorder) {
156
157 initialize(input);
158 }
159
160 AmbisonicRotator::~AmbisonicRotator(void) {
161
162 mInputPort->mUGen->removeOutput(this);
163 delete mInputPort;
164
165 delete [] mSinAngle;
166 delete [] mCosAngle;
167
168 delete [] mInputChannelIndex;
169 delete [] mChannelIndex;
170
171 for (int i = 0; i < mNumChannels; i++) { // one sample pointer for each ambisonic output channel
172 delete mOutPtr[i];
173 delete mInPtr[i];
174 }
175
176 free(mOutPtr);
177 free(mInPtr);
178
179 }
180
181
182 /*******************INITIALIZE METHOD*******************/
183
184 void AmbisonicRotator::initialize(UnitGenerator &input) {
185 mInputPort = new Port(&input);
186 input.addOutput(this); // be sure to add me as an output of the other guy
187
188 // are there enough encoded channels in the input to rotate the specified order?
189 if (orderToChannels(mOrder) > input.numChannels()) {
190 // reduce order if not enough channels in the input ugen
191 while (orderToChannels(mOrder) > input.numChannels()) {
192 if (mOrder.verticalOrder >= mOrder.horizontalOrder) {
193 --mOrder.verticalOrder;
194 if (orderToChannels(mOrder) <= input.numChannels())
195 break;
196 }
197
198 --mOrder.horizontalOrder;
199 }
200
201 logMsg(kLogWarning, "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",
202 input.numChannels(), mOrder.horizontalOrder, mOrder.verticalOrder);
203 }
204
205 // recalculate the greater order & the uniformity (overriding values calculated in inherited HOA_AmbisonicFramestream constructor)
206 mGreaterOrder = greaterOrder(mOrder);
207 mOrder.isUniform = (mOrder.horizontalOrder == mOrder.verticalOrder);
208 mNumChannels = orderToChannels(mOrder); // derive necessary output channels from ambisonic order
209 mNumChannelsGreaterOrder = orderToChannels(mGreaterOrder); // channels required for uniform greater order
210
211 mInputChannelIndex = new unsigned[mNumChannelsGreaterOrder];
212 channelIndexer(mInputChannelIndex); // get a channel indexer based on the input order
213
214 mChannelIndex = new unsigned[mNumChannelsGreaterOrder];
215 channelIndexer(mChannelIndex); // get a channel indexer based on the output order (after fixing the mOrder)
216
217 mShouldRotate = mShouldTurn = mShouldTilt = false; // default initialise to no rotational activity
218
219 mOutPtr = (sample**) malloc(mNumChannels * sizeof(sample*)); // Allocating memory for audio output pointers.
220 mInPtr = (sample**) malloc(mNumChannels * sizeof(sample*)); // Allocating memory for audio output pointers.
221 for (unsigned i = 0; i < mNumChannels; i++) { // one sample pointer for each ambisonic output channel
222 mOutPtr[i] = new sample;
223 mInPtr[i] = new sample;
224 }
225
226 mSinAngle = new sample[mGreaterOrder]; // create arrays of samples to store sin & cosines of tilt, tumble & rotate angles
227 mCosAngle = new sample[mGreaterOrder]; // used in the rotator DSP loop
228
229 #ifdef CSL_DEBUG
230 logMsg("Created AmbisonicRotator with horizontal order = %d, vertical order = %d\n", mOrder.horizontalOrder, mOrder.verticalOrder);
231 #endif
232 }
233
234
235 /*******************SET_N-TH_INPUT METHOD*******************/
236
237 // function to set framestream sources for rotation, til tand tumble
238 void AmbisonicRotator::setNthInput(float amount, Axes axis)
239 {
240 // for example, use setNthInput(sine, kTILT);
241 switch (axis) {
242 case kTILT:
243 setTilt(amount);
244 break;
245 case kTUMBLE:
246 setTumble(amount);
247 break;
248 case kROTATE:
249 setRotate(amount);
250 break;
251 default:
252 logMsg(kLogError, "AmbisonicRotator; Attempted to add control signal to invalid axis channel ???\n");
253 }
254 }
255
256 void AmbisonicRotator::setTilt(float amount) {
257 if (mOrder.isUniform) {
258 mTilt = amount;
259 if (amount)
260 mShouldTilt = true;
261 else
262 mShouldTilt = false; // if amount == 0 then don't even do the operations. Nothing to be changed!
263 } else {
264 logMsg(kLogWarning, "Attempt to set Tilt control for a hybrid order rotator failed.\n");
265 }
266 }
267 void AmbisonicRotator::setTumble(float amount) {
268 if (mOrder.isUniform) {
269 mTumble = amount;
270
271 if (amount)
272 mShouldTurn = true;
273 else
274 mShouldTurn = false; // if amount == 0 then don't even do the operations. Nothing to be changed!
275 } else {
276 logMsg(kLogWarning, "Attempt to set Tumble control for a hybrid order rotator failed.\n");
277 }
278 }
279 void AmbisonicRotator::setRotate(float amount) {
280
281 mRotate = amount;
282
283 if (amount)
284 mShouldRotate = true;
285 else
286 mShouldRotate = false; // if amount == 0 then don't even do the operations. Nothing to be changed!
287
288 }
289
290
291 /*******************NEXT_BUFFER METHOD*******************/
292
293 // Overriding Framestream::nextBuffer(); To get my next buffer, get a buffer from the input, and then "process" it...
294 void AmbisonicRotator::nextBuffer(Buffer &outputBuffer, unsigned outBufNum) throw(CException) {
295 unsigned numFrames = mNumFrames = outputBuffer.mNumFrames; // block size
296 unsigned numChannels = mNumChannels;
297 unsigned greaterOrder = mGreaterOrder;
298 unsigned i; // to be used as the index in several for loops.
299 Buffer *inputBuffer;
300
301
302 if (outputBuffer.mNumChannels < numChannels) { // make sure that enough output channels were provided
303 logMsg(kLogError, "AmbisonicRotator needs a buffer with %d channels, found only %d\n", numChannels, outputBuffer.mNumChannels);
304 throw RunTimeError("Insufficient number of channels in buffer passed.");
305 }
306
307 // inputBuffer = mInputPort->pullInput(numFrames);
308 // THE BLOCK BELOW WAS PULL_INPUT OF THE UGEN PORT...
309 inputBuffer = mInputPort->mBuffer; // else get the buffer
310 inputBuffer->mNumFrames = numFrames;
311 inputBuffer->mType = kSamples;
312 mInputPort->mUGen->nextBuffer(*inputBuffer); // and ask the UGen for nextBuffer()
313 inputBuffer->mIsPopulated = true;
314
315 // Get a buffer from the input signal
316 // _input->nextBuffer(outputBuffer);
317
318 // assign our in & out pointers to their respective channels of input & output buffers
319 for (i = 0; i < mNumChannelsGreaterOrder; i++) {
320 mOutPtr[i] = outputBuffer.buffer(mChannelIndex[i]);
321 }
322
323 for (i = 0; i < numChannels; i++) {
324 mInPtr[i] = inputBuffer->buffer(mInputChannelIndex[i]);
325 }
326
327 // First, copy the input encoded audio to the output, ready for rotations if required
328 for (i = 0; i < numFrames; i++) //
329 for (unsigned j = 0; j < numChannels; j++)
330 *mOutPtr[j]++ = *mInPtr[j]++;
331
332 // do tilt processing if required
333 if (mShouldTilt) {
334 // reset outPtr to the begining of the buffer
335 for (i = 0; i < numChannels; i++)
336 mOutPtr[i] = outputBuffer.buffer(mChannelIndex[i]);
337
338 // set the sin & cosines of i multiples of the tilt angles
339 for (i = 0; i < greaterOrder; i++) {
340 mSinAngle[i] = sinf((i+1) * mTilt);
341 mCosAngle[i] = cosf((i+1) * mTilt);
342 }
343
344 // do the actual calculations
345 tiltFirstOrder();
346
347 }
348
349 // do tumble processing if required
350 if (mShouldTurn) {
351 // reset outPtr to the begining of the buffer
352 for (i = 0; i < numChannels; i++)
353 mOutPtr[i] = outputBuffer.buffer(mChannelIndex[i]);
354
355 // set the sin & cosines of i multiples of the tilt angles
356 for (i = 0; i < greaterOrder; i++) {
357 mSinAngle[i] = sinf((i+1) * mTumble);
358 mCosAngle[i] = cosf((i+1) * mTumble);
359 }
360
361 // do the actual calculations
362 tumbleFirstOrder();
363
364 }
365
366 // do the rotation processing if required
367 if (mShouldRotate) {
368 // reset outPtr to the begining of the buffer
369 for (i = 0; i < numChannels; i++)
370 mOutPtr[i] = outputBuffer.buffer(mChannelIndex[i]);
371
372 // set the sin & cosines of i multiples of the tilt angles
373 for (i = 0; i < greaterOrder; i++) {
374 mSinAngle[i] = sinf((i+1) * mRotate);
375 mCosAngle[i] = cosf((i+1) * mRotate);
376 }
377
378 // do the actual calculations - checking for hyrbid orders:
379 if (mOrder.horizontalOrder)
380 rotateFirstOrderHorizontal();
381 // note that there is no 1st order vertical calculation necessary for rotation
382 if (mOrder.verticalOrder > 1)
383 rotateSecondOrderVertical();
384
385 }
386
387 // done
388 return;
389
390 }
391
392 /*******************TILT, TUMBLE, ROTATE METHODS*******************/
393
394 void AmbisonicRotator::tiltFirstOrder() {
395 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
396 sample temp1;
397
398 for (unsigned i = 0; i < mNumFrames; i++) {
399 // W channel (0) remains unchanged
400 // X channel (1) remains unchanged
401 temp1 = *outPtr[2]; // needed to avoid recursive error
402 *outPtr[2]++ = mCosAngle[0] * (*outPtr[2]) - mSinAngle[0] * (*outPtr[3]);
403 *outPtr[3]++ = mSinAngle[0] * temp1 + mCosAngle[0] * (*outPtr[3]);
404
405 }
406
407 }
408
409 void AmbisonicRotator::tiltSecondOrder()
410 {
411 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
412 sample temp1, temp2, temp3; // needed to avoid recursive error
413
414 for (unsigned i = 0; i < mNumFrames; i++) {
415 temp1 = *outPtr[4]; // needed to avoid recursive error
416 temp2 = *outPtr[5];
417 temp3 = *outPtr[7];
418
419 *outPtr[4]++ = 0.5 * (mCosAngle[1] - 1) * (*outPtr[8]) + 0.5 * mSinAngle[1] * (*outPtr[7]) + 0.25 * (mCosAngle[1] + 3) * (*outPtr[4]);
420 *outPtr[5]++ = -mSinAngle[0] * (*outPtr[6]) + mCosAngle[0] * (*outPtr[5]);
421 *outPtr[6]++ = mCosAngle[0] * (*outPtr[6]) + mSinAngle[0] * temp2;
422 *outPtr[7]++ = -mSinAngle[1] * (*outPtr[8]) + mCosAngle[1] * (*outPtr[7]) - (0.5 * mSinAngle[1]) * temp1;
423 *outPtr[8]++ = 0.25 * (1 + 3 * mCosAngle[1]) * (*outPtr[8]) + 0.75 * mSinAngle[1] * temp3 + 0.375 * (mCosAngle[1] - 1) * temp1;
424
425 }
426
427 if (mGreaterOrder > 2) tiltThirdOrder();
428
429 }
430
431 void AmbisonicRotator::tiltThirdOrder()
432 {
433 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
434 sample temp1, temp2, temp3, temp4, temp5;
435
436 for (unsigned i = 0; i < mNumFrames; i++) {
437 // temporaries needed to avoid recursive error
438 temp1 = *outPtr[9];
439 temp2 = *outPtr[10];
440 temp3 = *outPtr[11];
441 temp4 = *outPtr[12];
442 temp5 = *outPtr[14];
443
444 // The source for the 3rd order tilt matrix is Zmoelnig's thesis, taking into consideration his different channel naming convention!
445 *outPtr[9]++ = 0.125 * (3 * mCosAngle[1] + 5) * (*outPtr[9]) + 1.5 * mSinAngle[1] * (*outPtr[12]) + 0.515625 * (mCosAngle[1] - 1) * (*outPtr[13]); // P channel
446 *outPtr[10]++ = 0.0625 * (15 * mCosAngle[0] + mCosAngle[2]) * (*outPtr[10]) - 0.375 * (5 * mSinAngle[0] + mSinAngle[2]) * (*outPtr[11]) - 0.2578125 * (mCosAngle[0] - mCosAngle[2]) * (*outPtr[14]) + 0.25 * (3 * mSinAngle[0] - mSinAngle[2]) * (*outPtr[15]); // Q channel
447 *outPtr[11]++ = 0.0625 * (5 * mSinAngle[0] + mSinAngle[2]) * temp2 + 0.125 * (5 * mCosAngle[0] + 3 * mCosAngle[2]) * (*outPtr[11]) + 0.0859375 * (- mSinAngle[0] + 3 * mSinAngle[2]) * (*outPtr[14]) + 0.25 * (- mCosAngle[0] + mCosAngle[2]) * (*outPtr[15]); // N channel
448 *outPtr[12]++ = - 0.25 * mSinAngle[1] * temp1 + mCosAngle[1] * (*outPtr[12]) - 0.34375 * mSinAngle[1] * (*outPtr[13]); // O channel
449 *outPtr[13]++ = 0.454545454545 * (mCosAngle[1] - 1) * temp1 + + 1.818181818182 * mSinAngle[1] * temp4 + 0.125 * (3 + 5 * mCosAngle[1]) * (*outPtr[13]); // L channel
450 *outPtr[14]++ = 0.227272727273 * (- mCosAngle[0] + mCosAngle[2]) * temp2 + 0.454545454545 * (mSinAngle[0] - 3 * mSinAngle[2]) * temp3 + 0.0625 * (mCosAngle[0] + 15 * mCosAngle[2]) * (*outPtr[14]) - 0.181818181818 * (mSinAngle[0] + 5 * mSinAngle[2]) * (*outPtr[15]); // M channel
451 *outPtr[15]++ = 0.15625 * (- 3 * mSinAngle[0] + mSinAngle[2]) * temp2 + 0.9375 * (- mCosAngle[0] + mCosAngle[2]) * temp3 + 0.12890625 * (mSinAngle[0] + 5 * mSinAngle[2]) * temp5 + 0.125 * (3 * mCosAngle[0] + 5 * mCosAngle[2]) * (*outPtr[15]); // K channel
452
453 }
454 // if (mOrder.greaterOrder > 3) tiltFourthOrder();
455 }
456
457
458 void AmbisonicRotator::tumbleFirstOrder() {
459 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
460 sample temp1;
461
462 for (unsigned i = 0; i < mNumFrames; i++) {
463 temp1 = *outPtr[1]; // needed to avoid recursive error
464 // W channel (0) remains unchanged
465 *outPtr[1]++ = mCosAngle[0] * (*outPtr[1]) - mSinAngle[0] * (*outPtr[3]);
466 // Y channel (2) remains unchanged
467 *outPtr[3]++ = mSinAngle[0] * temp1 + mCosAngle[0] * (*outPtr[3]);
468 }
469 }
470
471 void AmbisonicRotator::tumbleSecondOrder() {
472
473 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
474 sample temp1, temp2, temp3;
475
476 for (unsigned i = 0; i < mNumFrames; i++) {
477 // temporaries needed to avoid recursive error
478 temp1 = *outPtr[4];
479 temp2 = *outPtr[5];
480 temp3 = *outPtr[6];
481
482 *outPtr[4]++ = 0.5 * (1 - mCosAngle[1]) * (*outPtr[8]) - 0.5 * mSinAngle[1] * (*outPtr[6]) + 0.25 * (3 + mCosAngle[1]) * (*outPtr[5]);
483 *outPtr[5]++ = -mSinAngle[0] * (*outPtr[7]) + mCosAngle[0] * (*outPtr[5]);
484 *outPtr[6]++ = -mSinAngle[1] * (*outPtr[8]) + mCosAngle[1] * (*outPtr[6]) + 0.5 * mSinAngle[1] * temp1;
485 *outPtr[7]++ = mCosAngle[0] * (*outPtr[7]) + mSinAngle[0] * temp2;
486 *outPtr[8]++ = 0.25 * (1 + 3 * mCosAngle[1]) * (*outPtr[8]) + 0.75 * mSinAngle[1] * temp3 + 0.375 * (1 - mCosAngle[1]) * temp1;
487 }
488 if (mGreaterOrder > 2) tumbleThirdOrder();
489 }
490
491 void AmbisonicRotator::tumbleThirdOrder()
492 {
493 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
494 sample temp1, temp2, temp3, temp4, temp5;
495
496 for (unsigned i = 0; i < mNumFrames; i++) {
497 // temporaries needed to avoid recursive error
498 temp1 = *outPtr[9];
499 temp2 = *outPtr[10];
500 temp3 = *outPtr[11];
501 temp4 = *outPtr[12];
502 temp5 = *outPtr[13];
503
504 // The source for the 3rd order tumble matrix is Zmoelnig's thesis, taking into consideration his different channel naming convention!
505 *outPtr[9]++ = 0.0625 * (15 * mCosAngle[0] + mCosAngle[2]) * (*outPtr[9]) - 0.375 * (5 * mSinAngle[0] + mSinAngle[2]) * (*outPtr[11]) + 0.2578125 * (mCosAngle[0] - mCosAngle[2]) * (*outPtr[13]) + 0.25 * (- 3 * mSinAngle[0] + mSinAngle[2]) * (*outPtr[15]); // P channel
506 *outPtr[10]++ = 0.125 * (5 + 3 * mCosAngle[1]) * (*outPtr[10]) - 1.5 * mSinAngle[1] * (*outPtr[12]) + 0.515625 * (1 - mCosAngle[1]) * (*outPtr[14]); // Q channel
507 *outPtr[11]++ = 0.0625 * (5 * mSinAngle[0] + mSinAngle[2]) * temp1 + 0.125 * (5 * mCosAngle[0] + 3 * mCosAngle[2]) * (*outPtr[11])+ 0.0859375 * (mSinAngle[0] - 3 * mSinAngle[2]) * (*outPtr[13]) + 0.25 * (mCosAngle[0] - mCosAngle[2]) * (*outPtr[15]); // N channel
508 *outPtr[12]++ = 0.25 * mSinAngle[1] * temp2 + mCosAngle[1] * (*outPtr[12]) - 0.34375 * mSinAngle[1] * (*outPtr[14]); // O channel
509 *outPtr[13]++ = 0.227272727273 * (mCosAngle[0] - mCosAngle[2]) * temp1 + 0.454545454545 * (- mSinAngle[0] + 3 * mSinAngle[2]) * temp3 + 0.0625 * (mCosAngle[0] + 15 * mCosAngle[2]) * (*outPtr[13]) - 0.181818181818 * (mSinAngle[0] + 5 * mSinAngle[2]) * (*outPtr[15]); // L channel
510 *outPtr[14]++ = 0.454545454545 * (1 - mCosAngle[1]) * temp2 + 1.818181818182 * mSinAngle[1] * temp4 + 0.125 * (3 + 5 * mCosAngle[1]) * (*outPtr[14]); // M channel
511 *outPtr[15]++ = 0.15625 * (3 * mSinAngle[0] - mSinAngle[2]) * temp1 + 0.9375 * (mCosAngle[0] - mCosAngle[2]) * temp3 + 0.12890625 * (mSinAngle[0] + 5 * mSinAngle[2]) * temp5 + 0.125 * (3 * mCosAngle[0] + 5 * mCosAngle[2]) * (*outPtr[15]); // K channel
512
513 }
514 //if (mOrder.greaterOrder > 3) tumbleFourthOrder();
515 }
516
517
518
519 void AmbisonicRotator::rotateFirstOrderHorizontal()
520 {
521 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
522 sample temp1;
523
524 for (unsigned i = 0; i < mNumFrames; i++) {
525 temp1 = *outPtr[1]; // needed to avoid recursive error
526 // W channel (0) remains unchanged
527 *outPtr[1]++ = mCosAngle[0] * (*outPtr[1]) - mSinAngle[0] * (*outPtr[2]);
528 *outPtr[2]++ = mSinAngle[0] * temp1 + mCosAngle[0] * (*outPtr[2]);
529 // Z channel (3) remains unchanged
530 }
531
532 if (mOrder.horizontalOrder > 1) rotateSecondOrderHorizontal();
533 }
534
535 // note that 1st order vertical for rotation is not required
536
537 void AmbisonicRotator::rotateSecondOrderHorizontal()
538 {
539 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
540 sample temp1;
541
542 for (unsigned i = 0; i < mNumFrames; i++) {
543 temp1 = *outPtr[4]; // needed to avoid recursive error
544 *outPtr[4]++ = mCosAngle[1] * (*outPtr[4]) - mSinAngle[1] * (*outPtr[5]); // U Channel
545 *outPtr[5]++ = mSinAngle[1] * temp1 + mCosAngle[1] * (*outPtr[5]); // V Channel
546 }
547 if (mOrder.horizontalOrder > 2) rotateThirdOrderHorizontal();
548 }
549
550 void AmbisonicRotator::rotateSecondOrderVertical()
551 {
552 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
553 sample temp1;
554
555 for (unsigned i = 0; i < mNumFrames; i++) {
556 temp1 = *outPtr[6]; // needed to avoid recursive error
557 *outPtr[6]++ = mCosAngle[0] * (*outPtr[6]) - mSinAngle[0] * (*outPtr[7]); // S channel
558 *outPtr[7]++ = mSinAngle[0] * temp1 + mCosAngle[0] * (*outPtr[7]); // T channel
559 // R channel (8) remains unchanged
560 }
561 if (mOrder.verticalOrder > 2) rotateThirdOrderVertical();
562 }
563
564 void AmbisonicRotator::rotateThirdOrderHorizontal()
565 {
566 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
567 sample temp1;
568
569 for (unsigned i = 0; i < mNumFrames; i++) {
570 temp1 = *outPtr[9]; // needed to avoid recursive error
571
572 // The source for the 3rd order tilt matrix is Zmoelnig's thesis, taking into consideration his different channel naming convention!
573 // This matrix seems to be different from the one that Malham specifies in his thesis!
574 *outPtr[9]++ = mCosAngle[2] * (*outPtr[9]) - mSinAngle[2] * (*outPtr[10]); // P channel
575 *outPtr[10]++ = mSinAngle[2] * temp1 + mCosAngle[2] * (*outPtr[10]); // Q channel
576 }
577 // if (mOrder.horizontalOrder > 3) rotateFourthOrderHorizontal();
578 }
579
580 void AmbisonicRotator::rotateThirdOrderVertical()
581 {
582 SampleBufferVector outPtr = mOutPtr; // create local copy of output pointer
583 sample temp2, temp3;
584
585 for (unsigned i = 0; i < mNumFrames; i++) {
586 temp2 = *outPtr[11];// needed to avoid recursive error
587 temp3 = *outPtr[13];
588
589 // The source for the 3rd order tilt matrix is Zmoelnig's thesis, taking into consideration his different channel naming convention!
590 // This matrix seems to be different from the one that Malham specifies in his thesis!
591 *outPtr[11]++ = mCosAngle[1] * (*outPtr[11]) - mSinAngle[1] * (*outPtr[12]); // N channel
592 *outPtr[12]++ = mSinAngle[1] * temp2 + mCosAngle[1] * (*mOutPtr[12]); // O channel
593 *outPtr[13]++ = mCosAngle[0] * (*outPtr[13]) - mSinAngle[0] * (*outPtr[14]); // L channel
594 *outPtr[14]++ = mSinAngle[0] * temp3 + mCosAngle[0] * (*outPtr[14]); // M channel
595 // K channel (15) remains unchanged
596 }
597 // if (mOrder.verticalOrder > 3) rotateFourthOrderVertical();
598 }
599
600
601
602 /*******************MATRIX INVERSION UTILITY FUNCTION*******************/
603
604 // code based on the matrix library found at: http://home1.gte.net/edwin2/Matrix/
605
606 /*
607 PYTHAG computes sqrt(a^{2} + b^{2}) without destructive overflow or underflow.
608 */
609 static double at, bt, ct;
610 #define PYTHAG(a, b) ((at = fabs(a)) > (bt = fabs(b)) ? \
611 (ct = bt/at, at*sqrt(1.0+ct*ct)): (bt ? (ct = at/bt, bt*sqrt(1.0+ct*ct)): 0.0))
612
613 static double maxarg1, maxarg2;
614 #define MAX(a, b) (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? \
615 (maxarg1) : (maxarg2))
616
617 #define SIGN(a, b) ((b) < 0.0 ? -fabs(a): fabs(a))
618
619 /** \relates AmbisonicUnitGenerator
620 Given a matrix a[m][n], this routine computes its singular value
621 decomposition, A = U*W*V^{T}. The matrix U replaces a on output.
622 The diagonal matrix of singular values W is output as a vector w[n].
623 The matrix V (not the transpose V^{T}) is output as v[n][n].
624 m must be greater or equal to n; if it is smaller, then a should be
625 filled up to square with zero rows.
626 */
627 void csl::singularValueDecomposition(sample** a, int m, int n, sample* w, sample** v)
628 {
629 int flag, i, its, j, jj, k, l, nm;
630 double c, f, h, s, x, y, z;
631 double anorm = 0.0, g = 0.0, scale = 0.0;
632
633 // if (m < n)
634 // error("SingularValueDecomposition: Matrix A must be augmented with extra rows of zeros.");
635 double* rv1 = new double [n];
636
637 /* Householder reduction to bidiagonal form. */
638 for (i = 0; i < n; i++) {
639 l = i + 1;
640 rv1[i] = scale*g;
641 g = s = scale = 0.0;
642 if (i < m) {
643 for (k = i; k < m; k++)
644 scale += fabs(a[k][i]);
645 if (scale) {
646 for (k = i; k < m; k++) {
647 a[k][i] /= scale;
648 s += a[k][i]*a[k][i];
649 };
650 f = a[i][i];
651 g = -SIGN(sqrt(s), f);
652 h = f*g - s;
653 a[i][i] = f - g;
654 if (i != n - 1) {
655 for (j = l; j < n; j++) {
656 for (s = 0.0, k = i; k < m; k++)
657 s += a[k][i]*a[k][j];
658 f = s/h;
659 for (k = i; k < m; k++)
660 a[k][j] += f*a[k][i];
661 };
662 };
663 for (k = i; k < m; k++)
664 a[k][i] *= scale;
665 };
666 };
667 w[i] = scale*g;
668 g = s= scale = 0.0;
669 if (i < m && i != n - 1) {
670 for (k = l; k < n; k++)
671 scale += fabs(a[i][k]);
672 if (scale) {
673 for (k = l; k < n; k++) {
674 a[i][k] /= scale;
675 s += a[i][k]*a[i][k];
676 };
677 f = a[i][l];
678 g = -SIGN(sqrt(s), f);
679 h = f*g - s;
680 a[i][l] = f - g;
681 for (k = l; k < n; k++)
682 rv1[k] = a[i][k]/h;
683 if (i != m - 1) {
684 for (j = l; j < m; j++) {
685 for (s = 0.0, k = l; k < n; k++)
686 s += a[j][k]*a[i][k];
687 for (k = l; k < n; k++)
688 a[j][k] += s*rv1[k];
689 };
690 };
691 for (k = l; k < n; k++)
692 a[i][k] *= scale;
693 };
694 };
695 anorm = MAX(anorm, (fabs(w[i]) + fabs(rv1[i])));
696 };
697 /* Accumulation of right-hand transformations. */
698 for (i = n - 1; 0 <= i; i--)
699 {
700 if (i < n - 1)
701 {
702 if (g)
703 {
704 for (j = l; j < n; j++)
705 v[j][i] = (a[i][j]/a[i][l])/g;
706 /* Double division to avoid possible underflow: */
707 for (j = l; j < n; j++)
708 {
709 for (s = 0.0, k = l; k < n; k++)
710 s += a[i][k]*v[k][j];
711 for (k = l; k < n; k++)
712 v[k][j] += s*v[k][i];
713 };
714 };
715 for (j = l; j < n; j++)
716 v[i][j] = v[j][i] = 0.0;
717 };
718 v[i][i] = 1.0;
719 g = rv1[i];
720 l = i;
721 };
722 /* Accumulation of left-hand transformations. */
723 for (i = n - 1; 0 <= i; i--)
724 {
725 l = i + 1;
726 g = w[i];
727 if (i < n - 1)
728 for (j = l; j < n; j++)
729 a[i][j] = 0.0;
730 if (g)
731 {
732 g = 1.0/g;
733 if (i != n - 1)
734 {
735 for (j = l; j < n; j++)
736 {
737 for (s = 0.0, k = l; k < m; k++)
738 s += a[k][i]*a[k][j];
739 f = (s/a[i][i])*g;
740 for (k = i; k < m; k++)
741 a[k][j] += f*a[k][i];
742 };
743 };
744 for (j = i; j < m; j++)
745 a[j][i] *= g;
746 }
747 else
748 for (j = i; j < m; j++)
749 a[j][i] = 0.0;
750 ++a[i][i];
751 };
752 /* Diagonalization of the bidiagonal form. */
753 for (k = n - 1; 0 <= k; k--) /* Loop over singular values. */
754 {
755 for (its = 0; its < 30; its++) /* Loop over allowed iterations.*/
756 {
757 flag = 1;
758 for (l = k; 0 <= l; l--) /* Test for splitting: */
759 {
760 nm = l - 1; /* Note that rv1[0] is always zero.*/
761 if (fabs(rv1[l]) + anorm == anorm)
762 {
763 flag = 0;
764 break;
765 };
766 if (fabs(w[nm]) + anorm == anorm)
767 break;
768 };
769 if (flag)
770 {
771 c = 0.0; /* Cancellation of rv1[l], if l>0:*/
772 s = 1.0;
773 for (i = l; i <= k; i++) {
774 f = s*rv1[i];
775 if (fabs(f) + anorm != anorm)
776 {
777 g = w[i];
778 h = PYTHAG(f, g);
779 w[i] = h;
780 h = 1.0/h;
781 c = g*h;
782 s = (-f*h);
783 for (j = 0; j < m; j++)
784 {
785 y = a[j][nm];
786 z = a[j][i];
787 a[j][nm] = y*c + z*s;
788 a[j][i] = z*c - y*s;
789 };
790 };
791 };
792 };
793 z = w[k];
794 if (l == k) /* Convergence. */
795 {
796 if (z < 0.0) /* Singular value is made non-negative. */
797 {
798 w[k] = -z;
799 for (j = 0; j < n; j++)
800 v[j][k] = (-v[j][k]);
801 };
802 break;
803 };
804 // if (its == 29)
805 // error("No convergence in 30 SingularValueDecomposition iterations.");
806 x = w[l]; /* Shift from bottom 2-by-2 minor. */
807 nm = k - 1;
808 y = w[nm];
809 g = rv1[nm];
810 h = rv1[k];
811 f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y);
812 g = PYTHAG(f, 1.0);
813 f = ((x - z)*(x + z) + h*((y/(f + SIGN(g, f))) - h))/x;
814 /* Next QR transformation: */
815 c = s = 1.0;
816 for (j = l; j <= nm; j++)
817 {
818 i = j + 1;
819 g = rv1[i];
820 y = w[i];
821 h = s*g;
822 g = c*g;
823 z = PYTHAG(f, h);
824 rv1[j] = z;
825 c = f/z;
826 s = h/z;
827 f = x*c + g*s;
828 g = g*c - x*s;
829 h = y*s;
830 y = y*c;
831 for (jj = 0; jj < n; jj++)
832 {
833 x = v[jj][j];
834 z = v[jj][i];
835 v[jj][j] = x*c + z*s;
836 v[jj][i] = z*c - x*s;
837 };
838 z = PYTHAG(f, h);
839 w[j] = z; /* Rotation can be arbitrary if z = 0. */
840 if (z)
841 {
842 z = 1.0/z;
843 c = f*z;
844 s = h*z;
845 };
846 f = (c*g) + (s*y);
847 x = (c*y) - (s*g);
848 for (jj = 0; jj < m; jj++)
849 {
850 y = a[jj][j];
851 z = a[jj][i];
852 a[jj][j] = y*c + z*s;
853 a[jj][i] = z*c - y*s;
854 };
855 };
856 rv1[l] = 0.0;
857 rv1[k] = f;
858 w[k] = x;
859 };
860 };
861 delete [] rv1;
862 }
863
864 /// \relates AmbisonicUnitGenerator
865 void csl::fumaEncodingWeights(SampleBuffer weights, const AmbisonicOrder &order, float azimuth, float elevation) {
866
867 float x,y,z,x2,y2,z2;
868 // assume default location of 0,0, i.e directly in front on the plane
869 x= x2 = 1.f;
870 y = z = y2 = z2 = 0.f;
871
872 unsigned h = order.horizontalOrder;
873 unsigned v = order.verticalOrder;
874 unsigned channel = 1;
875
876 // skip this line, do it in the initialization of weights instead, because it never changes
877 //weights[0] = AMBI_INVSQRT2; // W channel, shouldn't it be already defined?
878
879 if (h > 0) {
880 float cosel = cosf(elevation);
881 x = cosf(azimuth) * cosel;
882 y = sinf(azimuth) * cosel;
883
884 weights[channel++] = x; // X = cos(A)cos(E)
885 weights[channel++] = y; // Y = sin(A)cos(E)
886
887 }
888
889 if (v > 0) {
890 z = sinf(elevation);
891 weights[channel++] = z; // Z = sin(E)
892 }
893
894 if (h > 1) {
895 x2 = x*x;
896 y2 = y*y;
897
898 weights[channel++] = x2 - y2; // U = cos(2A)cos2(E) = xx-yy
899 weights[channel++] = 2.f * x * y; // V = sin(2A)cos2(E) = 2xy
900
901 }
902
903 if (v > 1) {
904 z2 = z*z;
905
906 weights[channel++] = 2.f * z * x; // S = cos(A)sin(2E) = 2zx
907 weights[channel++] = 2.f * z * y; // T = sin(A)sin(2E) = 2yz
908 weights[channel++] = (1.5f * z2) - 0.5f; // R = 1.5sin2(E)-0.5 = 1.5zz-0.5
909
910 }
911
912 if (h > 2) {
913 weights[channel++] = x * (x2 - 3.f*y2); // P = cos(3A)cos3(E) = X(X2-3Y2)
914 weights[channel++] = y * (y2 - 3.f*x2); // Q = sin(3A)cos3(E) = Y(3X2-Y2)
915 }
916
917 if (v > 2) {
918 float pre = 8.f/11.f;
919
920 weights[channel++] = z * (x2-y2) * 0.5f; // N = cos(2A)sin(E)cos2(E) = Z(X2-Y2)/2
921 weights[channel++] = x * y * z; // O = sin(2A)sin(E)cos2(E) = XYZ
922 weights[channel++] = pre * x * (5.f*z2 - 1.f); // L = 8cos(A)cos(E)(5sin2(E) - 1)/11 = 8X(5Z2-1)/11
923 weights[channel++] = pre * y * (5.f*z2 - 1.f); // M = 8sin(A)cos(E)(5sin2(E) - 1)/11 = 8Y(5Z2-1)/11
924 weights[channel++] = z * 0.5 * (5.f*z2 - 3.f); // K = sin(E)(5sin2(E) - 3)/2 = Z(5Z2-3)/2
925 }
926
927 }
928
929 /// \relates AmbisonicUnitGenerator
930 void csl::fumaIndexedEncodingWeights(SampleBuffer weights, const AmbisonicOrder &order, sample &azimuth, sample &elevation)
931 {
932 float x,y,z,x2,y2,z2;
933 // assume default location of 0,0, i.e directly in front on the plane
934 x= x2 = 1.f;
935 y= z = y2 = z2 = 0.f;
936
937
938 unsigned h = order.horizontalOrder;
939 unsigned v = order.verticalOrder;
940
941 // skip this line, do it in the initialization of weights instead, because it never changes
942 //weights[0] = AMBI_INVSQRT2; // W channel, shouldn't it be already defined?
943
944 if (h > 0) {
945 float cosel = cosf(elevation);
946 x = cosf(azimuth) * cosel;
947 y = sinf(azimuth) * cosel;
948
949 weights[1] = x; // X = cos(A)cos(E)
950 weights[2] = y; // Y = sin(A)cos(E)
951
952 }
953
954 if (v > 0) {
955 z = sinf(elevation);
956 weights[3] = z; // Z = sin(E)
957 }
958
959 if (h > 1) {
960 x2 = x*x;
961 y2 = y*y;
962
963 weights[4] = x2 - y2; // U = cos(2A)cos2(E) = xx-yy
964 weights[5] = 2.f * x * y; // V = sin(2A)cos2(E) = 2xy
965
966 if (h > 2) {
967 weights[9] = x * (x2 - 3.f*y2); // P = cos(3A)cos3(E) = X(X2-3Y2)
968 weights[10]= y * (y2 - 3.f*x2); // Q = sin(3A)cos3(E) = Y(3X2-Y2)
969 }
970 }
971
972 if (v > 1) {
973 z2 = z*z;
974
975 weights[6] = 2.f * z * x; // S = cos(A)sin(2E) = 2zx
976 weights[7] = 2.f * z * y; // T = sin(A)sin(2E) = 2yz
977 weights[8] = (1.5f * z2) - 0.5f; // R = 1.5sin2(E)-0.5 = 1.5zz-0.5
978
979 if (v > 2) {
980 float pre = 8.f/11.f;
981
982 weights[11] = z * (x2-y2) * 0.5f; // N = cos(2A)sin(E)cos2(E) = Z(X2-Y2)/2
983 weights[12] = x * y * z; // O = sin(2A)sin(E)cos2(E) = XYZ
984 weights[13] = pre * x * (5.f*z2 - 1.f); // L = 8cos(A)cos(E)(5sin2(E) - 1)/11 = 8X(5Z2-1)/11
985 weights[14] = pre * y * (5.f*z2 - 1.f); // M = 8sin(A)cos(E)(5sin2(E) - 1)/11 = 8Y(5Z2-1)/11
986 weights[15] = z * 0.5 * (5.f*z2 - 3.f); // K = sin(E)(5sin2(E) - 3)/2 = Z(5Z2-3)/2
987 }
988 }
989 }