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

5 files changed, 321 insertions, 7 deletions

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 <nmax> 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;n<x->l;n++) {
         x->phi[n]=(double) atom_getfloat(argv+n);
@@ -127,6 +142,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_size<size) 
+    post("mtx_circular_harmonics: sparse matrix not yet supported: use \"mtx_check\"");
+  else if ((rows!=1)||(columns<1))
+     post("mtx_circular_harmonics: 1 X L matrix expected with phi vector, but got more rows/no entries");
+  else {
+     if (x->l!=columns) {
+        deleteMTXChdata(x);
+        x->l=columns;
+        allocMTXChdata(x);
+     }
+     for (n=0;n<x->l;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;n<in_size; n++) 
+           SETFLOAT(x->list_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)
 {
   mtx_spherical_harmonics_class = class_new 
@@ -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;