computeLP.c 12.1 KB
Newer Older
johan's avatar
johan committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
/*
 computeLP.c

 Copyright (C) 2011 Belledonne Communications, Grenoble, France
 Author : Johan Pascal
 
 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
19
 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
johan's avatar
johan committed
20 21 22 23 24 25 26 27 28
 */
#include "typedef.h"
#include "codecParameters.h"
#include "basicOperationsMacros.h"
#include "codebooks.h"
#include "utils.h"

#include "computeLP.h"

29 30 31 32 33 34 35
/*****************************************************************************/
/* autoCorrelation2LP: convert autocorrelation coefficients to LP using      */
/*                     Levinson-Durbin algo described in spec 3.2.2          */
/*  parameters :                                                             */
/*    -(i) autoCorrelationCoefficients : 11 values in variable scale         */
/*         scale is not needed here as a division cancel it                  */
/*    -(o) LPCoefficientsQ12: 10 LP coefficients in Q12                      */
36 37
/*    -(o) reflectionCoefficient: 10 values Q31, k generated during Levinson */
/*         Durbin LP coefficient generation and needed for VAD and RFC3389   */
38 39 40 41
/*    -(o) residualEnergy : with the same scale factor as input              */
/*         autoCorrelationCoefficients, needed by DTX                        */
/*                                                                           */
/*****************************************************************************/
42
void autoCorrelation2LP(word32_t autoCorrelationCoefficients[], word16_t LPCoefficientsQ12[], word32_t reflectionCoefficients[], word32_t *residualEnergy) {
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
	/*********************************************************************************/
	/* Compute the LP Coefficient using Levinson-Durbin algo spec 3.2.2              */
	/*********************************************************************************/
	/* start a iteration i=2, init values as after iteration i=1 :                   */
	/*        a[0] = 1                                                               */
	/*        a[1] = -r1/r0                                                          */
	/*        E = r0(1 - a[1]^2)                                                     */
	/*                                                                               */
	/*  iterations i = 2..10                                                         */
	/*       sum = r[i] + ∑ a[j]*r[i-j] with j = 1..i-1 (a[0] is always 1)           */
	/*       a[i] = -sum/E                                                           */
	/*       iterations j = 1..i-1                                                   */
	/*            a[j] += a[i]*a{i-1}[i-j] use a{i-1}: from previous iteration       */
	/*       E *=(1-a[i]^2)                                                          */
	/*                                                                               */
	/*  r in Q31 (normalised) stored in array autoCorrelationCoefficients            */
	/*  E in Q31 (can't be > 1)                                                      */
	/*  sum in Q27 (sum can't be > 1 but intermediate accumulation can)              */
	/*  a in Q4.27 with full range possible                                          */
	/*      Note: during iteration, current a[i] is in Q31 (can't be >1) and is      */
	/*            set to Q27 at the end of current iteration                         */
	/*                                                                               */
	/*********************************************************************************/
	word32_t previousIterationLPCoefficients[NB_LSP_COEFF+1]; /* to compute a[]*/
	word32_t LPCoefficients[NB_LSP_COEFF+1]; /* in Q4.27 */

	word32_t E = 0; /* in Q31 */
	word32_t sum = 0; /* in Q27 */
	int i,j;

	/* init */
	LPCoefficients[0] = ONE_IN_Q27;
	LPCoefficients[1] = -DIV32_32_Q27(autoCorrelationCoefficients[1], autoCorrelationCoefficients[0]); /* result in Q27(but<1) */
76
	reflectionCoefficients[0] = SHL(LPCoefficients[1],4); /* k[0] is -r1/r0 in Q31 */
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
	/* E = r0(1 - a[1]^2) in Q31 */
	E = MULT32_32_Q31(autoCorrelationCoefficients[0], SUB32(ONE_IN_Q31, MULT32_32_Q23(LPCoefficients[1], LPCoefficients[1]))); /* LPCoefficient[1] is in Q27, using a Q23 operation will result in a Q31 variable */
	
	for (i=2; i<NB_LSP_COEFF+1; i++) {
		/* update the previousIterationLPCoefficients needed for this one */
		for (j=1; j<i; j++) {
			previousIterationLPCoefficients[j] = LPCoefficients[j];
		}

		/* sum = r[i] + ∑ a[j]*r[i-j] with j = 1..i-1 (a[0] is always 1)           */
		sum = 0;
		for (j=1; j<i; j++) {
			sum = MAC32_32_Q31(sum, LPCoefficients[j], autoCorrelationCoefficients[i-j]);/* LPCoefficients in Q27, autoCorrelation in Q31 -> result in Q27 -> sum in Q27 */
		}
		sum = ADD32(SHL(sum, 4), autoCorrelationCoefficients[i]); /* set sum in Q31 and add r[0] */
		
		/* a[i] = -sum/E                                                           */
		LPCoefficients[i] = -DIV32_32_Q31(sum,E); /* LPCoefficient of current iteration is in Q31 for now, it will be set to Q27 at the end of this iteration */
95
		reflectionCoefficients[i-1] = LPCoefficients[i]; /* k[1] is needed by VAD others by RFC3389 RTP payload for Comfort Noise spectral information encoding */
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115

		/* iterations j = 1..i-1                                                   */
		/*      a[j] += a[i]*a[i-j]                                                */
		for (j=1; j<i; j++) {
			LPCoefficients[j] = MAC32_32_Q31(LPCoefficients[j], LPCoefficients[i], previousIterationLPCoefficients[i-j]); /*LPCoefficients in Q27 except for LPCoefficients[i] in Q31 */
		}
		/* E *=(1-a[i]^2)                                                          */
		E = MULT32_32_Q31(E, SUB32(ONE_IN_Q31, MULT32_32_Q31(LPCoefficients[i], LPCoefficients[i]))); /* all in Q31 */
		/* set LPCoefficients[i] from Q31 to Q27 */
		LPCoefficients[i] = SHR(LPCoefficients[i], 4);
	}
	*residualEnergy = E;

	/* convert with rounding the LP Coefficients form Q27 to Q12, ignore first coefficient which is always 1 */
	for (i=0; i<NB_LSP_COEFF; i++) {
		LPCoefficientsQ12[i] = (word16_t)SATURATE(PSHR(LPCoefficients[i+1], 15), MAXINT16);
	}
}


johan's avatar
johan committed
116 117 118 119 120 121
/*****************************************************************************/
/* computeLP : As described in spec 3.2.1 and 3.2.2 : Windowing,             */
/*      Autocorrelation and Levinson-Durbin algorithm                        */
/*    parameters:                                                            */
/*      -(i) signal: 240 samples in Q0, the last 40 are from next frame      */
/*      -(o) LPCoefficientsQ12: 10 LP coefficients in Q12                    */
122 123
/*      -(o) reflectionCoefficient: 10 values Q31, k generated by Levinson   */
/*         Durbin LP coefficient generation and needed for VAD and RFC3389   */
124 125 126 127 128 129 130 131 132
/*      -(o) reflectionCoefficient: in Q31, k[1] generated during Levinson   */
/*           Durbin LP coefficient generation and needed for VAD             */
/*      -(o) autoCorrelationCoefficients : used internally but needed by VAD */
/*            scale is variable                                              */
/*      -(o) noLagautoCorrelationCoefficients : needed by DTX                */
/*            scale is variable                                              */
/*      -(o) autoCorrelationCoefficientsScale : scale factor of previous buf */
/*      -(i) autoCorrelationCoefficientsNumber number of coeff to be computed*/
/*           13 if we are using them for VAD, only 11 otherwise              */
johan's avatar
johan committed
133
/*****************************************************************************/
134
void computeLP(word16_t signal[], word16_t LPCoefficientsQ12[], word32_t reflectionCoefficients[], word32_t autoCorrelationCoefficients[], word32_t noLagAutocorrelationCoefficients[], int8_t *autoCorrelationCoefficientsScale, uint8_t autoCorrelationCoefficientsNumber)
johan's avatar
johan committed
135 136
{
	int i,j;
137 138 139
	word16_t windowedSignal[L_LP_ANALYSIS_WINDOW];
	word64_t acc64=0; /* acc on 64 bits */ 
	int rightShiftToNormalise=0;
140
	word32_t residualEnergy; /* interally compute by autoCorrelation2LP and extracted for DTX only, useless here */
johan's avatar
johan committed
141 142 143 144 145 146 147 148 149 150 151

	/*********************************************************************/
	/* Compute the windowed signal according to spec 3.2.1 eq4           */
	/*********************************************************************/
	for (i=0; i<L_LP_ANALYSIS_WINDOW; i++) {
		windowedSignal[i] = MULT16_16_P15(signal[i], wlp[i]); /* signal in Q0, wlp in Q0.15, windowedSignal in Q0 */
	}

	/*********************************************************************************/
	/* Compute the autoCorrelation coefficients r[0..10] according to spec 3.2.1 eq5 */
	/*********************************************************************************/
152 153 154
	/* Compute autoCorrelationCoefficients[0] first as it is the highest number and normalise it on 32 bits then apply the same normalisation to the other coefficients */
	/* autoCorrelationCoefficients are normalised on 32 bits and then considered as Q31 in range [-1,1[ */
	/* autoCorrelationCoefficients[0] is computed on 64 bits as it is likely to overflow 32 bits */
johan's avatar
johan committed
155 156 157 158 159 160 161 162 163 164 165 166
	for (i=0; i<L_LP_ANALYSIS_WINDOW; i++) {
		acc64 = MAC64(acc64, windowedSignal[i], windowedSignal[i]);
	}
	if (acc64==0) {
		acc64 = 1; /* spec 3.2.1: To avoid arithmetic problems for low-level input signals the value of r(0) has a lower boundary of r(0) = 1.0 */
	}
	/* normalise the acc64 on 32 bits */
	if (acc64>MAXINT32) {
		do {
			acc64 = SHR(acc64,1);
			rightShiftToNormalise++;
		} while (acc64>MAXINT32);
167
		autoCorrelationCoefficients[0] = acc64;
johan's avatar
johan committed
168 169
	} else {
		rightShiftToNormalise = -countLeadingZeros((word32_t)acc64);
170
		autoCorrelationCoefficients[0] = SHL((word32_t)acc64, -rightShiftToNormalise);
johan's avatar
johan committed
171 172
	}

173 174 175 176
	/* give current autoCorrelation coefficients scale to the output */
	*autoCorrelationCoefficientsScale = -rightShiftToNormalise;

	/* compute autoCorrelationCoefficients 1 to requested number (10 - no VAD - or 12 if VAD is enabled */
johan's avatar
johan committed
177
	if (rightShiftToNormalise>0) { /* acc64 was not fitting on 32 bits so compute the other sum on 64 bits too */
178
		for (i=1; i<autoCorrelationCoefficientsNumber; i++) {
johan's avatar
johan committed
179 180 181 182 183 184
			/* compute the sum in the 64 bits acc*/
			acc64=0;
			for (j=i; j<L_LP_ANALYSIS_WINDOW; j++) {
				acc64 = ADD64_32(acc64, MULT16_16(windowedSignal[j], windowedSignal[j-i]));
			}
			/* normalise it */
185
			autoCorrelationCoefficients[i] = SHR(acc64 ,rightShiftToNormalise);
johan's avatar
johan committed
186 187
		}
	} else { /* acc64 was fitting on 32 bits, compute the other sum on 32 bits only as it is faster */
188
		for (i=1; i<autoCorrelationCoefficientsNumber; i++) {
johan's avatar
johan committed
189 190 191 192 193 194
			/* compute the sum in the 64 bits acc*/
			word32_t acc32=0;
			for (j=i; j<L_LP_ANALYSIS_WINDOW; j++) {
				acc32 = MAC16_16(acc32, windowedSignal[j], windowedSignal[j-i]);
			}
			/* normalise it */
195
			autoCorrelationCoefficients[i] = SHL(acc32, -rightShiftToNormalise);
johan's avatar
johan committed
196 197 198
		}
	}

199 200 201
	/* save autocorrelation before applying lag window as they are requested like this for DTX */
	for (i=0; i<autoCorrelationCoefficientsNumber; i++) {
		noLagAutocorrelationCoefficients[i] = autoCorrelationCoefficients[i];
johan's avatar
johan committed
202 203
	}

204
	/* apply lag window on the autocorrelation coefficients spec 3.2.1 eq7 */
johan's avatar
johan committed
205 206 207 208 209
	/* this check shall be useless but it makes some compiler happy */
	if (autoCorrelationCoefficientsNumber>NB_LSP_COEFF+3) {
		autoCorrelationCoefficientsNumber = NB_LSP_COEFF+3;
	}

210 211 212
	for (i=1; i<autoCorrelationCoefficientsNumber; i++) {
		autoCorrelationCoefficients[i] = MULT16_32_P15(wlag[i], autoCorrelationCoefficients[i]); /* wlag in Q15 */
		//autoCorrelationCoefficients[i] = MULT32_32_Q31(wlag[i], autoCorrelationCoefficients[i]); /* wlag in Q31 */
johan's avatar
johan committed
213 214
	}

215
	/* convert to LP using Levinson-Durbin algo described in sepc 3.2.2 */
216
	autoCorrelation2LP(autoCorrelationCoefficients, LPCoefficientsQ12, reflectionCoefficients, &residualEnergy);
johan's avatar
johan committed
217 218 219

	return;
}