Disclaimer Information in this document is provided in connection with Intel products. No license, express or implied, by estoppel or otherwise, to any intellectual property rights is granted by this document. Except as provided in Intel's Terms and Conditions of Sale for such products, Intel assumes no liability whatsoever, and Intel disclaims any express or implied warranty, relating to sale and/or use of Intel products including liability or warranties relating to fitness for a particular purpose, merchantability, or infringement of any patent, copyright or other intellectual property right. Intel products are not intended for use in medical, life saving, or life sustaining applications. Intel may make changes to specifications and product descriptions at any time, without notice. Copyright © Intel Corporation (1996). Thirdparty brands and names are the property of their respective owners.

The Intel Architecture (IA) media extensions include singleinstruction, multidata (SIMD) instructions. This application note presents examples of code that exploit these instructions. Specifically, the levinsonmmx function presented here illustrates how to use the new MMX technology multiplyandadd instruction (PMADDWD) to perform matrix multiplication more efficiently. The performance improvement relative to traditional IA code is due to the ability to efficiently perform multiple multiplyand add operations in fewer cycles. To perform a signed Q15 multiplyandadd would use two IA instructions, IMUL and ADD, and would take as many as 13 cycles. Using the MMX technology PMADDWD instruction, you can perform four word multiplies and two doubleword adds in three cycles. Performance gain can also be attributed to the fact that the MMX instructions operate on packed 64bit values instead of 32bit values.
The amount of data which represents a human voice or sound is usually too large to store on a typical PC. Therefore, encoding the sound and storing only a partial set of the data would be more practical. Voice encoding is one of the applications in which the LevinsonDurbin algorithm is used. This algorithm generates both a set of prediction coefficients and reflection coefficients in a recursive manner as follows:
a_{m} = K_{m} =  r(m) + r^{b}_{m1}a_{m1}  r(0) + r^{b}_{m1}a^{b}_{m1 }a_{m}(k) = a_{m1}(k) + K_{m}a_{m1}(mk) k = 1, 2, ..., m  1 m = 1, 2, ..., p where a is the prediction coefficient array, K is the reflection coefficient array, and r is the input vector.
Note: The superscript b denotes the vector with elements taken in reverse order.
The integer algorithm consists of three sections:
The input and output of this version of the LevinsonDurbin algorithm is represented in fixedpoint notation. The input matrix r[ ] and the resultant reflection coefficients k[ ] are represented as Q15 fractions stored in an array of short integers of size m. The prediction coefficients a[ ] are represented as Q13 fractions also stored in an array of short integers of size m. If two (2) signed Q15 fractions are multiplied together, the result will be a Q30 signed fraction. If a Q15 signed fraction is multiplied by a Q13 signed fraction, the result is a Q28 signed fraction. To determine the number of bits assigned to the fractional part of the result of a multiplication, add the number of bits assigned to the fractional parts of each multiplicand. To determine the number of bits assigned to the fractional part of the result of a division, subtract the number of bits assigned to the denominator from the number of bits assigned to the numerator. Refer to Figure 1 for further clarification.
There are various areas in this algorithm where errors can occur. These errors are due to noise in the input data and performing integer instead of floatingpoint arithmetic. This version of the LevinsonDurbin algorithm performs two error correction techniques to minimize the overall amount of error found in the output results. The first technique minimizes the error introduced by using integer arithmetic. This is accomplished by rounding results to the nearest digit of precision prior to converting from one precision to another (i.e. Q28 to Q13). Refer to Example 2 for a graphical representation of fixed point conversion techniques.
Converting from a Q28 to a Q13 signed fraction: 0x0A234238 Q28 signed fraction + 0x00004000 rounding factor 0x0A238238 rounded result >> 15 0x00001447 conversion result
The second technique brings a set of input data back into a usable range for the LevinsonDurbin algorithm. Input data can be determined to be not stable if any of the reflection coefficients produced are greater than one. To stabilize the input data, a scaling factor was used to scale both the prediction and reflection coefficients. This number may vary across applications and may need to be modified. For the system which was used to verify the MMX code optimized LevinsonDurbin algorithm, the appropriate scaling factor was was 0x7ff8.
There were numerous areas within the LevinsonDurbin algorithm which could be easily adapted and optimized for use with the MMX instruction set. There was also one area which could not. The remainder of this section discusses each part of the LevinsonDurbin algorithm and the optimizations and MMX technology tricks used during programming.
The formulas used to calculate R_{n} and R_{d} are as follows:
R_{n} = r(mi)*a(i) for i=0 to m1 R_{d} = r(i)*a(i) for i=0 to m1 where r() is the input matrix and a() is the prediction coefficient matrix.
These formulas are simple vector multiplications, but when using strictly IA integer instructions may take many cycles to calculate. The MMX instruction set supplies a multiplyadd instruction (PMADDWD) which performs four word multiplies and two doubleword adds in three cycles. This instruction can be pipelined with any MMX instruction, including another PMADDWD instruction. The result of a PMADDWD instruction can only be used by other instructions issued three clocks later. If it is used prior to the three clock latency, a stall will occur until the result is ready.
The inner loop which calculates R_{n} and R_{d} uses the PMADDWD instruction to perform four word multiplies and two doubleword adds for R_{n} and R_{d} with each iteration. The result of each multiply and add instruction is two doubleword Q28 signed fixed point fractions stored in MM0 (R_{d}) and MM2 (R_{n}). Each of these values must be added to an accumulator, which keeps a running summation of the current result (MM6 is used for R_{d} and MM7 is used for R_{n}).
In order to use the PMADDWD instruction within the inner loop to perform four word multiplies and two doubleword adds, two assumptions were made. The first is that the prediction coefficient array is initialized to zero. This prevents incorrect results for calculations of R_{n} and R_{d} which include three or fewer input values. The second is that the size (in bytes) of both the input vector and the prediction coefficient vector must be divisible by eight. This is to prevent reading past the end of either vector since four input values are read for each iteration of the loop.
Example 4 contains the code for the inner loop of the LevinsonDurbin algorithm.
calcRnRd: 1. movq mm2, [edx][eax*2] ; get the next 4 rT() matrix entries 2. pmaddwd mm0, mm1 ; calculate (r(i) * a(i)) + (r(i+1) * ; a(i+1)) and (r(i+2) * a(i+2)) + ; (r(i+3) * a(i+3)) add eax, 3. ; increment the aMatrix address 4. pmaddwd mm2, mm1 ; calculate (rT(i) * a(i)) + (rT(i+1) * ; a(i+1)) and (rT(i+2) * a(i+2)) + ; (rT(i+3) * a(i+3)) 5. ; Two penalties occur during this clock 6. ; cycle. The first is a pmaddwd stall ; since the instruction on line 8 uses ; the result of the pmaddwd instruction ; issued on line 2. The second penalty ; occurs because eax is incremented on ; line 3 and used the next cycle. 7. movq mm1, [ecx][eax*2] ; get the next 4 a() matrix entries 8. paddd mm6, mm0 ; accumulate Rd 9. movq mm0, [esi][eax*2] ; get the next 4 r() matrix entries 10. paddd mm7, mm2 ; accumulate Rn 11. cmp eax, ebx ; compare i and m 12. jl calcRnRd ; if (i < m) then keep incrementally ; calculating Rn and Rd
Assuming all data is present in the data cache, the inner loop which calculates R_{n} and R_{d} varies in the number of cycles it takes to execute. If the data being read from the reverse input matrix is aligned, no penalty is incurred and the loop executes in six cycles per pass. If the data being read is unaligned, the read incurs a three cycle penalty, increasing the overall loop time from six to nine cycles per pass. Data alignment is extremely critical to the number of cycles in which the loop executes. If the input matrix and the prediction coefficient matrix are not aligned on quad word boundaries, accesses to these matrices will also incur an unaligned penalty of three clock cycles.
Within the inner loop is a one clock cycle penalty which is caused by two violations of the pairing rules. The first penalty occurs on line 7 where a quad word move is issued. This instruction uses the EAX register to generate an index into the a matrix. Exactly one cycle prior to issuing the MOVQ instruction, EAX was used as a destination register on line 3. Pairing rules state that any integer register can not be used to address memory for two clocks after it has been modified. This is referred to as an address generation interlock (AGI). Therefore a penalty of one clock cycle is incurred. The second penalty occurs on line 8 where a packed add is issued. This instruction uses MM0 as an operand. MM0 contains the result of the PMADDWD instruction issued on line 2. Pairing rules state that the result of a multiply instruction can only be used by other instructions issued three clocks later. Therefore, a penalty of one clock cycle is incurred. Since both penalties occur in the same cycle, only a one clock cycle stall occurs instead of two.
This loop was optimized to support efficient calculations of both the prediction and reflection coefficients to an order of 10. For orders of 16 or greater, this loop can be optimized even further by unrolling the loop one more time. This would be done by performing eight word multiplies and four doubleword adds for both R_{n} and R_{d} each time through the loop. It would eliminate half of the loop overhead and may also have the potential of hiding the multiplyadd and AGI pipe stall at instruction 8. See Example 5 for an example.
calcRnRd: 1. movq mm2, [edx][eax*2] ; get the next 4 rT() matrix entries 2. pmaddwd mm0, mm1 ; calculate (r(i) * a(i)) + (r(i+1) * ; a(i+1)) and (r(i+2) * a(i+2)) + ; (r(i+3) * a(i+3)) 3. movq mm3, 8[edx][eax*2] ; get the next 4 rT() matrix entries 4. pmaddwd mm2, mm1 ; calculate (rT(i) * a(i)) + (rT(i+1) * ; a(i+1)) and (rT(i+2) * a(i+2)) + ; (rT(i+3) * a(i+3)) 7. movq mm1, 8[ecx][eax*2] ; get the next 4 a() matrix entries 8. paddd mm6, mm0 ; accumulate Rd 9. movq mm4, 8[esi][eax*2] ; get the next 4 r() matrix entries 10. pmaddwd mm3, mm1 ; calculate (r(i+4) * a(i+4)) + (r(i+5)* ; a(i+5)) + (r(i+6) * a(i+6)) + (r(i+7)* ; a(i+7)) 11. movq mm0, [esi][eax*2] ; get the next 4 r() matrix elements 12. pmaddwd mm4, mm1 ; calculate (rT(i+4)*a(i+4)) + (rT(i+5) ; * a(i+5)) + (rT(i+6) * a(i+6)) + ; (rT(i+7) + a(i+7)) 13. movq mm1, [ecx][eax*2] ; get the next 4 a() matrix elements 14. paddd mm7, mm2 ; accumulate Rn 15. add eax, 8 ; add 8 elements to the eax index 16. paddd mm6, mm3 ; accumulate Rd 17. paddd mm7, mm4 ; accumulate Rn 18. cmp eax, ebx ; compare i and m 19. jl calcRnRd ; if (i < m) then keep incrementally ; calculating Rn and Rd
The formula to compute R_{n}_{ }uses the values of the input vector in the reverse order. Three options were considered when deciding how to obtain this vector:
Option 3 was chosen since this read/store operation could be paired with other nonpaired instructions prior to executing the inner loop. Therefore, the reverse of the input matrix was obtained "free." Both options 1 and 2 were investigated and discarded because of the extra cycles added to the overall performance of the code. Option 1 added extra instructions and loop overhead to the function. Option 2 added extra cycles within the inner loop to rearrange data appropriately for the calculation of R_{n}.
After the execution of the inner loop, R_{n} and R_{d} have their results accumulated in MM7 and MM6, respectively. These registers each have two doubleword Q28 signed fractions which need to be separated and added together to obtain the final result. A copy of each accumulator is made and then the packedshiftrightlogicalquad (PSRLQ) instruction is used to shift the high doubleword into the low doubleword position. This is added to the accumulator to obtain the final R_{n} and R_{d} result in the low doubleword position. Refer to Example 6 for a code example. In this example, there are also statements which get the data ready for the division in the next code segment. This includes negating R_{n} and clearing registers MM3 and MM0.
calcnew_A: 1. movq mm0, mm6 ; make a copy of the partial Rd 2. movq mm1, mm7 ; make a copy of the partial Rn 3. movd mm4, round_factor ; load the rounding factor for later 4. psrlq mm1, 32 ; put the high doubleword of Rn into ; the low doubleword position 5. paddd mm7, mm1 ; do the final Rn calculation 6. pxor mm3, mm3 ; zero out mm3 to negate Rn for the ; divide 7. psubd mm3, mm7 ; negate Rn prior to the divide 8. psrlq mm0, 32 ; put the high doubleword of Rd into ; the low doubleword position 9. paddd mm6, mm0 ; do the final Rd calculation 10. pxor mm1, mm1 ; zero out mm1 for unpacking the ; rounding factor later
The m^{th} order reflection and prediction coefficients are calculated by dividing R_{n} by R_{d}. This integer division can not be avoided since the denominator changes for each iteration of the main loop. It takes 43 clock cycles every time that it is executed. If the denominator value remained constant with each iteration of the outer loop, one divide could have been executed prior to entering the main loop to determine 1/R_{d}. This value could then be used in an integer multiplication using the IMUL instruction and would take 10 cycles to execute. The IMUL instruction would have to be used in this instance since the MMX instruction set does not supply a doubleword multiply instruction. A floating point divide could also be executed in place of the fixed point integer divide, but would incur a 100 clock cycle penalty of switching from floating point code to MMX code. Therefore it was not used.
This segment of code recalculates the prediction coefficients from 0 through m1 based on the newly calculated m^{th} coefficient. The pseudocode is as follows:
for ( i = 1; i < m; i++) b[i] = a[i] + (K[m] * a[mi]); for ( i = 1; i < m; i++) a[i] = b[i];
To avoid needing the temporary array b, the MMX code calculates both the new a[i] and a[mi] with each iteration of the loop. The MMX instruction set does not implicitly provide a doubleword multiply. Therefore, the PMADDWD instruction was used. Refer to Figure 2 for a graphical representation of the data setup. After the data has been setup, the instruction PMADDWD MM0, MM2 is issued to calculate K[m]*a[i] and K[m]*a[mi]. The result of this multiplyadd is a Q28 signed fraction. It is then added to MM3, which contains a[i] in the low doubleword and a[mi] in the high doubleword, each represented as Q28 signed fractions. This addition obtains the new a[i] and a[mi] values.
An approximate overall performance of this version of the LevinsonDurbin algorithm can be calculated based on the performance of three sections of the code: the calculation of R_{n} and R_{d}, the calculation of K[m], and the recalculation of the prediction coefficients. The following discussion about performance assumes all data is present in the data cache and that all data is properly aligned.
There are two cycle count numbers which are important to take into consideration when approximating the number of cycles it takes to calculate R_{n} and R_{d}. When accesses to the reverse of the input vector are aligned, it takes six clock cycles per iteration to multiply and accumulate four values for both R_{n} and R_{d}. When accesses to the reverse of the input vector are unaligned, it takes nine clock cycles per iteration. Since half of all accesses to the reverse of the input vector are aligned and half are unaligned, on average it takes 7.5 cycles to calculate four values for both R_{n} and R_{d} or .9375 clock cycles for one value.
To calculate a[m] and K[m] and to set up register values prior to recalculating the prediction coefficients takes approximately 64 clock cycles. The majority of these cycles can be directly attributed to the integer divide which executes in 43 clock cycles. As mentioned previously, this integer divide can not be avoided.
For the recalculation of the prediction coefficients, there are two cycle count numbers which must be taken into consideration. When accesses to the prediction coefficient array are aligned, it takes 14 clock cycles to recalculate two elements. When accesses to the prediction coefficient array are unaligned, it takes 20 clock cycles to recalculate two elements. Since half of all accesses to the prediction coefficient array are aligned and half are unaligned, on average it takes 17 cycles to recalculate two elements or 8.5 clock cycles to recalculate one.
Applying these cycle count numbers to an application which calculates pSize reflection coefficients using an input vector of size pSize+1, an approximate overall performance can be obtained as follows:
For the C language version of the LevinsonDurbin algorithm, refer to Appendix A. For the complete MMX code optimized version of the LevinsonDurbin algorithm, refer to Appendix B.
The MMX code optimized version is reentrant but the user may wish to dispose of the local variables and pass them in as arguments to the procedure if stack usage is an issue.
/ ********************************************************************************* * Description: * LevinsonDurbin is the scalar (not MMX code) version of the LevinsonDurbin * algorithm. It is used to calculate the reflection and prediction coefficients of * a given set of normal equations. * * Inputs: * r short int * a pointer to the first element * of the input 'r' matrix * a short int * a pointer to the output * prediction coefficients array * k short int * a pointer to the output * reflection coefficients array * p short int the number of reflection coefficients * to solve for (typically 10 or 16). * *********************************************************************************/ void Levinson_Durbin (short *r, short a*, short *k, short p) { // begin Levinson_Durbin() short i, // inner loop index m, // outer loop index b[11]; // temporary vector to store the prediction coefficients long Rn, // inner loop numerator accumulator Rd, // inner loop denominator accumulator temp; // temporary variable used for intermediate result calculations // initialize a[0] to 1/4 a[0] = 8192; // For each order, calculate a new prediction and reflection coefficient for (m = 1; m < p + 1; m++) { // Initialize the numerator and denominator accumulators to zero Rn = Rd = 0; // Calculate the numerator and denominator values for the integer division for (i = 0; i < m; i++) { Rn = Rn + (r[m1] * (long)a[i]); Rd = Rd + (r[i] * (long)a[i]); } // Calculate the reflection coefficient k[m]. Round the Q28 number prior // to converting it to a Q15 number. Also, scale it by the scaling factor // to help keep the input data in proper range. k[m] = Rn / ((Rd + 0x4000) >> 15); k[m] = (((long)k[m] * 0x7ff8) + 0x4000) >> 15; // Calculate the new prediction coefficient by converting k[m] from // a Q15 to a Q13 number b[m] = (k[m] + 0x2) >> 2; // Calculate the new prediction coefficients for the next iteration for (i = 1; i < m; i++) b[m] = (((long)a[i] << 15) + (k[m] * (long)a[mi]) + 0x4000) >> 15); // Copy the prediction coefficients from the temporary b[] array to a[] for ( i = 1; i < m+1; i++) a[i] = b[i]; } } // end Levinson_Durbin()
;***************************************************************************/ ;* Description: ;* The purpose of this file is to provide the MMX code for the ;* levinsondurbin algorithm as an instructional example to those who ;* are just beginning to code using MMX instructions. ;* ;* Assumptions: ;* 1. The set of normal equations given are stable ;* 2. The output matrix 'a' containing the prediction coefficients ;* is initialized to all 0's. ;* 3. The normal equations are represented by one matrix 'r' which ;* contains the coefficients gamma(0) through gamma(p) given as ;* Q15 signed short integers. ;* 4. The resultant reflection coefficients will be returned as ;* Q15 signed short integers. ;* 5. The resultant prediction coefficients will be returned as ;* Q13 signed short integers. ;* 6. The rTmatrix defined below is of size (p*2) + 16 bytes. ;****************************************************************************/ ;* TEXT SEGMENT ;**************************************************************************** _TEXT SEGMENT ; ; Declare levinson_mmx as a public routine to allow the 'C' code to ; call it. ; PUBLIC levinson_mmx ;***************************************************************************/ ;* _levinson_mmx is the MMX optimized code version of the LevinsonDurbin ;* algorithm. It is used to calculate the reflection coefficients of ;* a given set of normal equations. ;* ;* Inputs: ;* rPtr short int * a pointer to the first element ;* of the input 'r' matrix ;* aPtr short int * a pointer to the output ;* prediction coefficients array ;* kPtr short int * a pointer to the output ;* reflection coefficients array ;* pSize short int the number of reflection coefficients ;* to solve for (typically 10 or 16). ;* ;****************************************************************************/ levinson_mmx PROC C USES ebx ecx edx esi, rMatrix:PTR WORD, aMatrix:PTR WORD, kMatrix:PTR WORD, pSize:DWORD ; ; Declare local variables ; 1. rTmatrix  used to store the reverse of rMatrix ; 2. rTindex  used to store the rTmatrix index ; 3. mSave  used to store the loop counter 'm' ; 4. round_factor  contains the rounding factor 0x4000 ; 5. scale_factor  contains the scaling factor 0x7ff8 ; LOCAL rTmatrix[72]:WORD LOCAL rTindex:DWORD LOCAL mSave:DWORD LOCAL round_factor:DWORD LOCAL scale_factor:DWORD ; ; This section of code initializes the local variables round_factor and scale_factor ; and sets up some of the items which need to be done prior to executing the algorithm ; mov round_factor, 4000H ; initialize the rounding factor mov eax, pSize ; get the value of 'p' mov scale_factor, 7ff8H ; initialize the scaling factor mov esi, aMatrix ; get the pointer to the 'a' matrix mov WORD PTR [esi], 2000H ; store the scaling element into a[0] lea edx, rTmatrix[eax*2] ; get the first rTmatrix index address mov rTindex, edx ; save the first rTmatrix address mov ebx, 1 ; initialize the main loop counter 'm' ; ; This section of code sets up the pointer to ; the input array 'r' and rTmatrix used to calculate Rn and Rd. It also ; performs the test: ; while (p >= m) ; _levinson_outer_loop: mov esi, rMatrix ; initialize the rMatrix pointer cmp eax, ebx ; is (p >= m) ? mov ecx, [esi][ebx*2] ; get the next rMatrix element jl _levinson_done ; if (p < m) we are done! ; ; This section of code sets up registers used to calculate Rn and Rd. It also ; incrementally builds the rTmatrix with each iteration of the main levinson loop, ; adding one more rMatrix element to the front of the array. The reverse of the ; input matrix 'r' is needed when calculating Rn. Rn uses the elements 1 to m of ; rMatrix in reverse as one of the multiplicands with aMatrix elements as the other. ; _setup_RnRd_calc: mov eax, 0 ; initialize the loop counter to 0 mov [edx], cx ; store the next rTmatrix element mov ecx, aMatrix ; get the pointer to the 'a' matrix movq mm0, [esi] ; get the first quad word of 'r' pxor mm6, mm6 ; initialize Rd to 0 movq mm1, [ecx] ; get the first quad word of aMatrix pxor mm7, mm7 ; initialize Rn to 0 cmp eax, ebx ; is (i < m) ? jge _calc_new_A ; if not, calculate Rn and Rd ; ; This section of code calculates Rn and Rd. The calculation used is: ; while (i < m) { ; Rn = Rn + (rTmatrix[i]*aMatrix[i]); ; Rd = Rd + (rMatrix[i+1]*aMatrix[i]); ; } ; where Rn and Rd are long, Q28 signed fractions, rTmatrix, ; and rMatrix are short, Q15 signed fractions and aMatrix contains ; short, Q13 signed fractions. Up to 4 multiply/adds are done at a time ; for both Rn and Rd. ; _calc_RnRd: movq mm2, [edx][eax*2] ; get the next quad word of rTmatrix pmaddwd mm0, mm1 ; calculate 'r' * aMatrix quad multiply/add add eax, 4 ; increment the loop counter pmaddwd mm2, mm1 ; calculate rTmatrix * aMatrix quad mul/add movq mm1, [ecx][eax*2] ; get the next aMatrix quad word paddd mm6, mm0 ; Rd = Rd + (rMatrix * aMatrix) quad word movq mm0, [esi][eax*2] ; get the next rMatrix quad word paddd mm7, mm2 ; Rn = Rn + (rTmatrix * aMatrix) quad word cmp eax, ebx ; is (i < m)? jl _calc_RnRd ; if so, then keep incrementally calculating ; Rn and Rd ; ; Upon leaving the loop that calculates Rn and Rd, Rd (mm6) and Rn (mm7) contain two long, ; signed Q28 fractions which need to be added together to get the final Rn ; and Rd result. This code performs this calculation. ; _calc_new_A: movq mm0, mm6 ; make a copy of the partial Rd result movq mm1, mm7 ; make a copy of the partial Rn result movd mm4, round_factor ; load up the rounding factor psrlq mm1, 32 ; shift the Rn copy over by 32 for the add paddd mm7, mm1 ; do the last Rn calculation pxor mm3, mm3 ; clear this MMX register psubd mm3, mm7 ; negate Rn prior to doing the divide psrlq mm0, 32 ; shift the Rd copy over by 32 for the add ; ; This section of code performs the calculation of a[m] (the prediction coefficient) and k[m] ; (the reflection coefficient). The calculation is as follows: ; K[m] = Rn / ((Rd + 0x4000) >> 15); ; This step produces a Q15 signed result ; To obtain a short, Q15 signed reflection coefficient we must divide the ; Q28 signed numerator by a Q13 signed denominator. This ; will provide the desired result with the correct precision. The integer division ; can not be avoided since for every 'm' the denominator will be different. ; paddd mm6, mm0 ; do the last Rd calculation pxor mm1, mm1 ; clear mm1 movd eax, mm3 ; move the numerator to an integer register paddd mm6, mm4 ; round the denomintor prior to the shifting movd mm3, scale_factor ; setup mm3 with a scaling factor 0x7ff8 psrad mm6, 15 ; perform the right shift of the denominator mov esi, aMatrix ; get the aMatrix pointer mov edx, eax ; make edx the sign extension of eax sar edx, 31 ; complete the sign extension for the idiv mov mSave, ebx ; save the value of the main loop counter 'm' movd ecx, mm6 ; move the denominator to an integer register punpcklwd mm4, mm1 ; prepare mm4 to be used as the rounding ; register for recalculating the prediction ; coefficients idiv ecx ; perform the 32bit integer division ; ; This section of code scales the new reflection/prediction coefficient, obtains the correct ; precision for the prediction coefficient (Q13), and stores the results into the 'a' ; matrix and 'K' matrix. The calculations it performs are as follows: ; K[m] = ((K[m] * 0x7ff8) + 0x4000) >> 15; ; This step scales the reflection/prediction coefficient by .999x ; a[m] = (K[m] + 0x2) >> 2; ; This step produces a Q13 signed prediction coefficient. ; It also performs some setup for recalculating the prediction coefficients. ; movd mm2, eax ; move the idiv result to an MMX register pxor mm5, mm5 ; clear out mm5 sub ebx, 1 ; dec 'm' since the loop bounds are 0 and m1 pmaddwd mm2, mm3 ; scale K[m] by 0x7ff8 movd mm0, 2[esi] ; get a[1] for loop punpckldq mm4, mm4 ; complete the rounding register setup movd mm1, [esi][ebx*2] ; get a[m1] for loop preparation movq mm3, mm0 ; make a copy of a[1] through a[4] mov edx, kMatrix ; get the kMatrix pointer paddd mm2, mm4 ; round K[m] prior to shifting it mov ecx, 1 ; initialize the loop counter psrad mm2, 15 ; shift K[m] to obtain a Q15 fraction movd eax, mm2 ; move K[m] to an integer register psllq mm0, 32 ; move a[1] and a[2] to the upper doubleword mov 2[edx][ebx*2], ax ; store the idiv result into K[m] add eax, 2 ; round a[m] prior to shifting sar eax, 2 ; shift a[m] to obtain a Q13 fraction punpcklwd mm2, mm5 ; set an MMX register to 0, a[m], 0, a[m] as ; the a[m] multiplicand used below mov 2[esi][ebx*2], ax ; store the idiv result into a[m] punpckldq mm2, mm2 ; finish setting up the 0, a[m], 0, a[m] ; register cmp ecx, ebx ; is (i >= m) ? jge _done_alter_A ; if so, then jump past this loop ; ; This section of code calculates the new prediction coefficients to be used for the ; calculation of the (m+1)'th order reflection coefficient. The following algorithm is used: ; for (i = 0, j = m1; i < j; i++, j) ; { ; tempA = K[m] * a[j]; ; a[j] = (((long)a[j] << 15) + (K[m] * a[i]) + 0x4000) >> 15; ; a[i] = (((long)a[i] << 15) + tempA + 0x4000) >> 15; ; } ; If, upon exiting this loop, i == j, then we have one more element to do. This code can be ; found after the end of this loop. Two calculations are done for each loop iteration: ; a[i], a[j] ; _alter_A: por mm0, mm1 ; set mm0 to X, a[i], X, a[j] punpckldq mm3, mm1 ; set mm3 to X, a[j], X, a[i] pmaddwd mm0, mm2 ; multiply a[i] and a[j] by K[m] pslld mm3, 16 ; convert from Q13 to Q28 signed fraction sub ebx, 1 ; decrement j psrad mm3, 1 ; finish the conversion movd mm5, 2[esi][ecx*2] ; get the next a[i] element paddd mm3, mm4 ; add the rounding factor to a[i] and a[j] movd mm1, [esi][ebx*2] ; get the next a[j] element paddd mm0, mm3 ; add a[i] to a[j]*K[m] and a[j] to a[i]*K[m] movq mm3, mm5 ; make a copy of the next a[i] element psrad mm0, 15 ; convert the new a[i] and a[j] to Q13 ; signed fraction movd edx, mm0 ; move a[i] to an integer register psrlq mm0, 32 ; move a[j] to the low doubleword mov [esi][ecx*2], dx ; store a[i] psllq mm5, 32 ; move the next a[i] to the high doubleword movd edx, mm0 ; move a[j] to an integer register movq mm0, mm5 ; move the next a[i] to mm0 mov 2[esi][ebx*2], dx ; store a[j] inc ecx ; increment i cmp ecx, ebx ; is (i < j)? jl _alter_A ; if so, then keep looping! ; ; If, upon exiting the above loop, i == j, then perform the last calculation on a[i] by doing: ; a[i] = (((long)a[i] << 15) + (K[m]*a[i]) + 0x4000) >> 15; ; _done_alter_A: pmaddwd mm1, mm2 ; perform a[i]*K[m] pslld mm3, 16 ; convert from Q13 to Q28 signed fraction mov eax, rTindex ; get the rTmatrix current address psrad mm3, 1 ; complete the conversion mov ebx, mSave ; get the value of 'm', the main loop counter paddd mm1, mm4 ; perform a[i] + a[i]*K[m] paddd mm1, mm3 ; add the rounding factor before shifting psrad mm1, 15 ; convert the new a[i] to Q13 signed fraction jne _main_loop_done ; if (i != j) then skip this part movd edx, mm1 ; move the result to an integer register mov [esi][ecx*2], dx ; store a[i] ; ; We have now completed finding the (m+1)'th order prediction coefficients and the mth ; reflection coefficient. This section of code sets up registers to check if we need to ; continue with the next 'm'. ; _main_loop_done: add ebx, 1 ; increment 'm', the main loop counter sub eax, 2 ; since we need to add an rMatrix element on ; to the beginning of our rTmatrix, subtract ; 2 from the address mov rTindex, eax ; store this new address back to memory mov edx, eax ; setup the rTmatrix array pointer mov eax, pSize ; get the value of 'p' jmp _levinson_outer_loop ; go to the top of the main loop _levinson_done: ret 0 ; we are done! levinson_mmx ENDP _TEXT ENDS END