git.gag.com Git - fw/stlink/blob - exampleF4/CMSIS/DSP_Lib/Source/TransformFunctions/arm_rfft_f32.c

   1 /* ----------------------------------------------------------------------
   2 * Copyright (C) 2010 ARM Limited. All rights reserved.
   3 *
   4 * $Date:        15. July 2011
   5 * $Revision:    V1.0.10
   6 *
   7 * Project:          CMSIS DSP Library
   8 * Title:            arm_rfft_f32.c
   9 *
  10 * Description:  RFFT & RIFFT Floating point process function
  11 *
  12 * Target Processor: Cortex-M4/Cortex-M3/Cortex-M0
  13 *
  14 * Version 1.0.10 2011/7/15
  15 *    Big Endian support added and Merged M0 and M3/M4 Source code.
  16 *
  17 * Version 1.0.3 2010/11/29
  18 *    Re-organized the CMSIS folders and updated documentation.
  19 *
  20 * Version 1.0.2 2010/11/11
  21 *    Documentation updated.
  22 *
  23 * Version 1.0.1 2010/10/05
  24 *    Production release and review comments incorporated.
  25 *
  26 * Version 1.0.0 2010/09/20
  27 *    Production release and review comments incorporated.
  28 *
  29 * Version 0.0.7  2010/06/10
  30 *    Misra-C changes done
  31 * -------------------------------------------------------------------- */
  32
  33 #include "arm_math.h"
  34
  35 /**
  36  * @ingroup groupTransforms
  37  */
  38
  39 /**
  40  * @defgroup RFFT_RIFFT Real FFT Functions
  41  *
  42  * \par
  43  * Complex FFT/IFFT typically assumes complex input and output. However many applications use real valued data in time domain.
  44  * Real FFT/IFFT efficiently process real valued sequences with the advantage of requirement of low memory and with less complexity.
  45  *
  46  * \par
  47  * This set of functions implements Real Fast Fourier Transforms(RFFT) and Real Inverse Fast Fourier Transform(RIFFT)
  48  * for Q15, Q31, and floating-point data types.
  49  *
  50  *
  51  * \par Algorithm:
  52  *
  53  * <b>Real Fast Fourier Transform:</b>
  54  * \par
  55  * Real FFT of N-point is calculated using CFFT of N/2-point and Split RFFT process as shown below figure.
  56  * \par
  57  * \image html RFFT.gif "Real Fast Fourier Transform"
  58  * \par
  59  * The RFFT functions operate on blocks of input and output data and each call to the function processes
  60  * <code>fftLenR</code> samples through the transform.  <code>pSrc</code>  points to input array containing <code>fftLenR</code> values.
  61  * <code>pDst</code>  points to output array containing <code>2*fftLenR</code> values. \n
  62  * Input for real FFT is in the order of
  63  * <pre>{real[0], real[1], real[2], real[3], ..}</pre>
  64  * Output for real FFT is complex and are in the order of
  65  * <pre>{real(0), imag(0), real(1), imag(1), ...}</pre>
  66  *
  67  * <b>Real Inverse Fast Fourier Transform:</b>
  68  * \par
  69  * Real IFFT of N-point is calculated using Split RIFFT process and CFFT of N/2-point as shown below figure.
  70  * \par
  71  * \image html RIFFT.gif "Real Inverse Fast Fourier Transform"
  72  * \par
  73  * The RIFFT functions operate on blocks of input and output data and each call to the function processes
  74  * <code>2*fftLenR</code> samples through the transform.  <code>pSrc</code>  points to input array containing <code>2*fftLenR</code> values.
  75  * <code>pDst</code>  points to output array containing <code>fftLenR</code> values. \n
  76  * Input for real IFFT is complex and are in the order of
  77  * <pre>{real(0), imag(0), real(1), imag(1), ...}</pre>
  78  *  Output for real IFFT is real and in the order of
  79  * <pre>{real[0], real[1], real[2], real[3], ..}</pre>
  80  *
  81  * \par Lengths supported by the transform:
  82  * \par
  83  * Real FFT/IFFT supports the lengths [128, 512, 2048], as it internally uses CFFT/CIFFT.
  84  *
  85  * \par Instance Structure
  86  * A separate instance structure must be defined for each Instance but the twiddle factors can be reused.
  87  * There are separate instance structure declarations for each of the 3 supported data types.
  88  *
  89  * \par Initialization Functions
  90  * There is also an associated initialization function for each data type.
  91  * The initialization function performs the following operations:
  92  * - Sets the values of the internal structure fields.
  93  * - Initializes twiddle factor tables.
  94  * - Initializes CFFT data structure fields.
  95  * \par
  96  * Use of the initialization function is optional.
  97  * However, if the initialization function is used, then the instance structure cannot be placed into a const data section.
  98  * To place an instance structure into a const data section, the instance structure must be manually initialized.
  99  * Manually initialize the instance structure as follows:
 100  * <pre>
 101  *arm_rfft_instance_f32 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
 102  *arm_rfft_instance_q31 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
 103  *arm_rfft_instance_q15 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
 104  * </pre>
 105  * where <code>fftLenReal</code> length of RFFT/RIFFT; <code>fftLenBy2</code> length of CFFT/CIFFT.
 106  * <code>ifftFlagR</code> Flag for selection of RFFT or RIFFT(Set ifftFlagR to calculate RIFFT otherwise calculates RFFT);
 107  * <code>bitReverseFlagR</code> Flag for selection of output order(Set bitReverseFlagR to output in normal order otherwise output in bit reversed order);
 108  * <code>twidCoefRModifier</code> modifier for twiddle factor table which supports 128, 512, 2048 RFFT lengths with same table;
 109  * <code>pTwiddleAReal</code>points to A array of twiddle coefficients; <code>pTwiddleBReal</code>points to B array of twiddle coefficients;
 110  * <code>pCfft</code> points to the CFFT Instance structure. The CFFT structure also needs to be initialized, refer to arm_cfft_radix4_f32() for details regarding
 111  * static initialization of cfft structure.
 112  *
 113  * \par Fixed-Point Behavior
 114  * Care must be taken when using the fixed-point versions of the RFFT/RIFFT function.
 115  * Refer to the function specific documentation below for usage guidelines.
 116  */
 117
 118 /*--------------------------------------------------------------------
 119  *              Internal functions prototypes
 120  *--------------------------------------------------------------------*/
 121
 122 void arm_split_rfft_f32(
 123   float32_t * pSrc,
 124   uint32_t fftLen,
 125   float32_t * pATable,
 126   float32_t * pBTable,
 127   float32_t * pDst,
 128   uint32_t modifier);
 129 void arm_split_rifft_f32(
 130   float32_t * pSrc,
 131   uint32_t fftLen,
 132   float32_t * pATable,
 133   float32_t * pBTable,
 134   float32_t * pDst,
 135   uint32_t modifier);
 136
 137 /**
 138  * @addtogroup RFFT_RIFFT
 139  * @{
 140  */
 141
 142 /**
 143  * @brief Processing function for the floating-point RFFT/RIFFT.
 144  * @param[in]  *S    points to an instance of the floating-point RFFT/RIFFT structure.
 145  * @param[in]  *pSrc points to the input buffer.
 146  * @param[out] *pDst points to the output buffer.
 147  * @return none.
 148  */
 149
 150 void arm_rfft_f32(
 151   const arm_rfft_instance_f32 * S,
 152   float32_t * pSrc,
 153   float32_t * pDst)
 154 {
 155   const arm_cfft_radix4_instance_f32 *S_CFFT = S->pCfft;
 156
 157
 158   /* Calculation of Real IFFT of input */
 159   if(S->ifftFlagR == 1u)
 160   {
 161     /*  Real IFFT core process */
 162     arm_split_rifft_f32(pSrc, S->fftLenBy2, S->pTwiddleAReal,
 163                         S->pTwiddleBReal, pDst, S->twidCoefRModifier);
 164
 165
 166     /* Complex radix-4 IFFT process */
 167     arm_radix4_butterfly_inverse_f32(pDst, S_CFFT->fftLen,
 168                                      S_CFFT->pTwiddle,
 169                                      S_CFFT->twidCoefModifier,
 170                                      S_CFFT->onebyfftLen);
 171
 172     /* Bit reversal process */
 173     if(S->bitReverseFlagR == 1u)
 174     {
 175       arm_bitreversal_f32(pDst, S_CFFT->fftLen,
 176                           S_CFFT->bitRevFactor, S_CFFT->pBitRevTable);
 177     }
 178   }
 179   else
 180   {
 181
 182     /* Calculation of RFFT of input */
 183
 184     /* Complex radix-4 FFT process */
 185     arm_radix4_butterfly_f32(pSrc, S_CFFT->fftLen,
 186                              S_CFFT->pTwiddle, S_CFFT->twidCoefModifier);
 187
 188     /* Bit reversal process */
 189     if(S->bitReverseFlagR == 1u)
 190     {
 191       arm_bitreversal_f32(pSrc, S_CFFT->fftLen,
 192                           S_CFFT->bitRevFactor, S_CFFT->pBitRevTable);
 193     }
 194
 195
 196     /*  Real FFT core process */
 197     arm_split_rfft_f32(pSrc, S->fftLenBy2, S->pTwiddleAReal,
 198                        S->pTwiddleBReal, pDst, S->twidCoefRModifier);
 199   }
 200
 201 }
 202
 203 /**
 204    * @} end of RFFT_RIFFT group
 205    */
 206
 207 /**
 208  * @brief  Core Real FFT process
 209  * @param[in]   *pSrc                           points to the input buffer.
 210  * @param[in]   fftLen                          length of FFT.
 211  * @param[in]   *pATable                        points to the twiddle Coef A buffer.
 212  * @param[in]   *pBTable                        points to the twiddle Coef B buffer.
 213  * @param[out]  *pDst                           points to the output buffer.
 214  * @param[in]   modifier                twiddle coefficient modifier that supports different size FFTs with the same twiddle factor table.
 215  * @return none.
 216  */
 217
 218 void arm_split_rfft_f32(
 219   float32_t * pSrc,
 220   uint32_t fftLen,
 221   float32_t * pATable,
 222   float32_t * pBTable,
 223   float32_t * pDst,
 224   uint32_t modifier)
 225 {
 226   uint32_t i;                                    /* Loop Counter */
 227   float32_t outR, outI;                          /* Temporary variables for output */
 228   float32_t *pCoefA, *pCoefB;                    /* Temporary pointers for twiddle factors */
 229   float32_t CoefA1, CoefA2, CoefB1;              /* Temporary variables for twiddle coefficients */
 230   float32_t *pDst1 = &pDst[2], *pDst2 = &pDst[(4u * fftLen) - 1u];      /* temp pointers for output buffer */
 231   float32_t *pSrc1 = &pSrc[2], *pSrc2 = &pSrc[(2u * fftLen) - 1u];      /* temp pointers for input buffer */
 232
 233
 234   pSrc[2u * fftLen] = pSrc[0];
 235   pSrc[(2u * fftLen) + 1u] = pSrc[1];
 236
 237   /* Init coefficient pointers */
 238   pCoefA = &pATable[modifier * 2u];
 239   pCoefB = &pBTable[modifier * 2u];
 240
 241   i = fftLen - 1u;
 242
 243   while(i > 0u)
 244   {
 245     /*
 246        outR = (pSrc[2 * i] * pATable[2 * i] - pSrc[2 * i + 1] * pATable[2 * i + 1]
 247        + pSrc[2 * n - 2 * i] * pBTable[2 * i] +
 248        pSrc[2 * n - 2 * i + 1] * pBTable[2 * i + 1]);
 249      */
 250
 251     /* outI = (pIn[2 * i + 1] * pATable[2 * i] + pIn[2 * i] * pATable[2 * i + 1] +
 252        pIn[2 * n - 2 * i] * pBTable[2 * i + 1] -
 253        pIn[2 * n - 2 * i + 1] * pBTable[2 * i]); */
 254
 255     /* read pATable[2 * i] */
 256     CoefA1 = *pCoefA++;
 257     /* pATable[2 * i + 1] */
 258     CoefA2 = *pCoefA;
 259
 260     /* pSrc[2 * i] * pATable[2 * i] */
 261     outR = *pSrc1 * CoefA1;
 262     /* pSrc[2 * i] * CoefA2 */
 263     outI = *pSrc1++ * CoefA2;
 264
 265     /* (pSrc[2 * i + 1] + pSrc[2 * fftLen - 2 * i + 1]) * CoefA2 */
 266     outR -= (*pSrc1 + *pSrc2) * CoefA2;
 267     /* pSrc[2 * i + 1] * CoefA1 */
 268     outI += *pSrc1++ * CoefA1;
 269
 270     CoefB1 = *pCoefB;
 271
 272     /* pSrc[2 * fftLen - 2 * i + 1] * CoefB1 */
 273     outI -= *pSrc2-- * CoefB1;
 274     /* pSrc[2 * fftLen - 2 * i] * CoefA2 */
 275     outI -= *pSrc2 * CoefA2;
 276
 277     /* pSrc[2 * fftLen - 2 * i] * CoefB1 */
 278     outR += *pSrc2-- * CoefB1;
 279
 280     /* write output */
 281     *pDst1++ = outR;
 282     *pDst1++ = outI;
 283
 284     /* write complex conjugate output */
 285     *pDst2-- = -outI;
 286     *pDst2-- = outR;
 287
 288     /* update coefficient pointer */
 289     pCoefB = pCoefB + (modifier * 2u);
 290     pCoefA = pCoefA + ((modifier * 2u) - 1u);
 291
 292     i--;
 293
 294   }
 295
 296   pDst[2u * fftLen] = pSrc[0] - pSrc[1];
 297   pDst[(2u * fftLen) + 1u] = 0.0f;
 298
 299   pDst[0] = pSrc[0] + pSrc[1];
 300   pDst[1] = 0.0f;
 301
 302 }
 303
 304
 305 /**
 306  * @brief  Core Real IFFT process
 307  * @param[in]   *pSrc                           points to the input buffer.
 308  * @param[in]   fftLen                          length of FFT.
 309  * @param[in]   *pATable                        points to the twiddle Coef A buffer.
 310  * @param[in]   *pBTable                        points to the twiddle Coef B buffer.
 311  * @param[out]  *pDst                           points to the output buffer.
 312  * @param[in]   modifier                twiddle coefficient modifier that supports different size FFTs with the same twiddle factor table.
 313  * @return none.
 314  */
 315
 316 void arm_split_rifft_f32(
 317   float32_t * pSrc,
 318   uint32_t fftLen,
 319   float32_t * pATable,
 320   float32_t * pBTable,
 321   float32_t * pDst,
 322   uint32_t modifier)
 323 {
 324   float32_t outR, outI;                          /* Temporary variables for output */
 325   float32_t *pCoefA, *pCoefB;                    /* Temporary pointers for twiddle factors */
 326   float32_t CoefA1, CoefA2, CoefB1;              /* Temporary variables for twiddle coefficients */
 327   float32_t *pSrc1 = &pSrc[0], *pSrc2 = &pSrc[(2u * fftLen) + 1u];
 328
 329   pCoefA = &pATable[0];
 330   pCoefB = &pBTable[0];
 331
 332   while(fftLen > 0u)
 333   {
 334     /*
 335        outR = (pIn[2 * i] * pATable[2 * i] + pIn[2 * i + 1] * pATable[2 * i + 1] +
 336        pIn[2 * n - 2 * i] * pBTable[2 * i] -
 337        pIn[2 * n - 2 * i + 1] * pBTable[2 * i + 1]);
 338
 339        outI = (pIn[2 * i + 1] * pATable[2 * i] - pIn[2 * i] * pATable[2 * i + 1] -
 340        pIn[2 * n - 2 * i] * pBTable[2 * i + 1] -
 341        pIn[2 * n - 2 * i + 1] * pBTable[2 * i]);
 342
 343      */
 344
 345     CoefA1 = *pCoefA++;
 346     CoefA2 = *pCoefA;
 347
 348     /* outR = (pSrc[2 * i] * CoefA1 */
 349     outR = *pSrc1 * CoefA1;
 350
 351     /* - pSrc[2 * i] * CoefA2 */
 352     outI = -(*pSrc1++) * CoefA2;
 353
 354     /* (pSrc[2 * i + 1] + pSrc[2 * fftLen - 2 * i + 1]) * CoefA2 */
 355     outR += (*pSrc1 + *pSrc2) * CoefA2;
 356
 357     /* pSrc[2 * i + 1] * CoefA1 */
 358     outI += (*pSrc1++) * CoefA1;
 359
 360     CoefB1 = *pCoefB;
 361
 362     /* - pSrc[2 * fftLen - 2 * i + 1] * CoefB1 */
 363     outI -= *pSrc2-- * CoefB1;
 364
 365     /* pSrc[2 * fftLen - 2 * i] * CoefB1 */
 366     outR += *pSrc2 * CoefB1;
 367
 368     /* pSrc[2 * fftLen - 2 * i] * CoefA2 */
 369     outI += *pSrc2-- * CoefA2;
 370
 371     /* write output */
 372     *pDst++ = outR;
 373     *pDst++ = outI;
 374
 375     /* update coefficient pointer */
 376     pCoefB = pCoefB + (modifier * 2u);
 377     pCoefA = pCoefA + ((modifier * 2u) - 1u);
 378
 379     /* Decrement loop count */
 380     fftLen--;
 381   }
 382
 383 }