diff options
author | musil <tmusil@users.sourceforge.net> | 2006-03-09 16:00:50 +0000 |
---|---|---|
committer | musil <tmusil@users.sourceforge.net> | 2006-03-09 16:00:50 +0000 |
commit | d2a273ebc92f3d1336212f01891a070bf4a31b41 (patch) | |
tree | 1298a5648425bf96cf7d1c8c0827af436eaa50e3 /src | |
parent | 7457e39d10cb100c07f295211f07c45da148cc80 (diff) |
update to version JAN_06 with all sources
svn path=/trunk/externals/iem/iem_bin_ambi/; revision=4666
Diffstat (limited to 'src')
-rw-r--r-- | src/bin_ambi_calc_HRTF.c | 4 | ||||
-rw-r--r-- | src/bin_ambi_reduced_decode.c | 71 | ||||
-rw-r--r-- | src/bin_ambi_reduced_decode2.c | 1360 | ||||
-rw-r--r-- | src/bin_ambi_reduced_decode_fft.c | 1512 | ||||
-rw-r--r-- | src/bin_ambi_reduced_decode_fft2.c | 1388 | ||||
-rw-r--r-- | src/bin_ambi_reduced_decode_fir.c | 1385 | ||||
-rw-r--r-- | src/bin_ambi_reduced_decode_fir2.c | 1260 | ||||
-rw-r--r-- | src/iem_bin_ambi.c | 24 | ||||
-rw-r--r-- | src/iem_bin_ambi.h | 8 | ||||
-rw-r--r-- | src/iemlib.h | 45 | ||||
-rw-r--r-- | src/makefile_win | 82 |
11 files changed, 7061 insertions, 78 deletions
diff --git a/src/bin_ambi_calc_HRTF.c b/src/bin_ambi_calc_HRTF.c index 57460b3..d9e8ee3 100644 --- a/src/bin_ambi_calc_HRTF.c +++ b/src/bin_ambi_calc_HRTF.c @@ -1,7 +1,7 @@ /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution. -iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2003 */ +iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */ #ifdef NT #pragma warning( disable : 4244 ) @@ -69,7 +69,7 @@ static void bin_ambi_calc_HRTF_init_cos(t_bin_ambi_calc_HRTF *x) { f = g * (float)i; (*sincos).real = cos(f); - (*sincos).imag = sin(f);/*change*/ + (*sincos).imag = -sin(f);/*FFT*/ sincos++; } } diff --git a/src/bin_ambi_reduced_decode.c b/src/bin_ambi_reduced_decode.c index 8626473..bce4f83 100644 --- a/src/bin_ambi_reduced_decode.c +++ b/src/bin_ambi_reduced_decode.c @@ -1,7 +1,7 @@ /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution. -iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2003 */ +iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */ #ifdef NT #pragma warning( disable : 4244 ) @@ -88,7 +88,7 @@ static void bin_ambi_reduced_decode_init_cos(t_bin_ambi_reduced_decode *x) { f = g * (float)i; (*sincos).real = cos(f); - (*sincos).imag = sin(f);/*change*/ + (*sincos).imag = -sin(f);/*FFT*/ sincos++; } } @@ -741,7 +741,7 @@ static void bin_ambi_reduced_decode_check_arrays(t_bin_ambi_reduced_decode *x, f if(npoints < fftsize) { - post("warning: %s-array-size: %d", hrir->s_name, npoints); + post("bin_ambi_reduced_decode-WARNING: %s-array-size: %d < FFT-size: %d", hrir->s_name, npoints, fftsize); } vec = x->x_beg_hrir; vec += index * fftsize; @@ -858,16 +858,18 @@ static void bin_ambi_reduced_decode_calc_reduced(t_bin_ambi_reduced_decode *x, f val[i] = old1; } } - - vec_hrtf_re[0] = val[0].real; - vec_hrtf_im[0] = 0.0f; - for( i = 1; i < fs2; i++ ) + for(i = 0; i<fs2; i++) { - vec_hrtf_re[i] = 2.0f*val[i].real; - vec_hrtf_im[i] = 2.0f*val[i].imag; + vec_hrtf_re[i] = val[i].real; + vec_hrtf_im[i] = val[i].imag; } - vec_hrtf_re[fs2] = 0.0f; + vec_hrtf_re[fs2] = val[fs2].real; vec_hrtf_im[fs2] = 0.0f; + for(i = fs2+1; i < fftsize; i++) + { + vec_hrtf_re[i] = 0.0f; + vec_hrtf_im[i] = 0.0f; + } } static void bin_ambi_reduced_decode_decoding(t_bin_ambi_reduced_decode *x) @@ -1031,7 +1033,7 @@ static void bin_ambi_reduced_decode_ambi_weight(t_bin_ambi_reduced_decode *x, t_ } } else - post("bin_ambi_reduced_decode-ERROR: ambi_weight needs %d float weights", x->x_n_order); + post("bin_ambi_reduced_decode-ERROR: ambi_weight needs %d float weights", x->x_n_order+1); } static void bin_ambi_reduced_decode_sing_range(t_bin_ambi_reduced_decode *x, t_floatarg f) @@ -1081,12 +1083,12 @@ static void *bin_ambi_reduced_decode_new(t_symbol *s, int argc, t_atom *argv) { t_bin_ambi_reduced_decode *x = (t_bin_ambi_reduced_decode *)pd_new(bin_ambi_reduced_decode_class); char buf[400]; - int i, j, fftok; - int n_order, n_dim, n_ambi, fftsize, prefix; - t_symbol *s_hrir; - t_symbol *s_hrtf_re; - t_symbol *s_hrtf_im; - t_symbol *s_fade_out_hrir; + int i, j, fftok, ok=0; + int n_order=0, n_dim=0, n_ambi=0, fftsize=0, prefix=0; + t_symbol *s_hrir=gensym("L_HRIR"); + t_symbol *s_hrtf_re=gensym("HRTF_re"); + t_symbol *s_hrtf_im=gensym("HRTF_im"); + t_symbol *s_fade_out_hrir=gensym("HRIR_win"); if((argc >= 8) && IS_A_FLOAT(argv,0) && @@ -1109,6 +1111,34 @@ static void *bin_ambi_reduced_decode_new(t_symbol *s, int argc, t_atom *argv) n_dim = (int)atom_getintarg(6, argc, argv); fftsize = (int)atom_getintarg(7, argc, argv); + ok = 1; + } + else if((argc >= 8) && + 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) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = gensym("L_HRIR"); + s_hrtf_re = gensym("HRTF_re"); + s_hrtf_im = gensym("HRTF_im"); + s_fade_out_hrir = gensym("HRIR_win"); + + n_order = (int)atom_getintarg(5, argc, argv); + n_dim = (int)atom_getintarg(6, argc, argv); + fftsize = (int)atom_getintarg(7, argc, argv); + + ok = 1; + } + + if(ok) + { if(n_order < 1) n_order = 1; @@ -1141,7 +1171,8 @@ static void *bin_ambi_reduced_decode_new(t_symbol *s, int argc, t_atom *argv) if(!fftok) { fftsize = 512; - post("bin_ambi_reduced_decode-ERROR: fftsize not equal to 2 ^ n !!!"); + post("bin_ambi_reduced_decode-WARNING: fftsize not equal to 2 ^ n !!!"); + post(" fftsize set to %d", fftsize); } x->x_n_dim = n_dim; @@ -1211,8 +1242,8 @@ static void *bin_ambi_reduced_decode_new(t_symbol *s, int argc, t_atom *argv) else { post("bin_ambi_reduced_decode-ERROR: need 1 float + 4 symbols + 3 floats arguments:"); - post(" prefix + hrir_name + hrtf_re_name + hrtf_im_name + hrir_fade_out_name +"); - post(" ambi_order + ambi_dimension + fftsize"); + post(" prefix(unique-number) + hrir_name + hrtf_re_name + hrtf_im_name + hrir_fade_out_name +"); + post(" + ambi_order + ambi_dimension + fftsize"); return(0); } } diff --git a/src/bin_ambi_reduced_decode2.c b/src/bin_ambi_reduced_decode2.c new file mode 100644 index 0000000..3b694fd --- /dev/null +++ b/src/bin_ambi_reduced_decode2.c @@ -0,0 +1,1360 @@ +/* For information on usage and redistribution, and for a DISCLAIMER OF ALL +* WARRANTIES, see the file, "LICENSE.txt," in this distribution. + +iem_bin_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_bin_ambi.h" +#include <math.h> +#include <stdio.h> +#include <string.h> + + +/* -------------------------- bin_ambi_reduced_decode2 ------------------------------ */ +/* + ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren ** + ** Inputs: ls + Liste von 3 floats: Index [1 .. 16] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] ** + ** Inputs: calc_inv ** + ** Inputs: load_HRIR + float index1..16 ** + ** 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 1016_1_L_HRIR .. 1016_16_L_HRIR heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_re .. 1016_16_HRTF_re heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_im .. 1016_16_HRTF_im heissen und existieren ** + */ + +typedef struct _bin_ambi_reduced_decode2 +{ + t_object x_obj; + t_atom x_at[2];/*output filename and tablename of HRIR*/ + int x_n_dim;/*dimension of ambisonic system*/ + int x_n_ambi;/*number of ambi channels*/ + int x_n_order;/*order of ambisonic system*/ + int x_n_ls;/*number of loudspeakers*/ + int x_seq_ok; + int x_fftsize;/*fftsize of blockfilter, should be 2 times of FIR-length of HRIR*/ + double *x_inv_work1;/* n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_work2;/* 2*n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_buf2;/* 2*n_ambi buffer for matrix inverse */ + double *x_transp;/* n_ls*n_ambi transposed input-buffer of decoder-matrix before inversion */ + double *x_ls_encode;/* n_ambi*n_ls straight input-buffer of decoder-matrix before inversion */ + double *x_prod;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_ambi_channel_weight; + int *x_delta; + int *x_phi; + int *x_phi_sym; + int *x_sym_flag; + BIN_AMBI_COMPLEX *x_spec; + BIN_AMBI_COMPLEX *x_sin_cos; + float *x_beg_fade_out_hrir; + float *x_beg_hrir; + float **x_beg_hrtf_re; + float **x_beg_hrtf_im; + t_symbol **x_hrir_filename; + t_symbol **x_s_hrir; + t_symbol **x_s_hrtf_re; + t_symbol **x_s_hrtf_im; + t_symbol *x_s_fade_out_hrir; + void *x_out_sign_sum; + 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; + double x_sing_range; +} t_bin_ambi_reduced_decode2; + +static t_class *bin_ambi_reduced_decode2_class; + +static void bin_ambi_reduced_decode2_init_cos(t_bin_ambi_reduced_decode2 *x) +{/* initialize a whole sine and cosine wave over a fftsize array */ + int i, fftsize = x->x_fftsize; + float f, g; + BIN_AMBI_COMPLEX *sincos = x->x_sin_cos; + + g = 2.0f * 3.1415926538f / (float)fftsize; + for(i=0; i<fftsize; i++) + { + f = g * (float)i; + (*sincos).real = cos(f); + (*sincos).imag = -sin(f); /*FFT*/ + sincos++; + } +} + +static void bin_ambi_reduced_decode2_quant(t_bin_ambi_reduced_decode2 *x, double *delta_deg2rad, double *phi_deg2rad, int index) +{ + double q = 1.0; + double d = *delta_deg2rad; + double p = *phi_deg2rad; + int i; + + if(d < -40.0) + d = -40.0; + if(d > 90.0) + d = 90.0; + while(p < 0.0) + p += 360.0; + while(p >= 360.0) + p -= 360.0; + + if(d < -35.0) + { + *delta_deg2rad = -40.0; + q = 360.0 / 56.0; + } + else if(d < -25.0) + { + *delta_deg2rad = -30.0; + q = 6.0; + } + else if(d < -15.0) + { + *delta_deg2rad = -20.0; + q = 5.0; + } + else if(d < -5.0) + { + *delta_deg2rad = -10.0; + q = 5.0; + } + else if(d < 5.0) + { + *delta_deg2rad = 0.0; + q = 5.0; + } + else if(d < 15.0) + { + *delta_deg2rad = 10.0; + q = 5.0; + } + else if(d < 25.0) + { + *delta_deg2rad = 20.0; + q = 5.0; + } + else if(d < 35.0) + { + *delta_deg2rad = 30.0; + q = 6.0; + } + else if(d < 45.0) + { + *delta_deg2rad = 40.0; + q = 360.0 / 56.0; + } + else if(d < 55.0) + { + *delta_deg2rad = 50.0; + q = 8.0; + } + else if(d < 65.0) + { + *delta_deg2rad = 60.0; + q = 10.0; + } + else if(d < 75.0) + { + *delta_deg2rad = 70.0; + q = 15.0; + } + else if(d < 85.0) + { + *delta_deg2rad = 80.0; + q = 30.0; + } + else + { + *delta_deg2rad = 90.0; + q = 360.0; + } + + p /= q; + i = (int)(p + 0.499999); + p = (double)i * q; + i = (int)(p + 0.499999); + while(i >= 360) + i -= 360; + p = (double)i; + *phi_deg2rad = p; + + x->x_delta[index] = (int)(*delta_deg2rad); + x->x_phi[index] = i; + + *delta_deg2rad *= x->x_pi_over_180; + *phi_deg2rad *= x->x_pi_over_180; +} + +static void bin_ambi_reduced_decode2_do_2d(t_bin_ambi_reduced_decode2 *x, int argc, t_atom *argv) +{ + double delta=0.0, phi; + double *dw = x->x_transp; + int index; + int order=x->x_n_order; + int n_ls = x->x_n_ls; + + if(argc < 2) + { + post("bin_ambi_reduced_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(index >= n_ls) + index = n_ls - 1; + + bin_ambi_reduced_decode2_quant(x, &delta, &phi, index); + + 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 bin_ambi_reduced_decode2_do_3d(t_bin_ambi_reduced_decode2 *x, int argc, t_atom *argv) +{ + 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 order=x->x_n_order; + + if(argc < 3) + { + post("bin_ambi_reduced_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(index >= n_ls) + index = n_ls - 1; + + bin_ambi_reduced_decode2_quant(x, &delta, &phi, index); + + cd = cos(delta); + sd = sin(delta); + cp = cos(phi); + sp = sin(phi); + + dw += index * x->x_n_ambi; + + *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 bin_ambi_reduced_decode2_ls(t_bin_ambi_reduced_decode2 *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode2_do_2d(x, argc, argv); + else + bin_ambi_reduced_decode2_do_3d(x, argc, argv); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode2_copy_row2buf(t_bin_ambi_reduced_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 bin_ambi_reduced_decode2_copy_buf2row(t_bin_ambi_reduced_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 bin_ambi_reduced_decode2_copy_row2row(t_bin_ambi_reduced_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 bin_ambi_reduced_decode2_xch_rows(t_bin_ambi_reduced_decode2 *x, int row1, int row2) +{ + bin_ambi_reduced_decode2_copy_row2buf(x, row1); + bin_ambi_reduced_decode2_copy_row2row(x, row2, row1); + bin_ambi_reduced_decode2_copy_buf2row(x, row2); +} + +static void bin_ambi_reduced_decode2_mul_row(t_bin_ambi_reduced_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 bin_ambi_reduced_decode2_mul_col(t_bin_ambi_reduced_decode2 *x, int col, double mul) +{ + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int i; + double *dw=x->x_inv_work2; + + dw += col; + for(i=0; i<n_ambi; i++) + { + (*dw) *= mul; + dw += n_ambi2; + } +} + +static void bin_ambi_reduced_decode2_mul_buf_and_add2row(t_bin_ambi_reduced_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 bin_ambi_reduced_decode2_eval_which_element_of_col_not_zero(t_bin_ambi_reduced_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 bin_ambi_reduced_decode2_mul1(t_bin_ambi_reduced_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; + 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 bin_ambi_reduced_decode2_mul2(t_bin_ambi_reduced_decode2 *x) +{ + int n_ls=x->x_n_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 bin_ambi_reduced_decode2_transp_back(t_bin_ambi_reduced_decode2 *x) +{ + double *vec, *transp=x->x_transp; + double *straight=x->x_ls_encode; + int n_ls=x->x_n_ls; + int n_ambi=x->x_n_ambi; + int i, j; + + transp = x->x_transp; + for(j=0; j<n_ambi; j++) + { + vec = transp; + for(i=0; i<n_ls; i++) + { + *straight++ = *vec; + vec += n_ambi; + } + transp++; + } +} + +static void bin_ambi_reduced_decode2_inverse(t_bin_ambi_reduced_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 = bin_ambi_reduced_decode2_eval_which_element_of_col_not_zero(x, i, i); + if(nz < 0) + { + post("bin_ambi_reduced_decode2 ERROR: matrix not regular !!!!"); + x->x_seq_ok = 0; + return; + } + else + { + if(nz != i) + bin_ambi_reduced_decode2_xch_rows(x, i, nz); + dv = db + i*n_ambi2 + i; + rcp = 1.0 /(*dv); + bin_ambi_reduced_decode2_mul_row(x, i, rcp); + bin_ambi_reduced_decode2_copy_row2buf(x, i); + for(j=i+1; j<n_ambi; j++) + { + dv += n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_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; + bin_ambi_reduced_decode2_copy_row2buf(x, i); + for(j=i-1; j>=0; j--) + { + dv -= n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode2_mul_buf_and_add2row(x, j, rcp); + } + } + + post("matrix_inverse regular"); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode2_calc_pinv(t_bin_ambi_reduced_decode2 *x) +{ + t_garray *a; + int npoints; + t_float *fadevec; + + if((int)(x->x_beg_fade_out_hrir) == 0) + { + if (!(a = (t_garray *)pd_findbyclass(x->x_s_fade_out_hrir, garray_class))) + error("%s: no such array", x->x_s_fade_out_hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &fadevec)) + error("%s: bad template for bin_ambi_reduced_decode2", x->x_s_fade_out_hrir->s_name); + else if (npoints < x->x_fftsize) + error("%s: bad array-size: %d", x->x_s_fade_out_hrir->s_name, npoints); + else + x->x_beg_fade_out_hrir = fadevec; + } + + bin_ambi_reduced_decode2_transp_back(x); + bin_ambi_reduced_decode2_mul1(x); + bin_ambi_reduced_decode2_inverse(x); + bin_ambi_reduced_decode2_mul2(x); +} + +/* +x_prod: +n_ambi columns; +n_ambi rows; +*/ + +static void bin_ambi_reduced_decode2_load_HRIR(t_bin_ambi_reduced_decode2 *x, float findex) +{ + int index=(int)findex - 1; + int p; + char buf[60]; + + if(index < 0) + index = 0; + if(index >= x->x_n_ls) + index = x->x_n_ls - 1; + + p = x->x_phi[index]; + + if(p)/*change*/ + p = 360 - p; + + if(p < 10) + sprintf(buf, "L%de00%da.wav", x->x_delta[index], p); + else if(p < 100) + sprintf(buf, "L%de0%da.wav", x->x_delta[index], p); + else + sprintf(buf, "L%de%da.wav", x->x_delta[index], p); + x->x_hrir_filename[index] = gensym(buf); + + SETSYMBOL(x->x_at, x->x_hrir_filename[index]); + SETSYMBOL(x->x_at+1, x->x_s_hrir[index]); + outlet_list(x->x_obj.ob_outlet, &s_list, 2, x->x_at); +} + +static void bin_ambi_reduced_decode2_check_HRIR_arrays(t_bin_ambi_reduced_decode2 *x, float findex) +{ + int index=(int)findex - 1; + int j, k, n; + int fftsize = x->x_fftsize; + int fs2=fftsize/2; + t_garray *a; + int npoints; + t_symbol *hrir; + t_float *vec_hrir, *vec, *vec_fade_out_hrir; + float decr, sum; + + if(index < 0) + index = 0; + if(index >= x->x_n_ls) + index = x->x_n_ls - 1; + + hrir = x->x_s_hrir[index]; + if (!(a = (t_garray *)pd_findbyclass(hrir, garray_class))) + error("%s: no such array", hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrir)) + error("%s: bad template for bin_ambi_reduced_decode2", hrir->s_name); + else + { + if(npoints < fftsize) + { + post("bin_ambi_reduced_decode2-WARNING: %s-array-size: %d < FFT-size: %d", hrir->s_name, npoints, fftsize); + } + vec = x->x_beg_hrir; + vec += index * fftsize; + + if((int)(x->x_beg_fade_out_hrir)) + { + vec_fade_out_hrir = x->x_beg_fade_out_hrir; + for(j=0; j<fs2; j++) + vec[j] = vec_hrir[j]*vec_fade_out_hrir[j]; + } + else + { + post("no HRIR-fade-out-window found"); + n = fs2 * 3; + n /= 4; + for(j=0; j<n; j++) + vec[j] = vec_hrir[j]; + sum = 1.0f; + decr = 4.0f / (float)fs2; + for(j=n, k=0; j<fs2; j++, k++) + { + sum -= decr; + vec[j] = vec_hrir[j] * sum; + } + } + } +} + +static void bin_ambi_reduced_decode2_check_HRTF_arrays(t_bin_ambi_reduced_decode2 *x, float findex) +{ + int index=(int)findex - 1; + t_garray *a; + int npoints; + int fftsize = x->x_fftsize; + t_float *vec_hrtf_re, *vec_hrtf_im; + t_symbol *hrtf_re, *hrtf_im; + + if(index < 0) + index = 0; + if(index >= x->x_n_ambi) + index = x->x_n_ambi - 1; + + hrtf_re = x->x_s_hrtf_re[index]; + hrtf_im = x->x_s_hrtf_im[index]; + + if (!(a = (t_garray *)pd_findbyclass(hrtf_re, garray_class))) + error("%s: no such array", hrtf_re->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrtf_re)) + error("%s: bad template for bin_ambi_reduced_decode2", hrtf_re->s_name); + else if (npoints < fftsize) + error("%s: bad array-size: %d", hrtf_re->s_name, npoints); + else if (!(a = (t_garray *)pd_findbyclass(hrtf_im, garray_class))) + error("%s: no such array", hrtf_im->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrtf_im)) + error("%s: bad template for bin_ambi_reduced_decode2", hrtf_im->s_name); + else if (npoints < fftsize) + error("%s: bad array-size: %d", hrtf_im->s_name, npoints); + else + { + x->x_beg_hrtf_re[index] = vec_hrtf_re; + x->x_beg_hrtf_im[index] = vec_hrtf_im; + } +} + +static void bin_ambi_reduced_decode2_calc_reduced(t_bin_ambi_reduced_decode2 *x, float findex) +{ + int index=(int)findex - 1; + int i, j, k, w_index, w_inc, i_inc, v_index, fs1, fs2; + int fftsize = x->x_fftsize; + BIN_AMBI_COMPLEX old1, old2, w; + BIN_AMBI_COMPLEX *sincos = x->x_sin_cos; + BIN_AMBI_COMPLEX *val = x->x_spec;/*weighted decoder-matrix with n_ls lines and n_ambi columns*/ + t_float *vec_hrir, *vec_hrtf_re, *vec_hrtf_im; + double *dv, *db=x->x_prod; + int n_ambi = x->x_n_ambi; + int n_ls = x->x_n_ls; + float mul; + + if(x->x_seq_ok) + { + if(index < 0) + index = 0; + if(index >= n_ambi) + index = n_ambi - 1; + + vec_hrtf_re = x->x_beg_hrtf_re[index]; + vec_hrtf_im = x->x_beg_hrtf_im[index]; + + dv = db + index; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + for(k=0; k<fftsize; k++)/*first step of acumulating the HRIRs*/ + { + val[k].real = mul * vec_hrir[k]; + val[k].imag = 0.0f; + } + + for(j=1; j<n_ls; j++) + { + dv += n_ambi; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + vec_hrir += j * fftsize; + for(k=0; k<fftsize; k++) + { + val[k].real += mul * vec_hrir[k]; + } + } + + fs1 = fftsize - 1; + fs2 = fftsize / 2; + i_inc = fs2; + w_inc = 1; + for(i=1; i<fftsize; i<<=1) + { + v_index = 0; + for(j=0; j<i; j++) + { + w_index = 0; + for(k=0; k<i_inc; k++) + { + old1 = val[v_index]; + old2 = val[v_index+i_inc]; + w = sincos[w_index]; + val[v_index+i_inc].real = (old1.real-old2.real)*w.real - (old1.imag-old2.imag)*w.imag; + val[v_index+i_inc].imag = (old1.imag - old2.imag)*w.real + (old1.real-old2.real)*w.imag; + val[v_index].real = old1.real+old2.real; + val[v_index].imag = old1.imag+old2.imag; + w_index += w_inc; + v_index++; + } + v_index += i_inc; + } + w_inc <<= 1; + i_inc >>= 1; + } + + j = 0; + for(i=1;i<fs1;i++) + { + k = fs2; + while(k <= j) + { + j = j - k; + k >>= 1; + } + j = j + k; + if(i < j) + { + old1 = val[j]; + val[j] = val[i]; + val[i] = old1; + } + } + for(i = 0; i<fs2; i++) + { + vec_hrtf_re[i] = val[i].real; + vec_hrtf_im[i] = val[i].imag; + } + vec_hrtf_re[fs2] = val[fs2].real; + vec_hrtf_im[fs2] = 0.0f; + for(i = fs2+1; i < fftsize; i++) + { + vec_hrtf_re[i] = 0.0f; + vec_hrtf_im[i] = 0.0f; + } + } +} + +static void bin_ambi_reduced_decode2_decoding(t_bin_ambi_reduced_decode2 *x) +{ + int *sym=x->x_phi_sym; + int n_ls = x->x_n_ls; + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int *phi=x->x_phi; + int *delta=x->x_delta; + int *flag=x->x_sym_flag; + int i, j, d, p, regular=1; + double *dv, *db=x->x_prod; + double a, b, q; + int sym_max=0; + int pos_sym_counter, neg_sym_counter; + char plus_minus[100]; + char abc[100]; + char onenine[100]; + char ten[100]; + + if(x->x_seq_ok) + { + for(i=0; i<n_ls; i++) + sym[i] = -1; + + plus_minus[0] = 0; + if(x->x_n_dim == 2) + { + strcpy(abc, " WXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXY"); + strcpy(ten, " 000000000000000000011111111111111111111"); + strcpy(onenine, " 011223344556677889900112233445566778899"); + abc[n_ambi+22] = 0; + ten[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(ten); + post(onenine); + } + else + { + strcpy(abc, " AABCABCDEABCDEFGABCDEFGHIABCDEFGHIJKABCDEFGHIJKLM"); + strcpy(onenine, " 0111222223333333444444444555555555556666666666666"); + abc[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(onenine); + } + + for(i=0; i<n_ls; i++)/* looking for azimuth symmetries of loudspeakers */ + { + d = delta[i]; + p = phi[i]; + for(j=i+1; j<n_ls; j++) + { + if(d == delta[j]) + { + if((p + phi[j]) == 360) + { + if((sym[i] < 0) && (sym[j] < 0)) + { + sym[i] = j; + sym[j] = i; + j = n_ls + 1; + sym_max++; + } + } + } + } + } + + for(p=0; p<n_ambi; p++)/*each col*/ + { + pos_sym_counter = 0; + neg_sym_counter = 0; + for(i=0; i<n_ls; i++) + flag[i] = 1; + dv = db + p; + for(i=0; i<n_ls; i++) + { + if((sym[i] >= 0) && flag[i]) + { + j = sym[i]; + flag[i] = 0; + flag[j] = 0; + a = dv[i*n_ambi]; + b = dv[j*n_ambi]; + if((a < 5.0e-4)&&(a > -5.0e-4)&&(b < 5.0e-4)&&(b > -5.0e-4)) + { + pos_sym_counter++; + neg_sym_counter++; + } + else + { + q = a / b; + if((q < 1.005)&&(q > 0.995)) + pos_sym_counter++; + else if((q > -1.005)&&(q < -0.995)) + neg_sym_counter++; + } + } + } + if(pos_sym_counter == sym_max) + strcat(plus_minus, "+"); + else if(neg_sym_counter == sym_max) + strcat(plus_minus, "-"); + else + { + strcat(plus_minus, "?"); + regular = 0; + } + } + post("sum of right channel: %s", plus_minus); + + if(regular) + { + for(i=0; i<n_ambi; i++) + { + /*change*/ + if(plus_minus[i] == '+') + SETFLOAT(x->x_at, 1.0f); + else if(plus_minus[i] == '-') + SETFLOAT(x->x_at, 2.0f); + SETFLOAT(x->x_at+1, (float)(i+1)); + outlet_list(x->x_out_sign_sum, &s_list, 2, x->x_at); + } + } + } +} + +static void bin_ambi_reduced_decode2_ambi_weight(t_bin_ambi_reduced_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] = (double)atom_getfloat(argv++); + k++; + if(x->x_n_dim == 2) + { + for(i=1; i<=n; i++) + { + d = (double)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 = (double)atom_getfloat(argv++); + m = 2*i + 1; + for(j=0; j<m; j++) + { + x->x_ambi_channel_weight[k] = d; + k++; + } + } + } + } + else + post("bin_ambi_reduced_decode2-ERROR: ambi_weight needs %d float weights", x->x_n_order+1); +} + +static void bin_ambi_reduced_decode2_sing_range(t_bin_ambi_reduced_decode2 *x, t_floatarg f) +{ + if(f < 0.0f) + x->x_sing_range = -(double)f; + else + x->x_sing_range = (double)f; +} + +static void bin_ambi_reduced_decode2_free(t_bin_ambi_reduced_decode2 *x) +{ + freebytes(x->x_hrir_filename, x->x_n_ls * sizeof(t_symbol *)); + freebytes(x->x_s_hrir, x->x_n_ls * sizeof(t_symbol *)); + freebytes(x->x_s_hrtf_re, x->x_n_ambi * sizeof(t_symbol *)); + freebytes(x->x_s_hrtf_im, x->x_n_ambi * sizeof(t_symbol *)); + 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_ambi * sizeof(double)); + freebytes(x->x_ls_encode, x->x_n_ls * x->x_n_ambi * sizeof(double)); + freebytes(x->x_prod, x->x_n_ls * x->x_n_ambi * sizeof(double)); + freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double)); + + freebytes(x->x_delta, x->x_n_ls * sizeof(int)); + freebytes(x->x_phi, x->x_n_ls * sizeof(int)); + freebytes(x->x_phi_sym, x->x_n_ls * sizeof(int)); + freebytes(x->x_sym_flag, x->x_n_ls * sizeof(int)); + + freebytes(x->x_spec, x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + freebytes(x->x_sin_cos, x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + + freebytes(x->x_beg_hrir, x->x_fftsize * x->x_n_ls * sizeof(float)); + freebytes(x->x_beg_hrtf_re, x->x_n_ambi * sizeof(float *)); + freebytes(x->x_beg_hrtf_im, x->x_n_ambi * sizeof(float *)); +} + +/* + 1.arg_ int prefix; + 2.arg: t_symbol *hrir_name; + 3.arg: t_symbol *hrtf_re_name; + 4.arg: t_symbol *hrtf_im_name; + 5.arg: t_symbol *hrir_fade_out_name; + 6.arg: int ambi_order; + 7.arg: int dim; + 8.arg: int number_of_loudspeakers; + 9.arg: int fftsize; + */ + +static void *bin_ambi_reduced_decode2_new(t_symbol *s, int argc, t_atom *argv) +{ + t_bin_ambi_reduced_decode2 *x = (t_bin_ambi_reduced_decode2 *)pd_new(bin_ambi_reduced_decode2_class); + char buf[400]; + int i, j, fftok, ok=0; + int n_order=0, n_dim=0, n_ls=0, n_ambi=0, fftsize=0, prefix=0; + t_symbol *s_hrir=gensym("L_HRIR"); + t_symbol *s_hrtf_re=gensym("HRTF_re"); + t_symbol *s_hrtf_im=gensym("HRTF_im"); + t_symbol *s_fade_out_hrir=gensym("HRIR_win"); + + if((argc >= 9) && + IS_A_FLOAT(argv,0) && + IS_A_SYMBOL(argv,1) && + IS_A_SYMBOL(argv,2) && + IS_A_SYMBOL(argv,3) && + IS_A_SYMBOL(argv,4) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = (t_symbol *)atom_getsymbolarg(1, argc, argv); + s_hrtf_re = (t_symbol *)atom_getsymbolarg(2, argc, argv); + s_hrtf_im = (t_symbol *)atom_getsymbolarg(3, argc, argv); + s_fade_out_hrir = (t_symbol *)atom_getsymbolarg(4, argc, argv); + + n_order = (int)atom_getintarg(5, argc, argv); + n_dim = (int)atom_getintarg(6, argc, argv); + n_ls = (int)atom_getintarg(7, argc, argv); + fftsize = (int)atom_getintarg(8, argc, argv); + + ok = 1; + } + else if((argc >= 9) && + 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) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = gensym("L_HRIR"); + s_hrtf_re = gensym("HRTF_re"); + s_hrtf_im = gensym("HRTF_im"); + s_fade_out_hrir = gensym("HRIR_win"); + + n_order = (int)atom_getintarg(5, argc, argv); + n_dim = (int)atom_getintarg(6, argc, argv); + n_ls = (int)atom_getintarg(7, argc, argv); + fftsize = (int)atom_getintarg(8, argc, argv); + + ok = 1; + } + + if(ok) + { + if(n_order < 1) + n_order = 1; + + if(n_dim == 3) + { + if(n_order > 5) + n_order = 5; + n_ambi = (n_order + 1)*(n_order + 1); + } + else + { + if(n_order > 12) + n_order = 12; + n_dim = 2; + n_ambi = 2 * n_order + 1; + } + + if(n_ls < 1) + n_ls = 1; + if(n_ls < n_ambi) + { + n_ls = n_ambi; + post("bin_ambi_reduced_decode2-WARNING: Number of Loudspeakers < Number of Ambisonic-Channels !!!!"); + post(" Number of Loudspeakers set to %d", n_ambi); + } + + j = 2; + fftok = 0; + for(i=0; i<21; i++) + { + if(j == fftsize) + { + fftok = 1; + i = 22; + } + j *= 2; + } + + if(!fftok) + { + fftsize = 512; + post("bin_ambi_reduced_decode2-WARNING: fftsize not equal to 2 ^ n !!!"); + post(" fftsize set to %d", fftsize); + } + + x->x_n_dim = n_dim; + x->x_n_ambi = n_ambi; + x->x_n_ls = n_ls; + x->x_n_order = n_order; + x->x_fftsize = fftsize; + + x->x_hrir_filename = (t_symbol **)getbytes(x->x_n_ls * sizeof(t_symbol *)); + x->x_s_hrir = (t_symbol **)getbytes(x->x_n_ls * sizeof(t_symbol *)); + x->x_s_hrtf_re = (t_symbol **)getbytes(x->x_n_ambi * sizeof(t_symbol *)); + x->x_s_hrtf_im = (t_symbol **)getbytes(x->x_n_ambi * sizeof(t_symbol *)); + + for(i=0; i<n_ls; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrir->s_name); + x->x_s_hrir[i] = gensym(buf); + } + + for(i=0; i<n_ambi; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrtf_re->s_name); + x->x_s_hrtf_re[i] = gensym(buf); + + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrtf_im->s_name); + x->x_s_hrtf_im[i] = gensym(buf); + } + + sprintf(buf, "%d_%s", prefix, s_fade_out_hrir->s_name); + x->x_s_fade_out_hrir = gensym(buf); + + 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_ambi * sizeof(double)); + x->x_ls_encode = (double *)getbytes(x->x_n_ls * x->x_n_ambi * sizeof(double)); + x->x_prod = (double *)getbytes(x->x_n_ls * x->x_n_ambi * sizeof(double)); + x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double)); + + x->x_delta = (int *)getbytes(x->x_n_ls * sizeof(int)); + x->x_phi = (int *)getbytes(x->x_n_ls * sizeof(int)); + x->x_phi_sym = (int *)getbytes(x->x_n_ls * sizeof(int)); + x->x_sym_flag = (int *)getbytes(x->x_n_ls * sizeof(int)); + + x->x_spec = (BIN_AMBI_COMPLEX *)getbytes(x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + x->x_sin_cos = (BIN_AMBI_COMPLEX *)getbytes(x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + + x->x_beg_fade_out_hrir = (float *)0; + x->x_beg_hrir = (float *)getbytes(x->x_fftsize * x->x_n_ls * sizeof(float)); + x->x_beg_hrtf_re = (float **)getbytes(x->x_n_ambi * sizeof(float *)); + x->x_beg_hrtf_im = (float **)getbytes(x->x_n_ambi * sizeof(float *)); + + 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; + x->x_seq_ok = 1; + for(i=0; i<n_ambi; i++) + x->x_ambi_channel_weight[i] = 1.0; + + bin_ambi_reduced_decode2_init_cos(x); + + outlet_new(&x->x_obj, &s_list); + x->x_out_sign_sum = outlet_new(&x->x_obj, &s_list); + return(x); + } + else + { + post("bin_ambi_reduced_decode2-ERROR: need 1 float + 4 symbols + 4 floats arguments:"); + post(" prefix(unique-number) + hrir_name + hrtf_re_name + hrtf_im_name + hrir_fade_out_name +"); + post(" + ambi_order + ambi_dimension + number_of_loudspeakers + fftsize"); + return(0); + } +} + +void bin_ambi_reduced_decode2_setup(void) +{ + bin_ambi_reduced_decode2_class = class_new(gensym("bin_ambi_reduced_decode2"), (t_newmethod)bin_ambi_reduced_decode2_new, (t_method)bin_ambi_reduced_decode2_free, + sizeof(t_bin_ambi_reduced_decode2), 0, A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_ls, gensym("ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_calc_pinv, gensym("calc_pinv"), 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_load_HRIR, gensym("load_HRIR"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_check_HRIR_arrays, gensym("check_HRIR_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_check_HRTF_arrays, gensym("check_HRTF_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_calc_reduced, gensym("calc_reduced"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_decoding, gensym("decoding"), 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_ambi_weight, gensym("ambi_weight"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode2_class, (t_method)bin_ambi_reduced_decode2_sing_range, gensym("sing_range"), A_DEFFLOAT, 0); + class_sethelpsymbol(bin_ambi_reduced_decode2_class, gensym("iemhelp2/help-bin_ambi_reduced_decode2")); +} +/* +Reihenfolge: +n_ls x bin_ambi_reduced_decode2_ls + +bin_ambi_reduced_decode2_calc_pinv + +n_ls x bin_ambi_reduced_decode2_load_HRIR + +n_ls x bin_ambi_reduced_decode2_check_HRIR_arrays + +n_ambi x bin_ambi_reduced_decode2_check_HRTF_arrays + +n_ambi x bin_ambi_reduced_decode2_calc_reduced + +bin_ambi_reduced_decode2_decoding + +*/ diff --git a/src/bin_ambi_reduced_decode_fft.c b/src/bin_ambi_reduced_decode_fft.c new file mode 100644 index 0000000..092bbfc --- /dev/null +++ b/src/bin_ambi_reduced_decode_fft.c @@ -0,0 +1,1512 @@ +/* For information on usage and redistribution, and for a DISCLAIMER OF ALL +* WARRANTIES, see the file, "LICENSE.txt," in this distribution. + +iem_bin_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_bin_ambi.h" +#include <math.h> +#include <stdio.h> +#include <string.h> + + +/* -------------------------- bin_ambi_reduced_decode_fft ------------------------------ */ +/* + ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren ** + ** Inputs: ls + Liste von 3 floats: Index [1 .. 16] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] ** + ** Inputs: calc_inv ** + ** Inputs: load_HRIR + float index1..16 ** + ** 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 1016_1_L_HRIR .. 1016_16_L_HRIR heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_re .. 1016_16_HRTF_re heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_im .. 1016_16_HRTF_im heissen und existieren ** + */ + +typedef struct _bin_ambi_reduced_decode_fft +{ + t_object x_obj; + t_atom x_at[2];/*output filename and tablename of HRIR*/ + int x_n_dim;/*dimension of ambisonic system*/ + int x_n_ambi;/*number of ambi channels*/ + int x_n_order;/*order of ambisonic system*/ + int x_n_ind_ls;/*number of independent loudspeakers*/ + int x_n_mrg_mir_ls;/*number of mirrored and merged loudspeakers*/ + int x_n_ph_ls;/*number of phantom loudspeakers*/ + int x_seq_ok; + int x_fftsize;/*fftsize of blockfilter, should be 2 times of FIR-length of HRIR*/ + double *x_inv_work1;/* n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_work2;/* 2*n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_buf2;/* 2*n_ambi buffer for matrix inverse */ + double *x_transp;/* n_ls*n_ambi transposed input-buffer of decoder-matrix before inversion */ + double *x_ls_encode;/* n_ambi*n_ls straight input-buffer of decoder-matrix before inversion */ + double *x_prod2;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_prod3;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_ambi_channel_weight; + double x_mirror_weight; + int *x_delta; + int *x_phi; + int *x_phi_sym; + int *x_sym_flag; + BIN_AMBI_COMPLEX *x_spec; + BIN_AMBI_COMPLEX *x_sin_cos; + float *x_beg_fade_out_hrir; + float *x_beg_hrir; + float **x_beg_hrtf_re; + float **x_beg_hrtf_im; + t_symbol **x_hrir_filename; + t_symbol **x_s_hrir; + t_symbol **x_s_hrtf_re; + t_symbol **x_s_hrtf_im; + t_symbol *x_s_fade_out_hrir; + void *x_out_sign_sum; + 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; + double x_sing_range; +} t_bin_ambi_reduced_decode_fft; + +static t_class *bin_ambi_reduced_decode_fft_class; + +static void bin_ambi_reduced_decode_fft_init_cos(t_bin_ambi_reduced_decode_fft *x) +{/* initialize a whole sine and cosine wave over a fftsize array */ + int i, fftsize = x->x_fftsize; + float f, g; + BIN_AMBI_COMPLEX *sincos = x->x_sin_cos; + + g = 2.0f * 3.1415926538f / (float)fftsize; + for(i=0; i<fftsize; i++) + { + f = g * (float)i; + (*sincos).real = cos(f); + (*sincos).imag = -sin(f); /*FFT*/ + sincos++; + } +} + +static void bin_ambi_reduced_decode_fft_quant(t_bin_ambi_reduced_decode_fft *x, double *delta_deg2rad, double *phi_deg2rad, int index) +{ + double q = 1.0; + double d = *delta_deg2rad; + double p = *phi_deg2rad; + int i; + + if(d < -40.0) + d = -40.0; + if(d > 90.0) + d = 90.0; + while(p < 0.0) + p += 360.0; + while(p >= 360.0) + p -= 360.0; + + if(d < -35.0) + { + *delta_deg2rad = -40.0; + q = 360.0 / 56.0; + } + else if(d < -25.0) + { + *delta_deg2rad = -30.0; + q = 6.0; + } + else if(d < -15.0) + { + *delta_deg2rad = -20.0; + q = 5.0; + } + else if(d < -5.0) + { + *delta_deg2rad = -10.0; + q = 5.0; + } + else if(d < 5.0) + { + *delta_deg2rad = 0.0; + q = 5.0; + } + else if(d < 15.0) + { + *delta_deg2rad = 10.0; + q = 5.0; + } + else if(d < 25.0) + { + *delta_deg2rad = 20.0; + q = 5.0; + } + else if(d < 35.0) + { + *delta_deg2rad = 30.0; + q = 6.0; + } + else if(d < 45.0) + { + *delta_deg2rad = 40.0; + q = 360.0 / 56.0; + } + else if(d < 55.0) + { + *delta_deg2rad = 50.0; + q = 8.0; + } + else if(d < 65.0) + { + *delta_deg2rad = 60.0; + q = 10.0; + } + else if(d < 75.0) + { + *delta_deg2rad = 70.0; + q = 15.0; + } + else if(d < 85.0) + { + *delta_deg2rad = 80.0; + q = 30.0; + } + else + { + *delta_deg2rad = 90.0; + q = 360.0; + } + + p /= q; + i = (int)(p + 0.499999); + p = (double)i * q; + i = (int)(p + 0.499999); + while(i >= 360) + i -= 360; + p = (double)i; + *phi_deg2rad = p; + + x->x_delta[index] = (int)(*delta_deg2rad); + x->x_phi[index] = i; + + *delta_deg2rad *= x->x_pi_over_180; + *phi_deg2rad *= x->x_pi_over_180; +} + +static void bin_ambi_reduced_decode_fft_do_2d(t_bin_ambi_reduced_decode_fft *x, int argc, t_atom *argv, int mode) +{ + double delta=0.0, phi; + double *dw = x->x_transp; + int index; + int order=x->x_n_order; + + if(argc < 2) + { + post("bin_ambi_reduced_decode_fft 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 == BIN_AMBI_LS_IND) + { + if(index >= x->x_n_ind_ls) + index = x->x_n_ind_ls - 1; + } + else if(mode == BIN_AMBI_LS_MRG) + { + if(x->x_n_mrg_mir_ls) + { + if(index >= x->x_n_mrg_mir_ls) + index = x->x_n_mrg_mir_ls - 1; + index += x->x_n_ind_ls; + } + else + return; + } + else if(mode == BIN_AMBI_LS_MIR) + { + if(x->x_n_mrg_mir_ls) + { + if(index >= x->x_n_mrg_mir_ls) + index = x->x_n_mrg_mir_ls - 1; + index += x->x_n_ind_ls; + index += x->x_n_mrg_mir_ls; + } + else + return; + } + else if(mode == BIN_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_ind_ls; + index += 2*x->x_n_mrg_mir_ls; + } + else + return; + } + else + return; + + bin_ambi_reduced_decode_fft_quant(x, &delta, &phi, index); + + 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 bin_ambi_reduced_decode_fft_do_3d(t_bin_ambi_reduced_decode_fft *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("bin_ambi_reduced_decode_fft 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 == BIN_AMBI_LS_IND) + { + if(index >= x->x_n_ind_ls) + index = x->x_n_ind_ls - 1; + } + else if(mode == BIN_AMBI_LS_MRG) + { + if(x->x_n_mrg_mir_ls) + { + if(index >= x->x_n_mrg_mir_ls) + index = x->x_n_mrg_mir_ls - 1; + index += x->x_n_ind_ls; + } + else + return; + } + else if(mode == BIN_AMBI_LS_MIR) + { + if(x->x_n_mrg_mir_ls) + { + if(index >= x->x_n_mrg_mir_ls) + index = x->x_n_mrg_mir_ls - 1; + index += x->x_n_ind_ls; + index += x->x_n_mrg_mir_ls; + } + else + return; + } + else if(mode == BIN_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_ind_ls; + index += 2*x->x_n_mrg_mir_ls; + } + else + return; + } + else + return; + + bin_ambi_reduced_decode_fft_quant(x, &delta, &phi, index); + + cd = cos(delta); + sd = sin(delta); + cp = cos(phi); + sp = sin(phi); + + dw += index * x->x_n_ambi; + + *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 bin_ambi_reduced_decode_fft_ind_ls(t_bin_ambi_reduced_decode_fft *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fft_do_2d(x, argc, argv, BIN_AMBI_LS_IND); + else + bin_ambi_reduced_decode_fft_do_3d(x, argc, argv, BIN_AMBI_LS_IND); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode_fft_mrg_ls(t_bin_ambi_reduced_decode_fft *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fft_do_2d(x, argc, argv, BIN_AMBI_LS_MRG); + else + bin_ambi_reduced_decode_fft_do_3d(x, argc, argv, BIN_AMBI_LS_MRG); +} + +static void bin_ambi_reduced_decode_fft_mir_ls(t_bin_ambi_reduced_decode_fft *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fft_do_2d(x, argc, argv, BIN_AMBI_LS_MIR); + else + bin_ambi_reduced_decode_fft_do_3d(x, argc, argv, BIN_AMBI_LS_MIR); +} + +static void bin_ambi_reduced_decode_fft_pht_ls(t_bin_ambi_reduced_decode_fft *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fft_do_2d(x, argc, argv, BIN_AMBI_LS_PHT); + else + bin_ambi_reduced_decode_fft_do_3d(x, argc, argv, BIN_AMBI_LS_PHT); +} + +static void bin_ambi_reduced_decode_fft_copy_row2buf(t_bin_ambi_reduced_decode_fft *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 bin_ambi_reduced_decode_fft_copy_buf2row(t_bin_ambi_reduced_decode_fft *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 bin_ambi_reduced_decode_fft_copy_row2row(t_bin_ambi_reduced_decode_fft *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 bin_ambi_reduced_decode_fft_xch_rows(t_bin_ambi_reduced_decode_fft *x, int row1, int row2) +{ + bin_ambi_reduced_decode_fft_copy_row2buf(x, row1); + bin_ambi_reduced_decode_fft_copy_row2row(x, row2, row1); + bin_ambi_reduced_decode_fft_copy_buf2row(x, row2); +} + +static void bin_ambi_reduced_decode_fft_mul_row(t_bin_ambi_reduced_decode_fft *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 bin_ambi_reduced_decode_fft_mul_col(t_bin_ambi_reduced_decode_fft *x, int col, double mul) +{ + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int i; + double *dw=x->x_inv_work2; + + dw += col; + for(i=0; i<n_ambi; i++) + { + (*dw) *= mul; + dw += n_ambi2; + } +} + +static void bin_ambi_reduced_decode_fft_mul_buf_and_add2row(t_bin_ambi_reduced_decode_fft *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 bin_ambi_reduced_decode_fft_eval_which_element_of_col_not_zero(t_bin_ambi_reduced_decode_fft *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 bin_ambi_reduced_decode_fft_mul1(t_bin_ambi_reduced_decode_fft *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_ind_ls+2*x->x_n_mrg_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 bin_ambi_reduced_decode_fft_mul2(t_bin_ambi_reduced_decode_fft *x) +{ + int n_ls=x->x_n_ind_ls+2*x->x_n_mrg_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_prod2; + 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 bin_ambi_reduced_decode_fft_mul3(t_bin_ambi_reduced_decode_fft *x) +{ + int i, n; + double *dv3=x->x_prod3; + double *dv2=x->x_prod2; + double *dv1=x->x_prod2; + double mw=x->x_mirror_weight; + + n = x->x_n_ind_ls * x->x_n_ambi; + for(i=0; i<n; i++) + *dv3++ = *dv1++; + dv2 += n; + n = x->x_n_mrg_mir_ls * x->x_n_ambi; + dv2 += n; + for(i=0; i<n; i++) + *dv3++ = *dv1++ + (*dv2++)*mw; +} + +static void bin_ambi_reduced_decode_fft_transp_back(t_bin_ambi_reduced_decode_fft *x) +{ + double *vec, *transp=x->x_transp; + double *straight=x->x_ls_encode; + int n_ls=x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls; + int n_ambi=x->x_n_ambi; + int i, j; + + transp = x->x_transp; + for(j=0; j<n_ambi; j++) + { + vec = transp; + for(i=0; i<n_ls; i++) + { + *straight++ = *vec; + vec += n_ambi; + } + transp++; + } +} + +static void bin_ambi_reduced_decode_fft_inverse(t_bin_ambi_reduced_decode_fft *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 = bin_ambi_reduced_decode_fft_eval_which_element_of_col_not_zero(x, i, i); + if(nz < 0) + { + post("bin_ambi_reduced_decode_fft ERROR: matrix not regular !!!!"); + x->x_seq_ok = 0; + return; + } + else + { + if(nz != i) + bin_ambi_reduced_decode_fft_xch_rows(x, i, nz); + dv = db + i*n_ambi2 + i; + rcp = 1.0 /(*dv); + bin_ambi_reduced_decode_fft_mul_row(x, i, rcp); + bin_ambi_reduced_decode_fft_copy_row2buf(x, i); + for(j=i+1; j<n_ambi; j++) + { + dv += n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode_fft_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; + bin_ambi_reduced_decode_fft_copy_row2buf(x, i); + for(j=i-1; j>=0; j--) + { + dv -= n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode_fft_mul_buf_and_add2row(x, j, rcp); + } + } + + post("matrix_inverse regular"); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode_fft_calc_pinv(t_bin_ambi_reduced_decode_fft *x) +{ + t_garray *a; + int npoints; + t_float *fadevec; + + if((int)(x->x_beg_fade_out_hrir) == 0) + { + if (!(a = (t_garray *)pd_findbyclass(x->x_s_fade_out_hrir, garray_class))) + error("%s: no such array", x->x_s_fade_out_hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &fadevec)) + error("%s: bad template for bin_ambi_reduced_decode_fft", x->x_s_fade_out_hrir->s_name); + else if (npoints < x->x_fftsize) + error("%s: bad array-size: %d", x->x_s_fade_out_hrir->s_name, npoints); + else + x->x_beg_fade_out_hrir = fadevec; + } + + bin_ambi_reduced_decode_fft_transp_back(x); + bin_ambi_reduced_decode_fft_mul1(x); + bin_ambi_reduced_decode_fft_inverse(x); + bin_ambi_reduced_decode_fft_mul2(x); + bin_ambi_reduced_decode_fft_mul3(x); +} + +/* +x_prod: +n_ambi columns; +n_ambi rows; +*/ + +static void bin_ambi_reduced_decode_fft_load_HRIR(t_bin_ambi_reduced_decode_fft *x, float findex) +{ + int index=(int)findex - 1; + int p; + char buf[60]; + + if(index < 0) + index = 0; + if(index >= (x->x_n_ind_ls + x->x_n_mrg_mir_ls)) + index = x->x_n_ind_ls + x->x_n_mrg_mir_ls - 1; + + p = x->x_phi[index]; + + if(p)/*change*/ + p = 360 - p; + + if(p < 10) + sprintf(buf, "L%de00%da.wav", x->x_delta[index], p); + else if(p < 100) + sprintf(buf, "L%de0%da.wav", x->x_delta[index], p); + else + sprintf(buf, "L%de%da.wav", x->x_delta[index], p); + x->x_hrir_filename[index] = gensym(buf); + + SETSYMBOL(x->x_at, x->x_hrir_filename[index]); + SETSYMBOL(x->x_at+1, x->x_s_hrir[index]); + outlet_list(x->x_obj.ob_outlet, &s_list, 2, x->x_at); +} + +static void bin_ambi_reduced_decode_fft_check_HRIR_arrays(t_bin_ambi_reduced_decode_fft *x, float findex) +{ + int index=(int)findex - 1; + int j, k, n; + int fftsize = x->x_fftsize; + int fs2=fftsize/2; + t_garray *a; + int npoints; + t_symbol *hrir; + t_float *vec_hrir, *vec, *vec_fade_out_hrir; + float decr, sum; + + if(index < 0) + index = 0; + if(index >= (x->x_n_ind_ls + x->x_n_mrg_mir_ls)) + index = x->x_n_ind_ls + x->x_n_mrg_mir_ls - 1; + + hrir = x->x_s_hrir[index]; + if (!(a = (t_garray *)pd_findbyclass(hrir, garray_class))) + error("%s: no such array", hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrir)) + error("%s: bad template for bin_ambi_reduced_decode_fft", hrir->s_name); + else + { + if(npoints < fftsize) + { + post("bin_ambi_reduced_decode_fft-WARNING: %s-array-size: %d < FFT-size: %d", hrir->s_name, npoints, fftsize); + } + vec = x->x_beg_hrir; + vec += index * fftsize; + + if((int)(x->x_beg_fade_out_hrir)) + { + vec_fade_out_hrir = x->x_beg_fade_out_hrir; + for(j=0; j<fs2; j++) + vec[j] = vec_hrir[j]*vec_fade_out_hrir[j]; + } + else + { + post("no HRIR-fade-out-window found"); + n = fs2 * 3; + n /= 4; + for(j=0; j<n; j++) + vec[j] = vec_hrir[j]; + sum = 1.0f; + decr = 4.0f / (float)fs2; + for(j=n, k=0; j<fs2; j++, k++) + { + sum -= decr; + vec[j] = vec_hrir[j] * sum; + } + } + } +} + +static void bin_ambi_reduced_decode_fft_check_HRTF_arrays(t_bin_ambi_reduced_decode_fft *x, float findex) +{ + int index=(int)findex - 1; + t_garray *a; + int npoints; + int fftsize = x->x_fftsize; + t_float *vec_hrtf_re, *vec_hrtf_im; + t_symbol *hrtf_re, *hrtf_im; + + if(index < 0) + index = 0; + if(index >= x->x_n_ambi) + index = x->x_n_ambi - 1; + + hrtf_re = x->x_s_hrtf_re[index]; + hrtf_im = x->x_s_hrtf_im[index]; + + if (!(a = (t_garray *)pd_findbyclass(hrtf_re, garray_class))) + error("%s: no such array", hrtf_re->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrtf_re)) + error("%s: bad template for bin_ambi_reduced_decode_fft", hrtf_re->s_name); + else if (npoints < fftsize) + error("%s: bad array-size: %d", hrtf_re->s_name, npoints); + else if (!(a = (t_garray *)pd_findbyclass(hrtf_im, garray_class))) + error("%s: no such array", hrtf_im->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrtf_im)) + error("%s: bad template for bin_ambi_reduced_decode_fft", hrtf_im->s_name); + else if (npoints < fftsize) + error("%s: bad array-size: %d", hrtf_im->s_name, npoints); + else + { + x->x_beg_hrtf_re[index] = vec_hrtf_re; + x->x_beg_hrtf_im[index] = vec_hrtf_im; + } +} + +static void bin_ambi_reduced_decode_fft_calc_reduced(t_bin_ambi_reduced_decode_fft *x, float findex) +{ + int index=(int)findex - 1; + int i, j, k, w_index, w_inc, i_inc, v_index, fs1, fs2; + int fftsize = x->x_fftsize; + BIN_AMBI_COMPLEX old1, old2, w; + BIN_AMBI_COMPLEX *sincos = x->x_sin_cos; + BIN_AMBI_COMPLEX *val = x->x_spec;/*weighted decoder-matrix with n_ls lines and n_ambi columns*/ + t_float *vec_hrir, *vec_hrtf_re, *vec_hrtf_im; + double *dv, *db=x->x_prod3; + int n_ambi = x->x_n_ambi; + int n_ls = x->x_n_ind_ls + x->x_n_mrg_mir_ls; + float mul; + + if(x->x_seq_ok) + { + if(index < 0) + index = 0; + if(index >= n_ambi) + index = n_ambi - 1; + + vec_hrtf_re = x->x_beg_hrtf_re[index]; + vec_hrtf_im = x->x_beg_hrtf_im[index]; + + dv = db + index; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + for(k=0; k<fftsize; k++)/*first step of acumulating the HRIRs*/ + { + val[k].real = mul * vec_hrir[k]; + val[k].imag = 0.0f; + } + + for(j=1; j<n_ls; j++) + { + dv += n_ambi; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + vec_hrir += j * fftsize; + for(k=0; k<fftsize; k++) + { + val[k].real += mul * vec_hrir[k]; + } + } + + fs1 = fftsize - 1; + fs2 = fftsize / 2; + i_inc = fs2; + w_inc = 1; + for(i=1; i<fftsize; i<<=1) + { + v_index = 0; + for(j=0; j<i; j++) + { + w_index = 0; + for(k=0; k<i_inc; k++) + { + old1 = val[v_index]; + old2 = val[v_index+i_inc]; + w = sincos[w_index]; + val[v_index+i_inc].real = (old1.real-old2.real)*w.real - (old1.imag-old2.imag)*w.imag; + val[v_index+i_inc].imag = (old1.imag - old2.imag)*w.real + (old1.real-old2.real)*w.imag; + val[v_index].real = old1.real+old2.real; + val[v_index].imag = old1.imag+old2.imag; + w_index += w_inc; + v_index++; + } + v_index += i_inc; + } + w_inc <<= 1; + i_inc >>= 1; + } + + j = 0; + for(i=1;i<fs1;i++) + { + k = fs2; + while(k <= j) + { + j = j - k; + k >>= 1; + } + j = j + k; + if(i < j) + { + old1 = val[j]; + val[j] = val[i]; + val[i] = old1; + } + } + for(i = 0; i<fs2; i++) + { + vec_hrtf_re[i] = val[i].real; + vec_hrtf_im[i] = val[i].imag; + } + vec_hrtf_re[fs2] = val[fs2].real; + vec_hrtf_im[fs2] = 0.0f; + for(i = fs2+1; i < fftsize; i++) + { + vec_hrtf_re[i] = 0.0f; + vec_hrtf_im[i] = 0.0f; + } + } +} + +static void bin_ambi_reduced_decode_fft_calc_sym(t_bin_ambi_reduced_decode_fft *x) +{ + int *sym=x->x_phi_sym; + int n_ls = x->x_n_ind_ls + x->x_n_mrg_mir_ls; + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int *phi=x->x_phi; + int *delta=x->x_delta; + int *flag=x->x_sym_flag; + int i, j, d, p, regular=1; + double *dv, *db=x->x_prod3; + double a, b, q; + int sym_max=0; + int pos_sym_counter, neg_sym_counter; + char plus_minus[100]; + char abc[100]; + char onenine[100]; + char ten[100]; + + if(x->x_seq_ok) + { + for(i=0; i<n_ls; i++) + sym[i] = -1; + + plus_minus[0] = 0; + if(x->x_n_dim == 2) + { + strcpy(abc, " WXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXY"); + strcpy(ten, " 000000000000000000011111111111111111111"); + strcpy(onenine, " 011223344556677889900112233445566778899"); + abc[n_ambi+22] = 0; + ten[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(ten); + post(onenine); + } + else + { + strcpy(abc, " AABCABCDEABCDEFGABCDEFGHIABCDEFGHIJKABCDEFGHIJKLM"); + strcpy(onenine, " 0111222223333333444444444555555555556666666666666"); + abc[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(onenine); + } + + for(i=0; i<n_ls; i++)/* looking for azimuth symmetries of loudspeakers */ + { + d = delta[i]; + p = phi[i]; + for(j=i+1; j<n_ls; j++) + { + if(d == delta[j]) + { + if((p + phi[j]) == 360) + { + if((sym[i] < 0) && (sym[j] < 0)) + { + sym[i] = j; + sym[j] = i; + j = n_ls + 1; + sym_max++; + } + } + } + } + } + + for(p=0; p<n_ambi; p++)/*each col*/ + { + pos_sym_counter = 0; + neg_sym_counter = 0; + for(i=0; i<n_ls; i++) + flag[i] = 1; + dv = db + p; + for(i=0; i<n_ls; i++) + { + if((sym[i] >= 0) && flag[i]) + { + j = sym[i]; + flag[i] = 0; + flag[j] = 0; + a = dv[i*n_ambi]; + b = dv[j*n_ambi]; + if((a < 5.0e-4)&&(a > -5.0e-4)&&(b < 5.0e-4)&&(b > -5.0e-4)) + { + pos_sym_counter++; + neg_sym_counter++; + } + else + { + q = a / b; + if((q < 1.005)&&(q > 0.995)) + pos_sym_counter++; + else if((q > -1.005)&&(q < -0.995)) + neg_sym_counter++; + } + } + } + if(pos_sym_counter == sym_max) + strcat(plus_minus, "+"); + else if(neg_sym_counter == sym_max) + strcat(plus_minus, "-"); + else + { + strcat(plus_minus, "?"); + regular = 0; + } + } + post("sum of right channel: %s", plus_minus); + + if(regular) + { + for(i=0; i<n_ambi; i++) + { + /*change*/ + if(plus_minus[i] == '+') + SETFLOAT(x->x_at, 1.0f); + else if(plus_minus[i] == '-') + SETFLOAT(x->x_at, 2.0f); + SETFLOAT(x->x_at+1, (float)(i+1)); + outlet_list(x->x_out_sign_sum, &s_list, 2, x->x_at); + } + } + } +} + +static void bin_ambi_reduced_decode_fft_ambi_weight(t_bin_ambi_reduced_decode_fft *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] = (double)atom_getfloat(argv++); + k++; + if(x->x_n_dim == 2) + { + for(i=1; i<=n; i++) + { + d = (double)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 = (double)atom_getfloat(argv++); + m = 2*i + 1; + for(j=0; j<m; j++) + { + x->x_ambi_channel_weight[k] = d; + k++; + } + } + } + } + else + post("bin_ambi_reduced_decode_fft-ERROR: ambi_weight needs %d float weights", x->x_n_order+1); +} + +static void bin_ambi_reduced_decode_fft_mirror_weight(t_bin_ambi_reduced_decode_fft *x, t_floatarg f) +{ + x->x_mirror_weight = (double)f; +} + +static void bin_ambi_reduced_decode_fft_sing_range(t_bin_ambi_reduced_decode_fft *x, t_floatarg f) +{ + if(f < 0.0f) + x->x_sing_range = -(double)f; + else + x->x_sing_range = (double)f; +} + +static void bin_ambi_reduced_decode_fft_free(t_bin_ambi_reduced_decode_fft *x) +{ + freebytes(x->x_hrir_filename, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(t_symbol *)); + freebytes(x->x_s_hrir, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(t_symbol *)); + freebytes(x->x_s_hrtf_re, x->x_n_ambi * sizeof(t_symbol *)); + freebytes(x->x_s_hrtf_im, x->x_n_ambi * sizeof(t_symbol *)); + 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_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_ls_encode, (x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_prod2, (x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_prod3, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double)); + + freebytes(x->x_delta, (x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * sizeof(int)); + freebytes(x->x_phi, (x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * sizeof(int)); + freebytes(x->x_phi_sym, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(int)); + freebytes(x->x_sym_flag, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(int)); + + freebytes(x->x_spec, x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + freebytes(x->x_sin_cos, x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + + freebytes(x->x_beg_hrir, x->x_fftsize * (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(float)); + freebytes(x->x_beg_hrtf_re, x->x_n_ambi * sizeof(float *)); + freebytes(x->x_beg_hrtf_im, x->x_n_ambi * sizeof(float *)); +} + +/* + 1.arg_ int prefix; + 2.arg: t_symbol *hrir_name; + 3.arg: t_symbol *hrtf_re_name; + 4.arg: t_symbol *hrtf_im_name; + 5.arg: t_symbol *hrir_fade_out_name; + 6.arg: int ambi_order; + 7.arg: int dim; + 8.arg: int number_of_loudspeakers; + 9.arg: int fftsize; + */ + +static void *bin_ambi_reduced_decode_fft_new(t_symbol *s, int argc, t_atom *argv) +{ + t_bin_ambi_reduced_decode_fft *x = (t_bin_ambi_reduced_decode_fft *)pd_new(bin_ambi_reduced_decode_fft_class); + char buf[400]; + int i, j, fftok, ok=0; + int n_order=0, n_dim=0, n_ind_ls=0, n_mrg_mir_ls=0, n_ph_ls=0, n_ambi=0, fftsize=0, prefix=0; + t_symbol *s_hrir=gensym("L_HRIR"); + t_symbol *s_hrtf_re=gensym("HRTF_re"); + t_symbol *s_hrtf_im=gensym("HRTF_im"); + t_symbol *s_fade_out_hrir=gensym("HRIR_win"); + + if((argc >= 11) && + IS_A_FLOAT(argv,0) && + IS_A_SYMBOL(argv,1) && + IS_A_SYMBOL(argv,2) && + IS_A_SYMBOL(argv,3) && + IS_A_SYMBOL(argv,4) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8) && + IS_A_FLOAT(argv,9) && + IS_A_FLOAT(argv,10)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = (t_symbol *)atom_getsymbolarg(1, argc, argv); + s_hrtf_re = (t_symbol *)atom_getsymbolarg(2, argc, argv); + s_hrtf_im = (t_symbol *)atom_getsymbolarg(3, argc, argv); + s_fade_out_hrir = (t_symbol *)atom_getsymbolarg(4, argc, argv); + + n_order = (int)atom_getintarg(5, argc, argv); + n_dim = (int)atom_getintarg(6, argc, argv); + n_ind_ls = (int)atom_getintarg(7, argc, argv); + n_mrg_mir_ls = (int)atom_getintarg(8, argc, argv); + n_ph_ls = (int)atom_getintarg(9, argc, argv); + fftsize = (int)atom_getintarg(10, argc, argv); + + ok = 1; + } + else if((argc >= 11) && + 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) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8) && + IS_A_FLOAT(argv,9) && + IS_A_FLOAT(argv,10)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = gensym("L_HRIR"); + s_hrtf_re = gensym("HRTF_re"); + s_hrtf_im = gensym("HRTF_im"); + s_fade_out_hrir = gensym("HRIR_win"); + + n_order = (int)atom_getintarg(5, argc, argv); + n_dim = (int)atom_getintarg(6, argc, argv); + n_ind_ls = (int)atom_getintarg(7, argc, argv); + n_mrg_mir_ls = (int)atom_getintarg(8, argc, argv); + n_ph_ls = (int)atom_getintarg(9, argc, argv); + fftsize = (int)atom_getintarg(10, argc, argv); + + ok = 1; + } + + if(ok) + { + if(n_order < 1) + n_order = 1; + + if(n_dim == 3) + { + if(n_order > 5) + n_order = 5; + n_ambi = (n_order + 1)*(n_order + 1); + } + else + { + if(n_order > 12) + n_order = 12; + n_dim = 2; + n_ambi = 2 * n_order + 1; + } + + if(n_ind_ls < 1) + n_ind_ls = 1; + if(n_mrg_mir_ls < 1) + n_mrg_mir_ls = 1; + if(n_ph_ls < 0) + n_ph_ls = 0; + + if((n_ind_ls+2*n_mrg_mir_ls+n_ph_ls) < n_ambi) + { + post("bin_ambi_reduced_decode_fft-WARNING: Number of all Loudspeakers < Number of Ambisonic-Channels !!!!"); + } + + j = 2; + fftok = 0; + for(i=0; i<21; i++) + { + if(j == fftsize) + { + fftok = 1; + i = 22; + } + j *= 2; + } + + if(!fftok) + { + fftsize = 512; + post("bin_ambi_reduced_decode_fft-WARNING: fftsize not equal to 2 ^ n !!!"); + post(" fftsize set to %d", fftsize); + } + + x->x_n_dim = n_dim; + x->x_n_ambi = n_ambi; + x->x_n_ind_ls = n_ind_ls; + x->x_n_mrg_mir_ls = n_mrg_mir_ls; + x->x_n_ph_ls = n_ph_ls; + x->x_n_order = n_order; + x->x_fftsize = fftsize; + + x->x_hrir_filename = (t_symbol **)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(t_symbol *)); + x->x_s_hrir = (t_symbol **)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(t_symbol *)); + x->x_s_hrtf_re = (t_symbol **)getbytes(x->x_n_ambi * sizeof(t_symbol *)); + x->x_s_hrtf_im = (t_symbol **)getbytes(x->x_n_ambi * sizeof(t_symbol *)); + + j = x->x_n_ind_ls + x->x_n_mrg_mir_ls; + for(i=0; i<j; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrir->s_name); + x->x_s_hrir[i] = gensym(buf); + } + + for(i=0; i<n_ambi; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrtf_re->s_name); + x->x_s_hrtf_re[i] = gensym(buf); + + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrtf_im->s_name); + x->x_s_hrtf_im[i] = gensym(buf); + } + + sprintf(buf, "%d_%s", prefix, s_fade_out_hrir->s_name); + x->x_s_fade_out_hrir = gensym(buf); + + 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_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + x->x_ls_encode = (double *)getbytes((x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + x->x_prod2 = (double *)getbytes((x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + x->x_prod3 = (double *)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * x->x_n_ambi * sizeof(double)); + x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double)); + + x->x_delta = (int *)getbytes((x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * sizeof(int)); + x->x_phi = (int *)getbytes((x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * sizeof(int)); + x->x_phi_sym = (int *)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(int)); + x->x_sym_flag = (int *)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(int)); + + x->x_spec = (BIN_AMBI_COMPLEX *)getbytes(x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + x->x_sin_cos = (BIN_AMBI_COMPLEX *)getbytes(x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + + x->x_beg_fade_out_hrir = (float *)0; + x->x_beg_hrir = (float *)getbytes(x->x_fftsize * (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(float)); + x->x_beg_hrtf_re = (float **)getbytes(x->x_n_ambi * sizeof(float *)); + x->x_beg_hrtf_im = (float **)getbytes(x->x_n_ambi * sizeof(float *)); + + 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; + x->x_seq_ok = 1; + for(i=0; i<n_ambi; i++) + x->x_ambi_channel_weight[i] = 1.0; + x->x_mirror_weight = 0.0; + + bin_ambi_reduced_decode_fft_init_cos(x); + + outlet_new(&x->x_obj, &s_list); + x->x_out_sign_sum = outlet_new(&x->x_obj, &s_list); + return(x); + } + else + { + post("bin_ambi_reduced_decode_fft-ERROR: need 1 float + 4 symbols + 6 floats arguments:"); + post(" prefix(unique-number) + hrir_name + hrtf_re_name + hrtf_im_name + hrir_fade_out_name +"); + post(" + ambi_order + ambi_dimension + number_of_independent_loudspeakers + "); + post(" + number_of_mirrored_and_merged_loudspeakers + number_of_phantom_loudspeakers + fftsize"); + return(0); + } +} + +void bin_ambi_reduced_decode_fft_setup(void) +{ + bin_ambi_reduced_decode_fft_class = class_new(gensym("bin_ambi_reduced_decode_fft"), (t_newmethod)bin_ambi_reduced_decode_fft_new, (t_method)bin_ambi_reduced_decode_fft_free, + sizeof(t_bin_ambi_reduced_decode_fft), 0, A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_ind_ls, gensym("ind_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_mrg_ls, gensym("mrg_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_mir_ls, gensym("mir_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_pht_ls, gensym("pht_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_mirror_weight, gensym("mirror_weight"), A_DEFFLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_calc_pinv, gensym("calc_pinv"), 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_load_HRIR, gensym("load_HRIR"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_check_HRIR_arrays, gensym("check_HRIR_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_check_HRTF_arrays, gensym("check_HRTF_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_calc_reduced, gensym("calc_reduced"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_calc_sym, gensym("calc_sym"), 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_ambi_weight, gensym("ambi_weight"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft_class, (t_method)bin_ambi_reduced_decode_fft_sing_range, gensym("sing_range"), A_DEFFLOAT, 0); + class_sethelpsymbol(bin_ambi_reduced_decode_fft_class, gensym("iemhelp2/help-bin_ambi_reduced_decode_fft")); +} +/* +Reihenfolge: +n_ls x bin_ambi_reduced_decode_fft_ls + +bin_ambi_reduced_decode_fft_calc_pinv + +n_ls x bin_ambi_reduced_decode_fft_load_HRIR + +n_ls x bin_ambi_reduced_decode_fft_check_HRIR_arrays + +n_ambi x bin_ambi_reduced_decode_fft_check_HRTF_arrays + +n_ambi x bin_ambi_reduced_decode_fft_calc_reduced + +bin_ambi_reduced_decode_fft_calc_sym + +*/ diff --git a/src/bin_ambi_reduced_decode_fft2.c b/src/bin_ambi_reduced_decode_fft2.c new file mode 100644 index 0000000..7e22682 --- /dev/null +++ b/src/bin_ambi_reduced_decode_fft2.c @@ -0,0 +1,1388 @@ +/* For information on usage and redistribution, and for a DISCLAIMER OF ALL +* WARRANTIES, see the file, "LICENSE.txt," in this distribution. + +iem_bin_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_bin_ambi.h" +#include <math.h> +#include <stdio.h> +#include <string.h> + + +/* -------------------------- bin_ambi_reduced_decode_fft2 ------------------------------ */ +/* + ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren ** + ** Inputs: ls + Liste von 3 floats: Index [1 .. 16] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] ** + ** Inputs: calc_inv ** + ** Inputs: load_HRIR + float index1..16 ** + ** 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 1016_1_L_HRIR .. 1016_16_L_HRIR heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_re .. 1016_16_HRTF_re heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_im .. 1016_16_HRTF_im heissen und existieren ** + */ + +typedef struct _bin_ambi_reduced_decode_fft2 +{ + t_object x_obj; + t_atom x_at[2];/*output filename and tablename of HRIR*/ + int x_n_dim;/*dimension of ambisonic system*/ + int x_n_ambi;/*number of ambi channels*/ + int x_n_order;/*order of ambisonic system*/ + int x_n_real_ls;/*number of real loudspeakers*/ + int x_n_pht_ls;/*number of phantom loudspeakers*/ + int x_seq_ok; + int x_fftsize;/*fftsize of blockfilter, should be 2 times of FIR-length of HRIR*/ + double *x_inv_work1;/* n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_work2;/* 2*n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_buf2;/* 2*n_ambi buffer for matrix inverse */ + double *x_transp;/* n_ls*n_ambi transposed input-buffer of decoder-matrix before inversion */ + double *x_ls_encode;/* n_ambi*n_ls straight input-buffer of decoder-matrix before inversion */ + double *x_prod2;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_prod3;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_ambi_channel_weight; + int *x_delta; + int *x_phi; + int *x_phi_sym; + int *x_sym_flag; + BIN_AMBI_COMPLEX *x_spec; + BIN_AMBI_COMPLEX *x_sin_cos; + float *x_beg_fade_out_hrir; + float *x_beg_hrir; + float **x_beg_hrtf_re; + float **x_beg_hrtf_im; + t_symbol **x_hrir_filename; + t_symbol **x_s_hrir; + t_symbol **x_s_hrtf_re; + t_symbol **x_s_hrtf_im; + t_symbol *x_s_fade_out_hrir; + void *x_out_sign_sum; + 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; + double x_sing_range; +} t_bin_ambi_reduced_decode_fft2; + +static t_class *bin_ambi_reduced_decode_fft2_class; + +static void bin_ambi_reduced_decode_fft2_init_cos(t_bin_ambi_reduced_decode_fft2 *x) +{/* initialize a whole sine and cosine wave over a fftsize array */ + int i, fftsize = x->x_fftsize; + float f, g; + BIN_AMBI_COMPLEX *sincos = x->x_sin_cos; + + g = 2.0f * 3.1415926538f / (float)fftsize; + for(i=0; i<fftsize; i++) + { + f = g * (float)i; + (*sincos).real = cos(f); + (*sincos).imag = -sin(f); /*FFT*/ + sincos++; + } +} + +static void bin_ambi_reduced_decode_fft2_convert(t_bin_ambi_reduced_decode_fft2 *x, double *delta_deg2rad, double *phi_deg2rad, int index) +{ + double q = 1.0; + double d = *delta_deg2rad; + double p = *phi_deg2rad; + int i; + + if(d < -90.0) + d = -90.0; + if(d > 90.0) + d = 90.0; + while(p < 0.0) + p += 360.0; + while(p >= 360.0) + p -= 360.0; + + x->x_delta[index] = (int)(d); + x->x_phi[index] = (int)(p); + + *delta_deg2rad = d*x->x_pi_over_180; + *phi_deg2rad = p*x->x_pi_over_180; +} + +static void bin_ambi_reduced_decode_fft2_do_2d(t_bin_ambi_reduced_decode_fft2 *x, int argc, t_atom *argv, int mode) +{ + double delta=0.0, phi; + double *dw = x->x_transp; + int index; + int order=x->x_n_order; + + if(argc < 2) + { + post("bin_ambi_reduced_decode_fft2 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 == BIN_AMBI_LS_REAL) + { + if(index >= x->x_n_real_ls) + index = x->x_n_real_ls - 1; + } + else if(mode == BIN_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; + + bin_ambi_reduced_decode_fft2_convert(x, &delta, &phi, index); + + 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 bin_ambi_reduced_decode_fft2_do_3d(t_bin_ambi_reduced_decode_fft2 *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("bin_ambi_reduced_decode_fft2 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 == BIN_AMBI_LS_REAL) + { + if(index >= x->x_n_real_ls) + index = x->x_n_real_ls - 1; + } + else if(mode == BIN_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; + + bin_ambi_reduced_decode_fft2_convert(x, &delta, &phi, index); + + cd = cos(delta); + sd = sin(delta); + cp = cos(phi); + sp = sin(phi); + + dw += index * x->x_n_ambi; + + *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 bin_ambi_reduced_decode_fft2_real_ls(t_bin_ambi_reduced_decode_fft2 *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fft2_do_2d(x, argc, argv, BIN_AMBI_LS_REAL); + else + bin_ambi_reduced_decode_fft2_do_3d(x, argc, argv, BIN_AMBI_LS_REAL); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode_fft2_pht_ls(t_bin_ambi_reduced_decode_fft2 *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fft2_do_2d(x, argc, argv, BIN_AMBI_LS_PHT); + else + bin_ambi_reduced_decode_fft2_do_3d(x, argc, argv, BIN_AMBI_LS_PHT); +} + +static void bin_ambi_reduced_decode_fft2_copy_row2buf(t_bin_ambi_reduced_decode_fft2 *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 bin_ambi_reduced_decode_fft2_copy_buf2row(t_bin_ambi_reduced_decode_fft2 *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 bin_ambi_reduced_decode_fft2_copy_row2row(t_bin_ambi_reduced_decode_fft2 *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 bin_ambi_reduced_decode_fft2_xch_rows(t_bin_ambi_reduced_decode_fft2 *x, int row1, int row2) +{ + bin_ambi_reduced_decode_fft2_copy_row2buf(x, row1); + bin_ambi_reduced_decode_fft2_copy_row2row(x, row2, row1); + bin_ambi_reduced_decode_fft2_copy_buf2row(x, row2); +} + +static void bin_ambi_reduced_decode_fft2_mul_row(t_bin_ambi_reduced_decode_fft2 *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 bin_ambi_reduced_decode_fft2_mul_col(t_bin_ambi_reduced_decode_fft2 *x, int col, double mul) +{ + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int i; + double *dw=x->x_inv_work2; + + dw += col; + for(i=0; i<n_ambi; i++) + { + (*dw) *= mul; + dw += n_ambi2; + } +} + +static void bin_ambi_reduced_decode_fft2_mul_buf_and_add2row(t_bin_ambi_reduced_decode_fft2 *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 bin_ambi_reduced_decode_fft2_eval_which_element_of_col_not_zero(t_bin_ambi_reduced_decode_fft2 *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 bin_ambi_reduced_decode_fft2_mul1(t_bin_ambi_reduced_decode_fft2 *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 bin_ambi_reduced_decode_fft2_mul2(t_bin_ambi_reduced_decode_fft2 *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_prod2; + 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; + } +} + +/* wird ersetzt durch haendische spiegelung und reduzierung + +static void bin_ambi_reduced_decode_fft2_mul3(t_bin_ambi_reduced_decode_fft2 *x) +{ + int i, n; + double *dv3=x->x_prod3; + double *dv2=x->x_prod2; + double *dv1=x->x_prod2; + double mw=x->x_mirror_weight; + + n = x->x_n_real_ls * x->x_n_ambi; + for(i=0; i<n; i++) + *dv3++ = *dv1++; + dv2 += n; + n = x->x_n_mrg_mir_ls * x->x_n_ambi; + dv2 += n; + for(i=0; i<n; i++) + *dv3++ = *dv1++ + (*dv2++)*mw; +}*/ + +static void bin_ambi_reduced_decode_fft2_transp_back(t_bin_ambi_reduced_decode_fft2 *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; + + transp = x->x_transp; + for(j=0; j<n_ambi; j++) + { + vec = transp; + for(i=0; i<n_ls; i++) + { + *straight++ = *vec; + vec += n_ambi; + } + transp++; + } +} + +static void bin_ambi_reduced_decode_fft2_inverse(t_bin_ambi_reduced_decode_fft2 *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 = bin_ambi_reduced_decode_fft2_eval_which_element_of_col_not_zero(x, i, i); + if(nz < 0) + { + post("bin_ambi_reduced_decode_fft2 ERROR: matrix not regular !!!!"); + x->x_seq_ok = 0; + return; + } + else + { + if(nz != i) + bin_ambi_reduced_decode_fft2_xch_rows(x, i, nz); + dv = db + i*n_ambi2 + i; + rcp = 1.0 /(*dv); + bin_ambi_reduced_decode_fft2_mul_row(x, i, rcp); + bin_ambi_reduced_decode_fft2_copy_row2buf(x, i); + for(j=i+1; j<n_ambi; j++) + { + dv += n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode_fft2_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; + bin_ambi_reduced_decode_fft2_copy_row2buf(x, i); + for(j=i-1; j>=0; j--) + { + dv -= n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode_fft2_mul_buf_and_add2row(x, j, rcp); + } + } + + post("matrix_inverse regular"); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode_fft2_calc_pinv(t_bin_ambi_reduced_decode_fft2 *x) +{ + t_garray *a; + int npoints; + t_float *fadevec; + int i, n; + double *dv3=x->x_prod3; + double *dv2=x->x_prod2; + + if((int)(x->x_beg_fade_out_hrir) == 0) + { + if (!(a = (t_garray *)pd_findbyclass(x->x_s_fade_out_hrir, garray_class))) + error("%s: no such array", x->x_s_fade_out_hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &fadevec)) + error("%s: bad template for bin_ambi_reduced_decode_fft2", x->x_s_fade_out_hrir->s_name); + else if (npoints < x->x_fftsize) + error("%s: bad array-size: %d", x->x_s_fade_out_hrir->s_name, npoints); + else + x->x_beg_fade_out_hrir = fadevec; + } + + bin_ambi_reduced_decode_fft2_transp_back(x); + bin_ambi_reduced_decode_fft2_mul1(x); + bin_ambi_reduced_decode_fft2_inverse(x); + bin_ambi_reduced_decode_fft2_mul2(x); + + n = x->x_n_real_ls * x->x_n_ambi; + for(i=0; i<n; i++) + *dv3++ = *dv2++; +} + +/* +x_prod: +n_ambi columns; +n_ambi rows; +*/ + +static void bin_ambi_reduced_decode_fft2_ipht_ireal_muladd(t_bin_ambi_reduced_decode_fft2 *x, t_symbol *s, int argc, t_atom *argv) +{ + int n = x->x_n_ambi; + int i, pht_index, real_index; + double *dv3=x->x_prod3; + double *dv2=x->x_prod2; + float mw; + + if(argc < 3) + { + post("bin_ambi_reduced_decode_fft2 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; + + dv3 += real_index*x->x_n_ambi; + dv2 += (x->x_n_real_ls+pht_index)*x->x_n_ambi; + for(i=0; i<n; i++) + { + *dv3 += (*dv2++)*mw; + dv3++; + } +} + +static void bin_ambi_reduced_decode_fft2_load_HRIR(t_bin_ambi_reduced_decode_fft2 *x, t_symbol *s, int argc, t_atom *argv) +{ + int index; + t_symbol *hrirname; + + if(argc < 2) + { + post("bin_ambi_reduced_decode_fft2 ERROR: load_HRIR needs 1 index and 1 HRIR-wav"); + return; + } + index = (int)atom_getint(argv++) - 1; + hrirname = atom_getsymbol(argv); + if(index < 0) + index = 0; + if(index >= x->x_n_real_ls) + index = x->x_n_real_ls - 1; + x->x_hrir_filename[index] = hrirname; + + SETSYMBOL(x->x_at, x->x_hrir_filename[index]); + SETSYMBOL(x->x_at+1, x->x_s_hrir[index]); + outlet_list(x->x_obj.ob_outlet, &s_list, 2, x->x_at); +} + +static void bin_ambi_reduced_decode_fft2_check_HRIR_arrays(t_bin_ambi_reduced_decode_fft2 *x, float findex) +{ + int index=(int)findex - 1; + int j, k, n; + int fftsize = x->x_fftsize; + int fs2=fftsize/2; + t_garray *a; + int npoints; + t_symbol *hrir; + t_float *vec_hrir, *vec, *vec_fade_out_hrir; + float decr, sum; + + if(index < 0) + index = 0; + if(index >= x->x_n_real_ls) + index = x->x_n_real_ls - 1; + + hrir = x->x_s_hrir[index]; + if (!(a = (t_garray *)pd_findbyclass(hrir, garray_class))) + error("%s: no such array", hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrir)) + error("%s: bad template for bin_ambi_reduced_decode_fft2", hrir->s_name); + else + { + if(npoints < fftsize) + { + post("bin_ambi_reduced_decode_fft2-WARNING: %s-array-size: %d < FFT-size: %d", hrir->s_name, npoints, fftsize); + } + vec = x->x_beg_hrir; + vec += index * fftsize; + + if((int)(x->x_beg_fade_out_hrir)) + { + vec_fade_out_hrir = x->x_beg_fade_out_hrir; + for(j=0; j<fs2; j++) + vec[j] = vec_hrir[j]*vec_fade_out_hrir[j]; + } + else + { + post("no HRIR-fade-out-window found"); + n = fs2 * 3; + n /= 4; + for(j=0; j<n; j++) + vec[j] = vec_hrir[j]; + sum = 1.0f; + decr = 4.0f / (float)fs2; + for(j=n, k=0; j<fs2; j++, k++) + { + sum -= decr; + vec[j] = vec_hrir[j] * sum; + } + } + } +} + +static void bin_ambi_reduced_decode_fft2_check_HRTF_arrays(t_bin_ambi_reduced_decode_fft2 *x, float findex) +{ + int index=(int)findex - 1; + t_garray *a; + int npoints; + int fftsize = x->x_fftsize; + t_float *vec_hrtf_re, *vec_hrtf_im; + t_symbol *hrtf_re, *hrtf_im; + + if(index < 0) + index = 0; + if(index >= x->x_n_ambi) + index = x->x_n_ambi - 1; + + hrtf_re = x->x_s_hrtf_re[index]; + hrtf_im = x->x_s_hrtf_im[index]; + + if (!(a = (t_garray *)pd_findbyclass(hrtf_re, garray_class))) + error("%s: no such array", hrtf_re->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrtf_re)) + error("%s: bad template for bin_ambi_reduced_decode_fft2", hrtf_re->s_name); + else if (npoints < fftsize) + error("%s: bad array-size: %d", hrtf_re->s_name, npoints); + else if (!(a = (t_garray *)pd_findbyclass(hrtf_im, garray_class))) + error("%s: no such array", hrtf_im->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrtf_im)) + error("%s: bad template for bin_ambi_reduced_decode_fft2", hrtf_im->s_name); + else if (npoints < fftsize) + error("%s: bad array-size: %d", hrtf_im->s_name, npoints); + else + { + x->x_beg_hrtf_re[index] = vec_hrtf_re; + x->x_beg_hrtf_im[index] = vec_hrtf_im; + } +} + +static void bin_ambi_reduced_decode_fft2_calc_reduced(t_bin_ambi_reduced_decode_fft2 *x, float findex) +{ + int index=(int)findex - 1; + int i, j, k, w_index, w_inc, i_inc, v_index, fs1, fs2; + int fftsize = x->x_fftsize; + BIN_AMBI_COMPLEX old1, old2, w; + BIN_AMBI_COMPLEX *sincos = x->x_sin_cos; + BIN_AMBI_COMPLEX *val = x->x_spec;/*weighted decoder-matrix with n_ls lines and n_ambi columns*/ + t_float *vec_hrir, *vec_hrtf_re, *vec_hrtf_im; + double *dv, *db=x->x_prod3; + int n_ambi = x->x_n_ambi; + int n_ls = x->x_n_real_ls; + float mul; + + if(x->x_seq_ok) + { + if(index < 0) + index = 0; + if(index >= n_ambi) + index = n_ambi - 1; + + vec_hrtf_re = x->x_beg_hrtf_re[index]; + vec_hrtf_im = x->x_beg_hrtf_im[index]; + + dv = db + index; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + for(k=0; k<fftsize; k++)/*first step of acumulating the HRIRs*/ + { + val[k].real = mul * vec_hrir[k]; + val[k].imag = 0.0f; + } + + for(j=1; j<n_ls; j++) + { + dv += n_ambi; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + vec_hrir += j * fftsize; + for(k=0; k<fftsize; k++) + { + val[k].real += mul * vec_hrir[k]; + } + } + + fs1 = fftsize - 1; + fs2 = fftsize / 2; + i_inc = fs2; + w_inc = 1; + for(i=1; i<fftsize; i<<=1) + { + v_index = 0; + for(j=0; j<i; j++) + { + w_index = 0; + for(k=0; k<i_inc; k++) + { + old1 = val[v_index]; + old2 = val[v_index+i_inc]; + w = sincos[w_index]; + val[v_index+i_inc].real = (old1.real-old2.real)*w.real - (old1.imag-old2.imag)*w.imag; + val[v_index+i_inc].imag = (old1.imag - old2.imag)*w.real + (old1.real-old2.real)*w.imag; + val[v_index].real = old1.real+old2.real; + val[v_index].imag = old1.imag+old2.imag; + w_index += w_inc; + v_index++; + } + v_index += i_inc; + } + w_inc <<= 1; + i_inc >>= 1; + } + + j = 0; + for(i=1;i<fs1;i++) + { + k = fs2; + while(k <= j) + { + j = j - k; + k >>= 1; + } + j = j + k; + if(i < j) + { + old1 = val[j]; + val[j] = val[i]; + val[i] = old1; + } + } + for(i = 0; i<fs2; i++) + { + vec_hrtf_re[i] = val[i].real; + vec_hrtf_im[i] = val[i].imag; + } + vec_hrtf_re[fs2] = val[fs2].real; + vec_hrtf_im[fs2] = 0.0f; + for(i = fs2+1; i < fftsize; i++) + { + vec_hrtf_re[i] = 0.0f; + vec_hrtf_im[i] = 0.0f; + } + } +} + +static void bin_ambi_reduced_decode_fft2_calc_sym(t_bin_ambi_reduced_decode_fft2 *x) +{ + int *sym=x->x_phi_sym; + int n_ls = x->x_n_real_ls; + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int *phi=x->x_phi; + int *delta=x->x_delta; + int *flag=x->x_sym_flag; + int i, j, d, p, regular=1; + double *dv, *db=x->x_prod3; + double a, b, q; + int sym_max=0; + int pos_sym_counter, neg_sym_counter; + char plus_minus[100]; + char abc[100]; + char onenine[100]; + char ten[100]; + + if(x->x_seq_ok) + { + for(i=0; i<n_ls; i++) + sym[i] = -1; + + plus_minus[0] = 0; + if(x->x_n_dim == 2) + { + strcpy(abc, " WXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXY"); + strcpy(ten, " 000000000000000000011111111111111111111"); + strcpy(onenine, " 011223344556677889900112233445566778899"); + abc[n_ambi+22] = 0; + ten[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(ten); + post(onenine); + } + else + { + strcpy(abc, " AABCABCDEABCDEFGABCDEFGHIABCDEFGHIJKABCDEFGHIJKLM"); + strcpy(onenine, " 0111222223333333444444444555555555556666666666666"); + abc[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(onenine); + } + + for(i=0; i<n_ls; i++)/* looking for azimuth symmetries of loudspeakers */ + { + d = delta[i]; + p = phi[i]; + for(j=i+1; j<n_ls; j++) + { + if(d == delta[j]) + { + if((p + phi[j]) == 360) + { + if((sym[i] < 0) && (sym[j] < 0)) + { + sym[i] = j; + sym[j] = i; + j = n_ls + 1; + sym_max++; + } + } + } + } + } + + for(p=0; p<n_ambi; p++)/*each col*/ + { + pos_sym_counter = 0; + neg_sym_counter = 0; + for(i=0; i<n_ls; i++) + flag[i] = 1; + dv = db + p; + for(i=0; i<n_ls; i++) + { + if((sym[i] >= 0) && flag[i]) + { + j = sym[i]; + flag[i] = 0; + flag[j] = 0; + a = dv[i*n_ambi]; + b = dv[j*n_ambi]; + if((a < 5.0e-4)&&(a > -5.0e-4)&&(b < 5.0e-4)&&(b > -5.0e-4)) + { + pos_sym_counter++; + neg_sym_counter++; + } + else + { + q = a / b; + if((q < 1.005)&&(q > 0.995)) + pos_sym_counter++; + else if((q > -1.005)&&(q < -0.995)) + neg_sym_counter++; + } + } + } + if(pos_sym_counter == sym_max) + strcat(plus_minus, "+"); + else if(neg_sym_counter == sym_max) + strcat(plus_minus, "-"); + else + { + strcat(plus_minus, "?"); + regular = 0; + } + } + post("sum of right channel: %s", plus_minus); + + if(regular) + { + for(i=0; i<n_ambi; i++) + { + /*change*/ + if(plus_minus[i] == '+') + SETFLOAT(x->x_at, 1.0f); + else if(plus_minus[i] == '-') + SETFLOAT(x->x_at, 2.0f); + SETFLOAT(x->x_at+1, (float)(i+1)); + outlet_list(x->x_out_sign_sum, &s_list, 2, x->x_at); + } + } + } +} + +static void bin_ambi_reduced_decode_fft2_ambi_weight(t_bin_ambi_reduced_decode_fft2 *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] = (double)atom_getfloat(argv++); + k++; + if(x->x_n_dim == 2) + { + for(i=1; i<=n; i++) + { + d = (double)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 = (double)atom_getfloat(argv++); + m = 2*i + 1; + for(j=0; j<m; j++) + { + x->x_ambi_channel_weight[k] = d; + k++; + } + } + } + } + else + post("bin_ambi_reduced_decode_fft2-ERROR: ambi_weight needs %d float weights", x->x_n_order+1); +} + +static void bin_ambi_reduced_decode_fft2_sing_range(t_bin_ambi_reduced_decode_fft2 *x, t_floatarg f) +{ + if(f < 0.0f) + x->x_sing_range = -(double)f; + else + x->x_sing_range = (double)f; +} + +static void bin_ambi_reduced_decode_fft2_free(t_bin_ambi_reduced_decode_fft2 *x) +{ + freebytes(x->x_hrir_filename, x->x_n_real_ls * sizeof(t_symbol *)); + freebytes(x->x_s_hrir, x->x_n_real_ls * sizeof(t_symbol *)); + freebytes(x->x_s_hrtf_re, x->x_n_ambi * sizeof(t_symbol *)); + freebytes(x->x_s_hrtf_im, x->x_n_ambi * sizeof(t_symbol *)); + 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_prod2, (x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_prod3, x->x_n_real_ls * x->x_n_ambi * sizeof(double)); + freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double)); + + freebytes(x->x_delta, (x->x_n_real_ls+x->x_n_pht_ls) * sizeof(int)); + freebytes(x->x_phi, (x->x_n_real_ls+x->x_n_pht_ls) * sizeof(int)); + freebytes(x->x_phi_sym, x->x_n_real_ls * sizeof(int)); + freebytes(x->x_sym_flag, x->x_n_real_ls * sizeof(int)); + + freebytes(x->x_spec, x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + freebytes(x->x_sin_cos, x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + + freebytes(x->x_beg_hrir, x->x_fftsize * x->x_n_real_ls * sizeof(float)); + freebytes(x->x_beg_hrtf_re, x->x_n_ambi * sizeof(float *)); + freebytes(x->x_beg_hrtf_im, x->x_n_ambi * sizeof(float *)); +} + +/* + 1.arg_ int prefix; + 2.arg: t_symbol *hrir_name; + 3.arg: t_symbol *hrtf_re_name; + 4.arg: t_symbol *hrtf_im_name; + 5.arg: t_symbol *hrir_fade_out_name; + 6.arg: int ambi_order; + 7.arg: int dim; + 8.arg: int number_of_loudspeakers; + 9.arg: int fftsize; + */ + +static void *bin_ambi_reduced_decode_fft2_new(t_symbol *s, int argc, t_atom *argv) +{ + t_bin_ambi_reduced_decode_fft2 *x = (t_bin_ambi_reduced_decode_fft2 *)pd_new(bin_ambi_reduced_decode_fft2_class); + char buf[400]; + int i, j, fftok, ok=0; + int n_order=0, n_dim=0, n_real_ls=0, n_pht_ls=0, n_ambi=0, fftsize=0, prefix=0; + t_symbol *s_hrir=gensym("L_HRIR"); + t_symbol *s_hrtf_re=gensym("HRTF_re"); + t_symbol *s_hrtf_im=gensym("HRTF_im"); + t_symbol *s_fade_out_hrir=gensym("HRIR_win"); + + if((argc >= 10) && + IS_A_FLOAT(argv,0) && + IS_A_SYMBOL(argv,1) && + IS_A_SYMBOL(argv,2) && + IS_A_SYMBOL(argv,3) && + IS_A_SYMBOL(argv,4) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8) && + IS_A_FLOAT(argv,9)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = (t_symbol *)atom_getsymbolarg(1, argc, argv); + s_hrtf_re = (t_symbol *)atom_getsymbolarg(2, argc, argv); + s_hrtf_im = (t_symbol *)atom_getsymbolarg(3, argc, argv); + s_fade_out_hrir = (t_symbol *)atom_getsymbolarg(4, argc, argv); + + n_order = (int)atom_getintarg(5, argc, argv); + n_dim = (int)atom_getintarg(6, argc, argv); + n_real_ls = (int)atom_getintarg(7, argc, argv); + n_pht_ls = (int)atom_getintarg(8, argc, argv); + fftsize = (int)atom_getintarg(9, argc, argv); + + ok = 1; + } + else if((argc >= 10) && + 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) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8) && + IS_A_FLOAT(argv,9)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = gensym("L_HRIR"); + s_hrtf_re = gensym("HRTF_re"); + s_hrtf_im = gensym("HRTF_im"); + s_fade_out_hrir = gensym("HRIR_win"); + + n_order = (int)atom_getintarg(5, argc, argv); + n_dim = (int)atom_getintarg(6, argc, argv); + n_real_ls = (int)atom_getintarg(7, argc, argv); + n_pht_ls = (int)atom_getintarg(8, argc, argv); + fftsize = (int)atom_getintarg(9, argc, argv); + + ok = 1; + } + + if(ok) + { + if(n_order < 1) + n_order = 1; + + if(n_dim == 3) + { + if(n_order > 5) + n_order = 5; + n_ambi = (n_order + 1)*(n_order + 1); + } + else + { + if(n_order > 12) + n_order = 12; + n_dim = 2; + n_ambi = 2 * n_order + 1; + } + + 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) < n_ambi) + { + post("bin_ambi_reduced_decode_fft2-WARNING: Number of all Loudspeakers < Number of Ambisonic-Channels !!!!"); + } + + j = 2; + fftok = 0; + for(i=0; i<21; i++) + { + if(j == fftsize) + { + fftok = 1; + i = 22; + } + j *= 2; + } + + if(!fftok) + { + fftsize = 512; + post("bin_ambi_reduced_decode_fft2-WARNING: fftsize not equal to 2 ^ n !!!"); + post(" fftsize set to %d", fftsize); + } + + x->x_n_dim = n_dim; + x->x_n_ambi = n_ambi; + x->x_n_real_ls = n_real_ls; + x->x_n_pht_ls = n_pht_ls; + x->x_n_order = n_order; + x->x_fftsize = fftsize; + + x->x_hrir_filename = (t_symbol **)getbytes(x->x_n_real_ls * sizeof(t_symbol *)); + x->x_s_hrir = (t_symbol **)getbytes(x->x_n_real_ls * sizeof(t_symbol *)); + x->x_s_hrtf_re = (t_symbol **)getbytes(x->x_n_ambi * sizeof(t_symbol *)); + x->x_s_hrtf_im = (t_symbol **)getbytes(x->x_n_ambi * sizeof(t_symbol *)); + + j = x->x_n_real_ls; + for(i=0; i<j; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrir->s_name); + x->x_s_hrir[i] = gensym(buf); + } + + for(i=0; i<n_ambi; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrtf_re->s_name); + x->x_s_hrtf_re[i] = gensym(buf); + + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrtf_im->s_name); + x->x_s_hrtf_im[i] = gensym(buf); + } + + sprintf(buf, "%d_%s", prefix, s_fade_out_hrir->s_name); + x->x_s_fade_out_hrir = gensym(buf); + + 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_prod2 = (double *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double)); + x->x_prod3 = (double *)getbytes(x->x_n_real_ls * x->x_n_ambi * sizeof(double)); + x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double)); + + x->x_delta = (int *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * sizeof(int)); + x->x_phi = (int *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * sizeof(int)); + x->x_phi_sym = (int *)getbytes(x->x_n_real_ls * sizeof(int)); + x->x_sym_flag = (int *)getbytes(x->x_n_real_ls * sizeof(int)); + + x->x_spec = (BIN_AMBI_COMPLEX *)getbytes(x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + x->x_sin_cos = (BIN_AMBI_COMPLEX *)getbytes(x->x_fftsize * sizeof(BIN_AMBI_COMPLEX)); + + x->x_beg_fade_out_hrir = (float *)0; + x->x_beg_hrir = (float *)getbytes(x->x_fftsize * x->x_n_real_ls * sizeof(float)); + x->x_beg_hrtf_re = (float **)getbytes(x->x_n_ambi * sizeof(float *)); + x->x_beg_hrtf_im = (float **)getbytes(x->x_n_ambi * sizeof(float *)); + + 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; + x->x_seq_ok = 1; + for(i=0; i<n_ambi; i++) + x->x_ambi_channel_weight[i] = 1.0; + + bin_ambi_reduced_decode_fft2_init_cos(x); + + outlet_new(&x->x_obj, &s_list); + x->x_out_sign_sum = outlet_new(&x->x_obj, &s_list); + return(x); + } + else + { + post("bin_ambi_reduced_decode_fft2-ERROR: need 1 float + 4 symbols + 5 floats arguments:"); + post(" prefix(unique-number) + hrir_name + hrtf_re_name + hrtf_im_name + hrir_fade_out_name +"); + post(" + ambi_order + ambi_dimension + number_of_real_loudspeakers + "); + post(" + number_of_phantom_loudspeakers + fftsize"); + return(0); + } +} + +void bin_ambi_reduced_decode_fft2_setup(void) +{ + bin_ambi_reduced_decode_fft2_class = class_new(gensym("bin_ambi_reduced_decode_fft2"), (t_newmethod)bin_ambi_reduced_decode_fft2_new, (t_method)bin_ambi_reduced_decode_fft2_free, + sizeof(t_bin_ambi_reduced_decode_fft2), 0, A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_real_ls, gensym("real_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_pht_ls, gensym("pht_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_calc_pinv, gensym("calc_pinv"), 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_ipht_ireal_muladd, gensym("ipht_ireal_muladd"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_load_HRIR, gensym("load_HRIR"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_check_HRIR_arrays, gensym("check_HRIR_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_check_HRTF_arrays, gensym("check_HRTF_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_calc_reduced, gensym("calc_reduced"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_calc_sym, gensym("calc_sym"), 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_ambi_weight, gensym("ambi_weight"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fft2_class, (t_method)bin_ambi_reduced_decode_fft2_sing_range, gensym("sing_range"), A_DEFFLOAT, 0); + class_sethelpsymbol(bin_ambi_reduced_decode_fft2_class, gensym("iemhelp2/help-bin_ambi_reduced_decode_fft2")); +} +/* +Reihenfolge: +n_ls x bin_ambi_reduced_decode_fft2_ls + +bin_ambi_reduced_decode_fft2_calc_pinv + +n_ls x bin_ambi_reduced_decode_fft2_load_HRIR + +n_ls x bin_ambi_reduced_decode_fft2_check_HRIR_arrays + +n_ambi x bin_ambi_reduced_decode_fft2_check_HRTF_arrays + +n_ambi x bin_ambi_reduced_decode_fft2_calc_reduced + +bin_ambi_reduced_decode_fft2_calc_sym + +*/ diff --git a/src/bin_ambi_reduced_decode_fir.c b/src/bin_ambi_reduced_decode_fir.c new file mode 100644 index 0000000..e935f04 --- /dev/null +++ b/src/bin_ambi_reduced_decode_fir.c @@ -0,0 +1,1385 @@ +/* For information on usage and redistribution, and for a DISCLAIMER OF ALL +* WARRANTIES, see the file, "LICENSE.txt," in this distribution. + +iem_bin_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_bin_ambi.h" +#include <math.h> +#include <stdio.h> +#include <string.h> + + +/* -------------------------- bin_ambi_reduced_decode_fir ------------------------------ */ +/* + ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren ** + ** Inputs: ls + Liste von 3 floats: Index [1 .. 16] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] ** + ** Inputs: calc_inv ** + ** Inputs: load_HRIR + float index1..16 ** + ** 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 1016_1_L_HRIR .. 1016_16_L_HRIR heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_re .. 1016_16_HRTF_re heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_im .. 1016_16_HRTF_im heissen und existieren ** + */ + +typedef struct _bin_ambi_reduced_decode_fir +{ + t_object x_obj; + t_atom x_at[2];/*output filename and tablename of HRIR*/ + int x_n_dim;/*dimension of ambisonic system*/ + int x_n_ambi;/*number of ambi channels*/ + int x_n_order;/*order of ambisonic system*/ + int x_n_ind_ls;/*number of independent loudspeakers*/ + int x_n_mrg_mir_ls;/*number of mirrored and merged loudspeakers*/ + int x_n_ph_ls;/*number of phantom loudspeakers*/ + int x_seq_ok; + int x_firsize;/*size of FIR-filter*/ + double *x_inv_work1;/* n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_work2;/* 2*n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_buf2;/* 2*n_ambi buffer for matrix inverse */ + double *x_transp;/* n_ls*n_ambi transposed input-buffer of decoder-matrix before inversion */ + double *x_ls_encode;/* n_ambi*n_ls straight input-buffer of decoder-matrix before inversion */ + double *x_prod2;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_prod3;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_ambi_channel_weight; + double x_mirror_weight; + int *x_delta; + int *x_phi; + int *x_phi_sym; + int *x_sym_flag; + float *x_beg_fade_out_hrir; + float *x_beg_hrir; + float **x_beg_hrir_red; + t_symbol **x_hrir_filename; + t_symbol **x_s_hrir; + t_symbol **x_s_hrir_red; + t_symbol *x_s_fade_out_hrir; + void *x_out_sign_sum; + 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; + double x_sing_range; +} t_bin_ambi_reduced_decode_fir; + +static t_class *bin_ambi_reduced_decode_fir_class; + +static void bin_ambi_reduced_decode_fir_quant(t_bin_ambi_reduced_decode_fir *x, double *delta_deg2rad, double *phi_deg2rad, int index) +{ + double q = 1.0; + double d = *delta_deg2rad; + double p = *phi_deg2rad; + int i; + + if(d < -40.0) + d = -40.0; + if(d > 90.0) + d = 90.0; + while(p < 0.0) + p += 360.0; + while(p >= 360.0) + p -= 360.0; + + if(d < -35.0) + { + *delta_deg2rad = -40.0; + q = 360.0 / 56.0; + } + else if(d < -25.0) + { + *delta_deg2rad = -30.0; + q = 6.0; + } + else if(d < -15.0) + { + *delta_deg2rad = -20.0; + q = 5.0; + } + else if(d < -5.0) + { + *delta_deg2rad = -10.0; + q = 5.0; + } + else if(d < 5.0) + { + *delta_deg2rad = 0.0; + q = 5.0; + } + else if(d < 15.0) + { + *delta_deg2rad = 10.0; + q = 5.0; + } + else if(d < 25.0) + { + *delta_deg2rad = 20.0; + q = 5.0; + } + else if(d < 35.0) + { + *delta_deg2rad = 30.0; + q = 6.0; + } + else if(d < 45.0) + { + *delta_deg2rad = 40.0; + q = 360.0 / 56.0; + } + else if(d < 55.0) + { + *delta_deg2rad = 50.0; + q = 8.0; + } + else if(d < 65.0) + { + *delta_deg2rad = 60.0; + q = 10.0; + } + else if(d < 75.0) + { + *delta_deg2rad = 70.0; + q = 15.0; + } + else if(d < 85.0) + { + *delta_deg2rad = 80.0; + q = 30.0; + } + else + { + *delta_deg2rad = 90.0; + q = 360.0; + } + + p /= q; + i = (int)(p + 0.499999); + p = (double)i * q; + i = (int)(p + 0.499999); + while(i >= 360) + i -= 360; + p = (double)i; + *phi_deg2rad = p; + + x->x_delta[index] = (int)(*delta_deg2rad); + x->x_phi[index] = i; + + *delta_deg2rad *= x->x_pi_over_180; + *phi_deg2rad *= x->x_pi_over_180; +} + +static void bin_ambi_reduced_decode_fir_do_2d(t_bin_ambi_reduced_decode_fir *x, int argc, t_atom *argv, int mode) +{ + double delta=0.0, phi; + double *dw = x->x_transp; + int index; + int order=x->x_n_order; + + if(argc < 2) + { + post("bin_ambi_reduced_decode_fir 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 == BIN_AMBI_LS_IND) + { + if(index >= x->x_n_ind_ls) + index = x->x_n_ind_ls - 1; + } + else if(mode == BIN_AMBI_LS_MRG) + { + if(x->x_n_mrg_mir_ls) + { + if(index >= x->x_n_mrg_mir_ls) + index = x->x_n_mrg_mir_ls - 1; + index += x->x_n_ind_ls; + } + else + return; + } + else if(mode == BIN_AMBI_LS_MIR) + { + if(x->x_n_mrg_mir_ls) + { + if(index >= x->x_n_mrg_mir_ls) + index = x->x_n_mrg_mir_ls - 1; + index += x->x_n_ind_ls; + index += x->x_n_mrg_mir_ls; + } + else + return; + } + else if(mode == BIN_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_ind_ls; + index += 2*x->x_n_mrg_mir_ls; + } + else + return; + } + else + return; + + bin_ambi_reduced_decode_fir_quant(x, &delta, &phi, index); + + 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 bin_ambi_reduced_decode_fir_do_3d(t_bin_ambi_reduced_decode_fir *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("bin_ambi_reduced_decode_fir 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 == BIN_AMBI_LS_IND) + { + if(index >= x->x_n_ind_ls) + index = x->x_n_ind_ls - 1; + } + else if(mode == BIN_AMBI_LS_MRG) + { + if(x->x_n_mrg_mir_ls) + { + if(index >= x->x_n_mrg_mir_ls) + index = x->x_n_mrg_mir_ls - 1; + index += x->x_n_ind_ls; + } + else + return; + } + else if(mode == BIN_AMBI_LS_MIR) + { + if(x->x_n_mrg_mir_ls) + { + if(index >= x->x_n_mrg_mir_ls) + index = x->x_n_mrg_mir_ls - 1; + index += x->x_n_ind_ls; + index += x->x_n_mrg_mir_ls; + } + else + return; + } + else if(mode == BIN_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_ind_ls; + index += 2*x->x_n_mrg_mir_ls; + } + else + return; + } + else + return; + + bin_ambi_reduced_decode_fir_quant(x, &delta, &phi, index); + + cd = cos(delta); + sd = sin(delta); + cp = cos(phi); + sp = sin(phi); + + dw += index * x->x_n_ambi; + + *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 bin_ambi_reduced_decode_fir_ind_ls(t_bin_ambi_reduced_decode_fir *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fir_do_2d(x, argc, argv, BIN_AMBI_LS_IND); + else + bin_ambi_reduced_decode_fir_do_3d(x, argc, argv, BIN_AMBI_LS_IND); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode_fir_mrg_ls(t_bin_ambi_reduced_decode_fir *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fir_do_2d(x, argc, argv, BIN_AMBI_LS_MRG); + else + bin_ambi_reduced_decode_fir_do_3d(x, argc, argv, BIN_AMBI_LS_MRG); +} + +static void bin_ambi_reduced_decode_fir_mir_ls(t_bin_ambi_reduced_decode_fir *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fir_do_2d(x, argc, argv, BIN_AMBI_LS_MIR); + else + bin_ambi_reduced_decode_fir_do_3d(x, argc, argv, BIN_AMBI_LS_MIR); +} + +static void bin_ambi_reduced_decode_fir_pht_ls(t_bin_ambi_reduced_decode_fir *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fir_do_2d(x, argc, argv, BIN_AMBI_LS_PHT); + else + bin_ambi_reduced_decode_fir_do_3d(x, argc, argv, BIN_AMBI_LS_PHT); +} + +static void bin_ambi_reduced_decode_fir_copy_row2buf(t_bin_ambi_reduced_decode_fir *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 bin_ambi_reduced_decode_fir_copy_buf2row(t_bin_ambi_reduced_decode_fir *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 bin_ambi_reduced_decode_fir_copy_row2row(t_bin_ambi_reduced_decode_fir *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 bin_ambi_reduced_decode_fir_xch_rows(t_bin_ambi_reduced_decode_fir *x, int row1, int row2) +{ + bin_ambi_reduced_decode_fir_copy_row2buf(x, row1); + bin_ambi_reduced_decode_fir_copy_row2row(x, row2, row1); + bin_ambi_reduced_decode_fir_copy_buf2row(x, row2); +} + +static void bin_ambi_reduced_decode_fir_mul_row(t_bin_ambi_reduced_decode_fir *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 bin_ambi_reduced_decode_fir_mul_col(t_bin_ambi_reduced_decode_fir *x, int col, double mul) +{ + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int i; + double *dw=x->x_inv_work2; + + dw += col; + for(i=0; i<n_ambi; i++) + { + (*dw) *= mul; + dw += n_ambi2; + } +} + +static void bin_ambi_reduced_decode_fir_mul_buf_and_add2row(t_bin_ambi_reduced_decode_fir *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 bin_ambi_reduced_decode_fir_eval_which_element_of_col_not_zero(t_bin_ambi_reduced_decode_fir *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 bin_ambi_reduced_decode_fir_mul1(t_bin_ambi_reduced_decode_fir *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_ind_ls+2*x->x_n_mrg_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 bin_ambi_reduced_decode_fir_mul2(t_bin_ambi_reduced_decode_fir *x) +{ + int n_ls=x->x_n_ind_ls+2*x->x_n_mrg_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_prod2; + 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 bin_ambi_reduced_decode_fir_mul3(t_bin_ambi_reduced_decode_fir *x) +{ + int i, n; + double *dv3=x->x_prod3; + double *dv2=x->x_prod2; + double *dv1=x->x_prod2; + double mw=x->x_mirror_weight; + + n = x->x_n_ind_ls * x->x_n_ambi; + for(i=0; i<n; i++) + *dv3++ = *dv1++; + dv2 += n; + n = x->x_n_mrg_mir_ls * x->x_n_ambi; + dv2 += n; + for(i=0; i<n; i++) + *dv3++ = *dv1++ + (*dv2++)*mw; +} + +static void bin_ambi_reduced_decode_fir_transp_back(t_bin_ambi_reduced_decode_fir *x) +{ + double *vec, *transp=x->x_transp; + double *straight=x->x_ls_encode; + int n_ls=x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls; + int n_ambi=x->x_n_ambi; + int i, j; + + transp = x->x_transp; + for(j=0; j<n_ambi; j++) + { + vec = transp; + for(i=0; i<n_ls; i++) + { + *straight++ = *vec; + vec += n_ambi; + } + transp++; + } +} + +static void bin_ambi_reduced_decode_fir_inverse(t_bin_ambi_reduced_decode_fir *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 = bin_ambi_reduced_decode_fir_eval_which_element_of_col_not_zero(x, i, i); + if(nz < 0) + { + post("bin_ambi_reduced_decode_fir ERROR: matrix not regular !!!!"); + x->x_seq_ok = 0; + return; + } + else + { + if(nz != i) + bin_ambi_reduced_decode_fir_xch_rows(x, i, nz); + dv = db + i*n_ambi2 + i; + rcp = 1.0 /(*dv); + bin_ambi_reduced_decode_fir_mul_row(x, i, rcp); + bin_ambi_reduced_decode_fir_copy_row2buf(x, i); + for(j=i+1; j<n_ambi; j++) + { + dv += n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode_fir_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; + bin_ambi_reduced_decode_fir_copy_row2buf(x, i); + for(j=i-1; j>=0; j--) + { + dv -= n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode_fir_mul_buf_and_add2row(x, j, rcp); + } + } + + post("matrix_inverse regular"); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode_fir_calc_pinv(t_bin_ambi_reduced_decode_fir *x) +{ + t_garray *a; + int npoints; + t_float *fadevec; + + if((int)(x->x_beg_fade_out_hrir) == 0) + { + if (!(a = (t_garray *)pd_findbyclass(x->x_s_fade_out_hrir, garray_class))) + error("%s: no such array", x->x_s_fade_out_hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &fadevec)) + error("%s: bad template for bin_ambi_reduced_decode_fir", x->x_s_fade_out_hrir->s_name); + else if (npoints < x->x_firsize) + error("%s: bad array-size: %d", x->x_s_fade_out_hrir->s_name, npoints); + else + x->x_beg_fade_out_hrir = fadevec; + } + + bin_ambi_reduced_decode_fir_transp_back(x); + bin_ambi_reduced_decode_fir_mul1(x); + bin_ambi_reduced_decode_fir_inverse(x); + bin_ambi_reduced_decode_fir_mul2(x); + bin_ambi_reduced_decode_fir_mul3(x); +} + +/* +x_prod: +n_ambi columns; +n_ambi rows; +*/ + +static void bin_ambi_reduced_decode_fir_load_HRIR(t_bin_ambi_reduced_decode_fir *x, float findex) +{ + int index=(int)findex - 1; + int p; + char buf[60]; + + if(index < 0) + index = 0; + if(index >= (x->x_n_ind_ls + x->x_n_mrg_mir_ls)) + index = x->x_n_ind_ls + x->x_n_mrg_mir_ls - 1; + + p = x->x_phi[index]; + + if(p)/*change*/ + p = 360 - p; + + if(p < 10) + sprintf(buf, "L%de00%da.wav", x->x_delta[index], p); + else if(p < 100) + sprintf(buf, "L%de0%da.wav", x->x_delta[index], p); + else + sprintf(buf, "L%de%da.wav", x->x_delta[index], p); + x->x_hrir_filename[index] = gensym(buf); + + SETSYMBOL(x->x_at, x->x_hrir_filename[index]); + SETSYMBOL(x->x_at+1, x->x_s_hrir[index]); + outlet_list(x->x_obj.ob_outlet, &s_list, 2, x->x_at); +} + +static void bin_ambi_reduced_decode_fir_check_HRIR_LS_arrays(t_bin_ambi_reduced_decode_fir *x, float findex) +{ + int index=(int)findex - 1; + int j, k, n; + int firsize = x->x_firsize; + int fs2=firsize/2; + t_garray *a; + int npoints; + t_symbol *hrir; + t_float *vec_hrir, *vec, *vec_fade_out_hrir; + float decr, sum; + + if(index < 0) + index = 0; + if(index >= (x->x_n_ind_ls + x->x_n_mrg_mir_ls)) + index = x->x_n_ind_ls + x->x_n_mrg_mir_ls - 1; + + hrir = x->x_s_hrir[index]; + if (!(a = (t_garray *)pd_findbyclass(hrir, garray_class))) + error("%s: no such array", hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrir)) + error("%s: bad template for bin_ambi_reduced_decode_fir", hrir->s_name); + else + { + if(npoints < firsize) + { + post("bin_ambi_reduced_decode_fir-WARNING: %s-array-size: %d < FFT-size: %d", hrir->s_name, npoints, firsize); + } + vec = x->x_beg_hrir; + vec += index * firsize; + + if((int)(x->x_beg_fade_out_hrir)) + { + vec_fade_out_hrir = x->x_beg_fade_out_hrir; + for(j=0; j<fs2; j++) + vec[j] = vec_hrir[j]*vec_fade_out_hrir[j]; + } + else + { + post("no HRIR-fade-out-window found"); + n = fs2 * 3; + n /= 4; + for(j=0; j<n; j++) + vec[j] = vec_hrir[j]; + sum = 1.0f; + decr = 4.0f / (float)fs2; + for(j=n, k=0; j<fs2; j++, k++) + { + sum -= decr; + vec[j] = vec_hrir[j] * sum; + } + } + } +} + +static void bin_ambi_reduced_decode_fir_check_HRIR_RED_arrays(t_bin_ambi_reduced_decode_fir *x, float findex) +{ + int index=(int)findex - 1; + t_garray *a; + int npoints; + int firsize = x->x_firsize; + t_float *vec_hrir_red; + t_symbol *hrir_red; + + if(index < 0) + index = 0; + if(index >= x->x_n_ambi) + index = x->x_n_ambi - 1; + + hrir_red = x->x_s_hrir_red[index]; + + if (!(a = (t_garray *)pd_findbyclass(hrir_red, garray_class))) + error("%s: no such array", hrir_red->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrir_red)) + error("%s: bad template for bin_ambi_reduced_decode_fir", hrir_red->s_name); + else if (npoints < firsize) + error("%s: bad array-size: %d", hrir_red->s_name, npoints); + else + { + x->x_beg_hrir_red[index] = vec_hrir_red; + } +} + +static void bin_ambi_reduced_decode_fir_calc_reduced(t_bin_ambi_reduced_decode_fir *x, float findex) +{ + int index=(int)findex - 1; + int i, j; + int firsize = x->x_firsize; + t_float *vec_hrir, *vec_hrir_red; + double *dv; + int n_ambi = x->x_n_ambi; + int n_ls = x->x_n_ind_ls + x->x_n_mrg_mir_ls; + float mul; + + if(x->x_seq_ok) + { + if(index < 0) + index = 0; + if(index >= n_ambi) + index = n_ambi - 1; + + vec_hrir_red = x->x_beg_hrir_red[index]; + + dv = x->x_prod3 + index; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + for(i=0; i<firsize; i++)/*first step of acumulating the HRIRs*/ + { + vec_hrir_red[i] = mul * vec_hrir[i]; + } + + for(j=1; j<n_ls; j++) + { + dv += n_ambi; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + vec_hrir += j * firsize; + for(i=0; i<firsize; i++) + { + vec_hrir_red[i] += mul * vec_hrir[i]; + } + } + } +} + +static void bin_ambi_reduced_decode_fir_calc_sym(t_bin_ambi_reduced_decode_fir *x) +{ + int *sym=x->x_phi_sym; + int n_ls=x->x_n_ind_ls+x->x_n_mrg_mir_ls; + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int *phi=x->x_phi; + int *delta=x->x_delta; + int *flag=x->x_sym_flag; + int i, j, d, p, regular=1; + double *dv, *db=x->x_prod3; + double a, b, q; + int sym_max=0; + int pos_sym_counter, neg_sym_counter; + char plus_minus[100]; + char abc[100]; + char onenine[100]; + char ten[100]; + + if(x->x_seq_ok) + { + for(i=0; i<n_ls; i++) + sym[i] = -1; + + plus_minus[0] = 0; + if(x->x_n_dim == 2) + { + strcpy(abc, " WXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXY"); + strcpy(ten, " 000000000000000000011111111111111111111"); + strcpy(onenine, " 011223344556677889900112233445566778899"); + abc[n_ambi+22] = 0; + ten[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(ten); + post(onenine); + } + else + { + strcpy(abc, " AABCABCDEABCDEFGABCDEFGHIABCDEFGHIJKABCDEFGHIJKLM"); + strcpy(onenine, " 0111222223333333444444444555555555556666666666666"); + abc[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(onenine); + } + + for(i=0; i<n_ls; i++)/* looking for azimuth symmetries of loudspeakers */ + { + d = delta[i]; + p = phi[i]; + for(j=i+1; j<n_ls; j++) + { + if(d == delta[j]) + { + if((p + phi[j]) == 360) + { + if((sym[i] < 0) && (sym[j] < 0)) + { + sym[i] = j; + sym[j] = i; + j = n_ls + 1; + sym_max++; + } + } + } + } + } + + for(p=0; p<n_ambi; p++)/*each col*/ + { + pos_sym_counter = 0; + neg_sym_counter = 0; + for(i=0; i<n_ls; i++) + flag[i] = 1; + dv = db + p; + for(i=0; i<n_ls; i++) + { + if((sym[i] >= 0) && flag[i]) + { + j = sym[i]; + flag[i] = 0; + flag[j] = 0; + a = dv[i*n_ambi]; + b = dv[j*n_ambi]; + if((a < 5.0e-4)&&(a > -5.0e-4)&&(b < 5.0e-4)&&(b > -5.0e-4)) + { + pos_sym_counter++; + neg_sym_counter++; + } + else + { + q = a / b; + if((q < 1.005)&&(q > 0.995)) + pos_sym_counter++; + else if((q > -1.005)&&(q < -0.995)) + neg_sym_counter++; + } + } + } + if(pos_sym_counter == sym_max) + strcat(plus_minus, "+"); + else if(neg_sym_counter == sym_max) + strcat(plus_minus, "-"); + else + { + strcat(plus_minus, "?"); + regular = 0; + } + } + post("sum of right channel: %s", plus_minus); + + if(regular) + { + for(i=0; i<n_ambi; i++) + { + /*change*/ + if(plus_minus[i] == '+') + SETFLOAT(x->x_at, 1.0f); + else if(plus_minus[i] == '-') + SETFLOAT(x->x_at, 2.0f); + SETFLOAT(x->x_at+1, (float)(i+1)); + outlet_list(x->x_out_sign_sum, &s_list, 2, x->x_at); + } + } + } +} + +static void bin_ambi_reduced_decode_fir_ambi_weight(t_bin_ambi_reduced_decode_fir *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] = (double)atom_getfloat(argv++); + k++; + if(x->x_n_dim == 2) + { + for(i=1; i<=n; i++) + { + d = (double)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 = (double)atom_getfloat(argv++); + m = 2*i + 1; + for(j=0; j<m; j++) + { + x->x_ambi_channel_weight[k] = d; + k++; + } + } + } + } + else + post("bin_ambi_reduced_decode_fir-ERROR: ambi_weight needs %d float weights", x->x_n_order+1); +} + +static void bin_ambi_reduced_decode_fir_mirror_weight(t_bin_ambi_reduced_decode_fir *x, t_floatarg f) +{ + x->x_mirror_weight = (double)f; +} + +static void bin_ambi_reduced_decode_fir_sing_range(t_bin_ambi_reduced_decode_fir *x, t_floatarg f) +{ + if(f < 0.0f) + x->x_sing_range = -(double)f; + else + x->x_sing_range = (double)f; +} + +static void bin_ambi_reduced_decode_fir_free(t_bin_ambi_reduced_decode_fir *x) +{ + freebytes(x->x_hrir_filename, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(t_symbol *)); + freebytes(x->x_s_hrir, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(t_symbol *)); + freebytes(x->x_s_hrir_red, x->x_n_ambi * sizeof(t_symbol *)); + 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_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_ls_encode, (x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_prod2, (x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_prod3, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double)); + + freebytes(x->x_delta, (x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * sizeof(int)); + freebytes(x->x_phi, (x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * sizeof(int)); + freebytes(x->x_phi_sym, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(int)); + freebytes(x->x_sym_flag, (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(int)); + + freebytes(x->x_beg_hrir, x->x_firsize * (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(float)); + freebytes(x->x_beg_hrir_red, x->x_n_ambi * sizeof(float *)); +} + +/* + 1.arg_ int prefix; + 2.arg: t_symbol *hrir_name; + 3.arg: t_symbol *hrir_red_name; + 4.arg: t_symbol *hrtf_im_name; + 5.arg: t_symbol *hrir_fade_out_name; + 6.arg: int ambi_order; + 7.arg: int dim; + 8.arg: int number_of_loudspeakers; + 9.arg: int firsize; + */ + +static void *bin_ambi_reduced_decode_fir_new(t_symbol *s, int argc, t_atom *argv) +{ + t_bin_ambi_reduced_decode_fir *x = (t_bin_ambi_reduced_decode_fir *)pd_new(bin_ambi_reduced_decode_fir_class); + char buf[400]; + int i, j, ok=0; + int n_order=0, n_dim=0, n_ind_ls=0, n_mrg_mir_ls=0, n_ph_ls=0, n_ambi=0, firsize=0, prefix=0; + t_symbol *s_hrir=gensym("L_HRIR"); + t_symbol *s_hrir_red=gensym("HRIR_red"); + t_symbol *s_fade_out_hrir=gensym("HRIR_win"); + + if((argc >= 10) && + IS_A_FLOAT(argv,0) && + IS_A_SYMBOL(argv,1) && + IS_A_SYMBOL(argv,2) && + IS_A_SYMBOL(argv,3) && + IS_A_FLOAT(argv,4) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8) && + IS_A_FLOAT(argv,9)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = (t_symbol *)atom_getsymbolarg(1, argc, argv); + s_hrir_red = (t_symbol *)atom_getsymbolarg(2, argc, argv); + s_fade_out_hrir = (t_symbol *)atom_getsymbolarg(3, argc, argv); + + n_order = (int)atom_getintarg(4, argc, argv); + n_dim = (int)atom_getintarg(5, argc, argv); + n_ind_ls = (int)atom_getintarg(6, argc, argv); + n_mrg_mir_ls = (int)atom_getintarg(7, argc, argv); + n_ph_ls = (int)atom_getintarg(8, argc, argv); + firsize = (int)atom_getintarg(9, argc, argv); + + ok = 1; + } + else if((argc >= 10) && + 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) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8) && + IS_A_FLOAT(argv,9)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = gensym("L_HRIR"); + s_hrir_red = gensym("HRIR_red"); + s_fade_out_hrir = gensym("HRIR_win"); + + n_order = (int)atom_getintarg(4, argc, argv); + n_dim = (int)atom_getintarg(5, argc, argv); + n_ind_ls = (int)atom_getintarg(6, argc, argv); + n_mrg_mir_ls = (int)atom_getintarg(7, argc, argv); + n_ph_ls = (int)atom_getintarg(8, argc, argv); + firsize = (int)atom_getintarg(9, argc, argv); + + ok = 1; + } + + if(ok) + { + if(n_order < 1) + n_order = 1; + + if(n_dim == 3) + { + if(n_order > 5) + n_order = 5; + n_ambi = (n_order + 1)*(n_order + 1); + } + else + { + if(n_order > 12) + n_order = 12; + n_dim = 2; + n_ambi = 2 * n_order + 1; + } + + if(n_ind_ls < 1) + n_ind_ls = 1; + if(n_mrg_mir_ls < 1) + n_mrg_mir_ls = 1; + if(n_ph_ls < 0) + n_ph_ls = 0; + + if((n_ind_ls+2*n_mrg_mir_ls+n_ph_ls) < n_ambi) + { + post("bin_ambi_reduced_decode_fir-WARNING: Number of all Loudspeakers < Number of Ambisonic-Channels !!!!"); + } + + if(firsize < 32) + firsize = 32; + + x->x_n_dim = n_dim; + x->x_n_ambi = n_ambi; + x->x_n_ind_ls = n_ind_ls; + x->x_n_mrg_mir_ls = n_mrg_mir_ls; + x->x_n_ph_ls = n_ph_ls; + x->x_n_order = n_order; + x->x_firsize = firsize; + + x->x_hrir_filename = (t_symbol **)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(t_symbol *)); + x->x_s_hrir = (t_symbol **)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(t_symbol *)); + x->x_s_hrir_red = (t_symbol **)getbytes(x->x_n_ambi * sizeof(t_symbol *)); + + j = x->x_n_ind_ls + x->x_n_mrg_mir_ls; + for(i=0; i<j; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrir->s_name); + x->x_s_hrir[i] = gensym(buf); + } + + for(i=0; i<n_ambi; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrir_red->s_name); + x->x_s_hrir_red[i] = gensym(buf); + } + + sprintf(buf, "%d_%s", prefix, s_fade_out_hrir->s_name); + x->x_s_fade_out_hrir = gensym(buf); + + 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_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + x->x_ls_encode = (double *)getbytes((x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + x->x_prod2 = (double *)getbytes((x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * x->x_n_ambi * sizeof(double)); + x->x_prod3 = (double *)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * x->x_n_ambi * sizeof(double)); + x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double)); + + x->x_delta = (int *)getbytes((x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * sizeof(int)); + x->x_phi = (int *)getbytes((x->x_n_ind_ls+2*x->x_n_mrg_mir_ls+x->x_n_ph_ls) * sizeof(int)); + x->x_phi_sym = (int *)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(int)); + x->x_sym_flag = (int *)getbytes((x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(int)); + + x->x_beg_fade_out_hrir = (float *)0; + x->x_beg_hrir = (float *)getbytes(x->x_firsize * (x->x_n_ind_ls+x->x_n_mrg_mir_ls) * sizeof(float)); + x->x_beg_hrir_red = (float **)getbytes(x->x_n_ambi * sizeof(float *)); + + 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; + x->x_seq_ok = 1; + for(i=0; i<n_ambi; i++) + x->x_ambi_channel_weight[i] = 1.0; + x->x_mirror_weight = 0.0; + + outlet_new(&x->x_obj, &s_list); + x->x_out_sign_sum = outlet_new(&x->x_obj, &s_list); + return(x); + } + else + { + post("bin_ambi_reduced_decode_fir-ERROR: need 1 float + 3 symbols + 6 floats arguments:"); + post(" prefix(unique-number) + hrir_loudspeaker_name + hrir_redused_name + hrir_fade_out_name +"); + post(" + ambi_order + ambi_dimension + number_of_independent_loudspeakers + "); + post(" + number_of_mirrored_and_merged_loudspeakers + number_of_phantom_loudspeakers + firsize"); + return(0); + } +} + +void bin_ambi_reduced_decode_fir_setup(void) +{ + bin_ambi_reduced_decode_fir_class = class_new(gensym("bin_ambi_reduced_decode_fir"), (t_newmethod)bin_ambi_reduced_decode_fir_new, (t_method)bin_ambi_reduced_decode_fir_free, + sizeof(t_bin_ambi_reduced_decode_fir), 0, A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_ind_ls, gensym("ind_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_mrg_ls, gensym("mrg_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_mir_ls, gensym("mir_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_pht_ls, gensym("pht_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_mirror_weight, gensym("mirror_weight"), A_DEFFLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_calc_pinv, gensym("calc_pinv"), 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_load_HRIR, gensym("load_HRIR"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_check_HRIR_LS_arrays, gensym("check_HRIR_LS_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_check_HRIR_RED_arrays, gensym("check_HRIR_RED_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_calc_reduced, gensym("calc_reduced"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_calc_sym, gensym("calc_sym"), 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_ambi_weight, gensym("ambi_weight"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir_class, (t_method)bin_ambi_reduced_decode_fir_sing_range, gensym("sing_range"), A_DEFFLOAT, 0); + class_sethelpsymbol(bin_ambi_reduced_decode_fir_class, gensym("iemhelp2/help-bin_ambi_reduced_decode_fir")); +} +/* +Reihenfolge: +n_ls x bin_ambi_reduced_decode_fir_ls + +bin_ambi_reduced_decode_fir_calc_pinv + +n_ls x bin_ambi_reduced_decode_fir_load_HRIR + +n_ls x bin_ambi_reduced_decode_fir_check_HRIR_arrays + +n_ambi x bin_ambi_reduced_decode_fir_check_HRTF_arrays + +n_ambi x bin_ambi_reduced_decode_fir_calc_reduced + +bin_ambi_reduced_decode_fir_calc_sym + +*/ diff --git a/src/bin_ambi_reduced_decode_fir2.c b/src/bin_ambi_reduced_decode_fir2.c new file mode 100644 index 0000000..fb0ce92 --- /dev/null +++ b/src/bin_ambi_reduced_decode_fir2.c @@ -0,0 +1,1260 @@ +/* For information on usage and redistribution, and for a DISCLAIMER OF ALL +* WARRANTIES, see the file, "LICENSE.txt," in this distribution. + +iem_bin_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_bin_ambi.h" +#include <math.h> +#include <stdio.h> +#include <string.h> + + +/* -------------------------- bin_ambi_reduced_decode_fir2 ------------------------------ */ +/* + ** berechnet ein reduziertes Ambisonic-Decoder-Set in die HRTF-Spektren ** + ** Inputs: ls + Liste von 3 floats: Index [1 .. 16] + Elevation [-90 .. +90 degree] + Azimut [0 .. 360 degree] ** + ** Inputs: calc_inv ** + ** Inputs: load_HRIR + float index1..16 ** + ** 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 1016_1_L_HRIR .. 1016_16_L_HRIR heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_re .. 1016_16_HRTF_re heissen und existieren ** + ** setzt voraus , dass die HRTF-tabele-names von 1016_1_HRTF_im .. 1016_16_HRTF_im heissen und existieren ** + */ + + +typedef struct _bin_ambi_reduced_decode_fir2 +{ + t_object x_obj; + t_atom x_at[2];/*output filename and tablename of HRIR*/ + int x_n_dim;/*dimension of ambisonic system*/ + int x_n_ambi;/*number of ambi channels*/ + int x_n_order;/*order of ambisonic system*/ + int x_n_real_ls;/*number of real loudspeakers*/ + int x_n_pht_ls;/*number of phantom loudspeakers*/ + int x_seq_ok; + int x_firsize;/*size of FIR-filter*/ + double *x_inv_work1;/* n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_work2;/* 2*n_ambi*n_ambi buffer for matrix inverse */ + double *x_inv_buf2;/* 2*n_ambi buffer for matrix inverse */ + double *x_transp;/* n_ls*n_ambi transposed input-buffer of decoder-matrix before inversion */ + double *x_ls_encode;/* n_ambi*n_ls straight input-buffer of decoder-matrix before inversion */ + double *x_prod2;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_prod3;/* n_ls*n_ambi output-buffer of decoder-matrix after inversion */ + double *x_ambi_channel_weight; + int *x_delta; + int *x_phi; + int *x_phi_sym; + int *x_sym_flag; + float *x_beg_fade_out_hrir; + float *x_beg_hrir; + float **x_beg_hrir_red; + t_symbol **x_hrir_filename; + t_symbol **x_s_hrir; + t_symbol **x_s_hrir_red; + t_symbol *x_s_fade_out_hrir; + void *x_out_sign_sum; + 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; + double x_sing_range; +} t_bin_ambi_reduced_decode_fir2; + +static t_class *bin_ambi_reduced_decode_fir2_class; + +static void bin_ambi_reduced_decode_fir2_convert(t_bin_ambi_reduced_decode_fir2 *x, double *delta_deg2rad, double *phi_deg2rad, int index) +{ + double q = 1.0; + double d = *delta_deg2rad; + double p = *phi_deg2rad; + int i; + + if(d < -90.0) + d = -90.0; + if(d > 90.0) + d = 90.0; + while(p < 0.0) + p += 360.0; + while(p >= 360.0) + p -= 360.0; + + x->x_delta[index] = (int)(d); + x->x_phi[index] = (int)(p); + + *delta_deg2rad = d*x->x_pi_over_180; + *phi_deg2rad = p*x->x_pi_over_180; +} + +static void bin_ambi_reduced_decode_fir2_do_2d(t_bin_ambi_reduced_decode_fir2 *x, int argc, t_atom *argv, int mode) +{ + double delta=0.0, phi; + double *dw = x->x_transp; + int index; + int order=x->x_n_order; + + if(argc < 2) + { + post("bin_ambi_reduced_decode_fir2 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 == BIN_AMBI_LS_REAL) + { + if(index >= x->x_n_real_ls) + index = x->x_n_real_ls - 1; + } + else if(mode == BIN_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; + + bin_ambi_reduced_decode_fir2_convert(x, &delta, &phi, index); + + 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 bin_ambi_reduced_decode_fir2_do_3d(t_bin_ambi_reduced_decode_fir2 *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("bin_ambi_reduced_decode_fir2 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 == BIN_AMBI_LS_REAL) + { + if(index >= x->x_n_real_ls) + index = x->x_n_real_ls - 1; + } + else if(mode == BIN_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; + + bin_ambi_reduced_decode_fir2_convert(x, &delta, &phi, index); + + cd = cos(delta); + sd = sin(delta); + cp = cos(phi); + sp = sin(phi); + + dw += index * x->x_n_ambi; + + *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 bin_ambi_reduced_decode_fir2_real_ls(t_bin_ambi_reduced_decode_fir2 *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fir2_do_2d(x, argc, argv, BIN_AMBI_LS_REAL); + else + bin_ambi_reduced_decode_fir2_do_3d(x, argc, argv, BIN_AMBI_LS_REAL); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode_fir2_pht_ls(t_bin_ambi_reduced_decode_fir2 *x, t_symbol *s, int argc, t_atom *argv) +{ + if(x->x_n_dim == 2) + bin_ambi_reduced_decode_fir2_do_2d(x, argc, argv, BIN_AMBI_LS_PHT); + else + bin_ambi_reduced_decode_fir2_do_3d(x, argc, argv, BIN_AMBI_LS_PHT); +} + +static void bin_ambi_reduced_decode_fir2_copy_row2buf(t_bin_ambi_reduced_decode_fir2 *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 bin_ambi_reduced_decode_fir2_copy_buf2row(t_bin_ambi_reduced_decode_fir2 *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 bin_ambi_reduced_decode_fir2_copy_row2row(t_bin_ambi_reduced_decode_fir2 *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 bin_ambi_reduced_decode_fir2_xch_rows(t_bin_ambi_reduced_decode_fir2 *x, int row1, int row2) +{ + bin_ambi_reduced_decode_fir2_copy_row2buf(x, row1); + bin_ambi_reduced_decode_fir2_copy_row2row(x, row2, row1); + bin_ambi_reduced_decode_fir2_copy_buf2row(x, row2); +} + +static void bin_ambi_reduced_decode_fir2_mul_row(t_bin_ambi_reduced_decode_fir2 *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 bin_ambi_reduced_decode_fir2_mul_col(t_bin_ambi_reduced_decode_fir2 *x, int col, double mul) +{ + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int i; + double *dw=x->x_inv_work2; + + dw += col; + for(i=0; i<n_ambi; i++) + { + (*dw) *= mul; + dw += n_ambi2; + } +} + +static void bin_ambi_reduced_decode_fir2_mul_buf_and_add2row(t_bin_ambi_reduced_decode_fir2 *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 bin_ambi_reduced_decode_fir2_eval_which_element_of_col_not_zero(t_bin_ambi_reduced_decode_fir2 *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 bin_ambi_reduced_decode_fir2_mul1(t_bin_ambi_reduced_decode_fir2 *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 bin_ambi_reduced_decode_fir2_mul2(t_bin_ambi_reduced_decode_fir2 *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_prod2; + 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; + } +} + +/* wird ersetzt durch haendische spiegelung und reduzierung +static void bin_ambi_reduced_decode_fir2_mul3(t_bin_ambi_reduced_decode_fir2 *x) +{ + int i, n; + double *dv3=x->x_prod3; + double *dv2=x->x_prod2; + double *dv1=x->x_prod2; + double mw=x->x_mirror_weight; + + n = x->x_n_ind_ls * x->x_n_ambi; + for(i=0; i<n; i++) + *dv3++ = *dv1++; + dv2 += n; + n = x->x_n_mrg_mir_ls * x->x_n_ambi; + dv2 += n; + for(i=0; i<n; i++) + *dv3++ = *dv1++ + (*dv2++)*mw; +}*/ + +static void bin_ambi_reduced_decode_fir2_transp_back(t_bin_ambi_reduced_decode_fir2 *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; + + transp = x->x_transp; + for(j=0; j<n_ambi; j++) + { + vec = transp; + for(i=0; i<n_ls; i++) + { + *straight++ = *vec; + vec += n_ambi; + } + transp++; + } +} + +static void bin_ambi_reduced_decode_fir2_inverse(t_bin_ambi_reduced_decode_fir2 *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 = bin_ambi_reduced_decode_fir2_eval_which_element_of_col_not_zero(x, i, i); + if(nz < 0) + { + post("bin_ambi_reduced_decode_fir2 ERROR: matrix not regular !!!!"); + x->x_seq_ok = 0; + return; + } + else + { + if(nz != i) + bin_ambi_reduced_decode_fir2_xch_rows(x, i, nz); + dv = db + i*n_ambi2 + i; + rcp = 1.0 /(*dv); + bin_ambi_reduced_decode_fir2_mul_row(x, i, rcp); + bin_ambi_reduced_decode_fir2_copy_row2buf(x, i); + for(j=i+1; j<n_ambi; j++) + { + dv += n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode_fir2_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; + bin_ambi_reduced_decode_fir2_copy_row2buf(x, i); + for(j=i-1; j>=0; j--) + { + dv -= n_ambi2; + rcp = -(*dv); + bin_ambi_reduced_decode_fir2_mul_buf_and_add2row(x, j, rcp); + } + } + + post("matrix_inverse regular"); + x->x_seq_ok = 1; +} + +static void bin_ambi_reduced_decode_fir2_calc_pinv(t_bin_ambi_reduced_decode_fir2 *x) +{ + t_garray *a; + int npoints; + t_float *fadevec; + int i, n; + double *dv3=x->x_prod3; + double *dv2=x->x_prod2; + + if((int)(x->x_beg_fade_out_hrir) == 0) + { + if (!(a = (t_garray *)pd_findbyclass(x->x_s_fade_out_hrir, garray_class))) + error("%s: no such array", x->x_s_fade_out_hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &fadevec)) + error("%s: bad template for bin_ambi_reduced_decode_fir2", x->x_s_fade_out_hrir->s_name); + else if (npoints < x->x_firsize) + error("%s: bad array-size: %d", x->x_s_fade_out_hrir->s_name, npoints); + else + x->x_beg_fade_out_hrir = fadevec; + } + + bin_ambi_reduced_decode_fir2_transp_back(x); + bin_ambi_reduced_decode_fir2_mul1(x); + bin_ambi_reduced_decode_fir2_inverse(x); + bin_ambi_reduced_decode_fir2_mul2(x); + + n = x->x_n_real_ls * x->x_n_ambi; + for(i=0; i<n; i++) + *dv3++ = *dv2++; +} + +/* +x_prod: +n_ambi columns; +n_ambi rows; +*/ + +static void bin_ambi_reduced_decode_fir2_ipht_ireal_muladd(t_bin_ambi_reduced_decode_fir2 *x, t_symbol *s, int argc, t_atom *argv) +{ + int n = x->x_n_ambi; + int i, pht_index, real_index; + double *dv3=x->x_prod3; + double *dv2=x->x_prod2; + float mw; + + if(argc < 3) + { + post("bin_ambi_reduced_decode_fir2 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; + + dv3 += real_index*x->x_n_ambi; + dv2 += (x->x_n_real_ls+pht_index)*x->x_n_ambi; + for(i=0; i<n; i++) + { + *dv3 += (*dv2++)*mw; + dv3++; + } +} + +static void bin_ambi_reduced_decode_fir2_load_HRIR(t_bin_ambi_reduced_decode_fir2 *x, t_symbol *s, int argc, t_atom *argv) +{ + int index; + t_symbol *hrirname; + + if(argc < 2) + { + post("bin_ambi_reduced_decode_fir2 ERROR: load_HRIR needs 1 index and 1 HRIR-wav"); + return; + } + index = (int)atom_getint(argv++) - 1; + hrirname = atom_getsymbol(argv); + if(index < 0) + index = 0; + if(index >= x->x_n_real_ls) + index = x->x_n_real_ls - 1; + x->x_hrir_filename[index] = hrirname; + + SETSYMBOL(x->x_at, x->x_hrir_filename[index]); + SETSYMBOL(x->x_at+1, x->x_s_hrir[index]); + outlet_list(x->x_obj.ob_outlet, &s_list, 2, x->x_at); +} + +static void bin_ambi_reduced_decode_fir2_check_HRIR_arrays(t_bin_ambi_reduced_decode_fir2 *x, float findex) +{ + int index=(int)findex - 1; + int j, k, n; + int firsize = x->x_firsize; + t_garray *a; + int npoints; + t_symbol *hrir; + t_float *vec_hrir, *vec, *vec_fade_out_hrir; + float decr, sum; + + if(index < 0) + index = 0; + if(index >= x->x_n_real_ls) + index = x->x_n_real_ls - 1; + + hrir = x->x_s_hrir[index]; + if (!(a = (t_garray *)pd_findbyclass(hrir, garray_class))) + error("%s: no such array", hrir->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrir)) + error("%s: bad template for bin_ambi_reduced_decode_fir2", hrir->s_name); + else + { + if(npoints < firsize) + { + post("bin_ambi_reduced_decode_fir2-WARNING: %s-array-size: %d < FIR-size: %d", hrir->s_name, npoints, firsize); + } + vec = x->x_beg_hrir; + vec += index * firsize; + + if((int)(x->x_beg_fade_out_hrir)) + { + vec_fade_out_hrir = x->x_beg_fade_out_hrir; + for(j=0; j<firsize; j++) + vec[j] = vec_hrir[j]*vec_fade_out_hrir[j]; + } + else + { + post("no HRIR-fade-out-window found"); + n = firsize * 3; + n /= 4; + for(j=0; j<n; j++) + vec[j] = vec_hrir[j]; + sum = 1.0f; + decr = 4.0f / (float)firsize; + for(j=n, k=0; j<firsize; j++, k++) + { + sum -= decr; + vec[j] = vec_hrir[j] * sum; + } + } + } +} + +static void bin_ambi_reduced_decode_fir2_check_HRIR_RED_arrays(t_bin_ambi_reduced_decode_fir2 *x, float findex) +{ + int index=(int)findex - 1; + t_garray *a; + int npoints; + int firsize = x->x_firsize; + t_float *vec_hrir_red; + t_symbol *hrir_red; + + if(index < 0) + index = 0; + if(index >= x->x_n_ambi) + index = x->x_n_ambi - 1; + + hrir_red = x->x_s_hrir_red[index]; + + if (!(a = (t_garray *)pd_findbyclass(hrir_red, garray_class))) + error("%s: no such array", hrir_red->s_name); + else if (!garray_getfloatarray(a, &npoints, &vec_hrir_red)) + error("%s: bad template for bin_ambi_reduced_decode_fir2", hrir_red->s_name); + else if (npoints < firsize) + error("%s: bad array-size: %d", hrir_red->s_name, npoints); + else + { + x->x_beg_hrir_red[index] = vec_hrir_red; + } +} + +static void bin_ambi_reduced_decode_fir2_calc_reduced(t_bin_ambi_reduced_decode_fir2 *x, float findex) +{ + int index=(int)findex - 1; + int i, j; + int firsize = x->x_firsize; + t_float *vec_hrir, *vec_hrir_red; + double *dv; + int n_ambi = x->x_n_ambi; + int n_ls = x->x_n_real_ls; + float mul; + + if(x->x_seq_ok) + { + if(index < 0) + index = 0; + if(index >= n_ambi) + index = n_ambi - 1; + + vec_hrir_red = x->x_beg_hrir_red[index]; + + dv = x->x_prod3 + index; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + for(i=0; i<firsize; i++)/*first step of acumulating the HRIRs*/ + { + vec_hrir_red[i] = mul * vec_hrir[i]; + } + + for(j=1; j<n_ls; j++) + { + dv += n_ambi; + mul = (float)(*dv); + vec_hrir = x->x_beg_hrir; + vec_hrir += j * firsize; + for(i=0; i<firsize; i++) + { + vec_hrir_red[i] += mul * vec_hrir[i]; + } + } + } +} + +static void bin_ambi_reduced_decode_fir2_calc_sym(t_bin_ambi_reduced_decode_fir2 *x) +{ + int *sym=x->x_phi_sym; + int n_ls = x->x_n_real_ls; + int n_ambi = x->x_n_ambi; + int n_ambi2 = 2*n_ambi; + int *phi=x->x_phi; + int *delta=x->x_delta; + int *flag=x->x_sym_flag; + int i, j, d, p, regular=1; + double *dv, *db=x->x_prod3; + double a, b, q; + int sym_max=0; + int pos_sym_counter, neg_sym_counter; + char plus_minus[100]; + char abc[100]; + char onenine[100]; + char ten[100]; + + if(x->x_seq_ok) + { + for(i=0; i<n_ls; i++) + sym[i] = -1; + + plus_minus[0] = 0; + if(x->x_n_dim == 2) + { + strcpy(abc, " WXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXYXY"); + strcpy(ten, " 000000000000000000011111111111111111111"); + strcpy(onenine, " 011223344556677889900112233445566778899"); + abc[n_ambi+22] = 0; + ten[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(ten); + post(onenine); + } + else + { + strcpy(abc, " AABCABCDEABCDEFGABCDEFGHIABCDEFGHIJKABCDEFGHIJKLM"); + strcpy(onenine, " 0111222223333333444444444555555555556666666666666"); + abc[n_ambi+22] = 0; + onenine[n_ambi+22] = 0; + post(abc); + post(onenine); + } + + for(i=0; i<n_ls; i++)/* looking for azimuth symmetries of loudspeakers */ + { + d = delta[i]; + p = phi[i]; + for(j=i+1; j<n_ls; j++) + { + if(d == delta[j]) + { + if((p + phi[j]) == 360) + { + if((sym[i] < 0) && (sym[j] < 0)) + { + sym[i] = j; + sym[j] = i; + j = n_ls + 1; + sym_max++; + } + } + } + } + } + + for(p=0; p<n_ambi; p++)/*each col*/ + { + pos_sym_counter = 0; + neg_sym_counter = 0; + for(i=0; i<n_ls; i++) + flag[i] = 1; + dv = db + p; + for(i=0; i<n_ls; i++) + { + if((sym[i] >= 0) && flag[i]) + { + j = sym[i]; + flag[i] = 0; + flag[j] = 0; + a = dv[i*n_ambi]; + b = dv[j*n_ambi]; + if((a < 5.0e-4)&&(a > -5.0e-4)&&(b < 5.0e-4)&&(b > -5.0e-4)) + { + pos_sym_counter++; + neg_sym_counter++; + } + else + { + q = a / b; + if((q < 1.005)&&(q > 0.995)) + pos_sym_counter++; + else if((q > -1.005)&&(q < -0.995)) + neg_sym_counter++; + } + } + } + if(pos_sym_counter == sym_max) + strcat(plus_minus, "+"); + else if(neg_sym_counter == sym_max) + strcat(plus_minus, "-"); + else + { + strcat(plus_minus, "?"); + regular = 0; + } + } + post("sum of right channel: %s", plus_minus); + + if(regular) + { + for(i=0; i<n_ambi; i++) + { + /*change*/ + if(plus_minus[i] == '+') + SETFLOAT(x->x_at, 1.0f); + else if(plus_minus[i] == '-') + SETFLOAT(x->x_at, 2.0f); + SETFLOAT(x->x_at+1, (float)(i+1)); + outlet_list(x->x_out_sign_sum, &s_list, 2, x->x_at); + } + } + } +} + +static void bin_ambi_reduced_decode_fir2_ambi_weight(t_bin_ambi_reduced_decode_fir2 *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] = (double)atom_getfloat(argv++); + k++; + if(x->x_n_dim == 2) + { + for(i=1; i<=n; i++) + { + d = (double)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 = (double)atom_getfloat(argv++); + m = 2*i + 1; + for(j=0; j<m; j++) + { + x->x_ambi_channel_weight[k] = d; + k++; + } + } + } + } + else + post("bin_ambi_reduced_decode_fir2-ERROR: ambi_weight needs %d float weights", x->x_n_order+1); +} + +static void bin_ambi_reduced_decode_fir2_sing_range(t_bin_ambi_reduced_decode_fir2 *x, t_floatarg f) +{ + if(f < 0.0f) + x->x_sing_range = -(double)f; + else + x->x_sing_range = (double)f; +} + +static void bin_ambi_reduced_decode_fir2_free(t_bin_ambi_reduced_decode_fir2 *x) +{ + freebytes(x->x_hrir_filename, x->x_n_real_ls * sizeof(t_symbol *)); + freebytes(x->x_s_hrir, x->x_n_real_ls * sizeof(t_symbol *)); + freebytes(x->x_s_hrir_red, x->x_n_ambi * sizeof(t_symbol *)); + 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_prod2, (x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double)); + freebytes(x->x_prod3, x->x_n_real_ls * x->x_n_ambi * sizeof(double)); + freebytes(x->x_ambi_channel_weight, x->x_n_ambi * sizeof(double)); + + freebytes(x->x_delta, (x->x_n_real_ls+x->x_n_pht_ls)* sizeof(int)); + freebytes(x->x_phi, (x->x_n_real_ls+x->x_n_pht_ls) * sizeof(int)); + freebytes(x->x_phi_sym, x->x_n_real_ls * sizeof(int)); + freebytes(x->x_sym_flag, x->x_n_real_ls * sizeof(int)); + + freebytes(x->x_beg_hrir, x->x_firsize * x->x_n_real_ls * sizeof(float)); + freebytes(x->x_beg_hrir_red, x->x_n_ambi * sizeof(float *)); +} + +/* + 1.arg_ int prefix; + 2.arg: t_symbol *hrir_name; + 3.arg: t_symbol *hrir_red_name; + 4.arg: t_symbol *hrtf_im_name; + 5.arg: t_symbol *hrir_fade_out_name; + 6.arg: int ambi_order; + 7.arg: int dim; + 8.arg: int number_of_loudspeakers; + 9.arg: int firsize; + */ + +static void *bin_ambi_reduced_decode_fir2_new(t_symbol *s, int argc, t_atom *argv) +{ + t_bin_ambi_reduced_decode_fir2 *x = (t_bin_ambi_reduced_decode_fir2 *)pd_new(bin_ambi_reduced_decode_fir2_class); + char buf[400]; + int i, j, ok=0; + int n_order=0, n_dim=0, n_real_ls=0, n_pht_ls=0, n_ambi=0, firsize=0, prefix=0; + t_symbol *s_hrir=gensym("L_HRIR"); + t_symbol *s_hrir_red=gensym("HRIR_red"); + t_symbol *s_fade_out_hrir=gensym("HRIR_win"); + + if((argc >= 9) && + IS_A_FLOAT(argv,0) && + IS_A_SYMBOL(argv,1) && + IS_A_SYMBOL(argv,2) && + IS_A_SYMBOL(argv,3) && + IS_A_FLOAT(argv,4) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = (t_symbol *)atom_getsymbolarg(1, argc, argv); + s_hrir_red = (t_symbol *)atom_getsymbolarg(2, argc, argv); + s_fade_out_hrir = (t_symbol *)atom_getsymbolarg(3, argc, argv); + + n_order = (int)atom_getintarg(4, argc, argv); + n_dim = (int)atom_getintarg(5, argc, argv); + n_real_ls = (int)atom_getintarg(6, argc, argv); + n_pht_ls = (int)atom_getintarg(7, argc, argv); + firsize = (int)atom_getintarg(8, argc, argv); + + ok = 1; + } + else if((argc >= 9) && + 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) && + IS_A_FLOAT(argv,5) && + IS_A_FLOAT(argv,6) && + IS_A_FLOAT(argv,7) && + IS_A_FLOAT(argv,8)) + { + prefix = (int)atom_getintarg(0, argc, argv); + + s_hrir = gensym("L_HRIR"); + s_hrir_red = gensym("HRIR_red"); + s_fade_out_hrir = gensym("HRIR_win"); + + n_order = (int)atom_getintarg(4, argc, argv); + n_dim = (int)atom_getintarg(5, argc, argv); + n_real_ls = (int)atom_getintarg(6, argc, argv); + n_pht_ls = (int)atom_getintarg(7, argc, argv); + firsize = (int)atom_getintarg(8, argc, argv); + + ok = 1; + } + + if(ok) + { + if(n_order < 1) + n_order = 1; + + if(n_dim == 3) + { + if(n_order > 5) + n_order = 5; + n_ambi = (n_order + 1)*(n_order + 1); + } + else + { + if(n_order > 12) + n_order = 12; + n_dim = 2; + n_ambi = 2 * n_order + 1; + } + + 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) < n_ambi) + { + post("bin_ambi_reduced_decode_fir2-WARNING: Number of all Loudspeakers < Number of Ambisonic-Channels !!!!"); + } + + if(firsize < 32) + firsize = 32; + + x->x_n_dim = n_dim; + x->x_n_ambi = n_ambi; + x->x_n_real_ls = n_real_ls; + x->x_n_pht_ls = n_pht_ls; + x->x_n_order = n_order; + x->x_firsize = firsize; + + x->x_hrir_filename = (t_symbol **)getbytes(x->x_n_real_ls * sizeof(t_symbol *)); + x->x_s_hrir = (t_symbol **)getbytes(x->x_n_real_ls * sizeof(t_symbol *)); + x->x_s_hrir_red = (t_symbol **)getbytes(x->x_n_ambi * sizeof(t_symbol *)); + + j = x->x_n_real_ls; + for(i=0; i<j; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrir->s_name); + x->x_s_hrir[i] = gensym(buf); + } + + for(i=0; i<n_ambi; i++) + { + sprintf(buf, "%d_%d_%s", prefix, i+1, s_hrir_red->s_name); + x->x_s_hrir_red[i] = gensym(buf); + } + + sprintf(buf, "%d_%s", prefix, s_fade_out_hrir->s_name); + x->x_s_fade_out_hrir = gensym(buf); + + 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_prod2 = (double *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * x->x_n_ambi * sizeof(double)); + x->x_prod3 = (double *)getbytes(x->x_n_real_ls * x->x_n_ambi * sizeof(double)); + x->x_ambi_channel_weight = (double *)getbytes(x->x_n_ambi * sizeof(double)); + + x->x_delta = (int *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * sizeof(int)); + x->x_phi = (int *)getbytes((x->x_n_real_ls+x->x_n_pht_ls) * sizeof(int)); + x->x_phi_sym = (int *)getbytes(x->x_n_real_ls * sizeof(int)); + x->x_sym_flag = (int *)getbytes(x->x_n_real_ls * sizeof(int)); + + x->x_beg_fade_out_hrir = (float *)0; + x->x_beg_hrir = (float *)getbytes(x->x_firsize * x->x_n_real_ls * sizeof(float)); + x->x_beg_hrir_red = (float **)getbytes(x->x_n_ambi * sizeof(float *)); + + 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; + x->x_seq_ok = 1; + for(i=0; i<n_ambi; i++) + x->x_ambi_channel_weight[i] = 1.0; + + outlet_new(&x->x_obj, &s_list); + x->x_out_sign_sum = outlet_new(&x->x_obj, &s_list); + return(x); + } + else + { + post("bin_ambi_reduced_decode_fir2-ERROR: need 1 float + 3 symbols + 5 floats arguments:"); + post(" prefix(unique-number) + hrir_loudspeaker_name + hrir_redused_name + hrir_fade_out_name +"); + post(" + ambi_order + ambi_dimension + number_of_real_loudspeakers + "); + post(" + number_of_phantom_loudspeakers + firsize"); + return(0); + } +} + +void bin_ambi_reduced_decode_fir2_setup(void) +{ + bin_ambi_reduced_decode_fir2_class = class_new(gensym("bin_ambi_reduced_decode_fir2"), (t_newmethod)bin_ambi_reduced_decode_fir2_new, (t_method)bin_ambi_reduced_decode_fir2_free, + sizeof(t_bin_ambi_reduced_decode_fir2), 0, A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_real_ls, gensym("real_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_pht_ls, gensym("pht_ls"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_calc_pinv, gensym("calc_pinv"), 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_ipht_ireal_muladd, gensym("ipht_ireal_muladd"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_load_HRIR, gensym("load_HRIR"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_check_HRIR_arrays, gensym("check_HRIR_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_check_HRIR_RED_arrays, gensym("check_HRIR_RED_arrays"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_calc_reduced, gensym("calc_reduced"), A_FLOAT, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_calc_sym, gensym("calc_sym"), 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_ambi_weight, gensym("ambi_weight"), A_GIMME, 0); + class_addmethod(bin_ambi_reduced_decode_fir2_class, (t_method)bin_ambi_reduced_decode_fir2_sing_range, gensym("sing_range"), A_DEFFLOAT, 0); + class_sethelpsymbol(bin_ambi_reduced_decode_fir2_class, gensym("iemhelp2/help-bin_ambi_reduced_decode_fir2")); +} +/* +Reihenfolge: +n_ls x bin_ambi_reduced_decode_fir2_ls + +bin_ambi_reduced_decode_fir2_calc_pinv + +n_ls x bin_ambi_reduced_decode_fir2_load_HRIR + +n_ls x bin_ambi_reduced_decode_fir2_check_HRIR_arrays + +n_ambi x bin_ambi_reduced_decode_fir2_check_HRTF_arrays + +n_ambi x bin_ambi_reduced_decode_fir2_calc_reduced + +bin_ambi_reduced_decode_fir2_calc_sym + +*/ diff --git a/src/iem_bin_ambi.c b/src/iem_bin_ambi.c index a580f76..3af6646 100644 --- a/src/iem_bin_ambi.c +++ b/src/iem_bin_ambi.c @@ -1,7 +1,7 @@ /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution. -iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2003 */ +iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */ #ifdef NT #pragma warning( disable : 4244 ) @@ -21,12 +21,26 @@ static void *iem_bin_ambi_new(void) return (x); } +void bin_ambi_calc_HRTF_setup(void); +void bin_ambi_reduced_decode_setup(void); +void bin_ambi_reduced_decode2_setup(void); +void bin_ambi_reduced_decode_fft_setup(void); +void bin_ambi_reduced_decode_fir_setup(void); +void bin_ambi_reduced_decode_fft2_setup(void); +void bin_ambi_reduced_decode_fir2_setup(void); + /* ------------------------ setup routine ------------------------- */ -void iem_bin_ambi_sources_setup(void); void iem_bin_ambi_setup(void) { - iem_bin_ambi_sources_setup(); - - post("iem_bin_ambi (R-1.15) library loaded!"); + bin_ambi_calc_HRTF_setup(); + bin_ambi_reduced_decode_setup(); + bin_ambi_reduced_decode2_setup(); + bin_ambi_reduced_decode_fft_setup(); + bin_ambi_reduced_decode_fir_setup(); + bin_ambi_reduced_decode_fft2_setup(); + bin_ambi_reduced_decode_fir2_setup(); + + post("iem_bin_ambi (R-1.16) library loaded! (c) Thomas Musil 05.2005"); + post(" musil%ciem.at iem KUG Graz Austria", '@'); } diff --git a/src/iem_bin_ambi.h b/src/iem_bin_ambi.h index 1757d2e..a9291c7 100644 --- a/src/iem_bin_ambi.h +++ b/src/iem_bin_ambi.h @@ -1,11 +1,17 @@ /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution. -iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2003 */ +iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */ #ifndef __IEMBINAMBI_H__ #define __IEMBINAMBI_H__ +#define BIN_AMBI_LS_REAL 0 +#define BIN_AMBI_LS_IND 0 +#define BIN_AMBI_LS_MRG 1 +#define BIN_AMBI_LS_MIR 2 +#define BIN_AMBI_LS_PHT 3 + typedef struct { float real; diff --git a/src/iemlib.h b/src/iemlib.h index 48a02bd..c71b0ed 100644 --- a/src/iemlib.h +++ b/src/iemlib.h @@ -1,7 +1,7 @@ /* For information on usage and redistribution, and for a DISCLAIMER OF ALL * WARRANTIES, see the file, "LICENSE.txt," in this distribution. -iem_bin_ambi written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2003 */ +iemlib.h written by Thomas Musil, Copyright (c) IEM KUG Graz Austria 2000 - 2005 */ #ifndef __IEMLIB_H__ #define __IEMLIB_H__ @@ -33,28 +33,40 @@ extern int sys_noloadbang; #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! */ +/* 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 */ +/* 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 NT - /* little-endian; most significant byte is at highest address */ +#ifdef MSW +/* little-endian; most significant byte is at highest address */ #define HIOFFSET 1 #define LOWOFFSET 0 #define int32 long #else -#ifdef __linux__ +#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 @@ -66,14 +78,25 @@ extern int sys_noloadbang; #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 /* NT */ +#endif /* MSW */ #endif /* SGI */ union tabfudge { - double tf_d; - int32 tf_i[2]; + 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 diff --git a/src/makefile_win b/src/makefile_win index a7ba3db..80f6286 100644 --- a/src/makefile_win +++ b/src/makefile_win @@ -1,39 +1,43 @@ -
-all: iem_bin_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\pthreadVC.lib \
- $(PD_INST_PATH)\bin\pd.lib
-
-
-SRC = bin_ambi_calc_HRTF.c \
- bin_ambi_reduced_decode.c \
- iem_bin_ambi.c
-
-
-OBJ = $(SRC:.c=.obj)
-
-.c.obj:
- cl $(PD_WIN_C_FLAGS) $(PD_WIN_INCLUDE_PATH) /c $*.c
-
-iem_bin_ambi.dll: $(OBJ)
- link $(PD_WIN_L_FLAGS) /dll /export:iem_bin_ambi_setup \
- /out:iem_bin_ambi.dll $(OBJ) $(PD_WIN_LIB)
-
-clean:
- del *.obj
+ +all: iem_bin_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 = bin_ambi_calc_HRTF.c \ + bin_ambi_reduced_decode.c \ + bin_ambi_reduced_decode2.c \ + bin_ambi_reduced_decode_fft.c \ + bin_ambi_reduced_decode_fir.c \ + bin_ambi_reduced_decode_fft2.c \ + bin_ambi_reduced_decode_fir2.c \ + iem_bin_ambi.c + + +OBJ = $(SRC:.c=.obj) + +.c.obj: + cl $(PD_WIN_C_FLAGS) $(PD_WIN_INCLUDE_PATH) /c $*.c + +iem_bin_ambi.dll: $(OBJ) + link $(PD_WIN_L_FLAGS) /dll /export:iem_bin_ambi_setup \ + /out:iem_bin_ambi.dll $(OBJ) $(PD_WIN_LIB) + +clean: + del *.obj |