[PD-cvs] 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

[musil]

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 at 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)- at 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"));
}