Recurrence formula

[sekisushai]

Hello, I'm having some trouble to compile a matrix where the terms are computed with a recurrence formula. The code is attached at the end of this report

This code compute a rotation matrix of size (M+1)^2 x (M+1)^2 to rotate a HOA soundscape described at order M.

Each term of the matrix is computed with a recurrence formula starting with rotation matrix of order M=1.

The recurrence formula are a bit complex and taken from ref [1] [1] Ivanic, J., & Ruedenberg, K. (1996). Rotation matrices for real spherical harmonics. Direct determination by recursion. The Journal of Physical Chemistry, 100(15), 6342–6347.

In the code, a term of the rotation matrix at order M is given by:

rot(M,n1,n2) with -M <= n1,n2 <= M

So far, everything compiles up to order M=3. However as soon as the order is set as M=4 or higher, faust never finish the compilation.

In order to have more clues on this, I can ask as output only on term of the matrix, for example

process = rot(M,M,M) // bottom right term of the matrix at order M

In the same vein, I can ask a whole line of the matrix:

process = row(M,i) // where 0 <= i <= (M+1)^2-1

For M=4, all the terms compile separately, as well as all the line of the matrix. However, as soon as I ask the whole matrix, it doesn't compile anymore..

Any ideas ?

Thanks, Pierre

declare name "HOA Rotator";
declare version "1.0";
declare author      "Pierre Lecomte";
declare license     "GPL";
declare copyright   "(c) Pierre Lecomte 2016";

import("stdfaust.lib");

// Description: This tool rotates the HOA scene around the x-axis (roll angle), y-axis (pitch angle), and z-axis (yaw angle). Driven with OSC from head-tracking, this tool can compensate the head rotations. Implementation according to [1] with corrections.

// References: 
//[1] J. Ivanic and K. Ruedenberg, “Rotation matrices for real spherical harmonics. Direct determination by recursion,” J. Phys. Chem., vol. 100, no. 15, pp. 6342–6347, 1996.
// In TABLE 2 of this paper, from the corrections, the formulas for the term V(l,m,m') should be taken from the original paper...

// Inputs: (M+1)^2
// Outputs: (M+1)^2

// Bus with gain
buswg(c) = R(c) with {
  R((c,cl)) = R(c),R(cl);
  R(1)      = _;
  //R(0)      = !;
  R(0)      = !:0;
  R(float(0)) = R(0);
  R(float(1)) = R(1);
  R(c)      = *(c);
};

// Sliders
yaw     =   hslider("Yaw[osc:/yaw 0 360]", 0, 0, 360, 0.01)*ma.PI/180; // Slider with yaw rotation angle
pitch   =   hslider("Pitch[osc:/picth 0 360]", 0, 0, 360, 0.01)*ma.PI/180; // Slider with pitch rotation angle
roll    =   hslider("Roll[osc:/roll 0 360]", 0, 0, 360, 0.01)*ma.PI/180; // Slider with roll rotation angle

// Maximum required order
M   =   4;

ins =   (M+1)^2;

// Zero-th order
rot(0,m,n) = 1;

// First order rotation matrix (n1, n2)
rot(1,-1,-1)    =   cos(roll)*cos(yaw) - sin(pitch)*sin(roll)*sin(yaw);
rot(1,-1,0)     =   -1*cos(pitch)*sin(roll);
rot(1,-1,1)     =   cos(yaw)*sin(pitch)*sin(roll) + cos(roll)*sin(yaw);
rot(1,0,-1)     =   cos(yaw)*sin(roll) + cos(roll)*sin(pitch)*sin(yaw);
rot(1,0,0)      =   cos(pitch)*cos(roll);
rot(1,0,1)      =   -1*cos(roll)*cos(yaw)*sin(pitch) + sin(roll)*sin(yaw);
rot(1,1,-1)     =   -1*cos(pitch)*sin(yaw);
rot(1,1,0)      =   sin(pitch);
rot(1,1,1)      =   cos(pitch)*cos(yaw);
rot(1,m,n)      =   0; // other cases for 1st order.

// Recurrence computation for higher-orders
denom(m,n2)     =   ba.if(abs(n2)<m,(m+n2)*(m-n2),2*m*(2*m-1));
u(m,n1,n2)      =   sqrt((m+n1)*(m-n1)/denom(m,n2));
v(m,n1,n2)      =   1/2*sqrt((1+(n1==0))*(m+abs(n1)-1)*(m+abs(n1))/denom(m,n2))*(1-2*(n1==0));
w(m,n1,n2)      =   -1/2*sqrt((m-abs(n1)-1)*(m-abs(n1))/denom(m,n2))*(1-(n1==0));

U(m,n1,n2)      =   ba.if(n1==0,P(0,m,0,n2),P(0,m,n1,n2));

V(m,n1,n2)      =   case{
                    (1,0,0)     =>P(1,m,1,n2)+P(-1,m,-1,n2);
                    (0,1,0)     =>P(1,m,n1-1,n2)*sqrt(1+(n1==1))-P(-1,m,-n1+1,n2)*(1-(n1==1));
                    (0,0,1)     =>P(1,m,n1+1,n2)*(1-(n1==-1))+P(-1,m,-n1-1,n2)*sqrt(1+(n1==-1)); // sqrt(1+(n1==1)) is right, in the correction of the paper it's sqrt(1-(n1==1))
                    }(n1==0,n1>0,n1<0);

W(m,n1,n2)      =   case{
                    (1,0,0)     => 0; // Shouldn't be defined but covers some pattern matching cases.
                    (0,1,0)     => P(1,m,n1+1,n2)+P(-1,m,-n1-1,n2);
                    (0,0,1)     => P(1,m,n1-1,n2)-P(-1,m,-n1+1,n2);
                    }(n1==0,n1>0,n1<0);

P(i,m,mu,n2)    =   case{
                    (1,0,0)     => rot(1,i,0)*rot(m-1,mu,n2);
                    (0,1,0)     => rot(1,i,1)*rot(m-1,mu,m-1)-rot(1,i,-1)*rot(m-1,mu,-m+1);
                    (0,0,1)     => rot(1,i,1)*rot(m-1,mu,-m+1)+rot(1,i,-1)*rot(m-1,mu,m-1);
                    }(abs(n2)<m,n2==m,n2==-m);

// Other cases
rot(m,n1,n2)    =   u(m,n1,n2)*U(m,n1,n2)+v(m,n1,n2)*V(m,n1,n2)+w(m,n1,n2)*W(m,n1,n2);

// Main-matrix row          
row(M,i)    =   par(m,M+1,
              par(j,2*m+1,term 
                with{term = ba.if((i >= m^2) & (i< (m+1)^2),rot(m,int(i-m^2)-m,j-m),0);}
                )
                );

// Matrix multiplication
// n = number of inputs
// m = number of outputs
matrix(n,m) = par(i,n,_) <: par(i,m,buswg(row(M,i)):>_);

//process = rot(M,-2,-1);
//process =   row(M,0);
process = matrix(ins,ins);

[orlarey]

Hi Pierre,

if(c,t,e) is very slow compared to pattern matching. Therefore I suggest replacing all remaining if in your code with the equivalent using pattern matching.

Let me know if that solves the problem...

Cheers

Yann

[sekisushai]

Hi Yann, Thank you for the suggestion. That solved it up to M=4 ! I've replaced ba.if occurrences with:

case{
                    (1)  => t;
                    (0)  => f;
                    }(test)

Unfortunately, the compilation never stops at order M=5. My guess is that the recurrence formula brings very long trigonometric formula, if intermediate simplifications are not done by the compiler. For example, the term rot(2,1,1) gives:

rot(2,1,1) = cos(pitch)^2*cos(roll)*cos(yaw) +
 sin(pitch)*(-cos(roll)*cos(yaw)*sin(pitch) + sin(roll)*sin(yaw))

Which, after factorization gives rot(2,1,1) = cos(roll)*cos(yaw)*(cos(pitch)^2 - sin(pitch)^2)) + sin(pitch)*sin(roll)*sin(yaw) However, knowing that cos(pitch)^2 - sin(pitch)^2 = cos(2*pitch) the latter expression could be simplified with: rot(2,1,1) = cos(roll)*cos(yaw)*cos(2*pitch) + sin(pitch)*sin(roll)*sin(yaw) which is a simpler expression. This kind of simplification can be done for majority of terms in the rotation matrix. I'm not sure if the compiler implements trigonometric simplifications, but this would be a great feature !