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

4 files changed, 135 insertions, 7 deletions

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;