From 243b48abf7922957fc3e0876bf9f12f893a8327f Mon Sep 17 00:00:00 2001 From: Franz Zotter Date: Wed, 17 Jun 2009 08:39:06 +0000 Subject: modified the chebyshev12 routine to put out readily normalized circular harmonics and added [mtx_circular_harmonics] to the file mtx_spherical_harmonics, as well as a corresponding helpfile. why: useful for circular (1D) ambisonics, e.g. ... svn path=/trunk/externals/iem/iemmatrix/; revision=11792 --- doc/mtx_circular_harmonics-help.pd | 186 +++++++++++++++++++++ src/mtx_spherical_harmonics.c | 121 +++++++++++++- src/mtx_spherical_harmonics/chebyshev12.c | 10 +- src/mtx_spherical_harmonics/sharmonics.c | 2 +- .../sharmonics_normalization.c | 9 + 5 files changed, 321 insertions(+), 7 deletions(-) create mode 100644 doc/mtx_circular_harmonics-help.pd diff --git a/doc/mtx_circular_harmonics-help.pd b/doc/mtx_circular_harmonics-help.pd new file mode 100644 index 0000000..2be06b0 --- /dev/null +++ b/doc/mtx_circular_harmonics-help.pd @@ -0,0 +1,186 @@ +#N canvas 61 39 921 503 10; +#X text 516 377 Franz Zotter \, 2009; +#X text 129 226 for -n<=m<=n:; +#X text 620 241 for m>=0; +#X text 619 257 for m< 0; +#X text 288 118 [mtx_circular_harmonics] requires a numerical creation +argument specifyiing the maximum order 0<=n<=nmax.; +#X text 74 54 [mtx_circular_harmonics] circular harmonics evaluated +at a set of points given in phi coordinates.; +#X text 284 160 for an L points 2xL input matrix \, [mtx_spherical_harmonics] +evaluates the (2*nmax+2) circular harmonics at L points and delivers +an Lx(2*nmax+2) output matrix.; +#X text 167 242 PHI_m(phi) = sqrt((2-delta_m) / (2*pi)) * cos(m*phi) +; +#X text 167 258 PHI_m(phi) = 1/sqrt(pi) * sin(m*phi); +#X text 126 291 The order of the harmonics in the output columns is +specified by the linear index k=nmax+m+1.; +#X text 125 328 [mtx_circular_harmonics] uses fully normalized PHI_m. +; +#N canvas 0 0 450 300 (subpatch) 0; +#X array circularharmonic1 100 float 1; +#A 0 0 0.106888 0.209904 0.305318 0.389673 0.459913 0.513495 0.548478 +0.563594 0.558297 0.532777 0.48796 0.425469 0.347567 0.257076 0.157273 +0.0517736 -0.055601 -0.160962 -0.260492 -0.350588 -0.427984 -0.489879 +-0.53403 -0.558838 -0.563404 -0.547564 -0.51189 -0.457675 -0.386883 +-0.302078 -0.206331 -0.10311 0.003845 0.110661 0.213468 0.308544 0.392444 +0.462129 0.515076 0.549366 0.563758 0.55773 0.5315 0.486019 0.422934 +0.34453 0.253647 0.153577 0.0479438 -0.059426 -0.164643 -0.263897 -0.353592 +-0.43048 -0.491775 -0.535258 -0.559353 -0.563189 -0.546625 -0.510261 +-0.455416 -0.384076 -0.298823 -0.202747 -0.0993281 0.00768981 0.114428 +0.217022 0.311755 0.395197 0.464324 0.516633 0.550229 0.563895 0.557136 +0.530198 0.484055 0.42038 0.341478 0.250207 0.149874 0.044112 -0.0632482 +-0.168317 -0.267289 -0.35658 -0.432955 -0.493648 -0.536461 -0.559843 +-0.562947 -0.54566 -0.508609 -0.453136 -0.38125 -0.295555 -0.199155 +-0.0955408 0.0115335; +#X coords 0 1 99 -1 200 50 1; +#X restore 720 14 graph; +#N canvas 0 0 450 300 (subpatch) 0; +#X array circularharmonic2 100 float 1; +#A 0 0 0.0714992 0.141845 0.209904 0.274579 0.334825 0.389673 0.438236 +0.479734 0.513495 0.538976 0.555766 0.563594 0.562334 0.552006 0.532777 +0.504957 0.468995 0.425469 0.375083 0.318648 0.257076 0.191358 0.122554 +0.0517737 -0.0198414 -0.0911363 -0.160962 -0.228192 -0.291742 -0.350588 +-0.40378 -0.450461 -0.489879 -0.521397 -0.544508 -0.558838 -0.564157 +-0.560378 -0.547564 -0.52592 -0.495795 -0.457675 -0.412175 -0.360029 +-0.302078 -0.239255 -0.172574 -0.10311 -0.031984 0.0396582 0.110661 +0.179879 0.246196 0.308544 0.365916 0.417388 0.462129 0.499418 0.528655 +0.549366 0.561219 0.564022 0.55773 0.542444 0.518411 0.486019 0.44579 +0.398372 0.344531 0.285133 0.221138 0.153577 0.0835392 0.0121547 -0.0594257 +-0.130048 -0.198574 -0.263897 -0.324964 -0.380792 -0.43048 -0.473225 +-0.50834 -0.535258 -0.553544 -0.562905 -0.563189 -0.554391 -0.536653 +-0.510261 -0.475642 -0.433353 -0.384075 -0.328605 -0.267836 -0.202748 +-0.13439 -0.0638651 0.00768928; +#X coords 0 1 99 -1 200 50 1; +#X restore 718 80 graph; +#N canvas 0 0 450 300 (subpatch) 0; +#X array circularharmonic3 100 float 1; +#A 0 0 0.0358219 0.0714992 0.106888 0.141845 0.176231 0.209904 0.242731 +0.274579 0.305318 0.334825 0.362981 0.389673 0.414791 0.438236 0.459913 +0.479734 0.497618 0.513495 0.527299 0.538976 0.548478 0.555766 0.560812 +0.563594 0.564102 0.562334 0.558297 0.552006 0.543489 0.532777 0.519916 +0.504957 0.48796 0.468995 0.448136 0.425469 0.401085 0.375083 0.347567 +0.318648 0.288444 0.257076 0.22467 0.191358 0.157273 0.122554 0.0873399 +0.0517737 0.0159984 -0.0198414 -0.055601 -0.0911363 -0.126304 -0.160962 +-0.19497 -0.228192 -0.260492 -0.291742 -0.321814 -0.350588 -0.377946 +-0.40378 -0.427984 -0.450461 -0.471121 -0.489879 -0.50666 -0.521397 +-0.53403 -0.544508 -0.552788 -0.558838 -0.562633 -0.564157 -0.563404 +-0.560378 -0.555091 -0.547564 -0.537827 -0.52592 -0.51189 -0.495795 +-0.477699 -0.457675 -0.435805 -0.412175 -0.386883 -0.360029 -0.331723 +-0.302078 -0.271213 -0.239255 -0.206331 -0.172574 -0.138121 -0.10311 +-0.0676836 -0.031984 0.00384473 ; +#X coords 0 1 99 -1 200 50 1; +#X restore 719 146 graph; +#N canvas 0 0 450 300 (subpatch) 0; +#X array circularharmonic4 100 float 1; +#A 0 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 0.398942 +0.398942 0.398942 0.398942 0.398942 0.398942; +#X coords 0 1 99 -1 200 50 1; +#X restore 718 217 graph; +#N canvas 0 0 450 300 (subpatch) 0; +#X array circularharmonic5 100 float 1; +#A 0 0.56419 0.563051 0.559641 0.553972 0.546068 0.53596 0.523689 0.509305 +0.492866 0.474437 0.454095 0.43192 0.408001 0.382437 0.355329 0.326787 +0.296927 0.265868 0.233737 0.200662 0.166777 0.13222 0.097129 0.061646 +0.0259142 -0.00992222 -0.0457185 -0.0813303 -0.116614 -0.151427 -0.185629 +-0.219082 -0.251651 -0.283204 -0.313614 -0.342759 -0.370521 -0.396788 +-0.421453 -0.444418 -0.465589 -0.484881 -0.502217 -0.517526 -0.530747 +-0.541826 -0.550718 -0.557388 -0.561809 -0.563963 -0.563841 -0.561443 +-0.55678 -0.54987 -0.540741 -0.52943 -0.515983 -0.500453 -0.482904 +-0.463407 -0.442039 -0.418887 -0.394045 -0.367613 -0.339697 -0.310411 +-0.279872 -0.248204 -0.215534 -0.181994 -0.14772 -0.112849 -0.0775239 +-0.0418852 -0.00607772 0.0297543 0.0654664 0.100914 0.135955 0.170447 +0.204251 0.237231 0.269253 0.300189 0.329914 0.358307 0.385255 0.410648 +0.434383 0.456366 0.476507 0.494725 0.510947 0.525107 0.537148 0.547022 +0.554687 0.560115 0.563282 0.564176; +#X coords 0 1 99 -1 200 50 1; +#X restore 718 286 graph; +#X obj 75 109 loadbang; +#N canvas 624 434 600 460 send_to_tables 0; +#X obj 21 14 inlet; +#X obj 72 211 mtx; +#X obj 21 36 t a a; +#X obj 21 59 mtx_size; +#X obj 74 105 until; +#X obj 74 82 t f b; +#X msg 117 82 0; +#X obj 101 124 + 1; +#X obj 74 124 f; +#X obj 73 147 t f f; +#X obj 72 249 s; +#X msg 105 168 symbol circularharmonic\$1; +#X obj 72 230 list prepend 0; +#X msg 73 192 column \$1; +#X connect 0 0 2 0; +#X connect 1 0 12 0; +#X connect 2 0 3 0; +#X connect 2 1 1 1; +#X connect 3 1 5 0; +#X connect 4 0 8 0; +#X connect 5 0 4 0; +#X connect 5 1 6 0; +#X connect 6 0 8 1; +#X connect 7 0 8 1; +#X connect 7 0 9 0; +#X connect 8 0 7 0; +#X connect 9 0 13 0; +#X connect 9 1 11 0; +#X connect 11 0 10 1; +#X connect 12 0 10 0; +#X connect 13 0 1 0; +#X restore 75 179 pd send_to_tables; +#X obj 75 133 mtx_linspace 0 6.29 100; +#X obj 159 114 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144 +-1 -1; +#X obj 75 157 mtx_circular_harmonics 3; +#N canvas 0 0 450 300 (subpatch) 0; +#X array circularharmonic6 100 float 1; +#A 0 0.56419 0.559641 0.546068 0.523689 0.492866 0.454095 0.408001 +0.355329 0.296927 0.233737 0.166777 0.097129 0.0259142 -0.0457185 -0.116614 +-0.185629 -0.251651 -0.313614 -0.370521 -0.421453 -0.465589 -0.502217 +-0.530747 -0.550718 -0.561809 -0.563841 -0.55678 -0.540741 -0.515983 +-0.482904 -0.442039 -0.394045 -0.339697 -0.279872 -0.215534 -0.14772 +-0.0775239 -0.00607772 0.0654664 0.135955 0.204251 0.269253 0.329914 +0.385255 0.434383 0.476507 0.510947 0.537148 0.554687 0.563282 0.562794 +0.553231 0.534746 0.507639 0.472346 0.429436 0.379601 0.323646 0.262471 +0.197064 0.128479 0.0578227 -0.0137663 -0.0851333 -0.155127 -0.22262 +-0.286523 -0.345805 -0.399512 -0.446776 -0.486836 -0.519045 -0.542885 +-0.557971 -0.564059 -0.561051 -0.548997 -0.528089 -0.498667 -0.461203 +-0.416302 -0.364688 -0.307193 -0.244745 -0.17835 -0.10908 -0.0380498 +0.0335932 0.104695 0.174108 0.240714 0.303438 0.361269 0.413275 0.458616 +0.496562 0.526501 0.54795 0.560563 0.564137; +#X coords 0 1 99 -1 200 50 1; +#X restore 718 356 graph; +#N canvas 0 0 450 300 (subpatch) 0; +#X array circularharmonic7 100 float 1; +#A 0 0.56419 0.553972 0.523689 0.474437 0.408001 0.326787 0.233737 +0.13222 0.0259142 -0.0813303 -0.185629 -0.283204 -0.370521 -0.444418 +-0.502217 -0.541826 -0.561809 -0.561443 -0.540741 -0.500453 -0.442039 +-0.367613 -0.279872 -0.181994 -0.0775239 0.0297544 0.135955 0.237231 +0.329914 0.410647 0.476507 0.525107 0.554687 0.564176 0.553231 0.522246 +0.472346 0.405336 0.323645 0.230232 0.128479 0.0220729 -0.0851331 -0.189256 +-0.286523 -0.373412 -0.446776 -0.503957 -0.542885 -0.562149 -0.561051 +-0.539632 -0.498667 -0.439639 -0.364688 -0.276527 -0.17835 -0.0737135 +0.0335929 0.139683 0.240714 0.333025 0.413275 0.478555 0.526501 0.555377 +0.564137 0.552464 0.520779 0.470232 0.402653 0.320489 0.226717 0.124732 +0.0182304 -0.0889317 -0.192873 -0.289828 -0.376285 -0.449113 -0.505674 +-0.543919 -0.562462 -0.560633 -0.538497 -0.496857 -0.437219 -0.361746 +-0.273169 -0.174699 -0.0698998 0.0374303 0.143405 0.244186 0.336121 +0.415882 0.48058 0.527871 0.556041 0.564072; +#X coords 0 1 99 -1 200 50 1; +#X restore 716 422 graph; +#X connect 16 0 18 0; +#X connect 18 0 20 0; +#X connect 19 0 18 0; +#X connect 20 0 17 0; diff --git a/src/mtx_spherical_harmonics.c b/src/mtx_spherical_harmonics.c index b0c270a..8ecc9a4 100644 --- a/src/mtx_spherical_harmonics.c +++ b/src/mtx_spherical_harmonics.c @@ -35,7 +35,22 @@ struct _MTXSh_ size_t l; }; -static void allocMTXdata (MTXSh *x) +static t_class *mtx_circular_harmonics_class; + +typedef struct _MTXCh_ MTXCh; +struct _MTXCh_ +{ + t_object x_obj; + t_outlet *list_ch_out; + t_atom *list_ch; + double *phi; + Cheby12WorkSpace *wc; + size_t nmax; + size_t l; +}; + + +static void allocMTXShdata (MTXSh *x) { x->phi=(double*)calloc(x->l,sizeof(double)); x->theta=(double*)calloc(x->l,sizeof(double)); @@ -104,7 +119,7 @@ static void mTXShMatrix (MTXSh *x, t_symbol *s, if (x->l!=columns) { deleteMTXShdata(x); x->l=columns; - allocMTXdata(x); + allocMTXShdata(x); } for (n=0;nl;n++) { x->phi[n]=(double) atom_getfloat(argv+n); @@ -125,6 +140,92 @@ static void mTXShMatrix (MTXSh *x, t_symbol *s, } +} + +static void allocMTXChdata (MTXCh *x) +{ + x->phi=(double*)calloc(x->l,sizeof(double)); + x->wc=chebyshev12_alloc(x->nmax,x->l); + x->list_ch=(t_atom*)calloc(x->l*(2*x->nmax+1)+2,sizeof(t_atom)); +} + +static void deleteMTXChdata (MTXCh *x) +{ + if (x->phi!=0) + free(x->phi); + if (x->list_ch!=0) + free(x->list_ch); + chebyshev12_free(x->wc); + x->wc=0; + x->list_ch=0; + x->phi=0; +} + +static void *newMTXCh (t_symbol *s, int argc, t_atom *argv) +{ + int nmax; + MTXCh *x = (MTXCh *) pd_new (mtx_circular_harmonics_class); + x->list_ch_out = outlet_new (&x->x_obj, gensym("matrix")); + x->list_ch = 0; + x->phi = 0; + x->wc = 0; + x->l=0; + nmax=(int) atom_getfloat(argv); + if (nmax<0) + nmax=0; + x->nmax=nmax; + + return ((void *) x); +} + +static void mTXChBang (MTXCh *x) +{ + if (x->list_ch!=0) { + outlet_anything(x->list_ch_out, gensym("matrix"), x->l*(2*x->nmax+1)+2, x->list_ch); + } +} + +static void mTXChMatrix (MTXCh *x, t_symbol *s, + int argc, t_atom *argv) +{ + int rows = atom_getint (argv++); + int columns = atom_getint (argv++); + int size = rows * columns; + int in_size = argc-2; + int n; + + + /* size check */ + if (!size) + post("mtx_circular_harmonics: invalid dimensions"); + else if (in_sizel!=columns) { + deleteMTXChdata(x); + x->l=columns; + allocMTXChdata(x); + } + for (n=0;nl;n++) { + x->phi[n]=(double) atom_getfloat(argv+n); + } + + if (x->wc!=0) { + chebyshev12(x->phi, x->wc); + in_size=x->l*(2*x->nmax+1); + SETFLOAT(x->list_ch,(float)x->l); + SETFLOAT(x->list_ch+1,(float)(2*x->nmax+1)); + for (n=0;nlist_ch+n+2,(float)x->wc->t[n]); + mTXChBang(x); + } + else + post("mtx_circular_harmonics: memory error, no operation"); + } + + } void mtx_spherical_harmonics_setup (void) @@ -139,6 +240,22 @@ void mtx_spherical_harmonics_setup (void) class_addmethod (mtx_spherical_harmonics_class, (t_method) mTXShMatrix, gensym("matrix"), A_GIMME,0); } + +void mtx_circular_harmonics_setup (void) +{ + mtx_circular_harmonics_class = class_new + (gensym("mtx_circular_harmonics"), + (t_newmethod) newMTXCh, + (t_method) deleteMTXChdata, + sizeof (MTXCh), + CLASS_DEFAULT, A_GIMME, 0); + class_addbang (mtx_circular_harmonics_class, (t_method) mTXChBang); + class_addmethod (mtx_circular_harmonics_class, (t_method) mTXChMatrix, gensym("matrix"), A_GIMME,0); +} + void iemtx_spherical_harmonics_setup(void){ + mtx_circular_harmonics_setup(); mtx_spherical_harmonics_setup(); } + + diff --git a/src/mtx_spherical_harmonics/chebyshev12.c b/src/mtx_spherical_harmonics/chebyshev12.c index 462341b..64af232 100644 --- a/src/mtx_spherical_harmonics/chebyshev12.c +++ b/src/mtx_spherical_harmonics/chebyshev12.c @@ -1,5 +1,5 @@ /* - Evaluates all circular harmonics + Evaluates all fully normalized circular harmonics at the angles phi up to the order nmax. using the recurrence for the Chebyshev polynomials of the first and second kind @@ -42,6 +42,8 @@ void chebyshev12(double *phi, Cheby12WorkSpace *wc) { const int incr=2*wc->nmax+1; double *cosphi; double *sinphi; + const double oneoversqrt2pi=1.0/sqrt(2.0*M_PI); + const double oneoversqrtpi=1.0/sqrt(M_PI); // memory allocation if ((wc!=0)&&(phi!=0)) { if ((cosphi=(double*)calloc(wc->l,sizeof(double)))==0) { @@ -56,9 +58,9 @@ void chebyshev12(double *phi, Cheby12WorkSpace *wc) { cosphi[l]=cos(phi[l]); sinphi[l]=sin(phi[l]); // initial value T_0=1 - wc->t[l0]=1; - wc->t[l0+1]=cosphi[l]; - wc->t[l0-1]=sinphi[l]; + wc->t[l0]=oneoversqrt2pi; + wc->t[l0+1]=cosphi[l]*oneoversqrtpi; + wc->t[l0-1]=sinphi[l]*oneoversqrtpi; } // recurrence for n>1 for (n=2; n<=wc->nmax; n++) { diff --git a/src/mtx_spherical_harmonics/sharmonics.c b/src/mtx_spherical_harmonics/sharmonics.c index 5d6f853..dbfcb18 100644 --- a/src/mtx_spherical_harmonics/sharmonics.c +++ b/src/mtx_spherical_harmonics/sharmonics.c @@ -43,7 +43,7 @@ static void sharmonics_initlegendrenormlzd(SHWorkSpace *ws) { } } -// multiplying Chebyshev sin/cos to the preliminary result +// multiplying normalized Chebyshev sin/cos to the preliminary result // Y_n^m(phi,theta) = Y_n^m(theta) * T_m(phi) // ny0 and nt0 denote where the position (n,m)=(n,0) or m=0 is in the arrays // ly0 and lt0 denote the starting position for one vertex in the arrays diff --git a/src/mtx_spherical_harmonics/sharmonics_normalization.c b/src/mtx_spherical_harmonics/sharmonics_normalization.c index ca97292..d7fe00a 100644 --- a/src/mtx_spherical_harmonics/sharmonics_normalization.c +++ b/src/mtx_spherical_harmonics/sharmonics_normalization.c @@ -27,8 +27,14 @@ SHNorml *sharmonics_normalization_new (const size_t nmax) { wn=0; } else { + /* + deprecated: // computing N_n^m for m=0, wrongly normalized wn->n[0]=sqrt(1/(2*M_PI)); + */ + + // computing N_n^m for m=0, + wn->n[0]=oneoversqrt2; for (n=1,n0=1; n<=nmax; n++) { wn->n[n0]=wn->n[0] * sqrt(2*n+1); n0+=n+1; @@ -40,11 +46,14 @@ SHNorml *sharmonics_normalization_new (const size_t nmax) { } n0+=n+1; } + /* + deprecated: // correcting normalization of N_n^0 for (n=0,n0=0; n<=nmax; n++) { wn->n[n0]*=oneoversqrt2; n0+=n+1; } + */ } } return wn; -- cgit v1.2.1