path: root/media/libstagefright/codecs/aacdec/fft_rx4_short.cpp

/* ------------------------------------------------------------------
 * Copyright (C) 1998-2009 PacketVideo
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
 * express or implied.
 * See the License for the specific language governing permissions
 * and limitations under the License.
 * -------------------------------------------------------------------
 */
/*

 Pathname: ./src/fft_rx4_short.c
 Funtions: fft_rx4_short

------------------------------------------------------------------------------
 REVISION HISTORY

 Description:
            (1) Eliminated search for max in the main loop.
            (2) Simplified the function by eliminating different conditions
                for exp.
            (3) Reduced precision on w_64rx4 from Q15 to Q12, so now the
                input can be as high as 1.0 and saturation will not occurre
                because the accumulation times the new Q12 format will never
                exceed 31 bits.

 Description:
            (1) Added comment to explain max search elimination and
                Q format during multiplications
            (2) Increased down shift from 1 to 2, to ensure that 32-bit
                numbers will not overflow when 2 consecutive adds are done
                This was found during code review.

 Who:                                   Date:
 Description:

------------------------------------------------------------------------------
 INPUT AND OUTPUT DEFINITIONS

 Inputs:
    Data       =  Input complex vector, arranged in the following order:
                  real, imag, real, imag...
                  This is a complex vector whose elements (real and Imag) are
                  Int32.
                  type Int32 *

    peak_value =  Input,  peak value of the input vector
                  Output,  peak value of the resulting vector
                  type Int32 *

 Local Stores/Buffers/Pointers Needed:
    None

 Global Stores/Buffers/Pointers Needed:
    None

 Outputs:
    exponent returns a shift to compensate the scaling introduced by
    overflow protection

 Pointers and Buffers Modified:
    calculation are done in-place and returned in Data

 Local Stores Modified:
    None

 Global Stores Modified:
    None

------------------------------------------------------------------------------
 FUNCTION DESCRIPTION

    Fast Fourier Transform, radix 4 with Decimation in Frequency and block
    floating point arithmetic.
    The radix-4 FFT  simply divides the FFT into four smaller FFTs. Each of
    the smaller FFTs is then further divided into smaller ones and so on.
    It consists of log 4 N stages and each stage consists of N/4 dragonflies.

    An FFT is nothing but a bundle of multiplications and summations which
    may overflow during calculations.


    This routine uses a scheme to test and scale the result output from
    each FFT stage in order to fix the accumulation overflow.

    The Input Data should be in Q13 format to get the highest precision.
    At the end of each dragonfly calculation, a test for possible bit growth
    is made, if bit growth is possible the Data is scale down back to Q13.


------------------------------------------------------------------------------
 REQUIREMENTS

    This function should provide a fixed point FFT for an input array
    of size 64.

------------------------------------------------------------------------------
 REFERENCES

    [1] Advance Digital Signal Processing, J. Proakis, C. Rader, F. Ling,
        C. Nikias, Macmillan Pub. Co.

------------------------------------------------------------------------------
 PSEUDO-CODE


   MODIFY( x[] )
   RETURN( exponent )

------------------------------------------------------------------------------
 RESOURCES USED
   When the code is written for a specific target processor the
     the resources used should be documented below.

 STACK USAGE: [stack count for this module] + [variable to represent
          stack usage for each subroutine called]

     where: [stack usage variable] = stack usage for [subroutine
         name] (see [filename].ext)

 DATA MEMORY USED: x words

 PROGRAM MEMORY USED: x words

 CLOCK CYCLES: [cycle count equation for this module] + [variable
           used to represent cycle count for each subroutine
           called]

     where: [cycle count variable] = cycle count for [subroutine
        name] (see [filename].ext)

------------------------------------------------------------------------------
*/
/*----------------------------------------------------------------------------
; INCLUDES
----------------------------------------------------------------------------*/

#include "pv_audio_type_defs.h"
#include "fft_rx4.h"
#include "pv_normalize.h"
#include "fxp_mul32.h"

/*----------------------------------------------------------------------------
; MACROS
; Define module specific macros here
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
; DEFINES
; Include all pre-processor statements here. Include conditional
; compile variables also.
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
; LOCAL FUNCTION DEFINITIONS
; Function Prototype declaration
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
; LOCAL VARIABLE DEFINITIONS
; Variable declaration - defined here and used outside this module
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
; EXTERNAL FUNCTION REFERENCES
; Declare functions defined elsewhere and referenced in this module
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
; EXTERNAL VARIABLES REFERENCES
; Declare variables used in this module but defined elsewhere
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
; EXTERNAL GLOBAL STORE/BUFFER/POINTER REFERENCES
; Declare variables used in this module but defined elsewhere
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
; FUNCTION CODE
----------------------------------------------------------------------------*/

Int fft_rx4_short(
    Int32      Data[],
    Int32      *peak_value)

{
    Int     n1;
    Int     n2;
    Int     n3;
    Int     j;
    Int     k;
    Int     i;
    Int32   exp_jw1;
    Int32   exp_jw2;
    Int32   exp_jw3;


    Int32   t1;
    Int32   t2;
    Int32   r1;
    Int32   r2;
    Int32   r3;
    Int32   s1;
    Int32   s2;
    Int32   s3;

    Int32   *pData1;
    Int32   *pData2;
    Int32   *pData3;
    Int32   *pData4;
    const Int32  *pw;
    Int32   temp1;
    Int32   temp2;
    Int32   temp3;
    Int32   temp4;
    Int32   max;
    Int     exp;
    Int     exponent = 0;
    Int     shift;


    max = *peak_value;
    exp = 0;

    if (max > 0x008000)
    {
        exp = 8 - pv_normalize(max);   /* use 24 bits  */

        exponent = exp;        /* keeps track of # of shifts */

    }

    n2 = FFT_RX4_SHORT;

    pw = W_64rx4;


    /* shift down to avoid possible overflow in first pass of the loop */
    shift = 2;

    for (k = FFT_RX4_SHORT; k > 4; k >>= 2)
    {

        n1 = n2;
        n2 >>= 2;
        n3 = n1 >> 1;

        exp -= 2;

        for (i = 0; i < FFT_RX4_SHORT; i += n1)
        {
            pData1 = &Data[ i<<1];
            pData3 = pData1 + n3;
            pData2 = pData1 + n1;
            pData4 = pData3 + n1;

            temp1   = *(pData1);
            temp2   = *(pData2);
            temp1   >>= shift;
            temp2   >>= shift;

            r1      = temp1 + temp2;
            r2      = temp1 - temp2;

            temp3   = *(pData3++);
            temp4   = *(pData4++);
            temp3   >>= shift;
            temp4   >>= shift;

            t1      = temp3 + temp4;
            t2      = temp3 - temp4;

            *(pData1++) = (r1 + t1) >> exp;
            *(pData2++) = (r1 - t1) >> exp;

            temp1   = *pData1;
            temp2   = *pData2;
            temp1   >>= shift;
            temp2   >>= shift;

            s1      = temp1 + temp2;
            s2      = temp1 - temp2;

            temp3   = *pData3;
            temp4   = *pData4;
            temp3   >>= shift;
            temp4   >>= shift;

            t1      = temp3 + temp4;
            r1      = temp3 - temp4;

            *pData1   = (s1 + t1) >> exp;
            *pData2   = (s1 - t1) >> exp;

            *pData4--    = (s2 + t2) >> exp;
            *pData4      = (r2 - r1) >> exp;

            *pData3--    = (s2 - t2) >> exp;
            *pData3      = (r2 + r1) >> exp;


        }  /* i */

        for (j = 1; j < n2; j++)
        {
            exp_jw1 = *pw++;
            exp_jw2 = *pw++;
            exp_jw3 = *pw++;


            for (i = j; i < FFT_RX4_SHORT; i += n1)
            {
                pData1 = &Data[ i<<1];
                pData3 = pData1 + n3;
                pData2 = pData1 + n1;
                pData4 = pData3 + n1;

                temp1   = *(pData1);
                temp2   = *(pData2++);
                temp1   >>= shift;
                temp2   >>= shift;

                r1      = temp1 + temp2;
                r2      = temp1 - temp2;
                temp3   = *(pData3++);
                temp4   = *(pData4++);
                temp3   >>= shift;
                temp4   >>= shift;

                t1      = temp3 + temp4;
                t2      = temp3 - temp4;

                *(pData1++) = (r1 + t1) >> exp;
                r1          = (r1 - t1) >> exp;

                temp1   = *pData1;
                temp2   = *pData2;
                temp1   >>= shift;
                temp2   >>= shift;

                s1      = temp1 + temp2;
                s2      = temp1 - temp2;

                s3      = (s2 + t2) >> exp;
                s2      = (s2 - t2) >> exp;

                temp3   = *pData3;
                temp4   = *pData4 ;
                temp3   >>= shift;
                temp4   >>= shift;

                t1      = temp3 + temp4;
                t2      = temp3 - temp4;

                *pData1  = (s1 + t1) >> exp;
                s1       = (s1 - t1) >> exp;


                *pData2--  = cmplx_mul32_by_16(s1, -r1, exp_jw2) << 1;
                *pData2    = cmplx_mul32_by_16(r1,  s1, exp_jw2) << 1;

                r3       = ((r2 - t2) >> exp);
                r2       = ((r2 + t2) >> exp);

                *pData3--  = cmplx_mul32_by_16(s2, -r2, exp_jw1) << 1;
                *pData3    = cmplx_mul32_by_16(r2,  s2, exp_jw1) << 1;

                *pData4--  = cmplx_mul32_by_16(s3, -r3, exp_jw3) << 1;
                *pData4    = cmplx_mul32_by_16(r3,  s3, exp_jw3) << 1;

            }  /* i */

        }  /*  j */

        /*
         *  this will reset exp and shift to zero for the second pass of the
         *  loop
         */
        exp   = 2;
        shift = 0;

    } /* k */


    max = 0;

    pData1 = Data - 7;

    for (i = ONE_FOURTH_FFT_RX4_SHORT; i != 0 ; i--)
    {
        pData1 += 7;

        pData3 = pData1 + 2;
        pData2 = pData1 + 4;
        pData4 = pData1 + 6;

        temp1   = *pData1;
        temp2   = *pData2++;

        r1      = temp1 + temp2;
        r2      = temp1 - temp2;

        temp1   = *pData3++;
        temp2   = *pData4++;

        t1      = temp1 + temp2;
        t2      = temp1 - temp2;

        temp1       = (r1 + t1);
        r1          = (r1 - t1);
        *(pData1++) = temp1;
        max        |= (temp1 >> 31) ^ temp1;



        temp1   = *pData1;
        temp2   = *pData2;

        s1      = temp1 + temp2;
        s2      = temp1 - temp2;

        s3      = (s2 + t2);
        s2      = (s2 - t2);

        temp1   = *pData3;
        temp2   = *pData4;

        t1      = temp1 + temp2;
        t2      = temp1 - temp2;

        temp1      = (s1 + t1);
        temp2      = (s1 - t1);
        *pData1    = temp1;
        *pData2--  = temp2;
        max       |= (temp1 >> 31) ^ temp1;
        max       |= (temp2 >> 31) ^ temp2;

        *pData2    = r1;
        *pData3--  = s2;
        *pData4--  = s3;
        max       |= (r1 >> 31) ^ r1;
        max       |= (s2 >> 31) ^ s2;
        max       |= (s3 >> 31) ^ s3;

        temp1      = (r2 - t2);
        temp2      = (r2 + t2);
        *pData4    = temp1;
        *pData3    = temp2;
        max       |= (temp1 >> 31) ^ temp1;
        max       |= (temp2 >> 31) ^ temp2;

    }  /* i */

    *peak_value = max;


    return (exponent);

}