rateconv.c

AAC音频解码算法程序

/**********************************************************************
audio sample rate converter

$Id: rateconv.c,v 1.1 2000/01/17 21:53:04 menno Exp $

Source file: rateconv.c

Authors:
NM    Nikolaus Meine, Uni Hannover (c/o Heiko Purnhagen)
HP    Heiko Purnhagen, Uni Hannover <[e-mail]>

Changes:
xx-jun-98   NM      resamp.c
18-sep-98   HP      converted into module
03-aug-99   NM      now even faster ...
04-nov-99   NM/HP   double ratio, htaps1 /= 2
11-nov-99   HP      improved RateConvInit() debuglevel
**********************************************************************/

/*
 * Sample-rate converter
 *
 * Realized in three steps:
 *
 * 1. Upsampling by factor two while doing appropriate lowpass-filtering.
 *    This is done by using an FFT-based convolution algorithm with a multi-tap
 *    Kaiser-windowed lowpass filter.
 *    If the cotoff-frequency is less than 0.5, only lowpass-filtering without
 *    the upsampling is done.
 * 2. Upsampling by factor 128 using a 15 tap Kaiser-windowed lowpass filter
 *    (alpha=12.5) and conventional convolution algorithm.
 *    Two values (the next neighbours) are computed for every sample needed.
 * 3. Linear interpolation between the two bounding values.
 *
 * Stereo and mono data is supported.
 * Up- and downsampling of any ratio is possible.
 * Input and output file format is Sun-audio.
 *
 * Written by N.Meine, 1998
 *
 */

/* Multi channel data is interleaved: l0 r0 l1 r1 ... */
/* Total number of samples (over all channels) is used. */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#include "aacenc.h"
#include "rateconv.h"

/* ---------- declarations ---------- */

#ifndef min
#define min(a,b) ((a) < (b) ? (a) : (b))
#endif
#ifndef min
#define max(a,b) ((a) > (b) ? (a) : (b))
#endif

/* ---------- declarations (structures) ---------- */

#define real	double	/* float seems to be precise enough, but it	*/
			/* doesn't run much faster than double		*/

typedef struct {
	real	re,im;
} complex;


struct RCBufStruct		/* buffer handle */
{
  int numChannel;		/* number of channels */
  long currentSampleIn;		/* number of input samples */
  long currentSampleOut;	/* number of output samples */

  real *wtab;
  complex *x;
  real *tp1;
  real *tp2;
  complex *ibuf;

  long nfn;
  long adv;
  long istart;
  long numode;
  long fm;
  double p1;
  double d1;

  long outSize;
  float *out;
};


/* ---------- variables ---------- */

static int RCdebugLevel = 0;	/* debug level */


#define deftaps1	485

#define ov2exp	7		/* Oversampling exponent in the second step */
#define ov2fac	(1<<ov2exp)	/* Oversampling factor in the second step */
#define htaps2	7		/* Num. taps of the 2nd filter: 2*htaps2+1 */
#define alpha2	12.5		/* Kaiser-parameter for the second filter */
#define cutoff2	0.97		/* cutoff-frequency for the second filter */

#define cexp	(CACHE1-(sizeof(real)==4 ? 3 : 4 ))
#define cache	(1<<cexp)

#define pi	3.1415926535897932
#define zpi	6.2831853071795965


/* ---------- local functions ---------- */


static void mkeiwutab(real *st,long fn)
{
	register int	i,k;
	register double	phi;

	for (k=fn; k; k>>=1)
	{
		phi = zpi/((double) k);
		for (i=0; i<k; i++) st[i]=(real) sin(phi*((double) i));
		st+=k;
	}
}

static void r4fftstep(register real *x0,register long n,register real *t)
{
	register long		n4 = n/4;
	register real		*x1 = x0+2*n4;
	register real		*x2 = x1+2*n4;
	register real		*x3 = x2+2*n4;
	register real		*y;
	register real		wr,wi,ar,ai,br,bi,cr,ci,dr,di,hr,hi;
	register long		i;


	if ( n>16 )
	{
		y=t+6*n4;
		r4fftstep(x0,n4,y);
		r4fftstep(x1,n4,y);
		r4fftstep(x2,n4,y);
		r4fftstep(x3,n4,y);
	}
	else
	{
		if ( n==8 )
		{
			for (y=x0; y<x0+16; y+=4)
			{
				ar=y[0]; ai=y[1];
				br=y[2]; bi=y[3];
				y[0]=ar+br;
				y[1]=ai+bi;
				y[2]=ar-br;
				y[3]=ai-bi;
			}
		}
		else
		{
			for (y=x0; y<x0+32; y+=8)
			{
				ar=y[0]; ai=y[1];
				br=y[2]; bi=y[3];
				cr=y[4]; ci=y[5];
				dr=y[6]; di=y[7];
				hr=ar+cr; wr=ar-cr; hi=ai+ci; wi=ai-ci;
				ar=br+dr; cr=br-dr; ai=bi+di; ci=bi-di;
				y[0]=hr+ar; y[1]=hi+ai;
				y[2]=wr+ci; y[3]=wi-cr;
				y[4]=hr-ar; y[5]=hi-ai;
				y[6]=wr-ci; y[7]=wi+cr;
			}
		}
	}

	y=t+n4;
	for (i=0; i<n4; i++)
	{
		ar=x0[0]; ai=x0[1];
		wr=y[  i]; wi=t[  i]; hr=x1[0]; hi=x1[1]; br=wr*hr+wi*hi; bi=wr*hi-wi*hr;
		wr=y[2*i]; wi=t[2*i]; hr=x2[0]; hi=x2[1]; cr=wr*hr+wi*hi; ci=wr*hi-wi*hr;
		wr=y[3*i]; wi=t[3*i]; hr=x3[0]; hi=x3[1]; dr=wr*hr+wi*hi; di=wr*hi-wi*hr;
		hr=ar+cr; wr=ar-cr; hi=ai+ci; wi=ai-ci;
		ar=br+dr; cr=br-dr; ai=bi+di; ci=bi-di;
		x0[0]=hr+ar; x0[1]=hi+ai;
		x1[0]=wr+ci; x1[1]=wi-cr;
		x2[0]=hr-ar; x2[1]=hi-ai;
		x3[0]=wr-ci; x3[1]=wi+cr;
		x0+=2; x1+=2; x2+=2; x3+=2;
	}
}


static void r4ifftstep(register real *x0,register long n,register real *t)
{
	register long		n4 = n/4;
	register real		*x1 = x0+2*n4;
	register real		*x2 = x1+2*n4;
	register real		*x3 = x2+2*n4;
	register real		*y;
	register real		wr,wi,ar,ai,br,bi,cr,ci,dr,di,hr,hi;
	register long		i;

	if ( n>16 )
	{
		y=t+6*n4;
		r4ifftstep(x0,n4,y);
		r4ifftstep(x1,n4,y);
		r4ifftstep(x2,n4,y);
		r4ifftstep(x3,n4,y);
	}
	else
	{
		if ( n==8 )
		{
			for (y=x0; y<x0+16; y+=4)
			{
				ar=y[0]; ai=y[1];
				br=y[2]; bi=y[3];
				y[0]=ar+br;
				y[1]=ai+bi;
				y[2]=ar-br;
				y[3]=ai-bi;
			}
		}
		else
		{
			for (y=x0; y<x0+32; y+=8)
			{
				ar=y[0]; ai=y[1];
				br=y[2]; bi=y[3];
				cr=y[4]; ci=y[5];
				dr=y[6]; di=y[7];
				hr=ar+cr; wr=ar-cr; hi=ai+ci; wi=ai-ci;
				ar=br+dr; cr=br-dr; ai=bi+di; ci=bi-di;
				y[0]=hr+ar; y[1]=hi+ai;
				y[2]=wr-ci; y[3]=wi+cr;
				y[4]=hr-ar; y[5]=hi-ai;
				y[6]=wr+ci; y[7]=wi-cr;
			}
		}
	}

	y=t+n4;
	for (i=0; i<n4; i++)
	{
		ar=x0[0]; ai=x0[1];
		wr=y[  i]; wi=t[  i]; hr=x1[0]; hi=x1[1]; br=wr*hr-wi*hi; bi=wr*hi+wi*hr;
		wr=y[2*i]; wi=t[2*i]; hr=x2[0]; hi=x2[1]; cr=wr*hr-wi*hi; ci=wr*hi+wi*hr;
		wr=y[3*i]; wi=t[3*i]; hr=x3[0]; hi=x3[1]; dr=wr*hr-wi*hi; di=wr*hi+wi*hr;
		hr=ar+cr; wr=ar-cr; hi=ai+ci; wi=ai-ci;
		ar=br+dr; cr=br-dr; ai=bi+di; ci=bi-di;
		x0[0]=hr+ar; x0[1]=hi+ai;
		x1[0]=wr-ci; x1[1]=wi+cr;
		x2[0]=hr-ar; x2[1]=hi-ai;
		x3[0]=wr+ci; x3[1]=wi-cr;
		x0+=2; x1+=2; x2+=2; x3+=2;
	}
}


static void perm4(real *x,long p)
{
	register real	a,b,c,d;
	register long	i,j,k,l;

	l=1<<(p-1);
	for (i=j=0; i<(2<<p); i+=2)
	{
		if ( j>i )
		{
			a=x[i]; b=x[i+1]; c=x[j]; d=x[j+1]; x[j]=a; x[j+1]=b; x[i]=c; x[i+1]=d;	
		}
		j+=l;
		for (k=l<<2; j&k; k>>=2) j=(j^=k)+(k>>4);
	}
}

static void perm42(real *x,long p)
{
	register real	*y,*z;
	register long	i,n;

	n=1<<p;
	perm4(x  ,p-1);
	perm4(x+n,p-1);

	y=x+2*n;
	z=x+n;
	for (i=0; i<n; i+=2)
	{
		y[2*i  ]=x[i  ];
		y[2*i+1]=x[i+1];
		y[2*i+2]=z[i  ];
		y[2*i+3]=z[i+1];
		x[i  ]=y[i  ];
		x[i+1]=y[i+1];
	}
	for (; i<2*n; i+=2)
	{
		x[i  ]=y[i  ];
		x[i+1]=y[i+1];
	}
}

static void r4fft(real *t,real *x,long p)
{
	register long	n  = 1<<p;

	if ( p&1 ) perm42(x,p); else perm4(x,p);
	r4fftstep(x,n,t);
}

static void r4ifft(real *t,real *x,long p)
{
	register long	n  = 1<<p;

	if ( p&1 ) perm42(x,p); else perm4(x,p);
	r4ifftstep(x,n,t);
}


static double bessel(double x)
{
  register double	p,s,ds,k;

  x*=0.5;
  k=0.0;
  p=s=1.0;
  do
    {
      k+=1.0;
      p*=x/k;
      ds=p*p;
      s+=ds;
    }
  while ( ds>1.0e-17*s ); 

  return s;
}

static void mktp1(complex *y,double fg,double a,long fn,long htaps1)
{
  register long	i;
  register double	f,g,x,px;

  y[0].re=fg;
  y[0].im=0.0;

  f=1.0/bessel(a);
  g=1.0/(htaps1+0.5);
  g*=g;

  for (i=1; i<=htaps1; i++)
    {
      x=(double) i;
      px=pi*x;
      y[i].re=(real) f*bessel(a*sqrt(1.0-g*x*x))*sin(fg*px)/px;
      y[i].im=(real) 0.0;
      y[fn-i]=y[i];
    }
  for (; i<fn-htaps1; i++) y[i].re=y[i].im=0.0;
}

static void mktp2(real *y,double fg,double a)
{
  register long	i;
  register double	f,g,x,px;

  y[0]=fg;

  f=1.0/bessel(a);
  g=1.0/(((double) ((htaps2+1)*ov2fac))-0.5);
  g*=g;

  for (i=1; i<(htaps2+1)*ov2fac; i++)
    {
      x=(double) i;
      px=(pi/ov2fac)*x;
      y[i]=(real) f*bessel(a*sqrt(1.0-g*x*x))*sin(fg*px)/px;
    }
  for (; i<(htaps2+2)*ov2fac; i++) y[i]=0.0;
}

static void ctp2_m(real *x,double p,real *k,float **out)
{
  register long	n,i,j;
  register double	h,l,r;
  register real	*ka,*kb;
  register real	*xa,*xb;

  h=ov2fac*p;
  i=(long) h;
  h-=(double) i;

  j=i&(ov2fac-1); i>>=ov2exp;