externals/iem/iem_ambi/src Make.config.in, NONE, 1.1 Makefile, NONE, 1.1 ambi_decode.c, NONE, 1.1 ambi_decode2.c, NONE, 1.1 ambi_decode3.c, NONE, 1.1 ambi_decode_cube.c, NONE, 1.1 ambi_encode.c, NONE, 1.1 ambi_rot.c, NONE, 1.1 configure.ac, NONE, 1.1 iem_ambi.c, NONE, 1.1 iem_ambi.dsp, NONE, 1.1 iem_ambi.dsw, NONE, 1.1 iem_ambi.h, NONE, 1.1 iem_ambi_sources.h, NONE, 1.1 iem_bin_ambi_sources.c, NONE, 1.1 iemlib.h, NONE, 1.1 makefile_win, NONE, 1.1 makesource.sh, NONE, 1.1 - Pd-cvs

9 Mar 2006


      Update of /cvsroot/pure-data/externals/iem/iem_ambi/src
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv21232/iem/iem_ambi/src
Added Files:
    Make.config.in Makefile ambi_decode.c ambi_decode2.c 
    ambi_decode3.c ambi_decode_cube.c ambi_encode.c ambi_rot.c 
    configure.ac iem_ambi.c iem_ambi.dsp iem_ambi.dsw iem_ambi.h 
    iem_ambi_sources.h iem_bin_ambi_sources.c iemlib.h 
    makefile_win makesource.sh 
Log Message:
initial version JAN_06 with all sources
--- NEW FILE: iem_ambi.dsp ---
# Microsoft Developer Studio Project File - Name="iem_ambi" - Package Owner=<4>
# Microsoft Developer Studio Generated Build File, Format Version 6.00
# ** NICHT BEARBEITEN **
# TARGTYPE "Win32 (x86) External Target" 0x0106
CFG=iem_ambi - Win32 Debug
!MESSAGE Dies ist kein gültiges Makefile. Zum Erstellen dieses Projekts mit NMAKE
!MESSAGE verwenden Sie den Befehl "Makefile exportieren" und führen Sie den Befehl
!MESSAGE 
!MESSAGE NMAKE /f "iem_ambi.mak".
!MESSAGE 
!MESSAGE Sie können beim Ausführen von NMAKE eine Konfiguration angeben
!MESSAGE durch Definieren des Makros CFG in der Befehlszeile. Zum Beispiel:
!MESSAGE 
!MESSAGE NMAKE /f "iem_ambi.mak" CFG="iem_ambi - Win32 Debug"
!MESSAGE 
!MESSAGE Für die Konfiguration stehen zur Auswahl:
!MESSAGE 
!MESSAGE "iem_ambi - Win32 Release" (basierend auf  "Win32 (x86) External Target")
!MESSAGE "iem_ambi - Win32 Debug" (basierend auf  "Win32 (x86) External Target")
!MESSAGE
# Begin Project
# PROP AllowPerConfigDependencies 0
# PROP Scc_ProjName ""
# PROP Scc_LocalPath ""
!IF  "$(CFG)" == "iem_ambi - Win32 Release"
# PROP BASE Use_Debug_Libraries 0
# PROP BASE Output_Dir "Release"
# PROP BASE Intermediate_Dir "Release"
# PROP BASE Cmd_Line "NMAKE /f makefile_win"
# PROP BASE Rebuild_Opt "/a"
# PROP BASE Target_File "makefile_win.exe"
# PROP BASE Bsc_Name "makefile_win.bsc"
# PROP BASE Target_Dir ""
# PROP Use_Debug_Libraries 0
# PROP Output_Dir "Release"
# PROP Intermediate_Dir "Release"
# PROP Cmd_Line "NMAKE /f makefile_win"
# PROP Rebuild_Opt "/a"
# PROP Target_File "iem_ambi.exe"
# PROP Bsc_Name "iem_ambi.bsc"
# PROP Target_Dir ""
!ELSEIF  "$(CFG)" == "iem_ambi - Win32 Debug"
# PROP BASE Use_Debug_Libraries 1
# PROP BASE Output_Dir "Debug"
# PROP BASE Intermediate_Dir "Debug"
# PROP BASE Cmd_Line "NMAKE /f makefile_win"
# PROP BASE Rebuild_Opt "/a"
# PROP BASE Target_File "makefile_win.exe"
# PROP BASE Bsc_Name "makefile_win.bsc"
# PROP BASE Target_Dir ""
# PROP Use_Debug_Libraries 1
# PROP Output_Dir "Debug"
# PROP Intermediate_Dir "Debug"
# PROP Cmd_Line "NMAKE /f makefile_win"
# PROP Rebuild_Opt "/a"
# PROP Target_File "iem_ambi.exe"
# PROP Bsc_Name "iem_ambi.bsc"
# PROP Target_Dir ""
!ENDIF
# Begin Target
# Name "iem_ambi - Win32 Release"
# Name "iem_ambi - Win32 Debug"
!IF  "$(CFG)" == "iem_ambi - Win32 Release"
!ELSEIF  "$(CFG)" == "iem_ambi - Win32 Debug"
!ENDIF
# Begin Source File
SOURCE=.\makefile_win
# End Source File
# End Target
# End Project
--- NEW FILE: iem_ambi.c ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif
#include "m_pd.h"
#include "iemlib.h"
static t_class *iem_ambi_class;
static void *iem_ambi_new(void)
{
    t_object *x = (t_object *)pd_new(iem_ambi_class);
return (x);
}
void ambi_encode_setup(void);
void ambi_decode_setup(void);
void ambi_decode2_setup(void);
void ambi_decode3_setup(void);
void ambi_decode_cube_setup(void);
void ambi_rot_setup(void);
/* ------------------------ setup routine ------------------------- */
void iem_ambi_setup(void)
{
    ambi_encode_setup();
    ambi_decode_setup();
    ambi_decode2_setup();
    ambi_decode3_setup();
    ambi_decode_cube_setup();
    ambi_rot_setup();
post("iem_ambi (R-1.16) library loaded!   (c) Thomas Musil 05.2005");
    post("   musil%ciem.at iem KUG Graz Austria", '@');
}
--- NEW FILE: iem_ambi.h ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifndef __IEMAMBI_H__
#define __IEMAMBI_H__
#define AMBI_LS_REAL 0
#define AMBI_LS_IND 0
#define AMBI_LS_MRG 1
#define AMBI_LS_MIR 2
#define AMBI_LS_PHT 3
#endif
--- NEW FILE: ambi_encode.c ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif
#include "m_pd.h"
#include "iemlib.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
/* -------------------------- ambi_encode ------------------------------ */
typedef struct _ambi_encode
{
    t_object	x_obj;
    t_atom		*x_at;
    int				x_size;
    int				x_size2d;
    int				x_size3d;
    float			x_sqrt3;
    float			x_sqrt10_4;
    float			x_sqrt15;
    float			x_sqrt6_4;
    float			x_sqrt35_2;
    float			x_sqrt70_4;
    float			x_sqrt5_2;
    float			x_sqrt126_16;
    float			x_sqrt315_2;
    float			x_sqrt105_2;
    float			x_pi_over_180;
    float			*x_ambi_order_weight;
    int				x_colrow;
    int				x_n_order;
} t_ambi_encode;
static t_class *ambi_encode_class;
static void ambi_encode_ambi_weight(t_ambi_encode *x, t_symbol *s, int argc, t_atom *argv)
{
    if(argc > x->x_n_order)
    {
    	int i, n=x->x_n_order;
for(i=0; i<=n; i++)
    	{
    		x->x_ambi_order_weight[i] = atom_getfloat(argv++);
    	}
    }
    else
    	post("ambi_encode-ERROR: ambi_weight needs %d float weights", x->x_n_order+1);
}
static void ambi_encode_do_2d(t_ambi_encode *x, t_floatarg phi)
{
    float c, s, cc, ss, c2, s2, c3, s3, s4, c4, s5, c5, s6, c6;
    float *awght = x->x_ambi_order_weight;
    t_atom *at=x->x_at;
    
    phi *= x->x_pi_over_180;
    c = cos(phi);
    s = sin(phi);
    cc = c*c;
    ss = s*s;
SETFLOAT(at, (float)x->x_colrow);
    at++;
SETFLOAT(at, awght[0]);
    at++;
SETFLOAT(at, c*awght[1]);
    at++;
    SETFLOAT(at, s*awght[1]);
    at++;
if(x->x_n_order >= 2)
    {
    	c2 = cc - ss;
    	s2 = 2.0f*s*c;
    	SETFLOAT(at, c2*awght[2]);
    	at++;
    	SETFLOAT(at, s2*awght[2]);
    	at++;
if(x->x_n_order >= 3)
    	{
    		c3 = c*(4.0f*cc - 3.0f);
    		s3 = s*(3.0f - 4.0f*ss);
    		SETFLOAT(at, c3*awght[3]);
    		at++;
    		SETFLOAT(at, s3*awght[3]);
    		at++;
if(x->x_n_order >= 4)
    		{
    			c4 = 1.0f + 8.0f*cc*(cc - 1.0f);
    			s4 = 2.0f*s2*c2;
    			SETFLOAT(at, c4*awght[4]);
    			at++;
    			SETFLOAT(at, s4*awght[4]);
    			at++;
if(x->x_n_order >= 5)
    			{
    				c5 = c*(1.0f + 4.0f*ss*(ss - 3.0f*cc));
    				s5 = s*(1.0f + 4.0f*cc*(cc - 3.0f*ss));
    				SETFLOAT(at, c5*awght[5]);
    				at++;
    				SETFLOAT(at, s5*awght[5]);
    				at++;
if(x->x_n_order >= 6)
    				{
    					c6 = c3*c3 - s3*s3;
    					s6 = 2.0f*s3*c3;
    					SETFLOAT(at, c6*awght[6]);
    					at++;
    					SETFLOAT(at, s6*awght[6]);
    					at++;
if(x->x_n_order >= 7)
    					{
    						SETFLOAT(at, cos(7.0f*phi)*awght[7]);
    						at++;
    						SETFLOAT(at, sin(7.0f*phi)*awght[7]);
    						at++;
if(x->x_n_order >= 8)
    						{
    							SETFLOAT(at, (c4*c4 - s4*s4)*awght[8]);
    							at++;
    							SETFLOAT(at, 2.0f*s4*c4*awght[8]);
    							at++;
if(x->x_n_order >= 9)
    							{
    								SETFLOAT(at, cos(9.0f*phi)*awght[9]);
    								at++;
    								SETFLOAT(at, sin(9.0f*phi)*awght[9]);
    								at++;
if(x->x_n_order >= 10)
    								{
    									SETFLOAT(at, (c5*c5 - s5*s5)*awght[10]);
    									at++;
    									SETFLOAT(at, 2.0f*s5*c5*awght[10]);
    									at++;
if(x->x_n_order >= 11)
    									{
    										SETFLOAT(at, cos(11.0f*phi)*awght[11]);
    										at++;
    										SETFLOAT(at, sin(11.0f*phi)*awght[11]);
    										at++;
if(x->x_n_order >= 12)
    										{
    											SETFLOAT(at, (c6*c6 - s6*s6)*awght[12]);
    											at++;
    											SETFLOAT(at, 2.0f*s6*c6*awght[12]);
    										}
if(x->x_n_order >= 13)
    											post("ambi_encode-ERROR: do not support Ambisonic-Order greater than 12 in 2d !!!");
}
    								}
    							}
    						}
    					}
    				}
    			}
    		}
    	}
    }
}
static void ambi_encode_do_3d(t_ambi_encode *x, t_symbol *s, int argc, t_atom *argv)
{
    float delta, phi;
    float cd, x1, y, z, x2, y2, z2, x2my2, x2m3y2, p3x2my2, xy, xz, yz, m1p5z2, m1p7z2, m3p7z2;
    float *awght = x->x_ambi_order_weight;
    t_atom *at=x->x_at;
delta = atom_getfloat(argv++)*x->x_pi_over_180;
    phi = atom_getfloat(argv)*x->x_pi_over_180;
cd = cos(delta);
    x1 = cd * cos(phi);
    y = cd * sin(phi);
    z = sin(delta);
xy = x1*y;
    xz = x1*z;
    yz = y*z;
    x2 = x1*x1;
    y2 = y*y;
    z2 = z*z;
x2my2 = x2 - y2;
    x2m3y2 = x2my2 - 2.0f*y2;
    p3x2my2 = 2.0f*x2 + x2my2;
    m1p5z2 = 5.0f*z2 - 1.0f;
    m1p7z2 = 2.0f*z2 + m1p5z2;
    m3p7z2 = m1p7z2 - 2.0f;
SETFLOAT(at, (float)x->x_colrow);
    at++;
SETFLOAT(at, awght[0]);
    at++;
SETFLOAT(at, x1*awght[1]);
    at++;
    SETFLOAT(at, y*awght[1]);
    at++;
    SETFLOAT(at, z*awght[1]);
    at++;
if(x->x_n_order >= 2)
    {
    	SETFLOAT(at, 0.5f*x->x_sqrt3*x2my2*awght[2]);
    	at++;
    	SETFLOAT(at, x->x_sqrt3*xy*awght[2]);
    	at++;
    	SETFLOAT(at, x->x_sqrt3*xz*awght[2]);
    	at++;
    	SETFLOAT(at, x->x_sqrt3*yz*awght[2]);
    	at++;
    	SETFLOAT(at, 0.5f*(3.0f*z2 - 1.0f)*awght[2]);
    	at++;
if(x->x_n_order >= 3)
    	{
    		SETFLOAT(at, x->x_sqrt10_4*x1*x2m3y2*awght[3]);
    		at++;
    		SETFLOAT(at, x->x_sqrt10_4*y*p3x2my2*awght[3]);
    		at++;
    		SETFLOAT(at, 0.5f*x->x_sqrt15*z*x2my2*awght[3]);
    		at++;
    		SETFLOAT(at, x->x_sqrt15*xy*z*awght[3]);
    		at++;
    		SETFLOAT(at, x->x_sqrt6_4*x1*m1p5z2*awght[3]);
    		at++;
    		SETFLOAT(at, x->x_sqrt6_4*y*m1p5z2*awght[3]);
    		at++;
    		SETFLOAT(at, 0.5f*z*(m1p5z2 - 2.0f)*awght[3]);
    		at++;
if(x->x_n_order >= 4)
    		{
    			SETFLOAT(at, 0.25f*x->x_sqrt35_2*(x2my2*x2my2 - 4.0f*x2*y2)*awght[4]);
    			at++;
    			SETFLOAT(at, x->x_sqrt35_2*xy*x2my2*awght[4]);
    			at++;
    			SETFLOAT(at, x->x_sqrt70_4*xz*x2m3y2*awght[4]);
    			at++;
    			SETFLOAT(at, x->x_sqrt70_4*yz*p3x2my2*awght[4]);
    			at++;
    			SETFLOAT(at, 0.5f*x->x_sqrt5_2*x2my2*m1p7z2*awght[4]);
    			at++;
    			SETFLOAT(at, x->x_sqrt5_2*xy*m1p7z2*awght[4]);
    			at++;
    			SETFLOAT(at, x->x_sqrt10_4*xz*m3p7z2*awght[4]);
    			at++;
    			SETFLOAT(at, x->x_sqrt10_4*yz*m3p7z2*awght[4]);
    			at++;
    			SETFLOAT(at, 0.125f*(5.0f*(z2 - 1.0f)*(m1p7z2 + 2.0f) + 8.0f)*awght[4]);
    			at++;
if(x->x_n_order >= 5)
    			{
    				SETFLOAT(at, x->x_sqrt126_16*x1*(x2*(x2 - 10.0f*y2) + 5.0f*y2*y2)*awght[5]);
    				at++;
    				SETFLOAT(at, x->x_sqrt126_16*y*(y2*(y2 - 10.0f*x2) + 5.0f*x2*x2)*awght[5]);
    				at++;
    				SETFLOAT(at, 0.25f*x->x_sqrt315_2*z*(y2*(y2 - 6.0f*x2) + x2*x2)*awght[5]);
    				at++;
    				SETFLOAT(at, x->x_sqrt315_2*xy*z*x2my2*awght[5]);
    				at++;
    				SETFLOAT(at, 0.25f*x->x_sqrt70_4*x1*(9.0f*z2 - 1.0f)*x2m3y2*awght[5]);
    				at++;
    				SETFLOAT(at, 0.25f*x->x_sqrt70_4*y*(9.0f*z2 - 1.0f)*p3x2my2*awght[5]);
    				at++;
    				SETFLOAT(at, 0.5f*x->x_sqrt105_2*x2my2*z*(3.0f*z2 - 1.0f)*awght[5]);
    				at++;
    				SETFLOAT(at, x->x_sqrt105_2*xy*z*(3.0f*z2 - 1.0f)*awght[5]);
    				at++;
    				SETFLOAT(at, 0.125f*x->x_sqrt15*x1*(z2*(21.0f*z2 - 14.0f) + 1.0f)*awght[5]);
    				at++;
    				SETFLOAT(at, 0.125f*x->x_sqrt15*y*(z2*(21.0f*z2 - 14.0f) + 1.0f)*awght[5]);
    				at++;
    				SETFLOAT(at, 0.125f*z*(z2*(63.0f*z2 - 70.0f) + 15.0f)*awght[5]);
    			}
if(x->x_n_order > 5)
    				post("ambi_encode-ERROR: do not support Ambisonic-Order greater than 5 in 3d !!!");
    		}
    	}
    }
}
static void ambi_encode_float(t_ambi_encode *x, t_floatarg phi)
{
    x->x_colrow = -1;
    ambi_encode_do_2d(x, phi);
    outlet_list(x->x_obj.ob_outlet, &s_list, x->x_size2d, x->x_at+1);
}
static void ambi_encode_list(t_ambi_encode *x, t_symbol *s, int argc, t_atom *argv)
{
    if(argc <= 0)
    {
    	post("ambi_encode ERROR: list-input needs 2 angles: delta [rad] and phi [rad]");
    	return;
    }
    else if(argc == 1)
    {
    	ambi_encode_float(x, atom_getfloat(argv));
    }
    else
    {
    	x->x_colrow = -1;
    	ambi_encode_do_3d(x, &s_list, 2, argv);
    	outlet_list(x->x_obj.ob_outlet, &s_list, x->x_size3d, x->x_at+1);
    }
}
static void ambi_encode_row(t_ambi_encode *x, t_symbol *s, int argc, t_atom *argv)
{
    if(argc == 2)
    {
    	x->x_colrow = (int)atom_getint(argv++);
    	ambi_encode_do_2d(x, atom_getfloat(argv));
    	outlet_anything(x->x_obj.ob_outlet, s, x->x_size2d+1, x->x_at);
    }
    else if(argc >= 3)
    {
    	x->x_colrow = (int)atom_getint(argv++);
    	ambi_encode_do_3d(x, &s_list, 2, argv);
    	outlet_anything(x->x_obj.ob_outlet, s, x->x_size3d+1, x->x_at);
    }
    else
    {
    	post("ambi_encode-ERROR: row needs <float> row-index + <float> angle ( + <float> angle)");
    }
}
static void ambi_encode_col(t_ambi_encode *x, t_symbol *s, int argc, t_atom *argv)
{
    if(argc == 2)
    {
    	x->x_colrow = (int)atom_getint(argv++);
    	ambi_encode_do_2d(x, atom_getfloat(argv));
    	outlet_anything(x->x_obj.ob_outlet, s, x->x_size2d+1, x->x_at);
    }
    else if(argc >= 3)
    {
    	x->x_colrow = (int)atom_getint(argv++);
    	ambi_encode_do_3d(x, &s_list, 2, argv);
    	outlet_anything(x->x_obj.ob_outlet, s, x->x_size3d+1, x->x_at);
    }
    else
    {
    	post("ambi_encode-ERROR: col needs <float> col-index + <float> angle ( + <float> angle)");
    }
}
static void ambi_encode_free(t_ambi_encode *x)
{
    freebytes(x->x_ambi_order_weight, (x->x_n_order+1) * sizeof(float));
    freebytes(x->x_at, x->x_size * sizeof(t_atom));
}
static void *ambi_encode_new(t_floatarg forder)
{
    t_ambi_encode *x = (t_ambi_encode *)pd_new(ambi_encode_class);
    t_atom *at;
    int i=(int)forder;
if(i < 1)
    	i = 1;
    if(i > 12)
    	i = 12;
    x->x_n_order = i;
    x->x_size = 6*6 + 1;
    x->x_size2d = 2*i + 1;
    x->x_size3d = (i + 1)*(i + 1);
x->x_sqrt3			= (float)(sqrt(3.0));
    x->x_sqrt5_2		= (float)(sqrt(5.0) / 2.0);
    x->x_sqrt6_4		= (float)(sqrt(6.0) / 4.0);
    x->x_sqrt10_4		= (float)(sqrt(10.0) / 4.0);
    x->x_sqrt15			= (float)(sqrt(15.0));
    x->x_sqrt35_2		= (float)(sqrt(35.0) / 2.0);
    x->x_sqrt70_4		= (float)(sqrt(70.0) / 4.0);
    x->x_sqrt126_16	= (float)(sqrt(126.0) / 16.0);
    x->x_sqrt315_2	= (float)(sqrt(315.0) / 2.0);
    x->x_sqrt105_2	= (float)(sqrt(105.0) / 2.0);
    x->x_pi_over_180 = (float)(4.0 * atan(1.0)/180.0);
    x->x_colrow = 0;
    x->x_ambi_order_weight = (float *)getbytes((x->x_n_order+1) * sizeof(float));
    x->x_at = (t_atom *)getbytes(x->x_size * sizeof(t_atom));
    at=x->x_at;
    SETFLOAT(at, -1.0f);/*row index*/
    at++;
    SETFLOAT(at, 1.0f);/*W channel*/
for(i=0; i<=x->x_n_order; i++)
    		x->x_ambi_order_weight[i] = 1.0f;
outlet_new(&x->x_obj, &s_list);
    return (x);
}
void ambi_encode_setup(void)
{
    ambi_encode_class = class_new(gensym("ambi_encode"), (t_newmethod)ambi_encode_new, (t_method)ambi_encode_free,
    			   sizeof(t_ambi_encode), 0, A_DEFFLOAT, 0);
    class_addlist(ambi_encode_class, (t_method)ambi_encode_list);
    class_addfloat(ambi_encode_class, (t_method)ambi_encode_float);
    class_addmethod(ambi_encode_class, (t_method)ambi_encode_row, gensym("row"), A_GIMME, 0);
    class_addmethod(ambi_encode_class, (t_method)ambi_encode_col, gensym("col"), A_GIMME, 0);
    class_addmethod(ambi_encode_class, (t_method)ambi_encode_ambi_weight, gensym("ambi_weight"), A_GIMME, 0);
    class_sethelpsymbol(ambi_encode_class, gensym("iemhelp2/help-ambi_encode"));
}
--- NEW FILE: iem_ambi.dsw ---
(This appears to be a binary file; contents omitted.)
--- NEW FILE: iem_bin_ambi_sources.c ---
/* iem_bin_ambi-setup autogenerated setup-file
 * generated by "./makesource.sh"
 * !! DO NOT MANUALLY EDIT  !!
 */
#include "iem_bin_ambi_sources.h"
void iem_bin_ambi_sources_setup(void)
{
    bin_ambi_calc_HRTF_setup(); /* bin_ambi_calc_HRTF.c */
    bin_ambi_reduced_decode_setup(); /* bin_ambi_reduced_decode.c */
}
--- NEW FILE: makefile_win ---
all: iem_ambi.dll
VIS_CPP_PATH = "C:\Programme\Microsoft Visual Studio\Vc98"
PD_INST_PATH = "C:\Programme\pd"
PD_WIN_INCLUDE_PATH = /I. /I$(PD_INST_PATH)\src /I$(VIS_CPP_PATH)\include
PD_WIN_C_FLAGS = /nologo /W3 /WX /DMSW /DNT /DPD /DWIN32 /DWINDOWS /Ox -DPA_LITTLE_ENDIAN
PD_WIN_L_FLAGS = /nologo
PD_WIN_LIB = /NODEFAULTLIB:libc /NODEFAULTLIB:oldnames /NODEFAULTLIB:kernel /NODEFAULTLIB:uuid \
    $(VIS_CPP_PATH)\lib\libc.lib \
    $(VIS_CPP_PATH)\lib\oldnames.lib \
    $(VIS_CPP_PATH)\lib\kernel32.lib \
    $(VIS_CPP_PATH)\lib\wsock32.lib \
    $(VIS_CPP_PATH)\lib\winmm.lib \
    $(PD_INST_PATH)\bin\pd.lib
SRC =	ambi_decode.c \
    	ambi_decode2.c \
    	ambi_decode3.c \
    	ambi_decode_cube.c \
    	ambi_encode.c \
    	ambi_rot.c \
    	iem_ambi.c
OBJ = $(SRC:.c=.obj)
.c.obj:
    cl $(PD_WIN_C_FLAGS) $(PD_WIN_INCLUDE_PATH) /c $*.c
iem_ambi.dll: $(OBJ)
    link $(PD_WIN_L_FLAGS) /dll /export:iem_ambi_setup \
    /out:iem_ambi.dll $(OBJ) $(PD_WIN_LIB)
clean:
    del *.obj
--- NEW FILE: ambi_rot.c ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif
#include "m_pd.h"
#include "iemlib.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
[...1495 lines suppressed...]
    x->x_sqrt10_4		= (float)(sqrt(10.0) / 4.0);
    x->x_sqrt14_16	= (float)(sqrt(14.0) / 16.0);
    x->x_sqrt15_8		= (float)(sqrt(15.0) / 8.0);
    x->x_sqrt35_64	= (float)(sqrt(35.0) / 64.0);
    x->x_sqrt70_32	= (float)(sqrt(70.0) / 32.0);
x->x_pi_over_180	= (float)(4.0 * atan(1.0) / 180.0);
x->x_s_matrix = gensym("matrix");
    return (x);
}
void ambi_rot_setup(void)
{
    ambi_rot_class = class_new(gensym("ambi_rot"), (t_newmethod)ambi_rot_new, (t_method)ambi_rot_free,
    				 sizeof(t_ambi_rot), 0, A_DEFFLOAT, 0);
    class_addfloat(ambi_rot_class, (t_method)ambi_rot_float);
    	class_addlist(ambi_rot_class, (t_method)ambi_rot_list);
    class_sethelpsymbol(ambi_rot_class, gensym("iemhelp2/help-ambi_rot"));
}
--- NEW FILE: makesource.sh ---
#!/bin/sh
IEMAMBI_H=iem_ambi_sources.h
IEMAMBI_C=iem_ambi_sources.c
EGREP=egrep
SED=sed
LS=ls
#################################
## functions
function head_h() {
 echo "/* iem_ambi-setup autogenerated header-file"
 echo " * generated by "$0""
 echo " * !! DO NOT MANUALLY EDIT  !!"
 echo " */"
 echo
 echo "#ifndef IEMAMBI_SOURCES_H__"
 echo "#define IEMAMBI_SOURCES_H__"
}
function foot_h() {
 echo "#endif /* IEMAMBI_SOURCES_H__ */"
 echo ""
}
function head_c() {
 echo "/* iem_ambi-setup autogenerated setup-file"
 echo " * generated by "$0""
 echo " * !! DO NOT MANUALLY EDIT  !!"
 echo " */"
 echo 
 echo "#include "$IEMAMBI_H""
 echo
 echo "void iem_ambi_sources_setup(void)"
 echo "{"
}
function foot_c() {
 echo "}"
 echo
}
##################################
## body
head_h > $IEMAMBI_H
head_c > $IEMAMBI_C
for i in $(${LS} *.c | ${EGREP} -v "iem_bin_ambi.*.c")
do
## each c-file in iem_ambi needs to have an <file>_setup()-function
## that calls all needed setup-functions
## any non-alpha-numeric-character is replaced by "_"
## e.g. "bla~.c" -> "bla__setup()"
  SETUPNAME=$(echo ${i%.c} | ${SED} -e 's/[^[:alnum:]]/_/g')_setup
  echo "void ${SETUPNAME}(void); /* $i */" >> $IEMAMBI_H
  echo "	${SETUPNAME}(); /* $i */" >> $IEMAMBI_C
done
foot_h >> $IEMAMBI_H
foot_c >> $IEMAMBI_C
--- NEW FILE: Makefile ---
default: all
.PHONEY: default all everything dist \
    clean realclean distclean \
    install install-bin install-doc install-abs
SOURCES=$(sort $(filter %.c, $(wildcard *.c)))
ifeq (,$(findstring clean, $(MAKECMDGOALS)))
Make.config: Make.config.in configure
    ./configure
endif
iem_bin_ambi_sources.c iem_bin_ambi_sources.h:
    ./makesource.sh
configure: configure.ac
    autoconf
## uaehh, here comes some magic
##  1st we don't want depend and config-makefiles to be included on "clean"-targets
ifeq (,$(findstring clean, $(MAKECMDGOALS)))
-include $(SOURCES:.c=.d)
endif
-include Make.config
##  2nd only generate depend-files when we have Make.config included
##  and thus MAKEDEP_FLAGS defined
ifdef MAKEDEP_FLAGS
## dependencies: as proposed by the GNU-make documentation
## see http://www.gnu.org/software/make/manual/html_node/make_47.html#SEC51
%.d: %.c
    @set -e; rm -f $@; \
     $(CPP) $(MAKEDEP_FLAGS) $(Z_CFLAGS) $< > $@.$$$$; \
     sed 's,($*).o[ :]*,\1.o $@ : ,g' < $@.$$$$ > $@; \
     rm -f $@.$$$$
endif
.SUFFIXES: .$(EXT)
TARGETS = $(SOURCES:.c=.o)
all: $(LIBNAME)
    cp $(LIBNAME).$(EXT) ..
$(LIBNAME): $(TARGETS) iem_bin_ambi_sources.c iem_bin_ambi_sources.h
    $(LD) $(LFLAGS) -o $(LIBNAME).$(EXT) *.o $(LIBS)
    $(STRIP) $(STRIPFLAGS) $(LIBNAME).$(EXT)
$(TARGETS): %.o : %.c
    $(CC) $(Z_CFLAGS) -c -o $@ $*.c
clean:
    -rm -f *.$(EXT) *.o
realclean: clean
    -rm -f *~ _* config.*
    -rm -f *.d *.d.*
distclean: realclean
    -rm -f Make.config ../*.$(EXT)
    -rm -f *.exp  *.lib  *.ncb  \
    	*.opt  *.plg
    -rm -rf autom4te.cache/
install: install-bin install-doc install-abs
install-bin:
    -install -d $(INSTALL_BIN)
    -install -m 644 $(LIBNAME).$(EXT) $(INSTALL_BIN)
install-doc:
    -install -d $(INSTALL_DOC)
    -install -m 644 ../examples/*.pd $(INSTALL_DOC)
install-abs:
    -install -d $(INSTALL_BIN)
    -install -m 644 ../abs/*.pd $(INSTALL_BIN)
dist: all realclean
    (cd ../..;tar czvf $(TARNAME) $(LIBNAME))
everything: clean all install distclean
--- NEW FILE: configure.ac ---
dnl  Process this file with autoconf to produce a configure script.
AC_INIT(iem_bin_ambi.c)
dnl Checks for programs.
AC_PROG_CC
AC_SUBST(STK)
AC_SUBST(DFLAGS)
AC_SUBST(LFLAGS)
AC_SUBST(EXT)
AC_SUBST(LD)
AC_SUBST(STRIP)
AC_SUBST(STRIPFLAGS)
AC_SUBST(IEMBINAMBI_VERSION)
AC_SUBST(REFERENCEPATH)
AC_SUBST(PDLIBDIR)
AC_SUBST(INCLUDES)
AC_ARG_WITH(pdversion, [  --with-pdversion=<ver>  enforce a certain pd-version (e.g. 0.37)])
AC_ARG_WITH(version,   [  --with-version=<ver>    enforce a certain iem_bin_ambi-version (e.g. 0.1)])
AC_ARG_WITH(extension, [  --with-extension=<ext>  enforce a certain extension for the dynamic library (e.g. dll)])
AC_ARG_WITH(pdpath,    [  --with-pd=</path/to/pd> where to look for pd-headers and and -libs])
AC_ARG_ENABLE(PIC,     [  --disable-PIC           disable compilation with PIC-flag])
dnl Checks for libraries.
dnl Replace `main' with a function in -lc:
AC_CHECK_LIB(c, main)
AC_CHECK_LIB(crtdll, fclose)
dnl Replace `main' with a function in -lm:
AC_CHECK_LIB(m, main)
dnl Replace `main' with a function in -lpthread:
dnl AC_CHECK_LIB(pthread, main)
dnl Replace `main' with a function in -lstk:
dnl AC_CHECK_LIB(stk, main, STK=yes)
if test "x$with_pd" != "x"; then
 if test -d "${with_pd}/src"; then
  INCLUDES="-I${with_pd}/src ${INCLUDES}"
 fi
 if test -d "${with_pd}/bin"; then
  LIBS="-L${with_pd}/bin ${LIBS}"
 fi
fi
if test "x$includedir" != "x"; then
 for id in $includedir
 do
   if test -d $id; then INCLUDES="-I$id $INCLUDES"; fi
 done
fi
if test "x$libdir" != "x"; then
 for id in $libdir
 do
   if test -d $id; then LIBS="-L$id $LIBS"; fi
 done
fi
AC_CHECK_LIB(pd, nullfn)
dnl Checks for header files.
AC_HEADER_STDC
AC_CHECK_HEADERS(stdlib.h stdio.h string.h math.h time.h sys/time.h)
dnl Checks for typedefs, structures, and compiler characteristics.
AC_HEADER_TIME
dnl Checks for library functions.
AC_FUNC_MMAP
AC_CHECK_FUNCS(select socket strerror)
### make-depend flags
if test "x$ac_cv_c_compiler_gnu" = "xyes"; then
 AC_SUBST(MAKEDEP_FLAGS, "-MM")
else
 AC_SUBST(MAKEDEP_FLAGS, "-M")
fi
dnl check for "-mms-bitfields" cflag
dnl why is there no generic compiler-check for a given flag ?
dnl it would make things so easy: AC_CHECK_FLAG([-mms-bitfields],,)
AC_MSG_CHECKING("ms-bitfields")
cat > conftest.c << EOF
int main(){
  return 0;
}
EOF
if ${CC} ${INCLUDES} ${DFLAGS} -o conftest.o conftest.c ${CFLAGS} -mms-bitfields > /dev/null 2>&1
then
  echo "yes"
  CFLAGS="${CFLAGS} -mms-bitfields"
else
  echo "no"
fi
dnl isn't there a better way to check for good linker/stripper ?
dnl if we don't have $LD set, we set it to $(CC)
dnl LD=${LD:=$CC}
if test "x$LD" = "x"
then
 if test "x$host" != "x"
 then
  LD=${host}-ld
  if $(which ${LD} > /dev/null)
  then
    :
  else
    LD=""
  fi
 fi
fi
LD=${LD:=$CC}
dnl if we don't have $STRIP set, we set it to ${host}-strip or strip
AC_CHECK_TOOL([STRIP], [strip], [true])
AC_MSG_CHECKING([if strip is GNU strip])
if $STRIP -V 2>&1 | grep GNU > /dev/null
then
    AC_SUBST(STRIPFLAGS, "--strip-unneeded")
    AC_MSG_RESULT([yes])
else
    AC_SUBST(STRIPFLAGS,"-x")
    AC_MSG_RESULT([no])
fi
DFLAGS=""
if test "x$enable_PIC" != "xno"; then
AC_MSG_CHECKING("PIC")
cat > conftest.c << EOF
int main(){
  return 0;
}
EOF
if ${CC} ${INCLUDES} ${DFLAGS} -o conftest.o conftest.c ${CFLAGS} -fPIC > /dev/null 2>&1
then
  echo "yes"
  CFLAGS="${CFLAGS} -fPIC"
else
  echo "no"
fi
fi
dnl
dnl  OK, checks for machines are here now
dnl
if test `uname -s` = Linux; 
then
  LFLAGS="-export_dynamic  -shared"
  CFLAGS="$CFLAGS"
  EXT=pd_linux	
fi
dnl This should use '-bundle_loader /path/to/pd/bin/pd' instead of'-undefined suppress'
dnl then strip might do something
if test `uname -s` = Darwin; 
then
  LD=cc
  LFLAGS="-bundle -undefined suppress -flat_namespace"
  EXT=pd_darwin
fi
if test `uname | sed -e 's/^MINGW.*/NT/'` = NT; 
then
  LD=gcc
  INCLUDES="-I@prefix@/src"
  DFLAGS="-D__WIN32__"
  LFLAGS="-shared @prefix@/bin/pd.dll"
  EXT=dll
else
  PDLIBDIR="/lib/pd"
fi
if test `uname -s` = IRIX64; 
then
  LFLAGS="-n32 -DUNIX -DIRIX -DN32 -woff 1080,1064,1185 \
    -OPT:roundoff=3 -OPT:IEEE_arithmetic=3 -OPT:cray_ivdep=true \
        -shared -rdata_shared"
  EXT=pd_irix6
  dnl DFLAGS="-DUNIX -DIRIX6"
fi
if test `uname -s` = IRIX32; 
then
  LFLAGS="-o32 -DUNIX -DIRIX -O2 
         -shared -rdata_shared"
  EXT=pd_irix5
  dnl DFLAGS="-DUNIX -DIRIX5"
fi
if test "x$with_extension" != "x"
then
 EXT=$with_extension
fi
dnl Checks for pd-version, to set the correct help-path
AC_MSG_CHECKING("pd>=0.37")
if test "$with_pdversion" != ""
then
echo -n "($with_pdversion)... "
  PD_VERSION="$with_pdversion"
else
if test "x$cross_compiling" = "xno"
then
cat > conftest.c << EOF
#include <stdio.h>
#include "m_pd.h"
int main(){
  printf("%d.%d\n", PD_MAJOR_VERSION, PD_MINOR_VERSION);
  return 0;
}
EOF
if $CC $INCLUDES -o conftest.o conftest.c > /dev/null 2>&1
 then
  PD_VERSION=`./conftest.o`
 else
  PD_VERSION=""
 fi
  echo -n $PD_VERSION
else
dnl we are cross-compiling...
 echo -n "(X)..."
 PD_VERSION="0.38"
fi
fi
let PD_MAJORVERSION=`echo $PD_VERSION | cut -d"." -f1`+0
let PD_MINORVERSION=`echo $PD_VERSION | cut -d"." -f2`+0
if test "$PD_MAJORVERSION" -gt 0 || test "$PD_MINORVERSION" -ge 37
then
  REFERENCEPATH=extra/help-
  echo " yes"
else
  REFERENCEPATH=doc/5.reference/
  echo " no"
fi
dnl check for iem_bin_ambi-version (but why...)
AC_MSG_CHECKING("iem_bin_ambi-version")
if test "$with_version" != ""
then
  echo -n "($with_version)... "
  IEMBINAMBI_VERSION="$with_version"
else
if test "x$cross_compiling" = "xno"
then
cat > conftest.c << EOF
#include <stdio.h>
#include "iem_bin_ambi.h"
int main(){
  printf("%s\n", VERSION);
  return 0;
}
EOF
if $CC $INCLUDES -o conftest.o conftest.c > /dev/null 2>&1
then
  IEMBINAMBI_VERSION=`./conftest.o`
  echo "$IEMBINAMBI_VERSION"
else
  IEMBINAMBI_VERSION=""
  echo "(unknown)"
fi
else
  IEMBINAMBI_VERSION="X"
  echo "(X)"
fi
fi
AC_OUTPUT(Make.config)
rm -f conftest.*
--- NEW FILE: ambi_decode2.c ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif
#include "m_pd.h"
#include "iemlib.h"
#include "iem_ambi.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
/* -------------------------- ambi_decode2 ------------------------------ */
/*
 ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren **
 ** Inputs: ls + Liste von 3 floats: Index [1 .. 25] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] **
 ** Inputs: calc_inv **
 ** Inputs: load_HRIR + float index1..25 **
 ** Outputs: List of 2 symbols: left-HRIR-File-name + HRIR-table-name **
 ** Inputs: calc_reduced **
 ** "output" ...  writes the HRTF into tables **
 **  **
 **  **
 ** setzt voraus , dass die HRIR-tabele-names von LS1_L_HRIR .. LS25_L_HRIR heissen und existieren **
 ** setzt voraus , dass die HRTF-tabele-names von LS1_HRTF_re .. LS25_HRTF_re heissen und existieren **
 ** setzt voraus , dass die HRTF-tabele-names von LS1_HRTF_im .. LS25_HRTF_im heissen und existieren **
 */
typedef struct _ambi_decode2
{
    t_object	x_obj;
    t_atom		*x_at;
    double		*x_inv_work1;
    double		*x_inv_work2;
    double		*x_inv_buf2;
    double		*x_transp;
    double		*x_ls_encode;
    double		*x_prod;
    double		*x_ambi_channel_weight;
    double		x_mirror_weight;
    double		x_sing_range;
    int				x_n_ambi;
    int				x_n_order;
    int				x_n_ls;
    int				x_n_ph_ls;
    int				x_n_mir_ls;
    int				x_n_dim;
    t_symbol	*x_s_matrix;
    double		x_sqrt3;
    double		x_sqrt10_4;
    double		x_sqrt15_2;
    double		x_sqrt6_4;
    double		x_sqrt35_8;
    double		x_sqrt70_4;
    double		x_sqrt5_2;
    double		x_sqrt126_16;
    double		x_sqrt315_8;
    double		x_sqrt105_4;
    double		x_pi_over_180;
} t_ambi_decode2;
static t_class *ambi_decode2_class;
static void ambi_decode2_copy_row2buf(t_ambi_decode2 *x, int row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*db++ = *dw++;
}
static void ambi_decode2_copy_buf2row(t_ambi_decode2 *x, int row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*dw++ = *db++;
}
static void ambi_decode2_copy_row2row(t_ambi_decode2 *x, int src_row, int dst_row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw_src=x->x_inv_work2;
    double *dw_dst=x->x_inv_work2;
dw_src += src_row*n_ambi2;
    dw_dst += dst_row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*dw_dst++ = *dw_src++;
}
static void ambi_decode2_xch_rows(t_ambi_decode2 *x, int row1, int row2)
{
    ambi_decode2_copy_row2buf(x, row1);
    ambi_decode2_copy_row2row(x, row2, row1);
    ambi_decode2_copy_buf2row(x, row2);
}
static void ambi_decode2_mul_row(t_ambi_decode2 *x, int row, double mul)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    {
    	(*dw) *= mul;
    	dw++;
    }
}
static void ambi_decode2_mul_buf_and_add2row(t_ambi_decode2 *x, int row, double mul)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    {
    	*dw += (*db)*mul;
    	dw++;
    	db++;
    }
}
static int ambi_decode2_eval_which_element_of_col_not_zero(t_ambi_decode2 *x, int col, int start_row)
{
    int n_ambi = x->x_n_ambi;
    int n_ambi2 = 2*n_ambi;
    int i, j;
    double *dw=x->x_inv_work2;
    double singrange=x->x_sing_range;
    int ret=-1;
dw += start_row*n_ambi2 + col;
    j = 0;
    for(i=start_row; i<n_ambi; i++)
    {
    	if((*dw > singrange) || (*dw < -singrange))
    	{
    		ret = i;
    		i = n_ambi+1;
    	}
    	dw += n_ambi2;
    }
    return(ret);
}
static void ambi_decode2_mul1(t_ambi_decode2 *x)
{
    double *vec1, *beg1=x->x_ls_encode;
    double *vec2, *beg2=x->x_ls_encode;
    double *inv=x->x_inv_work1;
    double sum;
    int n_ls=x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls;
    int n_ambi=x->x_n_ambi;
    int i, j, k;
for(k=0; k<n_ambi; k++)
    {
    	beg2=x->x_ls_encode;
    	for(j=0; j<n_ambi; j++)
    	{
    		sum = 0.0;
    		vec1 = beg1;
    		vec2 = beg2;
    		for(i=0; i<n_ls; i++)
    		{
    			sum += *vec1++ * *vec2++;
    		}
    		beg2 += n_ls;
    		*inv++ = sum;
    	}
    	beg1 += n_ls;
    }
}
static void ambi_decode2_mul2(t_ambi_decode2 *x)
{
    int n_ls=x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls;
    int n_ambi=x->x_n_ambi;
    int n_ambi2=2*n_ambi;
    int i, j, k;
    double *vec1, *beg1=x->x_transp;
    double *vec2, *beg2=x->x_inv_work2+n_ambi;
    double *vec3=x->x_prod;
    double *acw_vec=x->x_ambi_channel_weight;
    double sum;
for(k=0; k<n_ls; k++)
    {
    	beg2=x->x_inv_work2+n_ambi;
    	for(j=0; j<n_ambi; j++)
    	{
    		sum = 0.0;
    		vec1 = beg1;
    		vec2 = beg2;
    		for(i=0; i<n_ambi; i++)
    		{
    			sum += *vec1++ * *vec2;
    			vec2 += n_ambi2;
    		}
    		beg2++;
    		*vec3++ = sum * acw_vec[j];
    	}
    	beg1 += n_ambi;
    }
}
static void ambi_decode2_transp_back(t_ambi_decode2 *x)
{
    double *vec, *transp=x->x_transp;
    double *straight=x->x_ls_encode;
    int n_ls=x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls;
    int n_ambi=x->x_n_ambi;
    int i, j;
for(j=0; j<n_ambi; j++)
    {
    	vec = transp;
    	for(i=0; i<n_ls; i++)
    	{
    		*straight++ = *vec;
    		vec += n_ambi;
    	}
    	transp++;
    }
}
static void ambi_decode2_inverse(t_ambi_decode2 *x)
{
    int n_ambi = x->x_n_ambi;
    int n_ambi2 = 2*n_ambi;
    int i, j, nz;
    int r,c;
    double *src=x->x_inv_work1;
    double *db=x->x_inv_work2;
    double rcp, *dv;
dv = db;
    for(i=0; i<n_ambi; i++) /* init */
    {
    	for(j=0; j<n_ambi; j++)
    	{
    		*dv++ = *src++;
    	}
    	for(j=0; j<n_ambi; j++)
    	{
    		if(j == i)
    			*dv++ = 1.0;
    		else
    			*dv++ = 0.0;
    	}
    }
/* make 1 in main-diagonale, and 0 below */
    for(i=0; i<n_ambi; i++)
    {
    	nz = ambi_decode2_eval_which_element_of_col_not_zero(x, i, i);
    	if(nz < 0)
    	{
    		post("ambi_decode2 ERROR: matrix not regular !!!!");
    		return;
    	}
    	else
    	{
    		if(nz != i)
    			ambi_decode2_xch_rows(x, i, nz);
    		dv = db + i*n_ambi2 + i;
    		rcp = 1.0 /(*dv);
    		ambi_decode2_mul_row(x, i, rcp);
    		ambi_decode2_copy_row2buf(x, i);
    		for(j=i+1; j<n_ambi; j++)
    		{
    			dv += n_ambi2;
    			rcp = -(*dv);
    			ambi_decode2_mul_buf_and_add2row(x, j, rcp);
    		}
    	}
    }
/* make 0 above the main diagonale */
    for(i=n_ambi-1; i>=0; i--)
    {
    	dv = db + i*n_ambi2 + i;
    	ambi_decode2_copy_row2buf(x, i);
    	for(j=i-1; j>=0; j--)
    	{
    		dv -= n_ambi2;
    		rcp = -(*dv);
    		ambi_decode2_mul_buf_and_add2row(x, j, rcp);
    	}
    }
post("matrix_inverse regular");
}
static void ambi_decode2_pseudo_inverse(t_ambi_decode2 *x, t_symbol *s, int argc, t_atom *argv)
{
    t_atom *at=x->x_at;
    int i, n=x->x_n_ls*x->x_n_ambi;
    double *dv1=x->x_prod;
    double *dv2=x->x_prod;
    double mw=x->x_mirror_weight;
ambi_decode2_transp_back(x);
    ambi_decode2_mul1(x);
    ambi_decode2_inverse(x);
    ambi_decode2_mul2(x);
    at += 2;
    for(i=0; i<n; i++)
    {
    	SETFLOAT(at, (float)(*dv1));
    	dv1++;
    	at++;
    }
dv2 += n;
    n=x->x_n_mir_ls*x->x_n_ambi;
    dv2 += n;
    for(i=0; i<n; i++)
    {
    	SETFLOAT(at, (float)(*dv1 + *dv2*mw));
    	dv1++;
    	dv2++;
    	at++;
    }
outlet_anything(x->x_obj.ob_outlet, x->x_s_matrix, x->x_n_ambi*(x->x_n_ls+x->x_n_mir_ls)+2, x->x_at);
}
static void ambi_decode2_encode_ls_2d(t_ambi_decode2 *x, int argc, t_atom *argv, int mode)
{
    double phi;
    double *dw = x->x_transp;
    int index;
    int order=x->x_n_order;
if(argc < 2)
    {
    	post("ambi_decode2 ERROR: ls-input needs 1 index and 1 angle: ls_index + phi [degree]");
    	return;
    }
    index = (int)atom_getint(argv++) - 1;
    phi = (double)atom_getfloat(argv);
if(index < 0)
    	index = 0;
if(mode == AMBI_LS_IND)
    {
    	if(index >= x->x_n_ls)
    		index = x->x_n_ls - 1;
    }
    else if(mode == AMBI_LS_MRG)
    {
    	if(x->x_n_mir_ls)
    	{
    		if(index >= x->x_n_mir_ls)
    			index = x->x_n_mir_ls - 1;
    		index += x->x_n_ls;
    	}
    	else
    		return;
    }
    else if(mode == AMBI_LS_MIR)
    {
    	if(x->x_n_mir_ls)
    	{
    		if(index >= x->x_n_mir_ls)
    			index = x->x_n_mir_ls - 1;
    		index += x->x_n_ls;
    		index += x->x_n_mir_ls;
    	}
    	else
    		return;
    }
    else if(mode == AMBI_LS_PHT)
    {
    	if(x->x_n_ph_ls)
    	{
    		if(index >= x->x_n_ph_ls)
    			index = x->x_n_ph_ls - 1;
    		index += x->x_n_ls;
    		index += 2*x->x_n_mir_ls;
    	}
    	else
    		return;
    }
    else
    	return;
    
    phi *= x->x_pi_over_180;
dw += index * x->x_n_ambi;
*dw++ = 1.0;
    *dw++ = cos(phi);
    *dw++ = sin(phi);
if(order >= 2)
    {
    	*dw++ = cos(2.0*phi);
    	*dw++ = sin(2.0*phi);
if(order >= 3)
    	{
    		*dw++ = cos(3.0*phi);
    		*dw++ = sin(3.0*phi);
    		if(order >= 4)
    		{
    			*dw++ = cos(4.0*phi);
    			*dw++ = sin(4.0*phi);
if(order >= 5)
    			{
    				*dw++ = cos(5.0*phi);
    				*dw++ = sin(5.0*phi);
if(order >= 6)
    				{
    					*dw++ = cos(6.0*phi);
    					*dw++ = sin(6.0*phi);
if(order >= 7)
    					{
    						*dw++ = cos(7.0*phi);
    						*dw++ = sin(7.0*phi);
if(order >= 8)
    						{
    							*dw++ = cos(8.0*phi);
    							*dw++ = sin(8.0*phi);
if(order >= 9)
    							{
    								*dw++ = cos(9.0*phi);
    								*dw++ = sin(9.0*phi);
if(order >= 10)
    								{
    									*dw++ = cos(10.0*phi);
    									*dw++ = sin(10.0*phi);
if(order >= 11)
    									{
    										*dw++ = cos(11.0*phi);
    										*dw++ = sin(11.0*phi);
if(order >= 12)
    										{
    											*dw++ = cos(12.0*phi);
    											*dw++ = sin(12.0*phi);
    										}
    									}
    								}
    							}
    						}
    					}
    				}
    			}
    		}
    	}
    }
}
static void ambi_decode2_encode_ls_3d(t_ambi_decode2 *x, int argc, t_atom *argv, int mode)
{
    double delta, phi;
    double cd, sd, cd2, cd3, sd2, csd, cp, sp, cp2, sp2, cp3, sp3, cp4, sp4;
    double *dw = x->x_transp;
    int index;
    int order=x->x_n_order;
if(argc < 3)
    {
    	post("ambi_decode2 ERROR: ls-input needs 1 index and 2 angles: ls index + delta [degree] + phi [degree]");
    	return;
    }
    index = (int)atom_getint(argv++) - 1;
    delta = atom_getfloat(argv++);
    phi = atom_getfloat(argv);
if(index < 0)
    	index = 0;
    
    if(mode == AMBI_LS_IND)
    {
    	if(index >= x->x_n_ls)
    		index = x->x_n_ls - 1;
    }
    else if(mode == AMBI_LS_MRG)
    {
    	if(x->x_n_mir_ls)
    	{
    		if(index >= x->x_n_mir_ls)
    			index = x->x_n_mir_ls - 1;
    		index += x->x_n_ls;
    	}
    	else
    		return;
    }
    else if(mode == AMBI_LS_MIR)
    {
    	if(x->x_n_mir_ls)
    	{
    		if(index >= x->x_n_mir_ls)
    			index = x->x_n_mir_ls - 1;
    		index += x->x_n_ls;
    		index += x->x_n_mir_ls;
    	}
    	else
    		return;
    }
    else if(mode == AMBI_LS_PHT)
    {
    	if(x->x_n_ph_ls)
    	{
    		if(index >= x->x_n_ph_ls)
    			index = x->x_n_ph_ls - 1;
    		index += x->x_n_ls;
    		index += 2*x->x_n_mir_ls;
    	}
    	else
    		return;
    }
    else
    	return;
delta *= x->x_pi_over_180;
    phi *= x->x_pi_over_180;
dw += index * x->x_n_ambi;
cd = cos(delta);
    sd = sin(delta);
    cp = cos(phi);
    sp = sin(phi);
*dw++ = 1.0;
    *dw++ = cd * cp;
    *dw++ = cd * sp;
    *dw++ = sd;
if(order >= 2)
    {
    	cp2 = cos(2.0*phi);
    	sp2 = sin(2.0*phi);
    	cd2 = cd * cd;
    	sd2 = sd * sd;
    	csd = cd * sd;
    	*dw++ = 0.5 * x->x_sqrt3 * cd2 * cp2;
    	*dw++ = 0.5 * x->x_sqrt3 * cd2 * sp2;
    	*dw++ = x->x_sqrt3 * csd * cp;
    	*dw++ = x->x_sqrt3 * csd * sp;
    	*dw++ = 0.5 * (3.0 * sd2 - 1.0);
if(order >= 3)
    	{
    		cp3 = cos(3.0*phi);
    		sp3 = sin(3.0*phi);
    		cd3 = cd2 * cd;
    		*dw++ = x->x_sqrt10_4 * cd3 * cp3;
    		*dw++ = x->x_sqrt10_4 * cd3 * sp3;
    		*dw++ = x->x_sqrt15_2 * cd * csd * cp2;
    		*dw++ = x->x_sqrt15_2 * cd * csd * sp2;
    		*dw++ = x->x_sqrt6_4 * cd * (5.0 * sd2 - 1.0) * cp;
    		*dw++ = x->x_sqrt6_4 * cd * (5.0 * sd2 - 1.0) * sp;
    		*dw++ = 0.5 * sd * (5.0 * sd2 - 3.0);
if(order >= 4)
    		{
    			cp4 = cos(4.0*phi);
    			sp4 = sin(4.0*phi);
    			*dw++ = x->x_sqrt35_8 * cd2 * cd2 * cp4;
    			*dw++ = x->x_sqrt35_8 * cd2 * cd2 * sp4;
    			*dw++ = x->x_sqrt70_4 * cd2 * csd * cp3;
    			*dw++ = x->x_sqrt70_4 * cd2 * csd * sp3;
    			*dw++ = 0.5 * x->x_sqrt5_2 * cd2 * (7.0 * sd2 - 1.0) * cp2;
    			*dw++ = 0.5 * x->x_sqrt5_2 * cd2 * (7.0 * sd2 - 1.0) * sp2;
    			*dw++ = x->x_sqrt10_4 * csd * (7.0 * sd2 - 3.0) * cp;
    			*dw++ = x->x_sqrt10_4 * csd * (7.0 * sd2 - 3.0) * sp;
    			*dw++ = 0.125 * (sd2 * (35.0 * sd2 - 30.0) + 3.0);
if(order >= 5)
    			{
    				*dw++ = x->x_sqrt126_16 * cd3 * cd2 * cos(5.0*phi);
    				*dw++ = x->x_sqrt126_16 * cd3 * cd2 * sin(5.0*phi);
    				*dw++ = x->x_sqrt315_8 * cd3 * csd * cp4;
    				*dw++ = x->x_sqrt315_8 * cd3 * csd * sp4;
    				*dw++ = 0.25 * x->x_sqrt70_4 * cd3 * (9.0 * sd2 - 1.0) * cp3;
    				*dw++ = 0.25 * x->x_sqrt70_4 * cd3 * (9.0 * sd2 - 1.0) * sp3;
    				*dw++ = x->x_sqrt105_4 * cd * csd * (3.0 * sd2 - 1.0) * cp2;
    				*dw++ = x->x_sqrt105_4 * cd * csd * (3.0 * sd2 - 1.0) * sp2;
    				*dw++ = 0.25 * x->x_sqrt15_2 * cd * (sd2 * (21.0 * sd2 - 14.0) + 1.0) * cp;
    				*dw++ = 0.25 * x->x_sqrt15_2 * cd * (sd2 * (21.0 * sd2 - 14.0) + 1.0) * sp;
    				*dw = 0.125 * sd * (sd2 * (63.0 * sd2 - 70.0) + 15.0);
    			}
    		}
    	}
    }
}
static void ambi_decode2_ind_ls(t_ambi_decode2 *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode2_encode_ls_2d(x, argc, argv, AMBI_LS_IND);
    else
    	ambi_decode2_encode_ls_3d(x, argc, argv, AMBI_LS_IND);
}
static void ambi_decode2_mrg_ls(t_ambi_decode2 *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode2_encode_ls_2d(x, argc, argv, AMBI_LS_MRG);
    else
    	ambi_decode2_encode_ls_3d(x, argc, argv, AMBI_LS_MRG);
}
static void ambi_decode2_mir_ls(t_ambi_decode2 *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode2_encode_ls_2d(x, argc, argv, AMBI_LS_MIR);
    else
    	ambi_decode2_encode_ls_3d(x, argc, argv, AMBI_LS_MIR);
}
static void ambi_decode2_pht_ls(t_ambi_decode2 *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode2_encode_ls_2d(x, argc, argv, AMBI_LS_PHT);
    else
    	ambi_decode2_encode_ls_3d(x, argc, argv, AMBI_LS_PHT);
}
static void ambi_decode2_ambi_weight(t_ambi_decode2 *x, t_symbol *s, int argc, t_atom *argv)
{
    if(argc > x->x_n_order)
    {
    	int i, k=0, n=x->x_n_order;
    	double d;
x->x_ambi_channel_weight[k] = atom_getfloat(argv++);
    	k++;
    	if(x->x_n_dim == 2)
    	{
    		for(i=1; i<=n; i++)
    		{
    			d = atom_getfloat(argv++);
    			x->x_ambi_channel_weight[k] = d;
    			k++;
    			x->x_ambi_channel_weight[k] = d;
    			k++;
    		}
    	}
    	else
    	{
    		int j, m;
for(i=1; i<=n; i++)
    		{
    			d = atom_getfloat(argv++);
    			m = 2*i + 1;
    			for(j=0; j<m; j++)
    			{
    				x->x_ambi_channel_weight[k] = d;
    				k++;
    			}
    		}
    	}
    }
    else
    	post("ambi_decode2-ERROR: ambi_weight needs %d float weights", x->x_n_order+1);
}
static void ambi_decode2_mirror_weight(t_ambi_decode2 *x, t_floatarg f)
{
    x->x_mirror_weight = (double)f;
}
static void ambi_decode2_sing_range(t_ambi_decode2 *x, t_floatarg f)
{
    if(f < 0.0f)
    	x->x_sing_range = -(double)f;
    else
    	x->x_sing_range = (double)f;
}
static void ambi_decode2_free(t_ambi_decode2 *x)
{
    freebytes(x->x_inv_work1, x->x_n_ambi * x->x_n_ambi * sizeof(double));
    freebytes(x->x_inv_work2, 2 * x->x_n_ambi * x->x_n_ambi * sizeof(double));
    freebytes(x->x_inv_buf2, 2 * x->x_n_ambi * sizeof(double));
    freebytes(x->x_transp, (x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_ls_encode, (x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_prod, (x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double));
    freebytes(x->x_at, ((x->x_n_ls+x->x_n_mir_ls) * x->x_n_ambi + 2) * sizeof(t_atom));
}
static void *ambi_decode2_new(t_symbol *s, int argc, t_atom *argv)
{
    t_ambi_decode2 *x = (t_ambi_decode2 *)pd_new(ambi_decode2_class);
    int order, dim, i;
    int n_ls=0;/* number of loudspeakers */
    int n_mir_ls=0;/* number of mirror_loudspeakers */
    int n_ph_ls=0;/* number of phantom_loudspeakers */
if((argc >= 5) &&
    	IS_A_FLOAT(argv,0) &&
    	IS_A_FLOAT(argv,1) &&
    	IS_A_FLOAT(argv,2) &&
    	IS_A_FLOAT(argv,3) &&
    	IS_A_FLOAT(argv,4))
    {
    	order = (int)atom_getint(argv++);
    	dim = (int)atom_getint(argv++);
    	n_ls = (int)atom_getint(argv++);
    	n_mir_ls = (int)atom_getint(argv++);
    	n_ph_ls = (int)atom_getint(argv);
if(order < 1)
    		order = 1;
    	if(dim != 3)
    	{
    		dim = 2;
    		if(order > 12)
    			order = 12;
    		x->x_n_ambi = 2*order + 1;
    	}
    	else
    	{
    		if(order > 5)
    			order = 5;
    		x->x_n_ambi = (order + 1)*(order + 1);
    	}
    	x->x_n_dim = dim;
    	x->x_n_order = order;
    	if(n_ls < 1)
    		n_ls = 1;
    	if(n_mir_ls < 0)
    		n_mir_ls = 0;
    	if(n_ph_ls < 0)
    		n_ph_ls = 0;
    	if((n_ls + 2*n_mir_ls + n_ph_ls) < x->x_n_ambi)
    		post("ambi_decode2-WARNING: Number of Loudspeakers < Number of Ambisonic-Channels !!!!");
    	
    	x->x_n_ls = n_ls;
    	x->x_n_mir_ls = n_mir_ls;
    	x->x_n_ph_ls = n_ph_ls;
    	x->x_inv_work1 = (double *)getbytes(x->x_n_ambi * x->x_n_ambi * sizeof(double));
    	x->x_inv_work2 = (double *)getbytes(2 * x->x_n_ambi * x->x_n_ambi * sizeof(double));
    	x->x_inv_buf2 = (double *)getbytes(2 * x->x_n_ambi * sizeof(double));
    	x->x_transp = (double *)getbytes((x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double));
    	x->x_ls_encode = (double *)getbytes((x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double));
    	x->x_prod = (double *)getbytes((x->x_n_ls+2*x->x_n_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double));
    	x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double));
    	x->x_at = (t_atom *)getbytes(((x->x_n_ls+x->x_n_mir_ls) * x->x_n_ambi + 2) * sizeof(t_atom));
    	x->x_s_matrix = gensym("matrix");
    	/*change*/
    	SETFLOAT(x->x_at, (float)(x->x_n_ls+x->x_n_mir_ls));
    	SETFLOAT(x->x_at+1, (float)x->x_n_ambi);
    	x->x_mirror_weight = 0.0;
x->x_sqrt3				= sqrt(3.0);
    	x->x_sqrt5_2			= sqrt(5.0) / 2.0;
    	x->x_sqrt6_4			= sqrt(6.0) / 4.0;
    	x->x_sqrt10_4			= sqrt(10.0) / 4.0;
    	x->x_sqrt15_2			= sqrt(15.0) / 2.0;
    	x->x_sqrt35_8			= sqrt(35.0) / 8.0;
    	x->x_sqrt70_4			= sqrt(70.0) / 4.0;
    	x->x_sqrt126_16		= sqrt(126.0) / 16.0;
    	x->x_sqrt315_8		= sqrt(315.0) / 8.0;
    	x->x_sqrt105_4		= sqrt(105.0) / 4.0;
    	x->x_pi_over_180	= 4.0 * atan(1.0) / 180.0;
    	x->x_sing_range = 1.0e-10;
    	for(i=0; i<x->x_n_ambi; i++)
    		x->x_ambi_channel_weight[i] = 1.0;
    	outlet_new(&x->x_obj, &s_list);
    	return (x);
    }
    else
    {
    	post("ambi_decode2-ERROR: need 5 float arguments: ambi_order dimension number_of_independent_loudspeakers number_of_merged_and_mirrored_speakers number_of_canceled_phantom_speakers");
    	return(0);
    }
}
void ambi_decode2_setup(void)
{
    ambi_decode2_class = class_new(gensym("ambi_decode2"), (t_newmethod)ambi_decode2_new, (t_method)ambi_decode2_free,
    			   sizeof(t_ambi_decode2), 0, A_GIMME, 0);
    class_addmethod(ambi_decode2_class, (t_method)ambi_decode2_ind_ls, gensym("ind_ls"), A_GIMME, 0);
    class_addmethod(ambi_decode2_class, (t_method)ambi_decode2_mrg_ls, gensym("mrg_ls"), A_GIMME, 0);
    class_addmethod(ambi_decode2_class, (t_method)ambi_decode2_mir_ls, gensym("mir_ls"), A_GIMME, 0);
    class_addmethod(ambi_decode2_class, (t_method)ambi_decode2_pht_ls, gensym("pht_ls"), A_GIMME, 0);
    class_addmethod(ambi_decode2_class, (t_method)ambi_decode2_mirror_weight, gensym("mirror_weight"), A_DEFFLOAT, 0);
    class_addmethod(ambi_decode2_class, (t_method)ambi_decode2_ambi_weight, gensym("ambi_weight"), A_GIMME, 0);
    class_addmethod(ambi_decode2_class, (t_method)ambi_decode2_sing_range, gensym("sing_range"), A_DEFFLOAT, 0);
    class_addmethod(ambi_decode2_class, (t_method)ambi_decode2_pseudo_inverse, gensym("pseudo_inverse"), A_GIMME, 0);
    class_sethelpsymbol(ambi_decode2_class, gensym("iemhelp2/help-ambi_decode2"));
}
--- NEW FILE: iemlib.h ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iemlib.h written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifndef __IEMLIB_H__
#define __IEMLIB_H__
#define IS_A_POINTER(atom,index) ((atom+index)->a_type == A_POINTER)
#define IS_A_FLOAT(atom,index) ((atom+index)->a_type == A_FLOAT)
#define IS_A_SYMBOL(atom,index) ((atom+index)->a_type == A_SYMBOL)
#define IS_A_DOLLAR(atom,index) ((atom+index)->a_type == A_DOLLAR)
#define IS_A_DOLLSYM(atom,index) ((atom+index)->a_type == A_DOLLSYM)
#define IS_A_SEMI(atom,index) ((atom+index)->a_type == A_SEMI)
#define IS_A_COMMA(atom,index) ((atom+index)->a_type == A_COMMA)
#ifdef NT
int sys_noloadbang;
//t_symbol *iemgui_key_sym=0;
#include <io.h>
#else
extern int sys_noloadbang;
//extern t_symbol *iemgui_key_sym;
#include <unistd.h>
#endif
#define DEFDELVS 64
#define XTRASAMPS 4
#define SAMPBLK 4
#define UNITBIT32 1572864.  /* 3*2^19; bit 32 has place value 1 */
/* machine-dependent definitions.  These ifdefs really
should have been by CPU type and not by operating system! */
#ifdef IRIX
/* big-endian.  Most significant byte is at low address in memory */
#define HIOFFSET 0    /* word offset to find MSB */
#define LOWOFFSET 1    /* word offset to find LSB */
#define int32 long  /* a data type that has 32 bits */
#else
#ifdef MSW
/* little-endian; most significant byte is at highest address */
#define HIOFFSET 1
#define LOWOFFSET 0
#define int32 long
#else
#ifdef __FreeBSD__
#include <machine/endian.h>
#if BYTE_ORDER == LITTLE_ENDIAN
#define HIOFFSET 1
#define LOWOFFSET 0
#else
#define HIOFFSET 0    /* word offset to find MSB */
#define LOWOFFSET 1    /* word offset to find LSB */
#endif /* BYTE_ORDER */
#include <sys/types.h>
#define int32 int32_t
#endif
#ifdef __linux__
#include <endian.h>
#if !defined(__BYTE_ORDER) || !defined(__LITTLE_ENDIAN)                         
#error No byte order defined                                                    
#endif
#if __BYTE_ORDER == __LITTLE_ENDIAN                                             
#define HIOFFSET 1                                                              
#define LOWOFFSET 0                                                             
#else                                                                           
#define HIOFFSET 0    /* word offset to find MSB */                             
#define LOWOFFSET 1    /* word offset to find LSB */                            
#endif /* __BYTE_ORDER */
#include <sys/types.h>
#define int32 int32_t
#else
#ifdef __APPLE__
#define HIOFFSET 0    /* word offset to find MSB */
#define LOWOFFSET 1    /* word offset to find LSB */
#define int32 int  /* a data type that has 32 bits */
#endif /* __APPLE__ */
#endif /* __linux__ */
#endif /* MSW */
#endif /* SGI */
union tabfudge
{
  double tf_d;
  int32 tf_i[2];
};
#define IEM_DENORMAL(f) ((((*(unsigned int*)&(f))&0x60000000)==0) || \
(((*(unsigned int*)&(f))&0x60000000)==0x60000000))
/* more stringent test: anything not between 1e-19 and 1e19 in absolute val */
#endif
--- NEW FILE: iem_ambi_sources.h ---
/* iem_bin_ambi-setup autogenerated header-file
 * generated by "./makesource.sh"
 * !! DO NOT MANUALLY EDIT  !!
 */
#ifndef IEMAMBI_SOURCES_H__
#define IEMAMBI_SOURCES_H__
void bin_ambi_calc_HRTF_setup(void); /* bin_ambi_calc_HRTF.c */
void bin_ambi_reduced_decode_setup(void); /* bin_ambi_reduced_decode.c */
#endif /* IEMBINAMBI_SOURCES_H__ */
--- NEW FILE: ambi_decode_cube.c ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif
#include "m_pd.h"
#include "iemlib.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
/* -------------------------- ambi_decode_cube ------------------------------ */
/*
 ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren **
 ** Inputs: ls + Liste von 3 floats: Index [1 .. 25] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] **
 ** Inputs: calc_inv **
 ** Inputs: load_HRIR + float index1..25 **
 ** Outputs: List of 2 symbols: left-HRIR-File-name + HRIR-table-name **
 ** Inputs: calc_reduced **
 ** "output" ...  writes the HRTF into tables **
 **  **
 **  **
 ** setzt voraus , dass die HRIR-tabele-names von LS1_L_HRIR .. LS25_L_HRIR heissen und existieren **
 ** setzt voraus , dass die HRTF-tabele-names von LS1_HRTF_re .. LS25_HRTF_re heissen und existieren **
 ** setzt voraus , dass die HRTF-tabele-names von LS1_HRTF_im .. LS25_HRTF_im heissen und existieren **
 */
typedef struct _ambi_decode_cube
{
    t_object	x_obj;
    t_atom		*x_at;
    double		*x_inv_work1;
    double		*x_inv_work2;
    double		*x_inv_buf2;
    double		*x_transp;
    double		*x_ls_encode;
    double		*x_prod;
    double		*x_ambi_channel_weight;
    double		x_mir_wght;
    int				x_n_ambi;
    int				x_n_order;
    int				x_n_ls;
    int				x_n_phls;
    int				x_n_dim;
    int				x_realsum_beg;
    int				x_realsum_end;
    int				x_mirrorsum_beg;
    int				x_mirrorsum_end;
    t_symbol	*x_s_matrix;
    double		x_sqrt3;
    double		x_sqrt10_4;
    double		x_sqrt15_2;
    double		x_sqrt6_4;
    double		x_sqrt35_8;
    double		x_sqrt70_4;
    double		x_sqrt5_2;
    double		x_sqrt126_16;
    double		x_sqrt315_8;
    double		x_sqrt105_4;
    double		x_pi_over_180;
} t_ambi_decode_cube;
static t_class *ambi_decode_cube_class;
static void ambi_decode_cube_copy_row2buf(t_ambi_decode_cube *x, int row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*db++ = *dw++;
}
static void ambi_decode_cube_copy_buf2row(t_ambi_decode_cube *x, int row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*dw++ = *db++;
}
static void ambi_decode_cube_copy_row2row(t_ambi_decode_cube *x, int src_row, int dst_row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw_src=x->x_inv_work2;
    double *dw_dst=x->x_inv_work2;
dw_src += src_row*n_ambi2;
    dw_dst += dst_row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*dw_dst++ = *dw_src++;
}
static void ambi_decode_cube_xch_rows(t_ambi_decode_cube *x, int row1, int row2)
{
    ambi_decode_cube_copy_row2buf(x, row1);
    ambi_decode_cube_copy_row2row(x, row2, row1);
    ambi_decode_cube_copy_buf2row(x, row2);
}
static void ambi_decode_cube_mul_row(t_ambi_decode_cube *x, int row, double mul)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    {
    	(*dw) *= mul;
    	dw++;
    }
}
static void ambi_decode_cube_mul_buf_and_add2row(t_ambi_decode_cube *x, int row, double mul)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    {
    	*dw += (*db)*mul;
    	dw++;
    	db++;
    }
}
static int ambi_decode_cube_eval_which_element_of_col_not_zero(t_ambi_decode_cube *x, int col, int start_row)
{
    int n_ambi = x->x_n_ambi;
    int n_ambi2 = 2*n_ambi;
    int i, j;
    double *dw=x->x_inv_work2;
    int ret=-1;
dw += start_row*n_ambi2 + col;
    j = 0;
    for(i=start_row; i<n_ambi; i++)
    {
    	if((*dw > 1.0e-10) || (*dw < -1.0e-10))
    	{
    		ret = i;
    		i = n_ambi+1;
    	}
    	dw += n_ambi2;
    }
    return(ret);
}
static void ambi_decode_cube_mul1(t_ambi_decode_cube *x)
{
    double *vec1, *beg1=x->x_ls_encode;
    double *vec2, *beg2=x->x_ls_encode;
    double *inv=x->x_inv_work1;
    double sum;
    int n_ls=x->x_n_ls+x->x_n_phls;
    int n_ambi=x->x_n_ambi;
    int i, j, k;
for(k=0; k<n_ambi; k++)
    {
    	beg2=x->x_ls_encode;
    	for(j=0; j<n_ambi; j++)
    	{
    		sum = 0.0;
    		vec1 = beg1;
    		vec2 = beg2;
    		for(i=0; i<n_ls; i++)
    		{
    			sum += *vec1++ * *vec2++;
    		}
    		beg2 += n_ls;
    		*inv++ = sum;
    	}
    	beg1 += n_ls;
    }
}
static void ambi_decode_cube_mul2(t_ambi_decode_cube *x)
{
    int n_ls=x->x_n_ls+x->x_n_phls;
    int n_ambi=x->x_n_ambi;
    int n_ambi2=2*n_ambi;
    int i, j, k;
    double *vec1, *beg1=x->x_transp;
    double *vec2, *beg2=x->x_inv_work2+n_ambi;
    double *vec3=x->x_prod;
    double *acw_vec=x->x_ambi_channel_weight;
    double sum;
for(k=0; k<n_ls; k++)
    {
    	beg2=x->x_inv_work2+n_ambi;
    	for(j=0; j<n_ambi; j++)
    	{
    		sum = 0.0;
    		vec1 = beg1;
    		vec2 = beg2;
    		for(i=0; i<n_ambi; i++)
    		{
    			sum += *vec1++ * *vec2;
    			vec2 += n_ambi2;
    		}
    		beg2++;
    		*vec3++ = sum * acw_vec[j];
    	}
    	beg1 += n_ambi;
    }
}
static void ambi_decode_cube_transp_back(t_ambi_decode_cube *x)
{
    double *vec, *transp=x->x_transp;
    double *straight=x->x_ls_encode;
    int n_ls=x->x_n_ls+x->x_n_phls;
    int n_ambi=x->x_n_ambi;
    int i, j;
for(j=0; j<n_ambi; j++)
    {
    	vec = transp;
    	for(i=0; i<n_ls; i++)
    	{
    		*straight++ = *vec;
    		vec += n_ambi;
    	}
    	transp++;
    }
}
static void ambi_decode_cube_inverse(t_ambi_decode_cube *x)
{
    int n_ambi = x->x_n_ambi;
    int n_ambi2 = 2*n_ambi;
    int i, j, nz;
    int r,c;
    double *src=x->x_inv_work1;
    double *db=x->x_inv_work2;
    double rcp, *dv;
dv = db;
    for(i=0; i<n_ambi; i++) /* init */
    {
    	for(j=0; j<n_ambi; j++)
    	{
    		*dv++ = *src++;
    	}
    	for(j=0; j<n_ambi; j++)
    	{
    		if(j == i)
    			*dv++ = 1.0;
    		else
    			*dv++ = 0.0;
    	}
    }
/* make 1 in main-diagonale, and 0 below */
    for(i=0; i<n_ambi; i++)
    {
    	nz = ambi_decode_cube_eval_which_element_of_col_not_zero(x, i, i);
    	if(nz < 0)
    	{
    		post("ambi_decode_cube ERROR: matrix not regular !!!!");
    		return;
    	}
    	else
    	{
    		if(nz != i)
    			ambi_decode_cube_xch_rows(x, i, nz);
    		dv = db + i*n_ambi2 + i;
    		rcp = 1.0 /(*dv);
    		ambi_decode_cube_mul_row(x, i, rcp);
    		ambi_decode_cube_copy_row2buf(x, i);
    		for(j=i+1; j<n_ambi; j++)
    		{
    			dv += n_ambi2;
    			rcp = -(*dv);
    			ambi_decode_cube_mul_buf_and_add2row(x, j, rcp);
    		}
    	}
    }
/* make 0 above the main diagonale */
    for(i=n_ambi-1; i>=0; i--)
    {
    	dv = db + i*n_ambi2 + i;
    	ambi_decode_cube_copy_row2buf(x, i);
    	for(j=i-1; j>=0; j--)
    	{
    		dv -= n_ambi2;
    		rcp = -(*dv);
    		ambi_decode_cube_mul_buf_and_add2row(x, j, rcp);
    	}
    }
post("matrix_inverse regular");
}
static void ambi_decode_cube_pinv(t_ambi_decode_cube *x)
{
    t_atom *at=x->x_at;
ambi_decode_cube_transp_back(x);
    ambi_decode_cube_mul1(x);
    ambi_decode_cube_inverse(x);
    ambi_decode_cube_mul2(x);
    if((x->x_mirrorsum_end > x->x_mirrorsum_beg)&&
    		(x->x_realsum_end > x->x_realsum_beg)&&
    			((x->x_mirrorsum_end - x->x_mirrorsum_beg) == (x->x_realsum_end - x->x_realsum_beg)))
    {
    	double *mir=x->x_prod+x->x_mirrorsum_beg*x->x_n_ambi;
    	double *real=x->x_prod+x->x_realsum_beg*x->x_n_ambi;
    	double mwght=x->x_mir_wght;
    	int i, n=(x->x_mirrorsum_end - x->x_mirrorsum_beg)*x->x_n_ambi;
//		post("mirror");
    	for(i=0; i<n; i++)
    		real[i] += mir[i]*mwght;
    	
    	n = x->x_mirrorsum_beg*x->x_n_ambi;
    	real=x->x_prod;
    	SETFLOAT(at, (float)x->x_n_ambi);
    	at++;
    	SETFLOAT(at, (float)x->x_mirrorsum_beg);
    	at++;
    	for(i=0; i<n; i++)
    	{
    		SETFLOAT(at, (float)(*real));
    		real++;
    		at++;
    	}
    	outlet_anything(x->x_obj.ob_outlet, x->x_s_matrix, n+2, x->x_at);
    }
    else
    {
    	int i, n=x->x_n_ls*x->x_n_ambi;
    	double *dv=x->x_prod;
//		post("real");
    	SETFLOAT(at, (float)x->x_n_ambi);
    	at++;
    	SETFLOAT(at, (float)x->x_n_ls);
    	at++;
    	for(i=0; i<n; i++)
    	{
    		SETFLOAT(at, (float)(*dv));
    		dv++;
    		at++;
    	}
    	outlet_anything(x->x_obj.ob_outlet, x->x_s_matrix, n+2, x->x_at);
    }
}
static void ambi_decode_cube_encode_ls_2d(t_ambi_decode_cube *x, int argc, t_atom *argv, int ls0_ph1)
{
    double phi;
    double *dw = x->x_transp;
    int index;
    int n_ls=x->x_n_ls;
    int n_phls=x->x_n_phls;
    int order=x->x_n_order;
if(argc < 2)
    {
    	post("ambi_decode_cube ERROR: ls-input needs 1 index and 1 angle: ls_index + phi [degree]");
    	return;
    }
    index = (int)atom_getint(argv++) - 1;
    phi = (double)atom_getfloat(argv);
if(index < 0)
    	index = 0;
    if(ls0_ph1)
    {
    	if(n_phls)
    	{
    		if(index >= n_phls)
    			index = n_phls - 1;
    		index += n_ls;
    	}
    	else
    		return;
    }
    else
    {
    	if(index >= n_ls)
    		index = n_ls - 1;
    }
    
    phi *= x->x_pi_over_180;
dw += index * x->x_n_ambi;
*dw++ = 1.0;
    *dw++ = cos(phi);
    *dw++ = sin(phi);
if(order >= 2)
    {
    	*dw++ = cos(2.0*phi);
    	*dw++ = sin(2.0*phi);
if(order >= 3)
    	{
    		*dw++ = cos(3.0*phi);
    		*dw++ = sin(3.0*phi);
if(order >= 4)
    		{
    			*dw++ = cos(4.0*phi);
    			*dw++ = sin(4.0*phi);
if(order >= 5)
    			{
    				*dw++ = cos(5.0*phi);
    				*dw++ = sin(5.0*phi);
if(order >= 6)
    				{
    					*dw++ = cos(6.0*phi);
    					*dw++ = sin(6.0*phi);
if(order >= 7)
    					{
    						*dw++ = cos(7.0*phi);
    						*dw++ = sin(7.0*phi);
if(order >= 8)
    						{
    							*dw++ = cos(8.0*phi);
    							*dw++ = sin(8.0*phi);
if(order >= 9)
    							{
    								*dw++ = cos(9.0*phi);
    								*dw++ = sin(9.0*phi);
if(order >= 10)
    								{
    									*dw++ = cos(10.0*phi);
    									*dw++ = sin(10.0*phi);
if(order >= 11)
    									{
    										*dw++ = cos(11.0*phi);
    										*dw++ = sin(11.0*phi);
if(order >= 12)
    										{
    											*dw++ = cos(12.0*phi);
    											*dw++ = sin(12.0*phi);
    										}
    									}
    								}
    							}
    						}
    					}
    				}
    			}
    		}
    	}
    }
}
static void ambi_decode_cube_encode_ls_3d(t_ambi_decode_cube *x, int argc, t_atom *argv, int ls0_ph1)
{
    double delta, phi;
    double cd, sd, cd2, cd3, sd2, csd, cp, sp, cp2, sp2, cp3, sp3, cp4, sp4;
    double *dw = x->x_transp;
    int index;
    int n_ls=x->x_n_ls;
    int n_phls=x->x_n_phls;
    int order=x->x_n_order;
if(argc < 3)
    {
    	post("ambi_decode_cube ERROR: ls-input needs 1 index and 2 angles: ls index + delta [degree] + phi [degree]");
    	return;
    }
    index = (int)atom_getint(argv++) - 1;
    delta = atom_getfloat(argv++);
    phi = atom_getfloat(argv);
if(index < 0)
    	index = 0;
    if(ls0_ph1)
    {
    	if(n_phls)
    	{
    		if(index >= n_phls)
    			index = n_phls - 1;
    		index += n_ls;
    	}
    	else
    		return;
    }
    else
    {
    	if(index >= n_ls)
    		index = n_ls - 1;
    }
delta *= x->x_pi_over_180;
    phi *= x->x_pi_over_180;
dw += index * x->x_n_ambi;
cd = cos(delta);
    sd = sin(delta);
    cp = cos(phi);
    sp = sin(phi);
    
    *dw++ = 1.0;
    *dw++ = cd * cp;
    *dw++ = cd * sp;
    *dw++ = sd;
if(order >= 2)
    {
    	cp2 = cos(2.0*phi);
    	sp2 = sin(2.0*phi);
    	cd2 = cd * cd;
    	sd2 = sd * sd;
    	csd = cd * sd;
    	*dw++ = 0.5 * x->x_sqrt3 * cd2 * cp2;
    	*dw++ = 0.5 * x->x_sqrt3 * cd2 * sp2;
    	*dw++ = x->x_sqrt3 * csd * cp;
    	*dw++ = x->x_sqrt3 * csd * sp;
    	*dw++ = 0.5 * (3.0 * sd2 - 1.0);
if(order >= 3)
    	{
    		cp3 = cos(3.0*phi);
    		sp3 = sin(3.0*phi);
    		cd3 = cd2 * cd;
    		*dw++ = x->x_sqrt10_4 * cd3 * cp3;
    		*dw++ = x->x_sqrt10_4 * cd3 * sp3;
    		*dw++ = x->x_sqrt15_2 * cd * csd * cp2;
    		*dw++ = x->x_sqrt15_2 * cd * csd * sp2;
    		*dw++ = x->x_sqrt6_4 * cd * (5.0 * sd2 - 1.0) * cp;
    		*dw++ = x->x_sqrt6_4 * cd * (5.0 * sd2 - 1.0) * sp;
    		*dw++ = 0.5 * sd * (5.0 * sd2 - 3.0);
if(order >= 4)
    		{
    			cp4 = cos(4.0*phi);
    			sp4 = sin(4.0*phi);
    			*dw++ = x->x_sqrt35_8 * cd2 * cd2 * cp4;
    			*dw++ = x->x_sqrt35_8 * cd2 * cd2 * sp4;
    			*dw++ = x->x_sqrt70_4 * cd2 * csd * cp3;
    			*dw++ = x->x_sqrt70_4 * cd2 * csd * sp3;
    			*dw++ = 0.5 * x->x_sqrt5_2 * cd2 * (7.0 * sd2 - 1.0) * cp2;
    			*dw++ = 0.5 * x->x_sqrt5_2 * cd2 * (7.0 * sd2 - 1.0) * sp2;
    			*dw++ = x->x_sqrt10_4 * csd * (7.0 * sd2 - 3.0) * cp;
    			*dw++ = x->x_sqrt10_4 * csd * (7.0 * sd2 - 3.0) * sp;
    			*dw++ = 0.125 * (sd2 * (35.0 * sd2 - 30.0) + 3.0);
if(order >= 5)
    			{
    				*dw++ = x->x_sqrt126_16 * cd3 * cd2 * cos(5.0*phi);
    				*dw++ = x->x_sqrt126_16 * cd3 * cd2 * sin(5.0*phi);
    				*dw++ = x->x_sqrt315_8 * cd3 * csd * cp4;
    				*dw++ = x->x_sqrt315_8 * cd3 * csd * sp4;
    				*dw++ = 0.25 * x->x_sqrt70_4 * cd3 * (9.0 * sd2 - 1.0) * cp3;
    				*dw++ = 0.25 * x->x_sqrt70_4 * cd3 * (9.0 * sd2 - 1.0) * sp3;
    				*dw++ = x->x_sqrt105_4 * cd * csd * (3.0 * sd2 - 1.0) * cp2;
    				*dw++ = x->x_sqrt105_4 * cd * csd * (3.0 * sd2 - 1.0) * sp2;
    				*dw++ = 0.25 * x->x_sqrt15_2 * cd * (sd2 * (21.0 * sd2 - 14.0) + 1.0) * cp;
    				*dw++ = 0.25 * x->x_sqrt15_2 * cd * (sd2 * (21.0 * sd2 - 14.0) + 1.0) * sp;
    				*dw = 0.125 * sd * (sd2 * (63.0 * sd2 - 70.0) + 15.0);
    			}
    		}
    	}
    }
}
static void ambi_decode_cube_ls(t_ambi_decode_cube *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode_cube_encode_ls_2d(x, argc, argv, 0);
    else
    	ambi_decode_cube_encode_ls_3d(x, argc, argv, 0);
}
static void ambi_decode_cube_phls(t_ambi_decode_cube *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode_cube_encode_ls_2d(x, argc, argv, 1);
    else
    	ambi_decode_cube_encode_ls_3d(x, argc, argv, 1);
}
static void ambi_decode_cube_ambi_weight(t_ambi_decode_cube *x, t_symbol *s, int argc, t_atom *argv)
{
    if(argc > x->x_n_order)
    {
    	int i, k=0, n=x->x_n_order;
    	double d;
x->x_ambi_channel_weight[k] = atom_getfloat(argv++);
    	k++;
    	if(x->x_n_dim == 2)
    	{
    		for(i=1; i<=n; i++)
    		{
    			d = atom_getfloat(argv++);
    			x->x_ambi_channel_weight[k] = d;
    			k++;
    			x->x_ambi_channel_weight[k] = d;
    			k++;
    		}
    	}
    	else
    	{
    		int j, m;
for(i=1; i<=n; i++)
    		{
    			d = atom_getfloat(argv++);
    			m = 2*i + 1;
    			for(j=0; j<m; j++)
    			{
    				x->x_ambi_channel_weight[k] = d;
    				k++;
    			}
    		}
    	}
    }
    else
    	post("ambi_decode_cube-ERROR: ambi_weight needs %d float weights", x->x_n_order+1);
}
static void ambi_decode_cube_mirror_weight(t_ambi_decode_cube *x, t_floatarg mwght)
{
    x->x_mir_wght = mwght;
}
static void ambi_decode_cube_mirror_range(t_ambi_decode_cube *x, t_floatarg beg, t_floatarg end)
{
    int b=(int)beg;
    int e=(int)end;
if(b < 0)
    	b = 0;
    if(b > x->x_n_ls)
    	b = x->x_n_ls;
    if(e < 0)
    	e = 0;
    if(e > x->x_n_ls)
    	e = x->x_n_ls;
    x->x_mirrorsum_beg = b;
    x->x_mirrorsum_end = e;
}
static void ambi_decode_cube_real_sum_range(t_ambi_decode_cube *x, t_floatarg beg, t_floatarg end)
{
    int b=(int)beg;
    int e=(int)end;
if(b < 0)
    	b = 0;
    if(b > x->x_n_ls)
    	b = x->x_n_ls;
    if(e < 0)
    	e = 0;
    if(e > x->x_n_ls)
    	e = x->x_n_ls;
    x->x_realsum_beg = b;
    x->x_realsum_end = e;
}
static void ambi_decode_cube_free(t_ambi_decode_cube *x)
{
    freebytes(x->x_inv_work1, x->x_n_ambi * x->x_n_ambi * sizeof(double));
    freebytes(x->x_inv_work2, 2 * x->x_n_ambi * x->x_n_ambi * sizeof(double));
    freebytes(x->x_inv_buf2, 2 * x->x_n_ambi * sizeof(double));
    freebytes(x->x_transp, (x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_ls_encode, (x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_prod, (x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double));
    freebytes(x->x_at, (x->x_n_ls * x->x_n_ambi + 2) * sizeof(t_atom));
}
static void *ambi_decode_cube_new(t_symbol *s, int argc, t_atom *argv)
{
    t_ambi_decode_cube *x = (t_ambi_decode_cube *)pd_new(ambi_decode_cube_class);
    int nls, order, dim, i;
    int nphls=0;/* phantom_loudspeaker */
if(argc < 3)
    {
    	post("ambi_decode_cube-ERROR: need following arguments: ambi_order dimension number_of_loudspeakers (number_of_phantom_speakers)");
    	return(0);
    }
    else
    {
    	order = (int)atom_getint(argv++);
    	dim = (int)atom_getint(argv++);
    	nls = (int)atom_getint(argv++);
    	if((argc > 3)&&IS_A_FLOAT(argv,0))
    		nphls=(int)atom_getint(argv);
if(order < 1)
    		order = 1;
    	if(dim != 3)
    	{
    		dim = 2;
    		if(order > 12)
    			order = 12;
    		x->x_n_ambi = 2*order + 1;
    	}
    	else
    	{
    		if(order > 5)
    			order = 5;
    		x->x_n_ambi = (order + 1)*(order + 1);
    	}
    	x->x_n_dim = dim;
    	x->x_n_order = order;
    	if(nls < 1)
    		nls = 1;
    	if(nphls < 0)
    		nphls = 0;
    	if(nls < x->x_n_ambi)
    		post("ambi_decode_cube-WARNING: Number of Loudspeakers < Number of Ambisonic-Channels !!!!");
    	if(nphls > nls)
    	{
    		post("ambi_decode_cube-WARNING: Number of Phantom-Loudspeakers > Number of Loudspeakers !!!!");
    		nphls = nls;
    	}
    	x->x_n_ls = nls;
    	x->x_n_phls = nphls;
    	x->x_inv_work1 = (double *)getbytes(x->x_n_ambi * x->x_n_ambi * sizeof(double));
    	x->x_inv_work2 = (double *)getbytes(2 * x->x_n_ambi * x->x_n_ambi * sizeof(double));
    	x->x_inv_buf2 = (double *)getbytes(2 * x->x_n_ambi * sizeof(double));
    	x->x_transp = (double *)getbytes((x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    	x->x_ls_encode = (double *)getbytes((x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    	x->x_prod = (double *)getbytes((x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    	x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double));
    	x->x_at = (t_atom *)getbytes((x->x_n_ls * x->x_n_ambi + 2) * sizeof(t_atom));
    	x->x_s_matrix = gensym("matrix");
    	/*change*/
    	SETFLOAT(x->x_at, (float)x->x_n_ls);
    	SETFLOAT(x->x_at+1, (float)x->x_n_ambi);
    	x->x_sqrt3				= sqrt(3.0);
    	x->x_sqrt5_2			= sqrt(5.0) / 2.0;
    	x->x_sqrt6_4			= sqrt(6.0) / 4.0;
    	x->x_sqrt10_4			= sqrt(10.0) / 4.0;
    	x->x_sqrt15_2			= sqrt(15.0) / 2.0;
    	x->x_sqrt35_8			= sqrt(35.0) / 8.0;
    	x->x_sqrt70_4			= sqrt(70.0) / 4.0;
    	x->x_sqrt126_16		= sqrt(126.0) / 16.0;
    	x->x_sqrt315_8		= sqrt(315.0) / 8.0;
    	x->x_sqrt105_4		= sqrt(105.0) / 4.0;
    	x->x_pi_over_180	= 4.0 * atan(1.0) / 180.0;
    	x->x_mirrorsum_beg = x->x_n_ls;
    	x->x_mirrorsum_end = x->x_n_ls;
    	x->x_realsum_beg = 0;
    	x->x_realsum_end = 0;
    	for(i=0; i<x->x_n_ambi; i++)
    		x->x_ambi_channel_weight[i] = 1.0;
    	x->x_mir_wght = 1.0;
    	outlet_new(&x->x_obj, &s_list);
    	return (x);
    }
}
void ambi_decode_cube_setup(void)
{
    ambi_decode_cube_class = class_new(gensym("ambi_decode_cube"), (t_newmethod)ambi_decode_cube_new, (t_method)ambi_decode_cube_free,
    			   sizeof(t_ambi_decode_cube), 0, A_GIMME, 0);
    class_addmethod(ambi_decode_cube_class, (t_method)ambi_decode_cube_ls, gensym("ls"), A_GIMME, 0);
    class_addmethod(ambi_decode_cube_class, (t_method)ambi_decode_cube_phls, gensym("phls"), A_GIMME, 0);
    class_addmethod(ambi_decode_cube_class, (t_method)ambi_decode_cube_ambi_weight, gensym("ambi_weight"), A_GIMME, 0);
    class_addmethod(ambi_decode_cube_class, (t_method)ambi_decode_cube_pinv, gensym("pinv"), 0);
    class_addmethod(ambi_decode_cube_class, (t_method)ambi_decode_cube_mirror_weight, gensym("mirror_weight"), A_DEFFLOAT, 0);
    class_addmethod(ambi_decode_cube_class, (t_method)ambi_decode_cube_mirror_range, gensym("mirror_range"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_addmethod(ambi_decode_cube_class, (t_method)ambi_decode_cube_real_sum_range, gensym("real_sum_range"), A_DEFFLOAT, A_DEFFLOAT, 0);
    class_sethelpsymbol(ambi_decode_cube_class, gensym("iemhelp/help-ambi_decode_cube"));
}
--- NEW FILE: Make.config.in ---
LIBNAME    =iem_bin_ambi
PREFIX     =@prefix@@PDLIBDIR@
INSTALL_BIN=$(PREFIX)/extra
INSTALL_DOC=$(PREFIX)/@REFERENCEPATH@$(LIBNAME)
EXT = @EXT@ 
DEFS = @DFLAGS@
IFLAGS = -I. @INCLUDES@
CC = @CC@
LD = @LD@
STRIP = @STRIP@
STRIPFLAGS= @STRIPFLAGS@
AFLAGS = 
LFLAGS = @LFLAGS@
WFLAGS =
TARNAME =  $(LIBNAME)-@IEMBINAMBI_VERSION@.tgz
# ICCFLAGS=-march=pentiumiii -axK
Z_CFLAGS = $(IFLAGS) $(DEFS) -DPD $(WFLAGS) @CFLAGS@ $(CFLAGS)
MAKEDEP_FLAGS = @MAKEDEP_FLAGS@
LIBS = @LIBS@
--- NEW FILE: ambi_decode.c ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif
#include "m_pd.h"
#include "iemlib.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
/* -------------------------- ambi_decode ------------------------------ */
/*
 ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren **
 ** Inputs: ls + Liste von 3 floats: Index [1 .. 25] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] **
 ** Inputs: calc_inv **
 ** Inputs: load_HRIR + float index1..25 **
 ** Outputs: List of 2 symbols: left-HRIR-File-name + HRIR-table-name **
 ** Inputs: calc_reduced **
 ** "output" ...  writes the HRTF into tables **
 **  **
 **  **
 ** setzt voraus , dass die HRIR-tabele-names von LS1_L_HRIR .. LS25_L_HRIR heissen und existieren **
 ** setzt voraus , dass die HRTF-tabele-names von LS1_HRTF_re .. LS25_HRTF_re heissen und existieren **
 ** setzt voraus , dass die HRTF-tabele-names von LS1_HRTF_im .. LS25_HRTF_im heissen und existieren **
 */
typedef struct _ambi_decode
{
    t_object	x_obj;
    t_atom		*x_at;
    double		*x_inv_work1;
    double		*x_inv_work2;
    double		*x_inv_buf2;
    double		*x_transp;
    double		*x_ls_encode;
    double		*x_prod;
    double		*x_ambi_channel_weight;
    double		x_sing_range;
    int				x_n_ambi;
    int				x_n_order;
    int				x_n_ls;
    int				x_n_phls;
    int				x_n_dim;
    t_symbol	*x_s_matrix;
    double		x_sqrt3;
    double		x_sqrt10_4;
    double		x_sqrt15_2;
    double		x_sqrt6_4;
    double		x_sqrt35_8;
    double		x_sqrt70_4;
    double		x_sqrt5_2;
    double		x_sqrt126_16;
    double		x_sqrt315_8;
    double		x_sqrt105_4;
    double		x_pi_over_180;
} t_ambi_decode;
static t_class *ambi_decode_class;
static void ambi_decode_copy_row2buf(t_ambi_decode *x, int row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*db++ = *dw++;
}
static void ambi_decode_copy_buf2row(t_ambi_decode *x, int row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*dw++ = *db++;
}
static void ambi_decode_copy_row2row(t_ambi_decode *x, int src_row, int dst_row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw_src=x->x_inv_work2;
    double *dw_dst=x->x_inv_work2;
dw_src += src_row*n_ambi2;
    dw_dst += dst_row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*dw_dst++ = *dw_src++;
}
static void ambi_decode_xch_rows(t_ambi_decode *x, int row1, int row2)
{
    ambi_decode_copy_row2buf(x, row1);
    ambi_decode_copy_row2row(x, row2, row1);
    ambi_decode_copy_buf2row(x, row2);
}
static void ambi_decode_mul_row(t_ambi_decode *x, int row, double mul)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    {
    	(*dw) *= mul;
    	dw++;
    }
}
static void ambi_decode_mul_buf_and_add2row(t_ambi_decode *x, int row, double mul)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    {
    	*dw += (*db)*mul;
    	dw++;
    	db++;
    }
}
static int ambi_decode_eval_which_element_of_col_not_zero(t_ambi_decode *x, int col, int start_row)
{
    int n_ambi = x->x_n_ambi;
    int n_ambi2 = 2*n_ambi;
    int i, j;
    double *dw=x->x_inv_work2;
    double singrange=x->x_sing_range;
    int ret=-1;
dw += start_row*n_ambi2 + col;
    j = 0;
    for(i=start_row; i<n_ambi; i++)
    {
    	if((*dw > singrange) || (*dw < -singrange))
    	{
    		ret = i;
    		i = n_ambi+1;
    	}
    	dw += n_ambi2;
    }
    return(ret);
}
static void ambi_decode_mul1(t_ambi_decode *x)
{
    double *vec1, *beg1=x->x_ls_encode;
    double *vec2, *beg2=x->x_ls_encode;
    double *inv=x->x_inv_work1;
    double sum;
    int n_ls=x->x_n_ls+x->x_n_phls;
    int n_ambi=x->x_n_ambi;
    int i, j, k;
for(k=0; k<n_ambi; k++)
    {
    	beg2=x->x_ls_encode;
    	for(j=0; j<n_ambi; j++)
    	{
    		sum = 0.0;
    		vec1 = beg1;
    		vec2 = beg2;
    		for(i=0; i<n_ls; i++)
    		{
    			sum += *vec1++ * *vec2++;
    		}
    		beg2 += n_ls;
    		*inv++ = sum;
    	}
    	beg1 += n_ls;
    }
}
static void ambi_decode_mul2(t_ambi_decode *x)
{
    int n_ls=x->x_n_ls+x->x_n_phls;
    int n_ambi=x->x_n_ambi;
    int n_ambi2=2*n_ambi;
    int i, j, k;
    double *vec1, *beg1=x->x_transp;
    double *vec2, *beg2=x->x_inv_work2+n_ambi;
    double *vec3=x->x_prod;
    double *acw_vec=x->x_ambi_channel_weight;
    double sum;
for(k=0; k<n_ls; k++)
    {
    	beg2=x->x_inv_work2+n_ambi;
    	for(j=0; j<n_ambi; j++)
    	{
    		sum = 0.0;
    		vec1 = beg1;
    		vec2 = beg2;
    		for(i=0; i<n_ambi; i++)
    		{
    			sum += *vec1++ * *vec2;
    			vec2 += n_ambi2;
    		}
    		beg2++;
    		*vec3++ = sum * acw_vec[j];
    	}
    	beg1 += n_ambi;
    }
}
static void ambi_decode_transp_back(t_ambi_decode *x)
{
    double *vec, *transp=x->x_transp;
    double *straight=x->x_ls_encode;
    int n_ls=x->x_n_ls+x->x_n_phls;
    int n_ambi=x->x_n_ambi;
    int i, j;
for(j=0; j<n_ambi; j++)
    {
    	vec = transp;
    	for(i=0; i<n_ls; i++)
    	{
    		*straight++ = *vec;
    		vec += n_ambi;
    	}
    	transp++;
    }
}
static void ambi_decode_inverse(t_ambi_decode *x)
{
    int n_ambi = x->x_n_ambi;
    int n_ambi2 = 2*n_ambi;
    int i, j, nz;
    int r,c;
    double *src=x->x_inv_work1;
    double *db=x->x_inv_work2;
    double rcp, *dv;
dv = db;
    for(i=0; i<n_ambi; i++) /* init */
    {
    	for(j=0; j<n_ambi; j++)
    	{
    		*dv++ = *src++;
    	}
    	for(j=0; j<n_ambi; j++)
    	{
    		if(j == i)
    			*dv++ = 1.0;
    		else
    			*dv++ = 0.0;
    	}
    }
/* make 1 in main-diagonale, and 0 below */
    for(i=0; i<n_ambi; i++)
    {
    	nz = ambi_decode_eval_which_element_of_col_not_zero(x, i, i);
    	if(nz < 0)
    	{
    		post("ambi_decode ERROR: matrix not regular !!!!");
    		return;
    	}
    	else
    	{
    		if(nz != i)
    			ambi_decode_xch_rows(x, i, nz);
    		dv = db + i*n_ambi2 + i;
    		rcp = 1.0 /(*dv);
    		ambi_decode_mul_row(x, i, rcp);
    		ambi_decode_copy_row2buf(x, i);
    		for(j=i+1; j<n_ambi; j++)
    		{
    			dv += n_ambi2;
    			rcp = -(*dv);
    			ambi_decode_mul_buf_and_add2row(x, j, rcp);
    		}
    	}
    }
/* make 0 above the main diagonale */
    for(i=n_ambi-1; i>=0; i--)
    {
    	dv = db + i*n_ambi2 + i;
    	ambi_decode_copy_row2buf(x, i);
    	for(j=i-1; j>=0; j--)
    	{
    		dv -= n_ambi2;
    		rcp = -(*dv);
    		ambi_decode_mul_buf_and_add2row(x, j, rcp);
    	}
    }
post("matrix_inverse regular");
}
static void ambi_decode_pinv(t_ambi_decode *x)
{
    t_atom *at=x->x_at;
    int i, n=x->x_n_ls*x->x_n_ambi;
    double *dv=x->x_prod;
ambi_decode_transp_back(x);
    ambi_decode_mul1(x);
    ambi_decode_inverse(x);
    ambi_decode_mul2(x);
    at += 2;
    for(i=0; i<n; i++)
    {
    	SETFLOAT(at, (float)(*dv));
    	dv++;
    	at++;
    }
    outlet_anything(x->x_obj.ob_outlet, x->x_s_matrix, n+2, x->x_at);
}
static void ambi_decode_encode_ls_2d(t_ambi_decode *x, int argc, t_atom *argv, int ls0_ph1)
{
    double phi;
    double *dw = x->x_transp;
    int index;
    int n_ls=x->x_n_ls;
    int n_phls=x->x_n_phls;
    int order=x->x_n_order;
if(argc < 2)
    {
    	post("ambi_decode ERROR: ls-input needs 1 index and 1 angle: ls_index + phi [degree]");
    	return;
    }
    index = (int)atom_getint(argv++) - 1;
    phi = (double)atom_getfloat(argv);
if(index < 0)
    	index = 0;
    if(ls0_ph1)
    {
    	if(n_phls)
    	{
    		if(index >= n_phls)
    			index = n_phls - 1;
    		index += n_ls;
    	}
    	else
    		return;
    }
    else
    {
    	if(index >= n_ls)
    		index = n_ls - 1;
    }
    
    phi *= x->x_pi_over_180;
dw += index * x->x_n_ambi;
*dw++ = 1.0;
    *dw++ = cos(phi);
    *dw++ = sin(phi);
if(order >= 2)
    {
    	*dw++ = cos(2.0*phi);
    	*dw++ = sin(2.0*phi);
if(order >= 3)
    	{
    		*dw++ = cos(3.0*phi);
    		*dw++ = sin(3.0*phi);
    		if(order >= 4)
    		{
    			*dw++ = cos(4.0*phi);
    			*dw++ = sin(4.0*phi);
if(order >= 5)
    			{
    				*dw++ = cos(5.0*phi);
    				*dw++ = sin(5.0*phi);
if(order >= 6)
    				{
    					*dw++ = cos(6.0*phi);
    					*dw++ = sin(6.0*phi);
if(order >= 7)
    					{
    						*dw++ = cos(7.0*phi);
    						*dw++ = sin(7.0*phi);
if(order >= 8)
    						{
    							*dw++ = cos(8.0*phi);
    							*dw++ = sin(8.0*phi);
if(order >= 9)
    							{
    								*dw++ = cos(9.0*phi);
    								*dw++ = sin(9.0*phi);
if(order >= 10)
    								{
    									*dw++ = cos(10.0*phi);
    									*dw++ = sin(10.0*phi);
if(order >= 11)
    									{
    										*dw++ = cos(11.0*phi);
    										*dw++ = sin(11.0*phi);
if(order >= 12)
    										{
    											*dw++ = cos(12.0*phi);
    											*dw++ = sin(12.0*phi);
    										}
    									}
    								}
    							}
    						}
    					}
    				}
    			}
    		}
    	}
    }
}
static void ambi_decode_encode_ls_3d(t_ambi_decode *x, int argc, t_atom *argv, int ls0_ph1)
{
    double delta, phi;
    double cd, sd, cd2, cd3, sd2, csd, cp, sp, cp2, sp2, cp3, sp3, cp4, sp4;
    double *dw = x->x_transp;
    int index;
    int n_ls=x->x_n_ls;
    int n_phls=x->x_n_phls;
    int order=x->x_n_order;
if(argc < 3)
    {
    	post("ambi_decode ERROR: ls-input needs 1 index and 2 angles: ls index + delta [degree] + phi [degree]");
    	return;
    }
    index = (int)atom_getint(argv++) - 1;
    delta = atom_getfloat(argv++);
    phi = atom_getfloat(argv);
if(index < 0)
    	index = 0;
    if(ls0_ph1)
    {
    	if(n_phls)
    	{
    		if(index >= n_phls)
    			index = n_phls - 1;
    		index += n_ls;
    	}
    	else
    		return;
    }
    else
    {
    	if(index >= n_ls)
    		index = n_ls - 1;
    }
delta *= x->x_pi_over_180;
    phi *= x->x_pi_over_180;
dw += index * x->x_n_ambi;
cd = cos(delta);
    sd = sin(delta);
    cp = cos(phi);
    sp = sin(phi);
*dw++ = 1.0;
    *dw++ = cd * cp;
    *dw++ = cd * sp;
    *dw++ = sd;
if(order >= 2)
    {
    	cp2 = cos(2.0*phi);
    	sp2 = sin(2.0*phi);
    	cd2 = cd * cd;
    	sd2 = sd * sd;
    	csd = cd * sd;
    	*dw++ = 0.5 * x->x_sqrt3 * cd2 * cp2;
    	*dw++ = 0.5 * x->x_sqrt3 * cd2 * sp2;
    	*dw++ = x->x_sqrt3 * csd * cp;
    	*dw++ = x->x_sqrt3 * csd * sp;
    	*dw++ = 0.5 * (3.0 * sd2 - 1.0);
if(order >= 3)
    	{
    		cp3 = cos(3.0*phi);
    		sp3 = sin(3.0*phi);
    		cd3 = cd2 * cd;
    		*dw++ = x->x_sqrt10_4 * cd3 * cp3;
    		*dw++ = x->x_sqrt10_4 * cd3 * sp3;
    		*dw++ = x->x_sqrt15_2 * cd * csd * cp2;
    		*dw++ = x->x_sqrt15_2 * cd * csd * sp2;
    		*dw++ = x->x_sqrt6_4 * cd * (5.0 * sd2 - 1.0) * cp;
    		*dw++ = x->x_sqrt6_4 * cd * (5.0 * sd2 - 1.0) * sp;
    		*dw++ = 0.5 * sd * (5.0 * sd2 - 3.0);
if(order >= 4)
    		{
    			cp4 = cos(4.0*phi);
    			sp4 = sin(4.0*phi);
    			*dw++ = x->x_sqrt35_8 * cd2 * cd2 * cp4;
    			*dw++ = x->x_sqrt35_8 * cd2 * cd2 * sp4;
    			*dw++ = x->x_sqrt70_4 * cd2 * csd * cp3;
    			*dw++ = x->x_sqrt70_4 * cd2 * csd * sp3;
    			*dw++ = 0.5 * x->x_sqrt5_2 * cd2 * (7.0 * sd2 - 1.0) * cp2;
    			*dw++ = 0.5 * x->x_sqrt5_2 * cd2 * (7.0 * sd2 - 1.0) * sp2;
    			*dw++ = x->x_sqrt10_4 * csd * (7.0 * sd2 - 3.0) * cp;
    			*dw++ = x->x_sqrt10_4 * csd * (7.0 * sd2 - 3.0) * sp;
    			*dw++ = 0.125 * (sd2 * (35.0 * sd2 - 30.0) + 3.0);
if(order >= 5)
    			{
    				*dw++ = x->x_sqrt126_16 * cd3 * cd2 * cos(5.0*phi);
    				*dw++ = x->x_sqrt126_16 * cd3 * cd2 * sin(5.0*phi);
    				*dw++ = x->x_sqrt315_8 * cd3 * csd * cp4;
    				*dw++ = x->x_sqrt315_8 * cd3 * csd * sp4;
    				*dw++ = 0.25 * x->x_sqrt70_4 * cd3 * (9.0 * sd2 - 1.0) * cp3;
    				*dw++ = 0.25 * x->x_sqrt70_4 * cd3 * (9.0 * sd2 - 1.0) * sp3;
    				*dw++ = x->x_sqrt105_4 * cd * csd * (3.0 * sd2 - 1.0) * cp2;
    				*dw++ = x->x_sqrt105_4 * cd * csd * (3.0 * sd2 - 1.0) * sp2;
    				*dw++ = 0.25 * x->x_sqrt15_2 * cd * (sd2 * (21.0 * sd2 - 14.0) + 1.0) * cp;
    				*dw++ = 0.25 * x->x_sqrt15_2 * cd * (sd2 * (21.0 * sd2 - 14.0) + 1.0) * sp;
    				*dw = 0.125 * sd * (sd2 * (63.0 * sd2 - 70.0) + 15.0);
    			}
    		}
    	}
    }
}
static void ambi_decode_ls(t_ambi_decode *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode_encode_ls_2d(x, argc, argv, 0);
    else
    	ambi_decode_encode_ls_3d(x, argc, argv, 0);
}
static void ambi_decode_phls(t_ambi_decode *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode_encode_ls_2d(x, argc, argv, 1);
    else
    	ambi_decode_encode_ls_3d(x, argc, argv, 1);
}
static void ambi_decode_ambi_weight(t_ambi_decode *x, t_symbol *s, int argc, t_atom *argv)
{
    if(argc > x->x_n_order)
    {
    	int i, k=0, n=x->x_n_order;
    	double d;
x->x_ambi_channel_weight[k] = atom_getfloat(argv++);
    	k++;
    	if(x->x_n_dim == 2)
    	{
    		for(i=1; i<=n; i++)
    		{
    			d = atom_getfloat(argv++);
    			x->x_ambi_channel_weight[k] = d;
    			k++;
    			x->x_ambi_channel_weight[k] = d;
    			k++;
    		}
    	}
    	else
    	{
    		int j, m;
for(i=1; i<=n; i++)
    		{
    			d = atom_getfloat(argv++);
    			m = 2*i + 1;
    			for(j=0; j<m; j++)
    			{
    				x->x_ambi_channel_weight[k] = d;
    				k++;
    			}
    		}
    	}
    }
    else
    	post("ambi_decode-ERROR: ambi_weight needs %d float weights", x->x_n_order+1);
}
static void ambi_decode_sing_range(t_ambi_decode *x, t_floatarg f)
{
    if(f < 0.0f)
    	x->x_sing_range = -(double)f;
    else
    	x->x_sing_range = (double)f;
}
static void ambi_decode_free(t_ambi_decode *x)
{
    freebytes(x->x_inv_work1, x->x_n_ambi * x->x_n_ambi * sizeof(double));
    freebytes(x->x_inv_work2, 2 * x->x_n_ambi * x->x_n_ambi * sizeof(double));
    freebytes(x->x_inv_buf2, 2 * x->x_n_ambi * sizeof(double));
    freebytes(x->x_transp, (x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_ls_encode, (x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_prod, (x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double));
    freebytes(x->x_at, (x->x_n_ls * x->x_n_ambi + 2) * sizeof(t_atom));
}
static void *ambi_decode_new(t_symbol *s, int argc, t_atom *argv)
{
    t_ambi_decode *x = (t_ambi_decode *)pd_new(ambi_decode_class);
    int nls, order, dim, i;
    int nphls=0;/* phantom_loudspeaker */
if(argc < 3)
    {
    	post("ambi_decode-ERROR: need following arguments: ambi_order dimension number_of_loudspeakers (number_of_phantom_speakers)");
    	return(0);
    }
    else
    {
    	order = (int)atom_getint(argv++);
    	dim = (int)atom_getint(argv++);
    	nls = (int)atom_getint(argv++);
    	if((argc > 3)&&IS_A_FLOAT(argv,0))
    		nphls=(int)atom_getint(argv);
if(order < 1)
    		order = 1;
    	if(dim != 3)
    	{
    		dim = 2;
    		if(order > 12)
    			order = 12;
    		x->x_n_ambi = 2*order + 1;
    	}
    	else
    	{
    		if(order > 5)
    			order = 5;
    		x->x_n_ambi = (order + 1)*(order + 1);
    	}
    	x->x_n_dim = dim;
    	x->x_n_order = order;
    	if(nls < 1)
    		nls = 1;
    	if(nphls < 0)
    		nphls = 0;
    	if(nls < x->x_n_ambi)
    		post("ambi_decode-WARNING: Number of Loudspeakers < Number of Ambisonic-Channels !!!!");
    	if(nphls > nls)
    	{
    		post("ambi_decode-WARNING: Number of Phantom-Loudspeakers > Number of Loudspeakers !!!!");
    		nphls = nls;
    	}
    	x->x_n_ls = nls;
    	x->x_n_phls = nphls;
    	x->x_inv_work1 = (double *)getbytes(x->x_n_ambi * x->x_n_ambi * sizeof(double));
    	x->x_inv_work2 = (double *)getbytes(2 * x->x_n_ambi * x->x_n_ambi * sizeof(double));
    	x->x_inv_buf2 = (double *)getbytes(2 * x->x_n_ambi * sizeof(double));
    	x->x_transp = (double *)getbytes((x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    	x->x_ls_encode = (double *)getbytes((x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    	x->x_prod = (double *)getbytes((x->x_n_ls+x->x_n_phls) * x->x_n_ambi * sizeof(double));
    	x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double));
    	x->x_at = (t_atom *)getbytes((x->x_n_ls * x->x_n_ambi + 2) * sizeof(t_atom));
    	x->x_s_matrix = gensym("matrix");
    	/*change*/
    	SETFLOAT(x->x_at, (float)x->x_n_ls);
    	SETFLOAT(x->x_at+1, (float)x->x_n_ambi);
    	x->x_sqrt3				= sqrt(3.0);
    	x->x_sqrt5_2			= sqrt(5.0) / 2.0;
    	x->x_sqrt6_4			= sqrt(6.0) / 4.0;
    	x->x_sqrt10_4			= sqrt(10.0) / 4.0;
    	x->x_sqrt15_2			= sqrt(15.0) / 2.0;
    	x->x_sqrt35_8			= sqrt(35.0) / 8.0;
    	x->x_sqrt70_4			= sqrt(70.0) / 4.0;
    	x->x_sqrt126_16		= sqrt(126.0) / 16.0;
    	x->x_sqrt315_8		= sqrt(315.0) / 8.0;
    	x->x_sqrt105_4		= sqrt(105.0) / 4.0;
    	x->x_pi_over_180	= 4.0 * atan(1.0) / 180.0;
    	x->x_sing_range = 1.0e-10;
    	for(i=0; i<x->x_n_ambi; i++)
    		x->x_ambi_channel_weight[i] = 1.0;
    	outlet_new(&x->x_obj, &s_list);
    	return (x);
    }
}
void ambi_decode_setup(void)
{
    ambi_decode_class = class_new(gensym("ambi_decode"), (t_newmethod)ambi_decode_new, (t_method)ambi_decode_free,
    			   sizeof(t_ambi_decode), 0, A_GIMME, 0);
    class_addmethod(ambi_decode_class, (t_method)ambi_decode_ls, gensym("ls"), A_GIMME, 0);
    class_addmethod(ambi_decode_class, (t_method)ambi_decode_phls, gensym("phls"), A_GIMME, 0);
    class_addmethod(ambi_decode_class, (t_method)ambi_decode_ambi_weight, gensym("ambi_weight"), A_GIMME, 0);
    class_addmethod(ambi_decode_class, (t_method)ambi_decode_sing_range, gensym("sing_range"), A_DEFFLOAT, 0);
    class_addmethod(ambi_decode_class, (t_method)ambi_decode_pinv, gensym("pinv"), 0);
    class_sethelpsymbol(ambi_decode_class, gensym("iemhelp2/help-ambi_decode"));
}
--- NEW FILE: ambi_decode3.c ---
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */
#ifdef NT
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif
#include "m_pd.h"
#include "iemlib.h"
#include "iem_ambi.h"
#include <math.h>
#include <stdio.h>
#include <string.h>
/* -------------------------- ambi_decode3 ------------------------------ */
/*
 ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren **
 ** Inputs: ls + Liste von 3 floats: Index [1 .. 25] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] **
 ** Inputs: calc_inv **
 ** Inputs: load_HRIR + float index1..25 **
 ** Outputs: List of 2 symbols: left-HRIR-File-name + HRIR-table-name **
 ** Inputs: calc_reduced **
 ** "output" ...  writes the HRTF into tables **
 **  **
 **  **
 ** setzt voraus , dass die HRIR-tabele-names von LS1_L_HRIR .. LS25_L_HRIR heissen und existieren **
 ** setzt voraus , dass die HRTF-tabele-names von LS1_HRTF_re .. LS25_HRTF_re heissen und existieren **
 ** setzt voraus , dass die HRTF-tabele-names von LS1_HRTF_im .. LS25_HRTF_im heissen und existieren **
 */
typedef struct _ambi_decode3
{
    t_object	x_obj;
    t_atom		*x_at;
    double		*x_inv_work1;
    double		*x_inv_work2;
    double		*x_inv_buf2;
    double		*x_transp;
    double		*x_ls_encode;
    double		*x_prod;
    double		*x_ambi_channel_weight;
    double		x_sing_range;
    int				x_n_ambi;
    int				x_n_order;
    int				x_n_real_ls;
    int				x_n_pht_ls;
    int				x_n_dim;
    t_symbol	*x_s_matrix;
    double		x_sqrt3;
    double		x_sqrt10_4;
    double		x_sqrt15_2;
    double		x_sqrt6_4;
    double		x_sqrt35_8;
    double		x_sqrt70_4;
    double		x_sqrt5_2;
    double		x_sqrt126_16;
    double		x_sqrt315_8;
    double		x_sqrt105_4;
    double		x_pi_over_180;
} t_ambi_decode3;
static t_class *ambi_decode3_class;
static void ambi_decode3_copy_row2buf(t_ambi_decode3 *x, int row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*db++ = *dw++;
}
static void ambi_decode3_copy_buf2row(t_ambi_decode3 *x, int row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*dw++ = *db++;
}
static void ambi_decode3_copy_row2row(t_ambi_decode3 *x, int src_row, int dst_row)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw_src=x->x_inv_work2;
    double *dw_dst=x->x_inv_work2;
dw_src += src_row*n_ambi2;
    dw_dst += dst_row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    	*dw_dst++ = *dw_src++;
}
static void ambi_decode3_xch_rows(t_ambi_decode3 *x, int row1, int row2)
{
    ambi_decode3_copy_row2buf(x, row1);
    ambi_decode3_copy_row2row(x, row2, row1);
    ambi_decode3_copy_buf2row(x, row2);
}
static void ambi_decode3_mul_row(t_ambi_decode3 *x, int row, double mul)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    {
    	(*dw) *= mul;
    	dw++;
    }
}
static void ambi_decode3_mul_buf_and_add2row(t_ambi_decode3 *x, int row, double mul)
{
    int n_ambi2 = 2*x->x_n_ambi;
    int i;
    double *dw=x->x_inv_work2;
    double *db=x->x_inv_buf2;
dw += row*n_ambi2;
    for(i=0; i<n_ambi2; i++)
    {
    	*dw += (*db)*mul;
    	dw++;
    	db++;
    }
}
static int ambi_decode3_eval_which_element_of_col_not_zero(t_ambi_decode3 *x, int col, int start_row)
{
    int n_ambi = x->x_n_ambi;
    int n_ambi2 = 2*n_ambi;
    int i, j;
    double *dw=x->x_inv_work2;
    double singrange=x->x_sing_range;
    int ret=-1;
dw += start_row*n_ambi2 + col;
    j = 0;
    for(i=start_row; i<n_ambi; i++)
    {
    	if((*dw > singrange) || (*dw < -singrange))
    	{
    		ret = i;
    		i = n_ambi+1;
    	}
    	dw += n_ambi2;
    }
    return(ret);
}
static void ambi_decode3_mul1(t_ambi_decode3 *x)
{
    double *vec1, *beg1=x->x_ls_encode;
    double *vec2, *beg2=x->x_ls_encode;
    double *inv=x->x_inv_work1;
    double sum;
    int n_ls=x->x_n_real_ls+x->x_n_pht_ls;
    int n_ambi=x->x_n_ambi;
    int i, j, k;
for(k=0; k<n_ambi; k++)
    {
    	beg2=x->x_ls_encode;
    	for(j=0; j<n_ambi; j++)
    	{
    		sum = 0.0;
    		vec1 = beg1;
    		vec2 = beg2;
    		for(i=0; i<n_ls; i++)
    		{
    			sum += *vec1++ * *vec2++;
    		}
    		beg2 += n_ls;
    		*inv++ = sum;
    	}
    	beg1 += n_ls;
    }
}
static void ambi_decode3_mul2(t_ambi_decode3 *x)
{
    int n_ls=x->x_n_real_ls+x->x_n_pht_ls;
    int n_ambi=x->x_n_ambi;
    int n_ambi2=2*n_ambi;
    int i, j, k;
    double *vec1, *beg1=x->x_transp;
    double *vec2, *beg2=x->x_inv_work2+n_ambi;
    double *vec3=x->x_prod;
    double *acw_vec=x->x_ambi_channel_weight;
    double sum;
for(k=0; k<n_ls; k++)
    {
    	beg2=x->x_inv_work2+n_ambi;
    	for(j=0; j<n_ambi; j++)
    	{
    		sum = 0.0;
    		vec1 = beg1;
    		vec2 = beg2;
    		for(i=0; i<n_ambi; i++)
    		{
    			sum += *vec1++ * *vec2;
    			vec2 += n_ambi2;
    		}
    		beg2++;
    		*vec3++ = sum * acw_vec[j];
    	}
    	beg1 += n_ambi;
    }
}
static void ambi_decode3_transp_back(t_ambi_decode3 *x)
{
    double *vec, *transp=x->x_transp;
    double *straight=x->x_ls_encode;
    int n_ls=x->x_n_real_ls+x->x_n_pht_ls;
    int n_ambi=x->x_n_ambi;
    int i, j;
for(j=0; j<n_ambi; j++)
    {
    	vec = transp;
    	for(i=0; i<n_ls; i++)
    	{
    		*straight++ = *vec;
    		vec += n_ambi;
    	}
    	transp++;
    }
}
static void ambi_decode3_inverse(t_ambi_decode3 *x)
{
    int n_ambi = x->x_n_ambi;
    int n_ambi2 = 2*n_ambi;
    int i, j, nz;
    int r,c;
    double *src=x->x_inv_work1;
    double *db=x->x_inv_work2;
    double rcp, *dv;
dv = db;
    for(i=0; i<n_ambi; i++) /* init */
    {
    	for(j=0; j<n_ambi; j++)
    	{
    		*dv++ = *src++;
    	}
    	for(j=0; j<n_ambi; j++)
    	{
    		if(j == i)
    			*dv++ = 1.0;
    		else
    			*dv++ = 0.0;
    	}
    }
/* make 1 in main-diagonale, and 0 below */
    for(i=0; i<n_ambi; i++)
    {
    	nz = ambi_decode3_eval_which_element_of_col_not_zero(x, i, i);
    	if(nz < 0)
    	{
    		post("ambi_decode3 ERROR: matrix not regular !!!!");
    		return;
    	}
    	else
    	{
    		if(nz != i)
    			ambi_decode3_xch_rows(x, i, nz);
    		dv = db + i*n_ambi2 + i;
    		rcp = 1.0 /(*dv);
    		ambi_decode3_mul_row(x, i, rcp);
    		ambi_decode3_copy_row2buf(x, i);
    		for(j=i+1; j<n_ambi; j++)
    		{
    			dv += n_ambi2;
    			rcp = -(*dv);
    			ambi_decode3_mul_buf_and_add2row(x, j, rcp);
    		}
    	}
    }
/* make 0 above the main diagonale */
    for(i=n_ambi-1; i>=0; i--)
    {
    	dv = db + i*n_ambi2 + i;
    	ambi_decode3_copy_row2buf(x, i);
    	for(j=i-1; j>=0; j--)
    	{
    		dv -= n_ambi2;
    		rcp = -(*dv);
    		ambi_decode3_mul_buf_and_add2row(x, j, rcp);
    	}
    }
post("matrix_inverse regular");
}
static void ambi_decode3_begin_pseudo_inverse(t_ambi_decode3 *x)
{
    t_atom *at=x->x_at;
    int i, n=x->x_n_real_ls*x->x_n_ambi;
    double *dv1=x->x_prod;
ambi_decode3_transp_back(x);
    ambi_decode3_mul1(x);
    ambi_decode3_inverse(x);
    ambi_decode3_mul2(x);
    at += 2;
    for(i=0; i<n; i++)
    {
    	SETFLOAT(at, (float)(*dv1));
    	dv1++;
    	at++;
    }
}
static void ambi_decode3_ipht_ireal_muladd(t_ambi_decode3 *x, t_symbol *s, int argc, t_atom *argv)
{
    t_atom *at=x->x_at;
    int i, n=x->x_n_ambi;
    int pht_index, real_index;
    double mw;
    float dat1;
    double *dv2=x->x_prod;
if(argc < 3)
    {
    	post("ambi_decode3 ERROR: ipht_ireal_muladd needs 2 index and 1 mirrorweight: pht_ls_index + real_ls_index + mirror_weight_element");
    	return;
    }
    pht_index = (int)atom_getint(argv++) - 1;
    real_index = (int)atom_getint(argv++) - 1;
    mw = (double)atom_getfloat(argv);
if(pht_index < 0)
    	pht_index = 0;
    if(real_index < 0)
    	real_index = 0;
    if(real_index >= x->x_n_real_ls)
    	real_index = x->x_n_real_ls - 1;
    if(pht_index >= x->x_n_pht_ls)
    	pht_index = x->x_n_pht_ls - 1;
at += 2 + (real_index)*x->x_n_ambi;
    dv2 += (x->x_n_real_ls+pht_index)*x->x_n_ambi;
    for(i=0; i<n; i++)
    {
    	dat1 = atom_getfloat(at);
    	SETFLOAT(at, dat1 + (float)(*dv2*mw));
    	dv2++;
    	at++;
    }
}
static void ambi_decode3_end_pseudo_inverse(t_ambi_decode3 *x)
{
    outlet_anything(x->x_obj.ob_outlet, x->x_s_matrix, x->x_n_ambi*x->x_n_real_ls+2, x->x_at);
}
static void ambi_decode3_encode_ls_2d(t_ambi_decode3 *x, int argc, t_atom *argv, int mode)
{
    double phi;
    double *dw = x->x_transp;
    int index;
    int order=x->x_n_order;
if(argc < 2)
    {
    	post("ambi_decode3 ERROR: ls-input needs 1 index and 1 angle: ls_index + phi [degree]");
    	return;
    }
    index = (int)atom_getint(argv++) - 1;
    phi = (double)atom_getfloat(argv);
if(index < 0)
    	index = 0;
if(mode == AMBI_LS_REAL)
    {
    	if(index >= x->x_n_real_ls)
    		index = x->x_n_real_ls - 1;
    }
    else if(mode == AMBI_LS_PHT)
    {
    	if(x->x_n_pht_ls)
    	{
    		if(index >= x->x_n_pht_ls)
    			index = x->x_n_pht_ls - 1;
    		index += x->x_n_real_ls;
    	}
    	else
    		return;
    }
    else
    	return;
    
    phi *= x->x_pi_over_180;
dw += index * x->x_n_ambi;
*dw++ = 1.0;
    *dw++ = cos(phi);
    *dw++ = sin(phi);
if(order >= 2)
    {
    	*dw++ = cos(2.0*phi);
    	*dw++ = sin(2.0*phi);
if(order >= 3)
    	{
    		*dw++ = cos(3.0*phi);
    		*dw++ = sin(3.0*phi);
    		if(order >= 4)
    		{
    			*dw++ = cos(4.0*phi);
    			*dw++ = sin(4.0*phi);
if(order >= 5)
    			{
    				*dw++ = cos(5.0*phi);
    				*dw++ = sin(5.0*phi);
if(order >= 6)
    				{
    					*dw++ = cos(6.0*phi);
    					*dw++ = sin(6.0*phi);
if(order >= 7)
    					{
    						*dw++ = cos(7.0*phi);
    						*dw++ = sin(7.0*phi);
if(order >= 8)
    						{
    							*dw++ = cos(8.0*phi);
    							*dw++ = sin(8.0*phi);
if(order >= 9)
    							{
    								*dw++ = cos(9.0*phi);
    								*dw++ = sin(9.0*phi);
if(order >= 10)
    								{
    									*dw++ = cos(10.0*phi);
    									*dw++ = sin(10.0*phi);
if(order >= 11)
    									{
    										*dw++ = cos(11.0*phi);
    										*dw++ = sin(11.0*phi);
if(order >= 12)
    										{
    											*dw++ = cos(12.0*phi);
    											*dw++ = sin(12.0*phi);
    										}
    									}
    								}
    							}
    						}
    					}
    				}
    			}
    		}
    	}
    }
}
static void ambi_decode3_encode_ls_3d(t_ambi_decode3 *x, int argc, t_atom *argv, int mode)
{
    double delta, phi;
    double cd, sd, cd2, cd3, sd2, csd, cp, sp, cp2, sp2, cp3, sp3, cp4, sp4;
    double *dw = x->x_transp;
    int index;
    int order=x->x_n_order;
if(argc < 3)
    {
    	post("ambi_decode3 ERROR: ls-input needs 1 index and 2 angles: ls index + delta [degree] + phi [degree]");
    	return;
    }
    index = (int)atom_getint(argv++) - 1;
    delta = atom_getfloat(argv++);
    phi = atom_getfloat(argv);
if(index < 0)
    	index = 0;
    
    if(mode == AMBI_LS_REAL)
    {
    	if(index >= x->x_n_real_ls)
    		index = x->x_n_real_ls - 1;
    }
    else if(mode == AMBI_LS_PHT)
    {
    	if(x->x_n_pht_ls)
    	{
    		if(index >= x->x_n_pht_ls)
    			index = x->x_n_pht_ls - 1;
    		index += x->x_n_real_ls;
    	}
    	else
    		return;
    }
    else
    	return;
delta *= x->x_pi_over_180;
    phi *= x->x_pi_over_180;
dw += index * x->x_n_ambi;
cd = cos(delta);
    sd = sin(delta);
    cp = cos(phi);
    sp = sin(phi);
*dw++ = 1.0;
    *dw++ = cd * cp;
    *dw++ = cd * sp;
    *dw++ = sd;
if(order >= 2)
    {
    	cp2 = cos(2.0*phi);
    	sp2 = sin(2.0*phi);
    	cd2 = cd * cd;
    	sd2 = sd * sd;
    	csd = cd * sd;
    	*dw++ = 0.5 * x->x_sqrt3 * cd2 * cp2;
    	*dw++ = 0.5 * x->x_sqrt3 * cd2 * sp2;
    	*dw++ = x->x_sqrt3 * csd * cp;
    	*dw++ = x->x_sqrt3 * csd * sp;
    	*dw++ = 0.5 * (3.0 * sd2 - 1.0);
if(order >= 3)
    	{
    		cp3 = cos(3.0*phi);
    		sp3 = sin(3.0*phi);
    		cd3 = cd2 * cd;
    		*dw++ = x->x_sqrt10_4 * cd3 * cp3;
    		*dw++ = x->x_sqrt10_4 * cd3 * sp3;
    		*dw++ = x->x_sqrt15_2 * cd * csd * cp2;
    		*dw++ = x->x_sqrt15_2 * cd * csd * sp2;
    		*dw++ = x->x_sqrt6_4 * cd * (5.0 * sd2 - 1.0) * cp;
    		*dw++ = x->x_sqrt6_4 * cd * (5.0 * sd2 - 1.0) * sp;
    		*dw++ = 0.5 * sd * (5.0 * sd2 - 3.0);
if(order >= 4)
    		{
    			cp4 = cos(4.0*phi);
    			sp4 = sin(4.0*phi);
    			*dw++ = x->x_sqrt35_8 * cd2 * cd2 * cp4;
    			*dw++ = x->x_sqrt35_8 * cd2 * cd2 * sp4;
    			*dw++ = x->x_sqrt70_4 * cd2 * csd * cp3;
    			*dw++ = x->x_sqrt70_4 * cd2 * csd * sp3;
    			*dw++ = 0.5 * x->x_sqrt5_2 * cd2 * (7.0 * sd2 - 1.0) * cp2;
    			*dw++ = 0.5 * x->x_sqrt5_2 * cd2 * (7.0 * sd2 - 1.0) * sp2;
    			*dw++ = x->x_sqrt10_4 * csd * (7.0 * sd2 - 3.0) * cp;
    			*dw++ = x->x_sqrt10_4 * csd * (7.0 * sd2 - 3.0) * sp;
    			*dw++ = 0.125 * (sd2 * (35.0 * sd2 - 30.0) + 3.0);
if(order >= 5)
    			{
    				*dw++ = x->x_sqrt126_16 * cd3 * cd2 * cos(5.0*phi);
    				*dw++ = x->x_sqrt126_16 * cd3 * cd2 * sin(5.0*phi);
    				*dw++ = x->x_sqrt315_8 * cd3 * csd * cp4;
    				*dw++ = x->x_sqrt315_8 * cd3 * csd * sp4;
    				*dw++ = 0.25 * x->x_sqrt70_4 * cd3 * (9.0 * sd2 - 1.0) * cp3;
    				*dw++ = 0.25 * x->x_sqrt70_4 * cd3 * (9.0 * sd2 - 1.0) * sp3;
    				*dw++ = x->x_sqrt105_4 * cd * csd * (3.0 * sd2 - 1.0) * cp2;
    				*dw++ = x->x_sqrt105_4 * cd * csd * (3.0 * sd2 - 1.0) * sp2;
    				*dw++ = 0.25 * x->x_sqrt15_2 * cd * (sd2 * (21.0 * sd2 - 14.0) + 1.0) * cp;
    				*dw++ = 0.25 * x->x_sqrt15_2 * cd * (sd2 * (21.0 * sd2 - 14.0) + 1.0) * sp;
    				*dw = 0.125 * sd * (sd2 * (63.0 * sd2 - 70.0) + 15.0);
    			}
    		}
    	}
    }
}
static void ambi_decode3_real_ls(t_ambi_decode3 *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode3_encode_ls_2d(x, argc, argv, AMBI_LS_REAL);
    else
    	ambi_decode3_encode_ls_3d(x, argc, argv, AMBI_LS_REAL);
}
static void ambi_decode3_pht_ls(t_ambi_decode3 *x, t_symbol *s, int argc, t_atom *argv)
{
    if(x->x_n_dim == 2)
    	ambi_decode3_encode_ls_2d(x, argc, argv, AMBI_LS_PHT);
    else
    	ambi_decode3_encode_ls_3d(x, argc, argv, AMBI_LS_PHT);
}
static void ambi_decode3_ambi_weight(t_ambi_decode3 *x, t_symbol *s, int argc, t_atom *argv)
{
    if(argc > x->x_n_order)
    {
    	int i, k=0, n=x->x_n_order;
    	double d;
x->x_ambi_channel_weight[k] = atom_getfloat(argv++);
    	k++;
    	if(x->x_n_dim == 2)
    	{
    		for(i=1; i<=n; i++)
    		{
    			d = atom_getfloat(argv++);
    			x->x_ambi_channel_weight[k] = d;
    			k++;
    			x->x_ambi_channel_weight[k] = d;
    			k++;
    		}
    	}
    	else
    	{
    		int j, m;
for(i=1; i<=n; i++)
    		{
    			d = atom_getfloat(argv++);
    			m = 2*i + 1;
    			for(j=0; j<m; j++)
    			{
    				x->x_ambi_channel_weight[k] = d;
    				k++;
    			}
    		}
    	}
    }
    else
    	post("ambi_decode3-ERROR: ambi_weight needs %d float weights", x->x_n_order+1);
}
static void ambi_decode3_sing_range(t_ambi_decode3 *x, t_floatarg f)
{
    if(f < 0.0f)
    	x->x_sing_range = -(double)f;
    else
    	x->x_sing_range = (double)f;
}
static void ambi_decode3_free(t_ambi_decode3 *x)
{
    freebytes(x->x_inv_work1, x->x_n_ambi * x->x_n_ambi * sizeof(double));
    freebytes(x->x_inv_work2, 2 * x->x_n_ambi * x->x_n_ambi * sizeof(double));
    freebytes(x->x_inv_buf2, 2 * x->x_n_ambi * sizeof(double));
    freebytes(x->x_transp, (x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_ls_encode, (x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_prod, (x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double));
    freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double));
    freebytes(x->x_at, (x->x_n_real_ls * x->x_n_ambi + 2) * sizeof(t_atom));
}
static void *ambi_decode3_new(t_symbol *s, int argc, t_atom *argv)
{
    t_ambi_decode3 *x = (t_ambi_decode3 *)pd_new(ambi_decode3_class);
    int order, dim, i;
    int n_real_ls=0;/* number of loudspeakers */
    int n_pht_ls=0;/* number of phantom_loudspeakers */
if((argc >= 4) &&
    	IS_A_FLOAT(argv,0) &&
    	IS_A_FLOAT(argv,1) &&
    	IS_A_FLOAT(argv,2) &&
    	IS_A_FLOAT(argv,3))
    {
    	order = (int)atom_getint(argv++);
    	dim = (int)atom_getint(argv++);
    	n_real_ls = (int)atom_getint(argv++);
    	n_pht_ls = (int)atom_getint(argv);
if(order < 1)
    		order = 1;
    	if(dim != 3)
    	{
    		dim = 2;
    		if(order > 12)
    			order = 12;
    		x->x_n_ambi = 2*order + 1;
    	}
    	else
    	{
    		if(order > 5)
    			order = 5;
    		x->x_n_ambi = (order + 1)*(order + 1);
    	}
    	x->x_n_dim = dim;
    	x->x_n_order = order;
    	if(n_real_ls < 1)
    		n_real_ls = 1;
    	if(n_pht_ls < 0)
    		n_pht_ls = 0;
    	if((n_real_ls + n_pht_ls) < x->x_n_ambi)
    		post("ambi_decode3-WARNING: Number of Loudspeakers < Number of Ambisonic-Channels !!!!");
    	
    	x->x_n_real_ls = n_real_ls;
    	x->x_n_pht_ls = n_pht_ls;
    	x->x_inv_work1 = (double *)getbytes(x->x_n_ambi * x->x_n_ambi * sizeof(double));
    	x->x_inv_work2 = (double *)getbytes(2 * x->x_n_ambi * x->x_n_ambi * sizeof(double));
    	x->x_inv_buf2 = (double *)getbytes(2 * x->x_n_ambi * sizeof(double));
    	x->x_transp = (double *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double));
    	x->x_ls_encode = (double *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double));
    	x->x_prod = (double *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double));
    	x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double));
    	x->x_at = (t_atom *)getbytes((x->x_n_real_ls * x->x_n_ambi + 2) * sizeof(t_atom));
    	x->x_s_matrix = gensym("matrix");
    	/*change*/
    	SETFLOAT(x->x_at, (float)x->x_n_real_ls);
    	SETFLOAT(x->x_at+1, (float)x->x_n_ambi);
x->x_sqrt3				= sqrt(3.0);
    	x->x_sqrt5_2			= sqrt(5.0) / 2.0;
    	x->x_sqrt6_4			= sqrt(6.0) / 4.0;
    	x->x_sqrt10_4			= sqrt(10.0) / 4.0;
    	x->x_sqrt15_2			= sqrt(15.0) / 2.0;
    	x->x_sqrt35_8			= sqrt(35.0) / 8.0;
    	x->x_sqrt70_4			= sqrt(70.0) / 4.0;
    	x->x_sqrt126_16		= sqrt(126.0) / 16.0;
    	x->x_sqrt315_8		= sqrt(315.0) / 8.0;
    	x->x_sqrt105_4		= sqrt(105.0) / 4.0;
    	x->x_pi_over_180	= 4.0 * atan(1.0) / 180.0;
    	x->x_sing_range = 1.0e-10;
    	for(i=0; i<x->x_n_ambi; i++)
    		x->x_ambi_channel_weight[i] = 1.0;
    	outlet_new(&x->x_obj, &s_list);
    	return (x);
    }
    else
    {
    	post("ambi_decode3-ERROR: need 4 float arguments: ambi_order dimension number_of_real_loudspeakers number_of_canceled_phantom_speakers");
    	return(0);
    }
}
void ambi_decode3_setup(void)
{
    ambi_decode3_class = class_new(gensym("ambi_decode3"), (t_newmethod)ambi_decode3_new, (t_method)ambi_decode3_free,
    			   sizeof(t_ambi_decode3), 0, A_GIMME, 0);
    class_addmethod(ambi_decode3_class, (t_method)ambi_decode3_real_ls, gensym("real_ls"), A_GIMME, 0);
    class_addmethod(ambi_decode3_class, (t_method)ambi_decode3_pht_ls, gensym("pht_ls"), A_GIMME, 0);
    class_addmethod(ambi_decode3_class, (t_method)ambi_decode3_ambi_weight, gensym("ambi_weight"), A_GIMME, 0);
    class_addmethod(ambi_decode3_class, (t_method)ambi_decode3_sing_range, gensym("sing_range"), A_DEFFLOAT, 0);
    class_addmethod(ambi_decode3_class, (t_method)ambi_decode3_begin_pseudo_inverse, gensym("begin_pseudo_inverse"), 0);
    class_addmethod(ambi_decode3_class, (t_method)ambi_decode3_ipht_ireal_muladd, gensym("ipht_ireal_muladd"), A_GIMME, 0);
    class_addmethod(ambi_decode3_class, (t_method)ambi_decode3_end_pseudo_inverse, gensym("end_pseudo_inverse"), 0);
    class_sethelpsymbol(ambi_decode3_class, gensym("iemhelp2/help-ambi_decode3"));
}