/*
 * Lpc port to Java.
 * Copyright (C) 1999  Christopher Edwards
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */
 
package	org.tritonus.lowlevel.gsm;
 
public class Lpc
{
	private int L_ACF[] = new int[9];



	public void Gsm_LPC_Analysis(
		short[]	so,	/* 0..159 signals       IN/OUT  */
		short[]	LARc)	/* 0..7   LARc's        OUT     */
	{
		Autocorrelation                   (so);
		Reflection_coefficients           (LARc);
		Transformation_to_Log_Area_Ratios (LARc);
		Quantization_and_coding           (LARc);
	}



	private void Autocorrelation(short[] so) /* [0..159]  IN/OUT  */
		throws IllegalArgumentException
	{
		int     i = 0, sp_index = 0;
		short   temp = 0, smax = 0, scalauto = 0;
		int     L_temp = 0;
		int     L_temp2 = 0;

		/*  Dynamic scaling of the array  s[0..159]
		 */

		/*  Search for the maximum.
		 */
		for (int k = 0; k <= 159; k++)
		{
			temp = Add.GSM_ADD( so[k], (short) 0 );
			if (temp > smax) smax = temp;
		}

		/*  Computation of the scaling factor.
		 */
		if (smax == 0)
		{
			scalauto = 0;
		}
		else
		{
			if ( ! (smax > 0) )
			{
				throw new IllegalArgumentException
					("Autocorrelation: smax = "
					 +smax+" should be > 0.");
			}
			scalauto = (short)
				(4 - Add.gsm_norm((int)(smax << 16))); /* sub(4,..) */
		}

		/*  Scaling of the array s[0...159]
		 */

		if (scalauto > 0)
		{
			if( ! (scalauto <= 4) )
			{
				throw new IllegalArgumentException
					("Autocorrelation: scalauto = "
					 +scalauto+" should be <= 4.");
			}
			switch (scalauto)
			{
			case 1:
				for (int k = 0; k <= 159; k++)
				{
					so[k] = Add.GSM_MULT_R( so[k], (short)16384 );
				}
				break;

			case 2:
				for (int k = 0; k <= 159; k++)
				{
					so[k] = Add.GSM_MULT_R( so[k], 
								(short) (16384 >> 1) );
				}
				break;

			case 3:
				for (int k = 0; k <= 159; k++)
				{
					so[k] = Add.GSM_MULT_R( so[k], 
								(short) (16384 >> 2) );
				}
				break;

			case 4:
				for (int k = 0; k <= 159; k++)
				{
					so[k] = Add.GSM_MULT_R( so[k], 
								(short) (16384 >> 3) );
				}
				break;
			}
		}

		/*  Compute the L_ACF[..].
		 */
		short[]  sp = so;
		short    sl = sp[sp_index];

		// Zero out L_ACF
		int[] temp_arr = {0, 0, 0, 0, 0, 0, 0, 0, 0};
		System.arraycopy(temp_arr , 0, L_ACF, 0, L_ACF.length);	    

		L_ACF[0] += (int)(sl * sp[ (sp_index - 0) ]);
        
		sl = sp[++sp_index];	    
		L_ACF[0] += (int)(sl * sp[ (sp_index - 0) ]);
		L_ACF[1] += (int)(sl * sp[ (sp_index - 1) ]);
	    
		sl = sp[++sp_index];	    
		L_ACF[0] += (int)(sl * sp[ (sp_index - 0) ]);
		L_ACF[1] += (int)(sl * sp[ (sp_index - 1) ]);
		L_ACF[2] += (int)(sl * sp[ (sp_index - 2) ]);
	    
		sl = sp[++sp_index];	    
		L_ACF[0] += (int)(sl * sp[ (sp_index - 0) ]);
		L_ACF[1] += (int)(sl * sp[ (sp_index - 1) ]);
		L_ACF[2] += (int)(sl * sp[ (sp_index - 2) ]);
		L_ACF[3] += (int)(sl * sp[ (sp_index - 3) ]);
	    
		sl = sp[++sp_index];
		L_ACF[0] += (int)(sl * sp[ (sp_index - 0) ]);
		L_ACF[1] += (int)(sl * sp[ (sp_index - 1) ]);
		L_ACF[2] += (int)(sl * sp[ (sp_index - 2) ]);
		L_ACF[3] += (int)(sl * sp[ (sp_index - 3) ]);
		L_ACF[4] += (int)(sl * sp[ (sp_index - 4) ]);
	
		sl = sp[++sp_index];
		L_ACF[0] += (int)(sl * sp[ (sp_index - 0) ]);
		L_ACF[1] += (int)(sl * sp[ (sp_index - 1) ]);
		L_ACF[2] += (int)(sl * sp[ (sp_index - 2) ]);
		L_ACF[3] += (int)(sl * sp[ (sp_index - 3) ]);
		L_ACF[4] += (int)(sl * sp[ (sp_index - 4) ]);
		L_ACF[5] += (int)(sl * sp[ (sp_index - 5) ]);
	
		sl = sp[++sp_index];
		L_ACF[0] += (int)(sl * sp[ (sp_index - 0) ]);
		L_ACF[1] += (int)(sl * sp[ (sp_index - 1) ]);
		L_ACF[2] += (int)(sl * sp[ (sp_index - 2) ]);
		L_ACF[3] += (int)(sl * sp[ (sp_index - 3) ]);
		L_ACF[4] += (int)(sl * sp[ (sp_index - 4) ]);
		L_ACF[5] += (int)(sl * sp[ (sp_index - 5) ]);
		L_ACF[6] += (int)(sl * sp[ (sp_index - 6) ]);
	    
		sl = sp[++sp_index];
		L_ACF[0] += (int)(sl * sp[ (sp_index - 0) ]);
		L_ACF[1] += (int)(sl * sp[ (sp_index - 1) ]);
		L_ACF[2] += (int)(sl * sp[ (sp_index - 2) ]);
		L_ACF[3] += (int)(sl * sp[ (sp_index - 3) ]);
		L_ACF[4] += (int)(sl * sp[ (sp_index - 4) ]);
		L_ACF[5] += (int)(sl * sp[ (sp_index - 5) ]);
		L_ACF[6] += (int)(sl * sp[ (sp_index - 6) ]);
		L_ACF[7] += (int)(sl * sp[ (sp_index - 7) ]);
	
		sl = sp[++sp_index];

		for (i = sp_index; i < 160; ++i)
		{

			sl = sp[i];

			L_ACF[0] += (int)(sl * sp[ (i - 0) ]);
			L_ACF[1] += (int)(sl * sp[ (i - 1) ]);
			L_ACF[2] += (int)(sl * sp[ (i - 2) ]);
			L_ACF[3] += (int)(sl * sp[ (i - 3) ]);
			L_ACF[4] += (int)(sl * sp[ (i - 4) ]);
			L_ACF[5] += (int)(sl * sp[ (i - 5) ]);
			L_ACF[6] += (int)(sl * sp[ (i - 6) ]);
			L_ACF[7] += (int)(sl * sp[ (i - 7) ]);
			L_ACF[8] += (int)(sl * sp[ (i - 8) ]);
		}

		for (int k = 0; k < 9; k++)
		{
			L_ACF[k] <<= 1;
		}

		/*   Rescaling of the array s[0..159]
		 */
		if (scalauto > 0)
		{
			if( ! (scalauto <= 4) )
			{
				throw new IllegalArgumentException
					("Autocorrelation: scalauto = "
					 +scalauto+" should be <= 4.");
			}

			for (int k = 0; k < 160; k++)
			{
				so[k] <<= scalauto;
			}
		}
	}



	private void Reflection_coefficients(short[] r  /* 0...7   OUT */)
		throws IllegalArgumentException
	{
		short   temp = 0;
		int 	ltmp = 0;
		int	r_index = 0;

		short[] ACF = new short[9]; /* 0..8 */
		short[] P   = new short[9]; /* 0..8 */
		short[] K   = new short[9]; /* 2..8 */

		/*  Schur recursion with 16 bits arithmetic.
		 */

		if (L_ACF[0] == 0)
		{    /* everything is the same. */
			for (int i = 0; i < 8; i++)
			{
				r[i] = 0;
			}
			return;
		}

		if( L_ACF[0] == 0 )
		{
			throw new IllegalArgumentException
				("Reflection_coefficients: L_ACF[0] = "
				 +L_ACF[0]+" should not = 0.");
		}

		temp = Add.gsm_norm( L_ACF[0] );

		if ( ! (temp >= 0 && temp < 32) )
		{
			throw new IllegalArgumentException
				("Reflection_coefficients: temp = "
				 +temp+" should be >= 0 and < 32.");
		}


		/* ? overflow ? */
		for (int i = 0; i <= 8; i++)
		{
			ACF[i] = Add.SASR( L_ACF[i] << temp, 16 );
		}

		/*   Initialize array P[..] and K[..] for the recursion.
		 */

		System.arraycopy(ACF , 0, K, 0, 7);
	    
		System.arraycopy(ACF , 0, P, 0, 8);

		/*   Compute reflection coefficients
		 */
		for (int n = 1; n <= 8; n++, r_index++)
		{

			temp = P[1];
			temp = Add.GSM_ABS(temp);
			if (P[0] < temp) {
				for (int i = n; i < 8; i++)
				{
					r[i] = 0;
				}
				return;
			}

			r[r_index] = Add.gsm_div( temp, P[0] );

			if ( ! (r[r_index] >= 0)  )
			{
				throw new IllegalArgumentException
					("Reflection_coefficients: r["+r_index+"] = "
					 +r[r_index]+" should be >= 0");
			}

			if (P[1] > 0)
			{
				/* r[n] = sub(0, r[n]) */
				r[r_index] = (short) (-(r[r_index])); 
		        }

			if ( r[r_index] == Gsm_Def.MIN_WORD )
			{
				throw new IllegalArgumentException
					("Reflection_coefficients: r["+r_index+"] = "
					 +r[r_index]+" should not be "+ Gsm_Def.MIN_WORD );
			}
			if (n == 8)
				return;

			/*  Schur recursion
			 */
			temp = Add.GSM_MULT_R( P[1], r[r_index] );
			P[0] = Add.GSM_ADD(    P[0], temp );

			for (int m = 1; m <= 8 - n; m++)
			{
				temp     = Add.GSM_MULT_R( K[ m   ], r[r_index] );
				P[m]     = Add.GSM_ADD(    P[ m+1 ], temp );

				temp     = Add.GSM_MULT_R( P[ m+1 ], r[r_index] );
				K[m]     = Add.GSM_ADD(    K[ m   ], temp );
			}
		}
	}



	/* 4.2.6 */

	/*
	 *  The following scaling for r[..] and LAR[..] has been used:
	 *
	 *  r[..]   = integer( real_r[..]*32768. ); -1 <= real_r < 1.
	 *  LAR[..] = integer( real_LAR[..] * 16384 );
	 *  with -1.625 <= real_LAR <= 1.625
	 */
	private void Transformation_to_Log_Area_Ratios(short[] r /* 0..7 IN/OUT */)
		throws IllegalArgumentException
	{
		short   temp;

		/* Computation of the LAR[0..7] from the r[0..7]
		 */
		for (int i = 0; i < 8; i++)
		{

			temp = r[i];
			temp = Add.GSM_ABS(temp);

			if( ! (temp >= 0) )
			{
				throw new IllegalArgumentException
					("Transformation_to_Log_Area_Ratios: temp = "
					 +temp+" should be >= 0 ");
			}

			if (temp < 22118)
			{
				temp >>= 1;
			}
			else if (temp < 31130)
			{

				if( ! (temp >= 11059) )
				{
			                throw new IllegalArgumentException
						("Transformation_to_Log_Area_Ratios: temp = "
						 +temp+" should be >= 11059 ");
				}

				temp -= 11059;
			}
			else
			{
				if( ! (temp >= 26112) )
				{
			                throw new IllegalArgumentException
						("Transformation_to_Log_Area_Ratios: temp = "
						 +temp+" should be >= 26112 ");
				}

				temp -= 26112;
				temp <<= 2;
			}

			r[i] = (short) (r[i] < 0 ? -temp : temp);

			if( r[i] == Gsm_Def.MIN_WORD )
			{
		                throw new IllegalArgumentException
					("Transformation_to_Log_Area_Ratios: r["+i+"] = "
					 +r[i]+" should not be = " +Gsm_Def.MIN_WORD);
			}
		}
	}



	/* 4.2.7 */

	/*  This procedure needs four tables; the following equations
	 *  give the optimum scaling for the constants:
	 *
	 *  A[0..7] = integer( real_A[0..7] * 1024 )
	 *  B[0..7] = integer( real_B[0..7] *  512 )
	 *  MAC[0..7] = maximum of the LARc[0..7]
	 *  MIC[0..7] = minimum of the LARc[0..7]
	 */
	private void Quantization_and_coding(short[] LAR /* [0..7] IN/OUT  */)
	{
		int     index = 0;

		STEP2(  20480,      0,  31, -32, LAR, index++ );
		STEP2(  20480,      0,  31, -32, LAR, index++ );
		STEP2(  20480,   2048,  15, -16, LAR, index++ );
		STEP2(  20480,  -2560,  15, -16, LAR, index++ );

		STEP2(  13964,     94,   7,  -8, LAR, index++ );
		STEP2(  15360,  -1792,   7,  -8, LAR, index++ );
		STEP2(   8534,   -341,   3,  -4, LAR, index++ );
		STEP2(   9036,  -1144,   3,  -4, LAR, index++ );

	}



	private void STEP2(int A, int B, int MAC, int MIC, 
			   short[] LAR, int index) 
	{
		short temp = 0;

		temp = Add.GSM_MULT( (short)A, LAR[index] );  
		temp = Add.GSM_ADD( temp, (short)B ); 
		temp = Add.GSM_ADD( temp, (short) 256 );
		temp = Add.SASR( temp, 9 );
		LAR[index] = (short) (temp > MAC ? MAC - MIC : 
				      (temp < MIC ? 0 : temp - MIC)); 
	}
}