lib/main/CMSIS/DSP/Source/TransformFunctions/arm_rfft_fast_f32.c

   1 /* ----------------------------------------------------------------------
   2  * Project:      CMSIS DSP Library
   3  * Title:        arm_rfft_f32.c
   4  * Description:  RFFT & RIFFT Floating point process function
   5  *
   6  * $Date:        27. January 2017
   7  * $Revision:    V.1.5.1
   8  *
   9  * Target Processor: Cortex-M cores
  10  * -------------------------------------------------------------------- */
  11 /*
  12  * Copyright (C) 2010-2017 ARM Limited or its affiliates. All rights reserved.
  13  *
  14  * SPDX-License-Identifier: Apache-2.0
  15  *
  16  * Licensed under the Apache License, Version 2.0 (the License); you may
  17  * not use this file except in compliance with the License.
  18  * You may obtain a copy of the License at
  19  *
  20  * www.apache.org/licenses/LICENSE-2.0
  21  *
  22  * Unless required by applicable law or agreed to in writing, software
  23  * distributed under the License is distributed on an AS IS BASIS, WITHOUT
  24  * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  25  * See the License for the specific language governing permissions and
  26  * limitations under the License.
  27  */
  28
  29 #include "arm_math.h"
  30
  31 void stage_rfft_f32(
  32   arm_rfft_fast_instance_f32 * S,
  33   float32_t * p, float32_t * pOut)
  34 {
  35    uint32_t  k;                                                            /* Loop Counter                     */
  36    float32_t twR, twI;                                             /* RFFT Twiddle coefficients        */
  37    float32_t * pCoeff = S->pTwiddleRFFT;  /* Points to RFFT Twiddle factors   */
  38    float32_t *pA = p;                                              /* increasing pointer               */
  39    float32_t *pB = p;                                              /* decreasing pointer               */
  40    float32_t xAR, xAI, xBR, xBI;                                /* temporary variables              */
  41    float32_t t1a, t1b;                                   /* temporary variables              */
  42    float32_t p0, p1, p2, p3;                               /* temporary variables              */
  43
  44
  45    k = (S->Sint).fftLen - 1;
  46
  47    /* Pack first and last sample of the frequency domain together */
  48
  49    xBR = pB[0];
  50    xBI = pB[1];
  51    xAR = pA[0];
  52    xAI = pA[1];
  53
  54    twR = *pCoeff++ ;
  55    twI = *pCoeff++ ;
  56
  57    // U1 = XA(1) + XB(1); % It is real
  58    t1a = xBR + xAR  ;
  59
  60    // U2 = XB(1) - XA(1); % It is imaginary
  61    t1b = xBI + xAI  ;
  62
  63    // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
  64    // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
  65    *pOut++ = 0.5f * ( t1a + t1b );
  66    *pOut++ = 0.5f * ( t1a - t1b );
  67
  68    // XA(1) = 1/2*( U1 - imag(U2) +  i*( U1 +imag(U2) ));
  69    pB  = p + 2*k;
  70    pA += 2;
  71
  72    do
  73    {
  74       /*
  75          function X = my_split_rfft(X, ifftFlag)
  76          % X is a series of real numbers
  77          L  = length(X);
  78          XC = X(1:2:end) +i*X(2:2:end);
  79          XA = fft(XC);
  80          XB = conj(XA([1 end:-1:2]));
  81          TW = i*exp(-2*pi*i*[0:L/2-1]/L).';
  82          for l = 2:L/2
  83             XA(l) = 1/2 * (XA(l) + XB(l) + TW(l) * (XB(l) - XA(l)));
  84          end
  85          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))));
  86          X = XA;
  87       */
  88
  89       xBI = pB[1];
  90       xBR = pB[0];
  91       xAR = pA[0];
  92       xAI = pA[1];
  93
  94       twR = *pCoeff++;
  95       twI = *pCoeff++;
  96
  97       t1a = xBR - xAR ;
  98       t1b = xBI + xAI ;
  99
 100       // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
 101       // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
 102       p0 = twR * t1a;
 103       p1 = twI * t1a;
 104       p2 = twR * t1b;
 105       p3 = twI * t1b;
 106
 107       *pOut++ = 0.5f * (xAR + xBR + p0 + p3 ); //xAR
 108       *pOut++ = 0.5f * (xAI - xBI + p1 - p2 ); //xAI
 109
 110       pA += 2;
 111       pB -= 2;
 112       k--;
 113    } while (k > 0U);
 114 }
 115
 116 /* Prepares data for inverse cfft */
 117 void merge_rfft_f32(
 118 arm_rfft_fast_instance_f32 * S,
 119 float32_t * p, float32_t * pOut)
 120 {
 121    uint32_t  k;                                                         /* Loop Counter                     */
 122    float32_t twR, twI;                                          /* RFFT Twiddle coefficients        */
 123    float32_t *pCoeff = S->pTwiddleRFFT;         /* Points to RFFT Twiddle factors   */
 124    float32_t *pA = p;                                           /* increasing pointer               */
 125    float32_t *pB = p;                                           /* decreasing pointer               */
 126    float32_t xAR, xAI, xBR, xBI;                        /* temporary variables              */
 127    float32_t t1a, t1b, r, s, t, u;                      /* temporary variables              */
 128
 129    k = (S->Sint).fftLen - 1;
 130
 131    xAR = pA[0];
 132    xAI = pA[1];
 133
 134    pCoeff += 2 ;
 135
 136    *pOut++ = 0.5f * ( xAR + xAI );
 137    *pOut++ = 0.5f * ( xAR - xAI );
 138
 139    pB  =  p + 2*k ;
 140    pA +=  2        ;
 141
 142    while (k > 0U)
 143    {
 144       /* G is half of the frequency complex spectrum */
 145       //for k = 2:N
 146       //    Xk(k) = 1/2 * (G(k) + conj(G(N-k+2)) + Tw(k)*( G(k) - conj(G(N-k+2))));
 147       xBI =   pB[1]    ;
 148       xBR =   pB[0]    ;
 149       xAR =  pA[0];
 150       xAI =  pA[1];
 151
 152       twR = *pCoeff++;
 153       twI = *pCoeff++;
 154
 155       t1a = xAR - xBR ;
 156       t1b = xAI + xBI ;
 157
 158       r = twR * t1a;
 159       s = twI * t1b;
 160       t = twI * t1a;
 161       u = twR * t1b;
 162
 163       // real(tw * (xA - xB)) = twR * (xAR - xBR) - twI * (xAI - xBI);
 164       // imag(tw * (xA - xB)) = twI * (xAR - xBR) + twR * (xAI - xBI);
 165       *pOut++ = 0.5f * (xAR + xBR - r - s ); //xAR
 166       *pOut++ = 0.5f * (xAI - xBI + t - u ); //xAI
 167
 168       pA += 2;
 169       pB -= 2;
 170       k--;
 171    }
 172
 173 }
 174
 175 /**
 176 * @ingroup groupTransforms
 177 */
 178
 179 /**
 180  * @defgroup RealFFT Real FFT Functions
 181  *
 182  * \par
 183  * The CMSIS DSP library includes specialized algorithms for computing the
 184  * FFT of real data sequences.  The FFT is defined over complex data but
 185  * in many applications the input is real.  Real FFT algorithms take advantage
 186  * of the symmetry properties of the FFT and have a speed advantage over complex
 187  * algorithms of the same length.
 188  * \par
 189  * The Fast RFFT algorith relays on the mixed radix CFFT that save processor usage.
 190  * \par
 191  * The real length N forward FFT of a sequence is computed using the steps shown below.
 192  * \par
 193  * \image html RFFT.gif "Real Fast Fourier Transform"
 194  * \par
 195  * The real sequence is initially treated as if it were complex to perform a CFFT.
 196  * Later, a processing stage reshapes the data to obtain half of the frequency spectrum
 197  * in complex format. Except the first complex number that contains the two real numbers
 198  * X[0] and X[N/2] all the data is complex. In other words, the first complex sample
 199  * contains two real values packed.
 200  * \par
 201  * The input for the inverse RFFT should keep the same format as the output of the
 202  * forward RFFT. A first processing stage pre-process the data to later perform an
 203  * inverse CFFT.
 204  * \par
 205  * \image html RIFFT.gif "Real Inverse Fast Fourier Transform"
 206  * \par
 207  * The algorithms for floating-point, Q15, and Q31 data are slightly different
 208  * and we describe each algorithm in turn.
 209  * \par Floating-point
 210  * The main functions are arm_rfft_fast_f32() and arm_rfft_fast_init_f32().
 211  * The older functions arm_rfft_f32() and arm_rfft_init_f32() have been
 212  * deprecated but are still documented.
 213  * \par
 214  * The FFT of a real N-point sequence has even symmetry in the frequency
 215  * domain. The second half of the data equals the conjugate of the first
 216  * half flipped in frequency. Looking at the data, we see that we can
 217  * uniquely represent the FFT using only N/2 complex numbers. These are
 218  * packed into the output array in alternating real and imaginary
 219  * components:
 220  * \par
 221  * X = { real[0], imag[0], real[1], imag[1], real[2], imag[2] ...
 222  * real[(N/2)-1], imag[(N/2)-1 }
 223  * \par
 224  * It happens that the first complex number (real[0], imag[0]) is actually
 225  * all real. real[0] represents the DC offset, and imag[0] should be 0.
 226  * (real[1], imag[1]) is the fundamental frequency, (real[2], imag[2]) is
 227  * the first harmonic and so on.
 228  * \par
 229  * The real FFT functions pack the frequency domain data in this fashion.
 230  * The forward transform outputs the data in this form and the inverse
 231  * transform expects input data in this form. The function always performs
 232  * the needed bitreversal so that the input and output data is always in
 233  * normal order. The functions support lengths of [32, 64, 128, ..., 4096]
 234  * samples.
 235  * \par Q15 and Q31
 236  * The real algorithms are defined in a similar manner and utilize N/2 complex
 237  * transforms behind the scenes.
 238  * \par
 239  * The complex transforms used internally include scaling to prevent fixed-point
 240  * overflows.  The overall scaling equals 1/(fftLen/2).
 241  * \par
 242  * A separate instance structure must be defined for each transform used but
 243  * twiddle factor and bit reversal tables can be reused.
 244  * \par
 245  * There is also an associated initialization function for each data type.
 246  * The initialization function performs the following operations:
 247  * - Sets the values of the internal structure fields.
 248  * - Initializes twiddle factor table and bit reversal table pointers.
 249  * - Initializes the internal complex FFT data structure.
 250  * \par
 251  * Use of the initialization function is optional.
 252  * However, if the initialization function is used, then the instance structure
 253  * cannot be placed into a const data section. To place an instance structure
 254  * into a const data section, the instance structure should be manually
 255  * initialized as follows:
 256  * <pre>
 257  *arm_rfft_instance_q31 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
 258  *arm_rfft_instance_q15 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
 259  * </pre>
 260  * where <code>fftLenReal</code> is the length of the real transform;
 261  * <code>fftLenBy2</code> length of  the internal complex transform.
 262  * <code>ifftFlagR</code> Selects forward (=0) or inverse (=1) transform.
 263  * <code>bitReverseFlagR</code> Selects bit reversed output (=0) or normal order
 264  * output (=1).
 265  * <code>twidCoefRModifier</code> stride modifier for the twiddle factor table.
 266  * The value is based on the FFT length;
 267  * <code>pTwiddleAReal</code>points to the A array of twiddle coefficients;
 268  * <code>pTwiddleBReal</code>points to the B array of twiddle coefficients;
 269  * <code>pCfft</code> points to the CFFT Instance structure. The CFFT structure
 270  * must also be initialized.  Refer to arm_cfft_radix4_f32() for details regarding
 271  * static initialization of the complex FFT instance structure.
 272  */
 273
 274 /**
 275 * @addtogroup RealFFT
 276 * @{
 277 */
 278
 279 /**
 280 * @brief Processing function for the floating-point real FFT.
 281 * @param[in]  *S              points to an arm_rfft_fast_instance_f32 structure.
 282 * @param[in]  *p              points to the input buffer.
 283 * @param[in]  *pOut           points to the output buffer.
 284 * @param[in]  ifftFlag        RFFT if flag is 0, RIFFT if flag is 1
 285 * @return none.
 286 */
 287
 288 void arm_rfft_fast_f32(
 289 arm_rfft_fast_instance_f32 * S,
 290 float32_t * p, float32_t * pOut,
 291 uint8_t ifftFlag)
 292 {
 293    arm_cfft_instance_f32 * Sint = &(S->Sint);
 294    Sint->fftLen = S->fftLenRFFT / 2;
 295
 296    /* Calculation of Real FFT */
 297    if (ifftFlag)
 298    {
 299       /*  Real FFT compression */
 300       merge_rfft_f32(S, p, pOut);
 301
 302       /* Complex radix-4 IFFT process */
 303       arm_cfft_f32( Sint, pOut, ifftFlag, 1);
 304    }
 305    else
 306    {
 307       /* Calculation of RFFT of input */
 308       arm_cfft_f32( Sint, p, ifftFlag, 1);
 309
 310       /*  Real FFT extraction */
 311       stage_rfft_f32(S, p, pOut);
 312    }
 313 }
 314
 315 /**
 316 * @} end of RealFFT group
 317 */