arm_dct4_f32.c

Go to the documentation of this file.

/* ----------------------------------------------------------------------      
* Copyright (C) 2011 ARM Limited. All rights reserved. 
*      
* $Date:        15. December 2011   
* $Revision:    V2.0.0  
*      
* Project:      Cortex-R DSP Library 
* Title:        arm_dct4_f32.c      
*      
* Description:  Processing function of DCT4 & IDCT4 F32.      
*      
* Target Processor:          Cortex-R4/R5
*
* Version 1.0.0 2011/03/08
*     Alpha release.
*
* Version 1.0.1 2011/09/30
*     Beta release.
*
* Version 2.0.0 2011/12/15
*     Final release. 
* 
* -------------------------------------------------------------------- */     
#include "arm_math.h"     
     
void arm_dct4_f32(     
  const arm_dct4_instance_f32 * S,     
  float32_t * pState,     
  float32_t * pInlineBuffer)     
{     
  uint32_t i;                                    /* Loop counter */     
  float32_t *weights = S->pTwiddle;              /* Pointer to the Weights table */     
  float32_t *cosFact = S->pCosFactor;            /* Pointer to the cos factors table */     
  float32_t *pS1, *pS2, *pbuff;                  /* Temporary pointers for input buffer and pState buffer */     
  float32_t in;                                  /* Temporary variable */     
  float32_t in1, in2, in3, in4;  
     
     
  /* DCT4 computation involves DCT2 (which is calculated using RFFT)      
   * along with some pre-processing and post-processing.      
   * Computational procedure is explained as follows:      
   * (a) Pre-processing involves multiplying input with cos factor,      
   *     r(n) = 2 * u(n) * cos(pi*(2*n+1)/(4*n))      
   *              where,      
   *                 r(n) -- output of preprocessing      
   *                 u(n) -- input to preprocessing(actual Source buffer)      
   * (b) Calculation of DCT2 using FFT is divided into three steps:      
   *                  Step1: Re-ordering of even and odd elements of input.      
   *                  Step2: Calculating FFT of the re-ordered input.      
   *                  Step3: Taking the real part of the product of FFT output and weights.      
   * (c) Post-processing - DCT4 can be obtained from DCT2 output using the following equation:      
   *                   Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)      
   *                        where,      
   *                           Y4 -- DCT4 output,   Y2 -- DCT2 output      
   * (d) Multiplying the output with the normalizing factor sqrt(2/N).      
   */     
     
        /*-------- Pre-processing ------------*/     
  /* Multiplying input with cos factor i.e. r(n) = 2 * x(n) * cos(pi*(2*n+1)/(4*n)) */     
  arm_scale_f32(pInlineBuffer, 2.0f, pInlineBuffer, S->N);     
  arm_mult_f32(pInlineBuffer, cosFact, pInlineBuffer, S->N);     
     
  /* ----------------------------------------------------------------      
   * Step1: Re-ordering of even and odd elements as,      
   *             pState[i] =  pInlineBuffer[2*i] and      
   *             pState[N-i-1] = pInlineBuffer[2*i+1] where i = 0 to N/2      
   ---------------------------------------------------------------------*/     
     
  /* pS1 initialized to pState */     
  pS1 = pState;     
     
  /* pS2 initialized to pState+N-1, so that it points to the end of the state buffer */     
  pS2 = pState + (S->N - 1u);     
     
  /* pbuff initialized to input buffer */     
  pbuff = pInlineBuffer;     
     
  /* Initializing the loop counter to N/2 >> 2 for loop unrolling by 4 */     
  i = (uint32_t) S->Nby2 >> 2u;     
     
  /* First part of the processing with loop unrolling.  Compute 4 outputs at a time.      
   ** a second loop below computes the remaining 1 to 3 samples. */     
  do     
  {     
    /* Re-ordering of even and odd elements */     
    pS2 -= 4u;    
    /* pState[i] =  pInlineBuffer[2*i] */     
    in1 = *pbuff++;  
    *pS1++ = in1;  
    /* pState[N-i-1] = pInlineBuffer[2*i+1] */     
    in2 = *pbuff++;  
    pS2[4] = in2;   
    in3 = *pbuff++;  
    *pS1++ = in3;  
    in4 = *pbuff++;  
    pS2[3] = in4;   
    in1 = *pbuff++;  
    *pS1++ = in1;  
    in2 = *pbuff++;  
    pS2[2] = in2;   
    in3 = *pbuff++;  
    *pS1++ = in3;  
    in4 = *pbuff++;  
    pS2[1] = in4;   
     
    /* Decrement the loop counter */     
    i--;     
  } while(i > 0u);     
     
  /* pbuff initialized to input buffer */     
  pbuff = pInlineBuffer;     
     
  /* pS1 initialized to pState */     
  pS1 = pState;     
     
  /* Initializing the loop counter to N/4 instead of N for loop unrolling */     
  i = (uint32_t) S->N >> 2u;     
     
  /* Processing with loop unrolling 4 times as N is always multiple of 4.      
   * Compute 4 outputs at a time */     
  do     
  {     
    /* Writing the re-ordered output back to inplace input buffer */     
    in1 = *pS1++;  
    *pbuff++ = in1;  
    in2 = *pS1++;  
    *pbuff++ = in2;  
    in3 = *pS1++;  
    *pbuff++ = in3;  
    in4 = *pS1++;    
    *pbuff++ = in4;  
     
    /* Decrement the loop counter */     
    i--;     
  } while(i > 0u);     
     
     
  /* ---------------------------------------------------------      
   *     Step2: Calculate RFFT for N-point input      
   * ---------------------------------------------------------- */     
  /* pInlineBuffer is real input of length N , pState is the complex output of length 2N */     
  arm_rfft_f32(S->pRfft, pInlineBuffer, pState);     
     
        /*----------------------------------------------------------------------      
     *  Step3: Multiply the FFT output with the weights.      
     *----------------------------------------------------------------------*/     
  arm_cmplx_mult_cmplx_f32(pState, weights, pState, S->N);     
     
  /* ----------- Post-processing ---------- */     
  /* DCT-IV can be obtained from DCT-II by the equation,      
   *       Y4(k) = Y2(k) - Y4(k-1) and Y4(-1) = Y4(0)      
   *       Hence, Y4(0) = Y2(0)/2  */     
  /* Getting only real part from the output and Converting to DCT-IV */     
     
  /* Initializing the loop counter to N >> 2 for loop unrolling by 4 */     
  i = ((uint32_t) S->N - 1u) >> 2u;     
     
  /* pbuff initialized to input buffer. */     
  pbuff = pInlineBuffer;     
     
  /* pS1 initialized to pState */     
  pS1 = pState;     
     
  /* Calculating Y4(0) from Y2(0) using Y4(0) = Y2(0)/2 */     
  in = *pS1++ * (float32_t) 0.5;     
  /* input buffer acts as inplace, so output values are stored in the input itself. */     
  *pbuff++ = in;     
     
  /* pState pointer is incremented twice as the real values are located alternatively in the array */     
  pS1++;     
     
  /* First part of the processing with loop unrolling.  Compute 4 outputs at a time.      
   ** a second loop below computes the remaining 1 to 3 samples. */     
  do     
  {     
    /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */     
    /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */     
    in1 = pS1[0];  
    in2 = pS1[2];  
  
    in1 = in1 - in;  
  
    in3 = pS1[4];  
  
    in2 = in2 - in1;  
  
    in4 = pS1[6];  
  
    *pbuff++ = in1;  
  
    in3 = in3 - in2;  
  
    *pbuff++ = in2;  
  
    in = in4 - in3;  
  
    *pbuff++ = in3;  
    *pbuff++ = in;  
  
    pS1 += 8u;       
     
    /* Decrement the loop counter */     
    i--;     
  } while(i > 0u);     
     
  /* If the blockSize is not a multiple of 4, compute any remaining output samples here.      
   ** No loop unrolling is used. */     
  i = ((uint32_t) S->N - 1u) % 0x4u;     
     
  while(i > 0u)     
  {     
    /* Calculating Y4(1) to Y4(N-1) from Y2 using equation Y4(k) = Y2(k) - Y4(k-1) */     
    /* pState pointer (pS1) is incremented twice as the real values are located alternatively in the array */     
    in = *pS1++ - in;     
    *pbuff++ = in;     
    /* points to the next real value */     
    pS1++;     
     
    /* Decrement the loop counter */     
    i--;     
  }     
     
     
        /*------------ Normalizing the output by multiplying with the normalizing factor ----------*/     
     
  /* Initializing the loop counter to N/4 instead of N for loop unrolling */     
  i = (uint32_t) S->N >> 2u;     
     
  /* pbuff initialized to the pInlineBuffer(now contains the output values) */     
  pbuff = pInlineBuffer;     
     
  /* Processing with loop unrolling 4 times as N is always multiple of 4.  Compute 4 outputs at a time */     
  do     
  {     
    /* Multiplying pInlineBuffer with the normalizing factor sqrt(2/N) */     
    in1 = pbuff[0];     
    in2 = pbuff[1];     
    in3 = pbuff[2];     
    in4 = pbuff[3];   
      
    in1 = in1 * S->normalize;    
    in2 = in2 * S->normalize;    
    in3 = in3 * S->normalize;    
    in4 = in4 * S->normalize;  
      
    pbuff[0] = in1;  
    pbuff[1] = in2;  
    pbuff[2] = in3;  
    pbuff[3] = in4;  
  
    pbuff += 4u;     
  
    /* Decrement the loop counter */     
    i--;     
  } while(i > 0u);     
     
}