feat: first implementation of sox-hoadecoder

Alexander Nguyen

yesterday 4bfa4ffc65f2d908c2341feba9d69bb9f78c4417

feat: first implementation of sox-hoadecoder

 sources/classes/sox-hoadecode.lisp

New file
@@ -0,0 +1,775 @@
;Authors: 
; A. Nguyen, 2025

(in-package :om)

; Requirements:
; qhull

; Known limitations:
; 1) Signal order must match decoder order.
; 2) It appears that SoX remix's resolution is limited to about 1e-5.

(ql:quickload :lla)
(ql:quickload :array-operations)
(ql:quickload :str)

;======================================================
; Utility functions
;======================================================
(defun epsilon-to-zero (val) 
"Returns 0 if val's absolute value is less or equal 2x 'machine epsilon', otherwise val.

Example:
(cos (/ pi 2)) ; => 6.123233995736766d-17
(epsilon-to-zero (cos (/ pi 2))) ; => 0.0d0"
    (if (<= (abs val) (* 2 DOUBLE-FLOAT-EPSILON))
        0d0
        val
    ))

(defun sph2cart (sph)
"Converts a list of (azimuth colatitude) coordinates in radians to cartesian coordinates (x y z).
    
Example:
(sph2cart (list 0 0))               ; => '(0.0d0 0.0d0 1.0d0)
(sph2cart (list 0 (/ pi 2)))        ; => '(1.0d0 0.0d0 0.0d0)
(sph2cart (list 0 pi))              ; => '(0.0d0 0.0d0 -1.0d0)
(sph2cart (list (/ pi 2) (/ pi 2))) ; => '(0.0d0 1.0d0 0.0d0)
"
    (let*
        (
            (azimuth    (first sph))
            (colatitude (second sph))
            (x (* (cos azimuth) (sin colatitude)))
            (y (* (sin azimuth) (sin colatitude)))
            (z (cos colatitude))
        )
        (mapcar #'epsilon-to-zero (list x y z))))

(defun cart2sph (cart)
"Converts a list of (x y z) coordinates to spherical coordinates (azimuth elevation) in radians (without radius).
    
Example:
(cart2sph '(0 0 0)) ; => (0.0d0 1.5707963267948966d0)
(cart2sph '(1 0 0)) ; => (0.0d0 1.5707963267948966d0)
(cart2sph '(0 1 0)) ; => (1.5707963267948966d0 1.5707963267948966d0)
(cart2sph '(0 0 1)) ; => (0.0d0 0.0d0)
"
    (let*
        (
            (x (coerce (nth 0 cart) 'double-float))
            (y (coerce (nth 1 cart) 'double-float))
            (z (coerce (nth 2 cart) 'double-float))
            (azimuth     (atan y x))
            (colatitude (- (/ pi 2) (atan z (sqrt (+ (expt x 2) (expt y 2))))))
        )
        (mapcar #'epsilon-to-zero (list azimuth colatitude))))

(defun nth-row (arr row)
"Returns n-th row of arr as vector.

Example:
(nth-row #2A((1 0 0) (0 1 0) (0 0 1)) 1) ; => #(0 1 0)"
    (make-array (array-dimension arr 1) 
      :displaced-to arr 
       :displaced-index-offset (* row (array-dimension arr 1))))

(defun list-to-2d-array (list)
"Converts a list of lists (interpreted as list of rows) to a 2d array (#2A).

Example:
(list-to-2d-array '((1 0 0)))                 ; => #2A((1 0 0))
(list-to-2d-array '((1 0 0) (0 1 0) (0 0 1))) ; => #2A((1 0 0) (0 1 0) (0 0 1))"
  (make-array (list (length list) (length (first list))) :initial-contents list))

(defun number-to-string (num)
"Returns the string representation of num (interpreted as single-float).
    
Example:
(number-to-string 1)  => 1.0
(number-to-string pi) => 3.1415927"
    (substitute #\e #\d (format nil "~d" (coerce num 'single-float))))

(defun unassociated-legendre-polynomial (order x)
"Evaluates x using the unassociated Legendre polynomial of some order. 
The current implementation supports max. order 7.

Example:
(unassociated-legendre-polynomial 0 pi) ; => 1
(unassociated-legendre-polynomial 1 pi) ; => 3.141592653589793d0
(unassociated-legendre-polynomial 2 pi) ; => 14.304406601634037d0
"
    (ecase order
        (0 1)
        (1 x)
        (2 (- (/ (* 3 (expt x 2)) 2) (/ 1 2)))
        (3 (+ (/ (* 3 x (- (expt x 2) 1)) 2) (expt x 3)))
        (4 (+ (* 3 (expt x 2) (- (expt x 2) 1)) (/ (* 3 (expt (- (expt x 2) 1) 2)) 8) (expt x 4)))
        (5 (+ (/ (* 15 x (expt (- (expt x 2) 1) 2)) 8) (* 5 (expt x 3) (- (expt x 2) 1)) (expt x 5)))
        (6 (+ (/ (* 15 (expt x 4) (- (expt x 2) 1)) 2) (/ (* 5 (expt (- (expt x 2) 1) 3)) 16) (/ (* 45 (expt x 2) (expt (- (expt x 2) 1) 2)) 8) (expt x 6)))
        (7 (+ (/ (* 35 x (expt (- (expt x 2) 1) 3)) 16) (/ (* 21 (expt x 5) (- (expt x 2) 1)) 2) (/ (* 105 (expt x 3) (expt (- (expt x 2) 1) 2)) 8) (expt x 7)))
    ))

;======================================================
; Internal functions
;======================================================

; VBAP
(defun vbap-convex-hull (points)
"Returns the 3D convex hull of a list of cartesian coordinates. 
External dependency: qhull
    
Example:

(vbap-convex-hull '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1))) 
=> ((0 1 4) (0 1 5) (0 3 4) (0 3 5) (1 2 4) (1 2 5) (2 3 4) (2 3 5))
"
    (let*
        (
            (qconvex-dimension 3)
            (qconvex-points-count (length points))
            (qconvex-points-list (map 'list #'(lambda (e) (map 'list #'(lambda (ee) (format nil "~15$" ee)) e)) points))
            (qconvex-points-list-str (with-output-to-string (stream) (dolist (point qconvex-points-list) (write-line (format nil "~{~A~^ ~}" point) stream))))
            (qconvex-cli-input-str (with-output-to-string (stream) 
                (progn 
                    (write-line (write-to-string qconvex-dimension) stream)
                    (write-line (write-to-string qconvex-points-count) stream)
                    (write-line qconvex-points-list-str stream)
                )))
             (qconvex-cli-output-str (with-output-to-string (out)
                (with-input-from-string (in qconvex-cli-input-str)
                    (uiop:run-program "qconvex Qt i" :output out :input in))))
        )
        (sort
            (map 'list 
                #'(lambda (e) 
                    (sort 
                        (map 'list 
                            (lambda (e) (nth-value 0 (parse-integer e))) 
                            (str:split " " (str:trim e))) 
                        #'<)) 
                (cdr (str:lines qconvex-cli-output-str)))
            #'<
            :key (lambda (e) (+ (* 1e6 (first e)) (* 1e3 (second e)) (* 1e0 (third e)))) ; (1 23 456) => 001023456,
        )))
        

(defun vbap-make-triangle_to_inverse_matrix (points)
"Computes the convex hull of given (cartesian) points, then returns a map of triangles (i.e., the index value of each point) to the triangle's inverse matrix (when interpreted as vector base).

Example:
    
(vbap-make-triangle_to_inverse_matrix '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1))) 
=> (((0 1 4) #2A((1.0d0 0.0d0 0.0d0) (-0.0d0 1.0d0 0.0d0) (-0.0d0 -0.0d0 1.0d0)))
    ((0 1 5) #2A((1.0d0 0.0d0 0.0d0) (-0.0d0 1.0d0 0.0d0) (0.0d0 0.0d0 -1.0d0)))
    ((0 3 4) #2A((1.0d0 0.0d0 0.0d0) (0.0d0 -1.0d0 0.0d0) (-0.0d0 -0.0d0 1.0d0)))
    ((0 3 5) #2A((1.0d0 0.0d0 0.0d0) (0.0d0 -1.0d0 0.0d0) (0.0d0 0.0d0 -1.0d0)))
    ((1 2 4) #2A((-0.0d0 1.0d0 0.0d0) (-1.0d0 0.0d0 0.0d0) (-0.0d0 -0.0d0 1.0d0)))
    ((1 2 5) #2A((-0.0d0 1.0d0 0.0d0) (-1.0d0 0.0d0 0.0d0) (0.0d0 0.0d0 -1.0d0)))
    ((2 3 4) #2A((-1.0d0 0.0d0 0.0d0) (0.0d0 -1.0d0 0.0d0) (-0.0d0 -0.0d0 1.0d0)))
    ((2 3 5) #2A((-1.0d0 0.0d0 0.0d0) (0.0d0 -1.0d0 0.0d0) (0.0d0 0.0d0 -1.0d0))))
"
    (let* 
        (
            (triangles (vbap-convex-hull points))
            (triangle_to_inverse_matrix 
                (map 
                    'list
                    (lambda (tri) 
                        (list 
                            tri 
                            (aops:permute '(1 0) (lla:invert (list-to-2d-array (map 'list (lambda (p) (nth p points)) tri)))))) 
                    triangles))
        )
        triangle_to_inverse_matrix
    ))

(defun vbap-triangulate (triangle_to_inverse_matrix real_positions_length real_positions_imaginary_idxs  azimuth colatitude)
"
Returns a vector of length real_positions_length with at most three non-zero elements.
It represents the gain values for use in VBAP (according to given azimuth and colatitude), i.e. its length is equal to 1 (except when imaginary positions are involved) and the n-th may be interpreted as the gain value for the n-th loudspeaker.
    
Example:
    
(vbap-triangulate (vbap-make-triangle_to_inverse_matrix '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1))) 6 '() 0 0)
=> #(0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0)
    
(vbap-triangulate (vbap-make-triangle_to_inverse_matrix '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1))) 6 '() 0 (/ pi 2))
=> #(1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0)

(vbap-triangulate (vbap-make-triangle_to_inverse_matrix '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1))) 6 '() (/ pi 4) (/ pi 2))
=> #(0.7071067811865476d0 0.7071067811865475d0 0.0d0 0.0d0 0.0d0 0.0d0)

(vbap-triangulate (vbap-make-triangle_to_inverse_matrix '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1))) 0 '() (/ pi 4) (/ pi 2))
=> ((0 1 4) (0.7071067811865476d0 0.7071067811865475d0 0.0d0))
"
    (let* 
        (
            (position (vector (* (cos azimuth) (sin colatitude)) (* (sin azimuth) (sin colatitude)) (cos colatitude)))
            (best_triangle_to_inverse_matrix 
                (loop for i from 0 to (1- (list-length triangle_to_inverse_matrix))
                    do (let* 
                        (
                            (triangle       (first  (nth i triangle_to_inverse_matrix)))
                            (inverse_matrix (second (nth i triangle_to_inverse_matrix)))
                            (gains (mapcar #'epsilon-to-zero (coerce (lla:mm inverse_matrix position) 'list)))
                            (min-gain (apply #'min gains))
                        )
                        (if (>= min-gain 0)
                            (return (list triangle gains)))
                    )
                )
            )
            (best_triangle_to_inverse_matrix_normalized 
                (list (first best_triangle_to_inverse_matrix)
                      (mapcar #'(lambda (v) (/ v (lla:nrm2 (coerce (second best_triangle_to_inverse_matrix) 'vector)))) (second best_triangle_to_inverse_matrix))
                ))
            (vec (make-array real_positions_length :initial-element 0.0d0))
        )
        (let*
            (
                (lsp_idxs (first  best_triangle_to_inverse_matrix_normalized))
                (gains    (second best_triangle_to_inverse_matrix_normalized))
                (imaginaryLspIdxs (intersection lsp_idxs real_positions_imaginary_idxs))
            )
            (progn
                ;(print (list "Triangulating" best_triangle_to_inverse_matrix_normalized))
                (if (> (list-length imaginaryLspIdxs) 0)
                    (let*
                        (
                            (imaginaryLspIdx (first imaginaryLspIdxs))
                            (imagGainIdx (position imaginaryLspIdx lsp_idxs))
                            (realGainIndex (remove imagGainIdx '(0 1 2)))
                            (connectedTriangles (loop for k from 0 to (1- (list-length triangle_to_inverse_matrix)) collect
                                (let*
                                    (
                                        (probe (first (nth k triangle_to_inverse_matrix)))   
                                    )
                                    (progn
                                        ;(print probe)
                                        (if (find imaginaryLspIdx probe)
                                            probe
                                        )                 
                                    )
                                )
                            ))
                            (connectedLsps (remove imaginaryLspIdx (remove-duplicates (flatten connectedTriangles))))
                            (kappa (getKappa (nth imagGainIdx gains) (nth (first realGainIndex) gains) (nth (second realGainIndex) gains) (list-length connectedLsps)))
                            (gainVector (make-array (list-length connectedLsps) :initial-element (* kappa (nth imagGainIdx gains))))
                        )
                        (progn
                            ;(print (list "Connected lsps:" connectedLsps))
                            (loop for j from 0 to 1 do
                                (let*
                                    (
                                        (con_idx (position (nth (nth j realGainIndex) lsp_idxs) connectedLsps))
                                    )
                                    (progn
                                        ;(print (list "Set gainVector" con_idx " = " (aref gainVector con_idx) " + " (nth (nth j realGainIndex) gains)))
                                        (setf (aref gainVector con_idx) (+ (aref gainVector con_idx) (nth (nth j realGainIndex) gains)))
                                    )
                                )
                            )
                            ;(print gainVector)
                            ;(print (list "Using kappa'd gains for " connectedLsps ":" gainVector))
                            (loop for n from 0 to (1- (list-length connectedLsps)) do
                                (progn
                                    (setf (aref vec (nth n connectedLsps)) (+ (aref vec (nth n connectedLsps)) (aref gainVector n))) 
                                )
                            )
                            vec
                        )
                    )
                    (progn
                        (loop for i from 0 to 2 do 
                            (setf 
                                (aref vec (nth i (first best_triangle_to_inverse_matrix))) 
                                ; (if (find (nth i (first best_triangle_to_inverse_matrix)) real_positions_imaginary_idxs)
                                ;     0
                                ;     (nth i (second best_triangle_to_inverse_matrix_normalized)))
                                (+ (aref vec (nth i (first best_triangle_to_inverse_matrix))) (nth i (second best_triangle_to_inverse_matrix_normalized)))
                            ))
                        vec
                    )
                )
            )
        )
    ))

;(print (vbap-triangulate (vbap-make-triangle_to_inverse_matrix (first (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs))) (list-length (first (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs))) (second (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs))  (/ pi 4) (/ pi 1)))


; AllRAD

(defun sox-hoadecode-t5200-spherical-positions ()
"Returns a list of 5200 (azimuth colatitude) coordinates of a t-design."
    (read-from-string (uiop:read-file-string "sources/classes/t5200.txt")))

(defun getKappa (gIm gRe1 gRe2 N)
"Example:
    
(getKappa 0.33175008679632517d0 0.8352087169013385d0 0.4385980838109653d0 12)
=> 0.11097524644830536d0
"
    (let* ((p (/ (* gIm (+ gRe1 gRe2)) (* N (expt gIm 2d0))))
           (q (/ (- (+ (expt gRe1 2d0) (expt gRe2 2d0)) 1)(* N (expt gIm 2d0)))))
        (+ (- p) (sqrt (max (- (expt p 2) q) 0)))
    ))

(defun allrad-vbap-matrix (virtual_spherical_positions real_cartesian_positions real_cartesian_positions_imaginary_idxs)
"Returns a L times J matrix to render J virtual_spherical_positions using VBAP and L real_cartesian_positions. 
Each of the J columns contains the VBAP gains for the j-th virtual position. Conversely, each of the L rows represents the linear combination of all J loudspeaker (signals), for which the l-th loudspeaker is responsible (in terms of VBAP).
    
Example:
    
(allrad-vbap-matrix (mapcar #'cart2sph '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1))) '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '())
=>  #2A((1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0)
        (0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0)
        (0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0)
        (0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0)
        (0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0)
        (0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0))

(allrad-vbap-matrix (mapcar #'cart2sph '((1 1 0))) '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '())
=>  #2A((0.7071067811865476d0)
        (0.7071067811865475d0)
        (0.0d0)
        (0.0d0)
        (0.0d0)
        (0.0d0))
    
(allrad-vbap-matrix (subseq (sox-hoadecode-t5200-spherical-positions) 0 1) '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '())
=>  #2A((0.0d0)
        (0.8575241661297169d0)
        (0.0072383074899637885d0)
        (0.0d0)
        (0.514392759871496d0)
        (0.0d0))
"
    (let*
        (
             (vbap-triangle-to-inverse-matrix (vbap-make-triangle_to_inverse_matrix real_cartesian_positions))
        )
        (aops:permute '(1 0) (list-to-2d-array (mapcar (lambda (j) (vbap-triangulate vbap-triangle-to-inverse-matrix (list-length real_cartesian_positions) real_cartesian_positions_imaginary_idxs (first j) (second j))) virtual_spherical_positions)))
    ))

(defun allrad-maxre-factor (order truncation_order)
"Returns a factor to multiply each of the channels of the ambisonics matrix with if using decoder flavor max-rE.
Values are taken from the IEM Plug-In Suite to be closer to their AllRAD implementation.
    
Example:
    
(allrad-maxre-factor 0 3)
=> 1.d0
    
(allrad-maxre-factor 1 3)
=> 0.8615507588765864d0
        
(allrad-maxre-factor 2 3)
=> 0.613404565181233d0
    
(allrad-maxre-factor 3 3)
=> 0.3064314417993654d0
"
    (let* 
        (
            (vals '((1.0d0)
                (1.0d0 5.7754104119288496d-01)
                (1.0d0 7.7520766107019334d-01 4.0142037667287966d-01)
                (1.0d0 8.6155075887658639d-01 6.1340456518123299d-01 3.0643144179936538d-01)
                (1.0d0 9.0644136637224459d-01 7.3245392600617265d-01 5.0224998490808703d-01 2.4736484001129033d-01)
                (1.0d0 9.3263709143129281d-01 8.0471791647013236d-01 6.2909156744472861d-01 4.2321128963220900d-01 2.0719132924646289d-01)
                (1.0d0 9.4921830632793713d-01 8.5152308960211620d-01 7.1432330396679700d-01 5.4794300713180655d-01 3.6475291657556469d-01 1.7813609450688817d-01)
                (1.0d0 9.6036452263662697d-01 8.8345002450861454d-01 7.7381375334313540d-01 6.3791321433685355d-01 4.8368159255186721d-01 3.2000849790781744d-01 1.5616185043093761d-01)))
        )
        (if (and (<= 0 order) (<= order truncation_order) (<= truncation_order 7))
            (nth order (nth truncation_order vals))
            nil
        )
    )
)


(defun allrad-maxre-energy-correction (truncation_order)
"Returns a correction factor to multiply the ambisonics matrix with if using decoder flavor max-rE.
    
Example:

(allrad-maxre-energy-correction 0)
=> 1.0d0
    
(allrad-maxre-energy-correction 1)
=> 1.061114703d0
"
    (if (and (<= 0 truncation_order) (<= truncation_order 7))
        (nth truncation_order '(1.0d0 1.061114703d0 1.083513331d0 1.095089593d0 1.102148983d0 1.106899961d0 1.110314372d0 1.112886124d0))
        nil)
)

(defun allrad-ambisonics-matrix (virtual_positions truncation_order normalization decoder_flavor)
"Returns a (list-length virtual_positions) times (expt (+ truncation_order 1) 2) matrix, which represents a higher-order Ambisonics sampling decoder.
In the case of AllRAD, the virtual_positons origin from a t-design. 

Parameters:
- truncation_order : 0, 1, ...
- normalization : :sn3d, :n3d
- decoder_flavor : :maxre, :basic (otherwise)
    
Example:
    
(allrad-ambisonics-matrix  '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) 3 :sn3d :basic)
=>  #2A((0.06249999999999999d0 0.0d0 0.10825317547305481d0 0.0d0 0.0d0 0.0d0
        0.13975424859373684d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.1653594569415369d0
        0.0d0 0.0d0 0.0d0)
        (0.06249999999999999d0 0.0d0 0.05848944032563985d0 0.09109190617389346d0
        0.0d0 0.0d0 -0.008680154186295676d0 0.11005293096391006d0
        0.0856986424020531d0 0.0d0 0.0d0 0.0d0 -0.06881135255279784d0
        0.03916471154585938d0 0.12250668358046454d0 0.07789085702121995d0)
        (0.06249999999999999d0 0.0d0 0.10825317547305481d0 0.0d0 0.0d0 0.0d0
        0.13975424859373684d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.1653594569415369d0
        0.0d0 0.0d0 0.0d0)
        (0.06249999999999999d0 0.0d0 0.058489440325639866d0 -0.09109190617389346d0
        0.0d0 0.0d0 -0.008680154186295676d0 -0.11005293096391006d0
        0.0856986424020531d0 0.0d0 0.0d0 0.0d0 -0.06881135255279783d0
        -0.03916471154585938d0 0.12250668358046451d0 -0.07789085702121995d0)
        (0.06249999999999999d0 0.0d0 0.10825317547305481d0 0.0d0 0.0d0 0.0d0
        0.13975424859373684d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.1653594569415369d0
        0.0d0 0.0d0 0.0d0)
        (0.06249999999999999d0 0.0d0 0.10825317547305481d0 0.0d0 0.0d0 0.0d0
        0.13975424859373684d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.1653594569415369d0
        0.0d0 0.0d0 0.0d0))
    
(allrad-ambisonics-matrix (subseq (sox-hoadecode-t5200-spherical-positions) 0 3) 0 :sn3d :basic)
=> #2A((0.9999999999999999d0) (0.9999999999999999d0) (0.9999999999999999d0))
    
(allrad-ambisonics-matrix (subseq (sox-hoadecode-t5200-spherical-positions) 0 1) 3 :sn3d :basic)
=>  #2A((0.06249999999999999d0 0.09282971402842526d0 0.055684649696437986d0
        -7.83569770838977d-4 -0.0015024811881936859d0 0.10677433172348148d0
        -0.01440868160163194d0 -9.012754107422361d-4 -0.08899332689825591d0
        -0.0824167168518417d0 -0.0020448097645174353d0 0.028047423091425824d0
        -0.07132258351652541d0 -2.3674653223258648d-4 -0.12111594158274905d0
        0.0020874195679599084d0))
"
    (let* 
        (
            (js (make-array (list (list-length virtual_positions) (expt (+ truncation_order 1) 2)) :initial-contents (loop for i from 0 to (1- (list-length virtual_positions))
                    collect (loop for acn from 0 to (- (expt (+ truncation_order 1) 2) 1) 
                        collect (let*
                            (
                                (order (isqrt acn))
                                (degree (- acn (expt order 2) order))
                                (j_azimuth    (first  (nth i virtual_positions)))
                                (j_colatitude (second (nth i virtual_positions)))
                                (azimuth_deg   (* -1 (sox-hoaencode-rad-to-deg j_azimuth)))
                                (elevation_deg (sox-hoaencode-rad-to-deg (- (/ pi 2) j_colatitude)))
                                ;(correction (if (eq normalization :sn3d) 1 (sqrt (+ (* 2 truncation_order) 1d0))))
                                (correction (* (sqrt (* 4 pi)) (/ 1 (expt (+ truncation_order 1) 2))))
                                (flavor-correction (if (eq decoder_flavor :maxre) 
                                    (* (allrad-maxre-energy-correction truncation_order) (allrad-maxre-factor order truncation_order)) 
                                    1))
                            )
                            (progn 
                                ;(print (list i acn order degree j_azimuth j_colatitude azimuth_deg elevation_deg correction (sox-hoaencode-gain-single-component order degree azimuth_deg elevation_deg :n3d)))
                                (epsilon-to-zero (* flavor-correction correction (sox-hoaencode-gain-single-component order degree azimuth_deg elevation_deg :n3d))) 
                            )
                        ))
                ))
            )
        )
        js
    ))

(defun allrad-decoder-matrix (truncation_order loudspeakers loudspeakers_imaginary_idxs normalization decoder-flavor)
"Returns the (list-length loudspeakers) times (expt (+ truncation_order 1) 2) matrix to decode a truncation_order-th Ambisonics signal using AllRAD and the given loudspeakers layout.
    
Example:

(allrad-decoder-matrix 0 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :basic)
=>  #2A((0.4082521539974241d0)
        (0.4082497204809079d0)
        (0.408252153997083d0)
        (0.4082497204819567d0)
        (0.40824299685796434d0)
        (0.4082429968576299d0))

(allrad-decoder-matrix 1 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :n3d :basic)
=>  #2A((0.26726180898030116d0 -4.121735356202743d-6 3.1176420689701896d-6 0.30860515455826876d0)
        (0.26726021588147736d0 0.30860530309758993d0 -5.841208957330852d-6 -4.5668569247884075d-6)
        (0.2672618089800778d0 -4.12173563713413d-6 3.1176426062467998d-6 -0.3086053865522143d0)
        (0.26726021588216414d0 -0.3086052380121168d0 -5.8412091936397d-6 -4.5668572057260935d-6)
        (0.2672558142693533d0 3.0372719711827937d-6 0.30860537041237485d0 -1.5744514926730053d-6)
        (0.26725581426913414d0 3.0372717348746766d-6 -0.30860517069778515d0 -1.5744509553960298d-6))

(allrad-decoder-matrix 1 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :n3d :maxre) 
=>  #2A((0.33964966990487044d0 -3.025221988304797d-6 2.2882496142125288d-6 0.22650631799273105d0)
        (0.33964764531523056d0 0.2265064270158442d0 -4.287260643603554d-6 -3.3519269900538543d-6)
        (0.3396496699045866d0 -3.025222194500096d-6 2.28825000855551d-6 -0.22650648826887057d0)
        (0.33964764531610314d0 -0.22650637924518757d0 -4.287260817050707d-6 -3.351927196251641d-6)
        (0.3396420515264632d0 2.2292605316972416d-6 0.22650647642274438d0 -1.155597063743201d-6)
        (0.33964205152618476d0 2.229260358254026d-6 -0.22650632983861987d0 -1.1555966693972261d-6))
    
(allrad-decoder-matrix 3 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :basic)
=>  #2A((0.23132635066908433d0 -3.567535526350099d-6 2.698450501623429d-6
        0.26711075732820044d0 -4.678820125933257d-8 -9.537977397150688d-7
        -0.06466023863439527d0 -4.0231181660137245d-8 0.1119858835087529d0
        5.279006872057261d-6 -2.2151444156856834d-8 7.072629025934007d-6
        1.0091958482011735d-6 -3.633471917485323d-7 -2.410565206410929d-6
        9.971727398691361d-9)
        (0.23132497177496347d0 0.2671108858952382d0 -5.055812339034126d-6
        -3.952806989031822d-6 -9.956844379916084d-8 9.411227751817857d-8
        -0.06466012701108666d0 -6.234428562016533d-6 -0.11198772931832059d0
        2.641464136965804d-7 -1.561448094516174d-7 -2.3553869555088308d-7
        4.481201793886413d-6 -3.842131666028257d-6 -9.45568891970116d-6
        -2.8381531449148024d-6)
        (0.2313263506688911d0 -3.567535769509013d-6 2.698450966658684d-6
        -0.2671109581287301d0 -4.678708503022157d-8 -9.537977516266857d-7
        -0.06466023863445543d0 -4.023117363675782d-8 0.11198588350869225d0
        5.27900697040363d-6 -2.2151467699660293d-8 7.0726289406049545d-6
        1.009195501156908d-6 -3.6334705642359096d-7 -2.410564548297249d-6
        9.971351015858804d-9)
        (0.23132497177555786d0 -0.2671108295610204d0 -5.055812543567884d-6
        -3.952807232196366d-6 -9.956859800096157d-8 9.411186784898843d-8
        -0.06466012701137755d0 -6.234428573927192d-6 -0.1119877293181395d0
        2.6414607633245387d-7 -1.56144648514084d-7 -2.355388105676657d-7
        4.481201777403412d-6 -3.842131751356731d-6 -9.4556887140023d-6
        -2.8381526208161127d-6)
        (0.2313211619943082d0 2.6288867974199045d-6 0.2671109441590194d0
        -1.3627540706050697d-6 -4.443157119897495d-6 -8.873940212880686d-9
        0.12931854808913654d0 2.655927037473874d-7 -1.3268000096613416d-6
        8.756408995091749d-7 3.3791429901405865d-7 -2.464371510682579d-6
        -1.9828624970594008d-7 1.2768740309038784d-6 2.749060859701345d-8
        9.98925916570816d-6)
        (0.2313211619941186d0 2.628886592885831d-6 -0.2671107712976314d0
        -1.3627536055677794d-6 -4.443157131807048d-6 -8.873508852416418d-9
        0.12931854808906065d0 2.6559238474805657d-7 -1.3267998008155014d-6
        8.756409738692105d-7 3.379139033873876d-7 -2.46437191352441d-6
        -1.9828635214785303d-7 1.276874475867459d-6 2.7490369427357343d-8
        9.989259368678423d-6))
"
    (let*
        (
            (vbap-matrix        (allrad-vbap-matrix       (sox-hoadecode-t5200-spherical-positions) loudspeakers loudspeakers_imaginary_idxs))
            (ambisonics-matrix  (allrad-ambisonics-matrix (sox-hoadecode-t5200-spherical-positions) truncation_order normalization decoder-flavor))
            (allrad-matrix      (lla:mm vbap-matrix ambisonics-matrix))
            (max-gains (loop for i from 0 to (1- (list-length loudspeakers))
                collect (if (find i loudspeakers_imaginary_idxs)
                    0
                    (let* 
                        (   
                            (l_azimuth    (first  (cart2sph (nth i loudspeakers))))
                            (l_colatitude (second (cart2sph (nth i loudspeakers))))
                            (azimuth_deg   (* -1d0 (sox-hoaencode-rad-to-deg l_azimuth)))
                            (elevation_deg (sox-hoaencode-rad-to-deg (- (/ pi 2) l_colatitude)))
                            (sh_raw (sox-hoaencode-gains-up-to-order truncation_order azimuth_deg elevation_deg :n3d))
                            (sh (loop for acn from 0 to (1- (list-length sh_raw)) collect
                                    (let* 
                                        (
                                            (order (isqrt acn))
                                            (correction 
                                                (if (eq decoder-flavor :maxre)
                                                    (* (sqrt (* 4 pi)) (/ 1d0 (allrad-maxre-energy-correction truncation_order)) (/ 1d0 (allrad-maxre-factor order truncation_order)))
                                                    (sqrt (* 4 pi))
                                                ))
                                            (e (nth acn sh_raw))
                                        )
                                        (epsilon-to-zero (* correction e))
                                    )
                                )
                            )
                            (sumOfSquares 
                                (sqrt (reduce #'+ (mapcar (lambda (e) (expt e 2d0)) (loop for m from 0 to (1- (aops:nrow allrad-matrix)) collect
                                    ;(print 
                                        (reduce #'+ (loop for n from 0 to (1- (aops:ncol allrad-matrix)) collect
                                            (progn
                                                ;(print 
                                                    (list i (nth i loudspeakers) m n (nth n sh) (aref allrad-matrix m n))
                                                ;) 
                                                    (* (nth n sh) (aref allrad-matrix m n))
                                            )
                                        ))
                                    ;)
                                ))))
                            )
                        )
                        sumOfSquares
                    )
                )))
            (max-gain (apply #'max max-gains))
        )
        (progn
            ;(print max-gains)
            ;(print max-gain)
            (aops:each-index (i j)
               (/ (aref allrad-matrix i j) max-gain))
        )
    ))

(defun allrad-decode-signal-remix-parameters (signal_normalization truncation_order loudspeakers loudspeakers_imaginary_idxs normalization decoder-flavor)
"Returns the parameters for SoX's remix effect.
    
Example:

(allrad-decode-signal-remix-parameters :sn3d 0 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :basic)
=> 1v0.40825215 1v0.4082497 1v0.40825215 1v0.4082497 1v0.408243 1v0.408243
    
(allrad-decode-signal-remix-parameters :n3d 0 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :n3d :basic)
=> 1v0.40825215 1v0.4082497 1v0.40825215 1v0.4082497 1v0.408243 1v0.408243

(allrad-decode-signal-remix-parameters :sn3d 1 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :basic)
=> 1v0.2672618,2v-7.139055e-6,3v5.3999142e-6,4v0.5345198 1v0.26726022,2v0.5345201,3v-1.0117271e-5,4v-7.910028e-6 1v0.2672618,2v-7.1390555e-6,3v5.3999156e-6,4v-0.5345202 1v0.26726022,2v-0.53452,3v-1.0117271e-5,4v-7.910029e-6 1v0.2672558,2v5.2607093e-6,3v0.5345202,4v-2.72703e-6 1v0.2672558,2v5.260709e-6,3v-0.53451985,4v-2.7270291e-6
    
(allrad-decode-signal-remix-parameters :n3d 1 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :n3d :basic)
=> 1v0.33964968,2v-5.239838e-6,3v3.9633646e-6,4v0.39232045 1v0.33964765,2v0.39232063,3v-7.425753e-6,4v-5.805708e-6 1v0.33964968,2v-5.2398386e-6,3v3.963365e-6,4v-0.39232075 1v0.33964765,2v-0.39232054,3v-7.4257537e-6,4v-5.805708e-6 1v0.33964205,2v3.8611925e-6,3v0.39232072,4v-2.0015527e-6 1v0.33964205,2v3.861192e-6,3v-0.39232048,4v-2.001552e-6
   
(allrad-decode-signal-remix-parameters :sn3d 1 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :maxre)
=> 1v0.33964968,2v-5.239838e-6,3v3.9633646e-6,4v0.39232045 1v0.33964765,2v0.39232063,3v-7.425753e-6,4v-5.805708e-6 1v0.33964968,2v-5.2398386e-6,3v3.963365e-6,4v-0.39232075 1v0.33964765,2v-0.39232054,3v-7.4257537e-6,4v-5.805708e-6 1v0.33964205,2v3.8611925e-6,3v0.39232072,4v-2.0015527e-6 1v0.33964205,2v3.861192e-6,3v-0.39232048,4v-2.001552e-6
    
(allrad-decode-signal-remix-parameters :n3d 1 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :n3d :maxre)
=>  1v0.33964968,2v-3.025222e-6,3v2.2882496e-6,4v0.22650632 1v0.33964765,2v0.22650643,3v-4.2872607e-6,4v-3.351927e-6 1v0.33964968,2v-3.0252222e-6,3v2.28825e-6,4v-0.22650649 1v0.33964765,2v-0.22650638,3v-4.2872607e-6,4v-3.3519273e-6 1v0.33964205,2v2.2292606e-6,3v0.22650647,4v-1.155597e-6 1v0.33964205,2v2.2292604e-6,3v-0.22650632,4v-1.1555967e-6

"
    (let* 
        (
            (decoder-matrix (allrad-decoder-matrix truncation_order loudspeakers loudspeakers_imaginary_idxs normalization decoder-flavor))
            (signal-normalization-correction 
                (let*
                    (
                        (components-count (expt (+ truncation_order 1) 2))
                        (matr (make-array (list components-count components-count)))
                    )
                    (progn
                        (loop for acn from 0 to (1- components-count) do
                            (let*
                                (
                                    (order (isqrt acn))
                                )
                                (if (eq signal_normalization :sn3d) 
                                    (setf (aref matr acn acn) (sqrt (+ (* 2 order) 1d0)))
                                    (setf (aref matr acn acn) 1d0)
                                )
                            )
                        )
                        matr
                    )
            ))
            (signal-decoder-matrix (lla:mm decoder-matrix signal-normalization-correction))
        )
        (let*
            (
                (str "")
            )
            (progn 
                (loop for loudspeaker from 0 to (1- (aops:nrow signal-decoder-matrix)) do
                    (progn
                        (let*
                            (
                                (channel-to-gain-list (loop for channel from 1 to (expt (+ truncation_order 1) 2) collect
                                    (let* 
                                        (
                                            (gain (aref signal-decoder-matrix loudspeaker (1- channel)))
                                        )
                                        (format nil "~av~d" channel (number-to-string gain))
                                    )
                                ))
                            )
                            (setf str (str:concat str (format nil "~{~A~^,~}" channel-to-gain-list)))
                        )
                    )
                )
                str
            )
        )
    )
)

; Others

(defun sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs ()
"Returns a list of two lists: the first containing cartesian coordinates of the 61 real loudspeakers (of the 3D Sound Lab (room MUT 210) of the ComputerStudio at Karlsruhe University of Music), 
plus 1 virtual loudspeaker, the second list containing the index value of the virtual loudspeaker."
    (list 
        (mapcar #'sph2cart
            '((0.0d0 1.5707963267948966d0)
            (-0.39269908169872414d0 1.5707963267948966d0)
            (-0.7853981633974483d0 1.5707963267948966d0)
            (-1.1780972450961724d0 1.5707963267948966d0)
            (-1.5707963267948966d0 1.5707963267948966d0)
            (-1.9634954084936207d0 1.5707963267948966d0)
            (-2.356194490192345d0 1.5707963267948966d0)
            (-2.748893571891069d0 1.5707963267948966d0)
            (-3.141592653589793d0 1.5707963267948966d0)
            (-3.5342917352885173d0 1.5707963267948966d0)
            (-3.9269908169872414d0 1.5707963267948966d0)
            (-4.319689898685966d0 1.5707963267948966d0)
            (-4.71238898038469d0 1.5707963267948966d0)
            (-5.105088062083414d0 1.5707963267948966d0)
            (-5.497787143782138d0 1.5707963267948966d0)
            (-5.8904862254808625d0 1.5707963267948966d0)
            (0.0d0 1.2566370614359172d0)
            (-0.5235987755982988d0 1.2566370614359172d0)
            (-1.0471975511965976d0 1.2566370614359172d0)
            (-1.5707963267948966d0 1.2566370614359172d0)
            (-2.0943951023931953d0 1.2566370614359172d0)
            (-2.6179938779914944d0 1.2566370614359172d0)
            (-3.141592653589793d0 1.2566370614359172d0)
            (-3.6651914291880923d0 1.2566370614359172d0)
            (-4.1887902047863905d0 1.2566370614359172d0)
            (-4.71238898038469d0 1.2566370614359172d0)
            (-5.235987755982989d0 1.2566370614359172d0)
            (-5.759586531581287d0 1.2566370614359172d0)
            (0.0d0 0.9773843811168245d0)
            (-0.7504915783575618d0 0.9773843811168245d0)
            (-1.5707963267948966d0 0.9250245035569946d0)
            (-2.3736477827122884d0 0.9773843811168245d0)
            (-3.141592653589793d0 0.9250245035569946d0)
            (-3.961897402027128d0 0.9773843811168245d0)
            (-4.71238898038469d0 0.9599310885968813d0)
            (-5.550147021341968d0 0.9773843811168245d0)
            (-0.47123889803846897d0 0.7155849933176751d0)
            (-1.0821041362364843d0 0.7155849933176751d0)
            (-2.0420352248333655d0 0.7155849933176751d0)
            (-2.652900463031381d0 0.7330382858376183d0)
            (-3.6302848441482056d0 0.7155849933176751d0)
            (-4.241150082346221d0 0.7155849933176751d0)
            (-5.201081170943102d0 0.7155849933176751d0)
            (-5.794493116621174d0 0.6981317007977318d0)
            (0.0d0 0.5235987755982989d0)
            (-1.5707963267948966d0 0.5235987755982989d0)
            (-3.141592653589793d0 0.5235987755982989d0)
            (-4.71238898038469d0 0.5235987755982989d0)
            (0.0d0 0.0d0)
            (0.0d0 2.0594885173533086d0)
            (-0.5235987755982988d0 2.0594885173533086d0)
            (-1.0471975511965976d0 2.0594885173533086d0)
            (-1.5707963267948966d0 2.0594885173533086d0)
            (-2.0943951023931953d0 2.0594885173533086d0)
            (-2.6179938779914944d0 2.0594885173533086d0)
            (-3.141592653589793d0 2.0594885173533086d0)
            (-3.6651914291880923d0 2.0594885173533086d0)
            (-4.1887902047863905d0 2.0594885173533086d0)
            (-4.71238898038469d0 2.0594885173533086d0)
            (-5.235987755982989d0 2.0594885173533086d0)
            (-5.759586531581287d0 2.0594885173533086d0)
            (0.0d0 3.141592653589793d0)
        ))
        '(61)))

; TESTS
; (allrad-vbap-matrix (subseq (sox-hoadecode-t5200-spherical-positions) 0 4) (first (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) (second (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)))
; (allrad-ambisonics-matrix (subseq (sox-hoadecode-t5200-spherical-positions) 0 4) 3 :sn3d :basic)
; (allrad-decoder-matrix 3 (first (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) (second (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) :sn3d :basic)

; 3rd-basic-mut210        : (print (allrad-decode-signal-remix-parameters :n3d 3 (first (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) (second (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) :sn3d :basic))
; 3rd-basic-octahedron    : (print (allrad-decode-signal-remix-parameters :n3d 3 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :basic))
; 3rd-maxre-mut210        : (print (allrad-decode-signal-remix-parameters :n3d 3 (first (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) (second (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) :sn3d :maxre))
; 3rd-maxre-octahedron    : (print (allrad-decode-signal-remix-parameters :n3d 3 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :maxre))

; TEST FILES
; ... using !!! SN3D !!! input ... 
; Use N3D if you want to reproduce IEM's matrix.
; 3rd-basic-mut210        : (print (allrad-decode-signal-remix-parameters :sn3d 3 (first (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) (second (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) :sn3d :basic))
; 3rd-basic-octahedron    : (print (allrad-decode-signal-remix-parameters :sn3d 3 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :basic))
; 3rd-maxre-mut210        : (print (allrad-decode-signal-remix-parameters :sn3d 3 (first (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) (second (sox-hoadecode-mut210-cartesian-positions-and-imaginary-idxs)) :sn3d :maxre))
; 3rd-maxre-octahedron    : (print (allrad-decode-signal-remix-parameters :sn3d 3 '((1 0 0) (0 1 0) (-1 0 0) (0 -1 0) (0 0 1) (0 0 -1)) '() :sn3d :maxre))


 sources/classes/t5200.txt

New file
Diff too large