Subversion Repositories DashDisplay

/* ----------------------------------------------------------------------

* Copyright (C) 2010-2014 ARM Limited. All rights reserved.

*

* $Date:        19. March 2015

* $Revision:    V.1.4.5

*

* Project:          CMSIS DSP Library

* Title:            arm_rfft_f32.c

*

* Description:  RFFT & RIFFT Floating point process function

*

* Target Processor: Cortex-M4/Cortex-M3/Cortex-M0

*

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions

* are met:

*   - Redistributions of source code must retain the above copyright

*     notice, this list of conditions and the following disclaimer.

*   - Redistributions in binary form must reproduce the above copyright

*     notice, this list of conditions and the following disclaimer in

*     the documentation and/or other materials provided with the

*     distribution.

*   - Neither the name of ARM LIMITED nor the names of its contributors

*     may be used to endorse or promote products derived from this

*     software without specific prior written permission.

*

* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS

* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE

* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,

* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,

* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;

* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER

* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT

* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN

* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

* POSSIBILITY OF SUCH DAMAGE.

* -------------------------------------------------------------------- */

#include "arm_math.h"

void stage_rfft_f32(

  arm_rfft_fast_instance_f32 * S,

  float32_t * p, float32_t * pOut)

{

   uint32_t  k;                                                            /* Loop Counter                     */

   float32_t twR, twI;                                             /* RFFT Twiddle coefficients        */

   float32_t * pCoeff = S->pTwiddleRFFT;  /* Points to RFFT Twiddle factors   */

   float32_t *pA = p;                                              /* increasing pointer               */

   float32_t *pB = p;                                              /* decreasing pointer               */

   float32_t xAR, xAI, xBR, xBI;                                /* temporary variables              */

   float32_t t1a, t1b;                                   /* temporary variables              */

   float32_t p0, p1, p2, p3;                               /* temporary variables              */

   k = (S->Sint).fftLen - 1;                                   

   /* Pack first and last sample of the frequency domain together */

   xBR = pB[0];

   xBI = pB[1];

   xAR = pA[0];

   xAI = pA[1];

   twR = *pCoeff++ ;

   twI = *pCoeff++ ;

   // U1 = XA(1) + XB(1); % It is real

   t1a = xBR + xAR  ;

   // U2 = XB(1) - XA(1); % It is imaginary

   t1b = xBI + xAI  ;

   // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);

   // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);

   *pOut++ = 0.5f * ( t1a + t1b );

   *pOut++ = 0.5f * ( t1a - t1b );

   // XA(1) = 1/2*( U1 - imag(U2) +  i*( U1 +imag(U2) ));

   pB  = p + 2*k;

   pA += 2;

   do

   {

      /*

         function X = my_split_rfft(X, ifftFlag)

         % X is a series of real numbers

         L  = length(X);

         XC = X(1:2:end) +i*X(2:2:end);

         XA = fft(XC);

         XB = conj(XA([1 end:-1:2]));

         TW = i*exp(-2*pi*i*[0:L/2-1]/L).';

         for l = 2:L/2

            XA(l) = 1/2 * (XA(l) + XB(l) + TW(l) * (XB(l) - XA(l)));

         end

         XA(1) = 1/2* (XA(1) + XB(1) + TW(1) * (XB(1) - XA(1))) + i*( 1/2*( XA(1) + XB(1) + i*( XA(1) - XB(1))));

         X = XA;

      */

      xBI = pB[1];

      xBR = pB[0];

      xAR = pA[0];

      xAI = pA[1];

      twR = *pCoeff++;

      twI = *pCoeff++;

      t1a = xBR - xAR ;

      t1b = xBI + xAI ;

      // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);

      // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);

      p0 = twR * t1a;

      p1 = twI * t1a;

      p2 = twR * t1b;

      p3 = twI * t1b;

      *pOut++ = 0.5f * (xAR + xBR + p0 + p3 ); //xAR

      *pOut++ = 0.5f * (xAI - xBI + p1 - p2 ); //xAI

      pA += 2;

      pB -= 2;

      k--;

   } while(k > 0u);

}

/* Prepares data for inverse cfft */

void merge_rfft_f32(

arm_rfft_fast_instance_f32 * S,

float32_t * p, float32_t * pOut)

{

   uint32_t  k;                                                         /* Loop Counter                     */

   float32_t twR, twI;                                          /* RFFT Twiddle coefficients        */

   float32_t *pCoeff = S->pTwiddleRFFT;         /* Points to RFFT Twiddle factors   */

   float32_t *pA = p;                                           /* increasing pointer               */

   float32_t *pB = p;                                           /* decreasing pointer               */

   float32_t xAR, xAI, xBR, xBI;                        /* temporary variables              */

   float32_t t1a, t1b, r, s, t, u;                      /* temporary variables              */

   k = (S->Sint).fftLen - 1;                                   

   xAR = pA[0];

   xAI = pA[1];

   pCoeff += 2 ;

   *pOut++ = 0.5f * ( xAR + xAI );

   *pOut++ = 0.5f * ( xAR - xAI );

   pB  =  p + 2*k ;

   pA +=  2        ;

   while(k > 0u)

   {

      /* G is half of the frequency complex spectrum */

      //for k = 2:N

      //    Xk(k) = 1/2 * (G(k) + conj(G(N-k+2)) + Tw(k)*( G(k) - conj(G(N-k+2))));

      xBI =   pB[1]    ;

      xBR =   pB[0]    ;

      xAR =  pA[0];

      xAI =  pA[1];

      twR = *pCoeff++;

      twI = *pCoeff++;

      t1a = xAR - xBR ;

      t1b = xAI + xBI ;

      r = twR * t1a;

      s = twI * t1b;

      t = twI * t1a;

      u = twR * t1b;

      // real(tw * (xA - xB)) = twR * (xAR - xBR) - twI * (xAI - xBI);

      // imag(tw * (xA - xB)) = twI * (xAR - xBR) + twR * (xAI - xBI);

      *pOut++ = 0.5f * (xAR + xBR - r - s ); //xAR

      *pOut++ = 0.5f * (xAI - xBI + t - u ); //xAI

      pA += 2;

      pB -= 2;

      k--;

   }

}

/**

* @ingroup groupTransforms

*/

/**

 * @defgroup Fast Real FFT Functions

 *

 * \par

 * The CMSIS DSP library includes specialized algorithms for computing the

 * FFT of real data sequences.  The FFT is defined over complex data but

 * in many applications the input is real.  Real FFT algorithms take advantage

 * of the symmetry properties of the FFT and have a speed advantage over complex

 * algorithms of the same length.

 * \par

 * The Fast RFFT algorith relays on the mixed radix CFFT that save processor usage.

 * \par

 * The real length N forward FFT of a sequence is computed using the steps shown below.

 * \par

 * \image html RFFT.gif "Real Fast Fourier Transform"

 * \par

 * The real sequence is initially treated as if it were complex to perform a CFFT.

 * Later, a processing stage reshapes the data to obtain half of the frequency spectrum

 * in complex format. Except the first complex number that contains the two real numbers

 * X[0] and X[N/2] all the data is complex. In other words, the first complex sample

 * contains two real values packed.

 * \par

 * The input for the inverse RFFT should keep the same format as the output of the

 * forward RFFT. A first processing stage pre-process the data to later perform an

 * inverse CFFT.

 * \par    

 * \image html RIFFT.gif "Real Inverse Fast Fourier Transform"    

 * \par    

 * The algorithms for floating-point, Q15, and Q31 data are slightly different

 * and we describe each algorithm in turn.

 * \par Floating-point

 * The main functions are <code>arm_rfft_fast_f32()</code>

 * and <code>arm_rfft_fast_init_f32()</code>.  The older functions

 * <code>arm_rfft_f32()</code> and <code>arm_rfft_init_f32()</code> have been

 * deprecated but are still documented.

 * \par

 * The FFT of a real N-point sequence has even symmetry in the frequency

 * domain.  The second half of the data equals the conjugate of the first half

 * flipped in frequency:

 * <pre>

 *X[0] - real data

 *X[1] - complex data

 *X[2] - complex data

 *...

 *X[fftLen/2-1] - complex data

 *X[fftLen/2] - real data

 *X[fftLen/2+1] - conjugate of X[fftLen/2-1]

 *X[fftLen/2+2] - conjugate of X[fftLen/2-2]

 *...

 *X[fftLen-1] - conjugate of X[1]

 * </pre>

 * Looking at the data, we see that we can uniquely represent the FFT using only

 * <pre>

 *N/2+1 samples:

 *X[0] - real data

 *X[1] - complex data

 *X[2] - complex data

 *...

 *X[fftLen/2-1] - complex data

 *X[fftLen/2] - real data

 * </pre>

 * Looking more closely we see that the first and last samples are real valued.

 * They can be packed together and we can thus represent the FFT of an N-point

 * real sequence by N/2 complex values:

 * <pre>

 *X[0],X[N/2] - packed real data: X[0] + jX[N/2]

 *X[1] - complex data

 *X[2] - complex data

 *...

 *X[fftLen/2-1] - complex data

 * </pre>

 * The real FFT functions pack the frequency domain data in this fashion.  The

 * forward transform outputs the data in this form and the inverse transform

 * expects input data in this form.  The function always performs the needed

 * bitreversal so that the input and output data is always in normal order.  The

 * functions support lengths of [32, 64, 128, ..., 4096] samples.

 * \par

 * The forward and inverse real FFT functions apply the standard FFT scaling; no

 * scaling on the forward transform and 1/fftLen scaling on the inverse

 * transform.

 * \par Q15 and Q31

 * The real algorithms are defined in a similar manner and utilize N/2 complex

 * transforms behind the scenes.  

 * \par

 * The complex transforms used internally include scaling to prevent fixed-point

 * overflows.  The overall scaling equals 1/(fftLen/2).

 * \par

 * A separate instance structure must be defined for each transform used but

 * twiddle factor and bit reversal tables can be reused.

 * \par

 * There is also an associated initialization function for each data type.

 * The initialization function performs the following operations:

 * - Sets the values of the internal structure fields.  

 * - Initializes twiddle factor table and bit reversal table pointers.

 * - Initializes the internal complex FFT data structure.

 * \par  

 * Use of the initialization function is optional.  

 * However, if the initialization function is used, then the instance structure

 * cannot be placed into a const data section. To place an instance structure

 * into a const data section, the instance structure should be manually

 * initialized as follows:

 * <pre>

 *arm_rfft_instance_q31 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};    

 *arm_rfft_instance_q15 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};    

 * </pre>

 * where <code>fftLenReal</code> is the length of the real transform;

 * <code>fftLenBy2</code> length of  the internal complex transform.

 * <code>ifftFlagR</code> Selects forward (=0) or inverse (=1) transform.

 * <code>bitReverseFlagR</code> Selects bit reversed output (=0) or normal order

 * output (=1).

 * <code>twidCoefRModifier</code> stride modifier for the twiddle factor table.

 * The value is based on the FFT length;

 * <code>pTwiddleAReal</code>points to the A array of twiddle coefficients;

 * <code>pTwiddleBReal</code>points to the B array of twiddle coefficients;    

 * <code>pCfft</code> points to the CFFT Instance structure. The CFFT structure

 * must also be initialized.  Refer to arm_cfft_radix4_f32() for details regarding    

 * static initialization of the complex FFT instance structure.    

 */

/**

* @addtogroup RealFFT

* @{

*/

/**

* @brief Processing function for the floating-point real FFT.

* @param[in]  *S              points to an arm_rfft_fast_instance_f32 structure.

* @param[in]  *p              points to the input buffer.

* @param[in]  *pOut           points to the output buffer.

* @param[in]  ifftFlag        RFFT if flag is 0, RIFFT if flag is 1

* @return none.

*/

void arm_rfft_fast_f32(

arm_rfft_fast_instance_f32 * S,

float32_t * p, float32_t * pOut,

uint8_t ifftFlag)

{

   arm_cfft_instance_f32 * Sint = &(S->Sint);

   Sint->fftLen = S->fftLenRFFT / 2;

   /* Calculation of Real FFT */

   if(ifftFlag)

   {

      /*  Real FFT compression */

      merge_rfft_f32(S, p, pOut);

      /* Complex radix-4 IFFT process */

      arm_cfft_f32( Sint, pOut, ifftFlag, 1);

   }

   else

   {

      /* Calculation of RFFT of input */

      arm_cfft_f32( Sint, p, ifftFlag, 1);

      /*  Real FFT extraction */

      stage_rfft_f32(S, p, pOut);

   }

}

/**    

* @} end of RealFFT group    

*/