zjf.h

这个程序是我写的语音编码G.729 以编译通过

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

#define B0		 0.46363718
#define B1		-0.92724705
#define B2	 	 0.46363718
#define A1		 1.9059465
#define A2		-0.9114024
#define pi 3.14159
double b30[31] = {
 0.898517,
 0.769271,
 0.448635,
 0.095915,
-0.134333,
-0.178528,
-0.084919,
 0.036952,
 0.095533,
 0.068936,
-0.000000,
-0.050404,
-0.050835,
-0.014169,
 0.023083,
 0.033543,
 0.016774,
-0.007466,
-0.019340,
-0.013755,
0.000000,
 0.009400,
 0.009029,
 0.002381,
-0.003658,
-0.005027,
-0.002405,
 0.001050,
 0.002780,
 0.002145,
 0.000000
};
void preprocessing(double *y, double *x, int start, int length)
{
	int i;
	int n0, n1, n2;

	for (i=start; i<(start+length) ; i++)
	{
		n0 = (i+240)   % 240;
		n1 = (i-1+240) % 240;
		n2 = (i-2+240) % 240;
		y[n0] = (B0*x[n0] + B1*x[n1] + B2*x[n2]+
				 A1*y[n1] + A2*y[n2]);
	}
}
void response(double a_h[10], double r[40], double h[40], double x[40], double ew[10])
{int n, i;
	double gamma;
	double ir, ts;
    h[0]=1.0;
	x[0] =r[0];
	for (n=1; n<40; n++)
	{ir = 0;
	 ts = 0;
	 gamma = .75;
	 for (i=1; i<=10; i++)
	 { if (n-i>=0)
	   {ir += a_h[i-1]*gamma*h[n-i];
		ts += a_h[i-1]*gamma*x[n-i];
		}	
		else
		{   /*  ir????*/
		 ts += a_h[i-1]*gamma*ew[10-i];
		}
		gamma *= .75;
	 }
	h[n] = -ir;
	x[n] = r[n]-ts;
	}
}
 void correlation(double x[40], double h[40], double x_b[40])
{
	int n,l;

	for (n=0; n<40; n++)	
		for (l=n, x_b[n]=0; l<=39; l++)
			x_b[n] += x[l]*h[l-n];
}
void pitchdelay(int tmin, int tmax, double x_b[40], double u[200], int *T, int *frac, int subframe1, double r[40])
{
	int n, k, t;
	int i, j;
	int zero_index = 160;
	double R, Rmax;	
	double u_kt[40];

	for (n=0; n<40; n++)
		u[zero_index+n] = r[n];

	*T = tmin;
	*frac = 0;
	Rmax = 0;
    for (n=0; n<40; n++)	
		Rmax += x_b[n]*u[n-tmin+zero_index];	

	for (j=tmin+1; j<=tmax; j++)
	{	// Equation A.7
		for (n=0, R=0; n<40; n++)			
			R += x_b[n]*u[n-j+zero_index];

		if (R > Rmax)
		{
			*T = j;
			Rmax = R;
		}
	}                                      
	
	// Check for interpolation conditions
	if ((subframe1 & (*T<85) ) | !subframe1)
	{		
		/* If this is the 1st subframe and the search index < 85, or if this is the 2nd subframe,
		   then the 3 fractional parts T-1/3, T and T+1/3 are interpolated and the maximum
		   correlation term is picked.
			T-1/3: t=2, k=T-1
			T    : t=0, k=T
			T+1/3: t=1, k=T			
		  */
		for (t=0; t<=2; t++)
		{
			k = *T - (t==2);

			// Equation (A.8)
			for (n=0; n<40; n++)
			{
				for (i=0, u_kt[n]=0; i<=9; i++)
				{
					u_kt[n] += u[n-k+i+zero_index]*b30[t+3*i];
					u_kt[n] += u[n-k+1+i+zero_index]*b30[3-t+3*i];
				}
			}
			
			for (n=0, R=0; n<40; n++)
				R += x_b[n]*u_kt[n];
			
			if (t==0)
			{
				*frac = t;
				Rmax = R;
			}			
			else if (R>Rmax)
			{							
				*frac = t;
				Rmax = R;
			}
		}
				
		if (*frac>1)
		{
			*frac = *frac - 3;
			*T = *T + 1;
		}
	}
}
void adaptivecodebookvector(double u[200], int T1, int frac1,double v[40])
{
	int n, i;
	int k, t;
	int zero_index;

	if (frac1 == -1)
	{
		t = 2;
		k = T1 - 1;
	}
	else
	{
		t = frac1;
		k = T1;
	}
	zero_index = 160 - k;
	for (n=0; n<40; n++)
	{
		for (i=0, v[n]=0; i<=9; i++)
			v[n] += u[n+zero_index+i]*b30[t+3*i];
		for (i=0; i<=9; i++)
			v[n] += u[n+zero_index+1+i]*b30[3-t+3*i];
	}

}
void filteredadaptivecodebookvector(double v[40], double h[40], double y[40])
{
	int n, i;

	for (n=0; n<40; n++)
		for (i=0, y[n]=0; i<=n; i++)
			y[n] += v[i]*h[n-i];
}

double adaptive0gain(double y[40],double x[40])
{ double gainp,gainmax=1.2;
  double sum0=0,sum1=0;
  int n;
  for (n=0;n<40;n++)
  {sum0=sum0+x[n]*y[n];
   sum1=sum1+y[n]*y[n];
  }
  if (sum1==0)
	 gainp=(2*(sum0>0)-1)*gainmax;
 else
    gainp=sum0/sum1;
 if (gainp>gainmax)
      gainp=gainmax;
  if(gainp<0)
	 gainp=0;
 return gainp;	
}
double sign(double x)
{if(x>0)
  return(1.0);
 else
  return(-1.0);
}
void fixedcodebooksearchsubframe( double x[40],double h[40],double y[40], double g_p, double beta, int T,double c[40])
 {
	double x_p[40], d[40];	
	double phi[40][40], phi_p[40][40];
	double CC,C0,C1,C2,E,E0,E1,E2;
	double C2perE, maxC2perE;	
	int i,j,n;
	int m0,m1,m2,m3;
	int optm0,optm1,optm2,optm3;


	// Equation (49)
	if (T<40)
		for (n=T; n<40; n++)
			h[n] += beta * h[n-T];  /* beta ???*/
   //for (n=0;n<40;n++)
   // printf("%f\n",h[n]);
	// Equation (50)
	for (n=0; n<40; n++)
		x_p[n] = x[n] - g_p * y[n];
   
	// Equation (51)
	for (i=0; i<40; i++)
		for (j=i; j<40; j++)
			for (n=j; n<40; n++)			
				phi[i][j] = h[n-i]*h[n-j];

	for (i=0; i<40; i++)
		for (j=0; j<i; j++)
			phi[i][j] = phi[j][i];
    
	// Equation (52)
	correlation(x_p, h, d);
	

	// Equation (56)
	for (i=0; i<40; i++)
		for (j=0; j<40; j++)
			phi_p[i][j] = sign(d[i])*sign(d[j])*phi[i][j];
	
	// Equation (57)
	for (i=0; i<40; i++)
		phi_p[i][i] *= 0.5;
					
	optm0 = m0 = 0;
	optm1 = m1 = 1;
	optm2 = m2 = 2;
	optm3 = m3 = 3;
    
	CC = fabs(d[m0]) + fabs(d[m1]) + fabs(d[m2]) + fabs(d[m3]);
	E = phi_p[m0][m0] + phi_p[m1][m1] + phi_p[m0][m1] + phi_p[m2][m2] + phi_p[m0][m2] 
	  + phi_p[m1][m2] + phi_p[m3][m3] + phi_p[m0][m3] + phi_p[m1][m3] + phi_p[m2][m3];
	maxC2perE = CC*CC/E;

	for (m0=0; m0<=35; m0+=5)
	{		
		C0 = fabs(d[m0]);
		E0 = phi_p[m0][m0];
		for (m1=1; m1<=36; m1+=5)
		{
			C1 = C0 + fabs(d[m1]);
			E1 = E0 + phi_p[m1][m1] + phi_p[m0][m1];
			for (m2=2; m2<=37; m2+=5)
			{
				C2 = C1 + fabs(d[m2]);
				E2 = E1 + phi_p[m2][m2] + phi_p[m0][m2] + phi_p[m1][m2];
				for (m3=3; m3<=38; m3+=5)
				{
					CC = C2 + fabs(d[m3]);
					E = E2 + phi_p[m3][m3] + phi_p[m0][m3] + phi_p[m1][m3] + phi_p[m2][m3];
					// For equation (59), E should be multiplied by 2, but this operation is saved
					//   because the scale factor is the same for all C^2/E
					C2perE = CC*CC/E;
					if (C2perE > maxC2perE)
					{
						optm0 = m0;
						optm1 = m1;
						optm2 = m2;
						optm3 = m3;	
						maxC2perE = C2perE;
					}
				}

				for (m3=4; m3<=39; m3+=5)
				{
					CC = C2 + fabs(d[m3]);
					E = E2 + phi_p[m3][m3] + phi_p[m0][m3] + phi_p[m1][m3] + phi_p[m2][m3];
					// For equation (59), E should be multiplied by 2, but this operation is saved
					//   because the scale factor is the same for all C^2/E
					C2perE = CC*CC/E;
					if (C2perE > maxC2perE)
					{
						optm0 = m0;
						optm1 = m1;
						optm2 = m2;
						optm3 = m3;
						maxC2perE = C2perE;
					}
				}
			}
		}
	}
/*if optm=3,8,--,38
	                          jx=0;
	                      else if optm=4,9,14,--39
	                               jx=1;
	                      so i think jx=((ptm+1)%5==0);*/
	                                      
   
 
	for (n=0; n<40; n++)
	     c[n] = 0;

	c[optm0] = sign(d[optm0]);
	c[optm1] = sign(d[optm1]);
	c[optm2] = sign(d[optm2]);
	c[optm3] = sign(d[optm3]);

	if (T<40)
		for (n=T; n<40; n++)
			c[n] += beta*c[n-T];
   
}

void fixedcodebookgain(double x[40],double y[40], double z[40],double c[40], double h[40], double g_p, double *g_c)
{
	double A, B, C;
    double e1=0, e2=0, e3=0, e4=0, e5=0, e6=0;
	int n, i;
    *g_c=0;
	for (n=0; n<40; n++)
		for (i=0, z[n]=0; i<=n; i++)
			z[n] += c[i]*h[n-i];

	for (n=0; n<40; n++)
	{
		e1 += x[n]*x[n];
		e2 += y[n]*y[n];
		e3 += z[n]*z[n];
		e4 += x[n]*y[n];
		e5 += x[n]*z[n];
		e6 += y[n]*z[n];
	}

	A = e3; 
	B = 2*(g_p*e6 - e5);
	C = e1 + g_p*g_p*e2 - 2*g_p*e4;
	
	if (A !=0) 
		*g_c = -B/2/A;
	else if (B != 0)
		*g_c = -C/B;
}
void memoryupdate(double u[200], double g_p, double v[40], double g_c,double c[40], double x[40], double y[40], double z[40], double ew[10])
{
	int n;
	int zero_index = 160;

	for (n=0; n<40; n++)		
	    u[zero_index+n] = g_p*v[n] + g_c*c[n];
   	 
	for (n=0; n<10; n++)
		ew[n] = x[n+30] - g_p*y[n+30] - g_c*z[n+30];

}