view release on metacpan or  search on metacpan

Data/fft.c  view on Meta::CPAN

           sqrt(SQR(xReal(x,i))+SQR(xImag(x,i))),
           180*(atan2(xImag(x,i),xReal(x,i)))/M_PI);
  }
}

void
Audio_difference(int n, float *a, float *b)
{
 int i;
 for (i=0; i < n; i++)
  b[i] = a[i+1] - a[i];
}


#ifdef ORIGINAL
void
Audio_autocorrelation(int N, float *x,unsigned p,float *r)
{
 int k;
 for (k=0; k <= p; k++)
  {int n;
   r[k] = 0.0;
   for (n=0; n < N-k; n++)
    r[k] += x[n]*x[n+k];
   /* What the heck is this doing ? ! 
      - well it is normalizing each entry which makes a kind 
        of sense when using it to spot period
	but makes result useless for LPC via Audio_durbin() below.
    */
   r[k] /= (float) (N-k);  
  }
}

#else 

void
Audio_autocorrelation(int N, float *x,unsigned p,float *r)
{
 unsigned i; 
 for (i = 0; i <= p; i++) 
  {
   unsigned j;
   float sum = 0.0;
   for(j = 0; j < N - i; j++) 
    sum += x[j] * x[j + i];
   r[i]= sum;
  }
}

#endif

#ifdef __GNUC__

void
Audio_durbin(int NUM_POLES, float *R,float *aa)
{double E[NUM_POLES+1];
 double k[NUM_POLES+1];
 double a[NUM_POLES+1][NUM_POLES+1];
 double G = R[0];
 int i;
 /* Set everything to NaN */
 memset(a,-1,sizeof(a));
 memset(k,-1,sizeof(k));
 memset(E,-1,sizeof(E));
 E[0] = R[0];
 for (i=1; i <= NUM_POLES; i++)
  {int j;
   k[i] = 0.0;
   for (j=1; j < i; j++)
    k[i] += a[j][i-1]*R[i-j];
   k[i] -= R[i];
   k[i] /= E[i-1];
   a[i][i] = -k[i];
   for (j=1; j < i; j++)
    {
     a[j][i] = a[j][i-1]+k[i]*a[i-j][i-1];
    }
   E[i] = (1.0 - k[i]*k[i])*E[i-1];
  }
 for (i=1; i <= NUM_POLES; i++)
  {
   aa[i] = a[i][NUM_POLES];
   G -= aa[i]*R[i];
  }
 /* Return gain as element zero */
 if (G < 0.0)
  G = -G;
 aa[0] = sqrt(G);
}


#else
#define a(I,J) *(A+N*(I)+(J))
void
Audio_durbin(int NUM_POLES, float *R,float *aa)
{
 int N = NUM_POLES+1;
 double *E;
 double *k;
 double *A;
 double G = R[0];
 int i;
 Newz('D',E,N,double);
 Newz('D',k,N,double);
 Newz('D',A,N*N,double);
 E[0] = R[0];
 for (i=1; i <= NUM_POLES; i++)
  {int j;
   k[i] = 0.0;
   for (j=1; j < i; j++)
    k[i] += a(j,i-1)*R[i-j];
   k[i] -= R[i];
   k[i] /= E[i-1];
   a(i,i) = -k[i];
   for (j=1; j < i; j++)
    {
     a(j,i) = a(j,i-1)+k[i]*a(i-j,i-1);
    }
   E[i] = (1.0 - k[i]*k[i])*E[i-1];
  }
 for (i=1; i <= NUM_POLES; i++)