arm_rfft_fast_f32.c

00001 /* ----------------------------------------------------------------------
00002  * Project:      CMSIS DSP Library
00003  * Title:        arm_rfft_f32.c
00004  * Description:  RFFT & RIFFT Floating point process function
00005  *
00006  * $Date:        27. January 2017
00007  * $Revision:    V.1.5.1
00008  *
00009  * Target Processor: Cortex-M cores
00010  * -------------------------------------------------------------------- */
00011 /*
00012  * Copyright (C) 2010-2017 ARM Limited or its affiliates. All rights reserved.
00013  *
00014  * SPDX-License-Identifier: Apache-2.0
00015  *
00016  * Licensed under the Apache License, Version 2.0 (the License); you may
00017  * not use this file except in compliance with the License.
00018  * You may obtain a copy of the License at
00019  *
00020  * www.apache.org/licenses/LICENSE-2.0
00021  *
00022  * Unless required by applicable law or agreed to in writing, software
00023  * distributed under the License is distributed on an AS IS BASIS, WITHOUT
00024  * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00025  * See the License for the specific language governing permissions and
00026  * limitations under the License.
00027  */
00028 
00029 #include "arm_math.h"
00030 
00031 void stage_rfft_f32(
00032   arm_rfft_fast_instance_f32 * S,
00033   float32_t * p, float32_t * pOut)
00034 {
00035    uint32_t  k;                                /* Loop Counter                     */
00036    float32_t twR, twI;                         /* RFFT Twiddle coefficients        */
00037    float32_t * pCoeff = S->pTwiddleRFFT;  /* Points to RFFT Twiddle factors   */
00038    float32_t *pA = p;                          /* increasing pointer               */
00039    float32_t *pB = p;                          /* decreasing pointer               */
00040    float32_t xAR, xAI, xBR, xBI;                /* temporary variables              */
00041    float32_t t1a, t1b;                       /* temporary variables              */
00042    float32_t p0, p1, p2, p3;                   /* temporary variables              */
00043 
00044 
00045    k = (S->Sint).fftLen - 1;
00046 
00047    /* Pack first and last sample of the frequency domain together */
00048 
00049    xBR = pB[0];
00050    xBI = pB[1];
00051    xAR = pA[0];
00052    xAI = pA[1];
00053 
00054    twR = *pCoeff++ ;
00055    twI = *pCoeff++ ;
00056 
00057    // U1 = XA(1) + XB(1); % It is real
00058    t1a = xBR + xAR  ;
00059 
00060    // U2 = XB(1) - XA(1); % It is imaginary
00061    t1b = xBI + xAI  ;
00062 
00063    // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
00064    // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
00065    *pOut++ = 0.5f * ( t1a + t1b );
00066    *pOut++ = 0.5f * ( t1a - t1b );
00067 
00068    // XA(1) = 1/2*( U1 - imag(U2) +  i*( U1 +imag(U2) ));
00069    pB  = p + 2*k;
00070    pA += 2;
00071 
00072    do
00073    {
00074       /*
00075          function X = my_split_rfft(X, ifftFlag)
00076          % X is a series of real numbers
00077          L  = length(X);
00078          XC = X(1:2:end) +i*X(2:2:end);
00079          XA = fft(XC);
00080          XB = conj(XA([1 end:-1:2]));
00081          TW = i*exp(-2*pi*i*[0:L/2-1]/L).';
00082          for l = 2:L/2
00083             XA(l) = 1/2 * (XA(l) + XB(l) + TW(l) * (XB(l) - XA(l)));
00084          end
00085          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))));
00086          X = XA;
00087       */
00088 
00089       xBI = pB[1];
00090       xBR = pB[0];
00091       xAR = pA[0];
00092       xAI = pA[1];
00093 
00094       twR = *pCoeff++;
00095       twI = *pCoeff++;
00096 
00097       t1a = xBR - xAR ;
00098       t1b = xBI + xAI ;
00099 
00100       // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
00101       // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
00102       p0 = twR * t1a;
00103       p1 = twI * t1a;
00104       p2 = twR * t1b;
00105       p3 = twI * t1b;
00106 
00107       *pOut++ = 0.5f * (xAR + xBR + p0 + p3 ); //xAR
00108       *pOut++ = 0.5f * (xAI - xBI + p1 - p2 ); //xAI
00109 
00110       pA += 2;
00111       pB -= 2;
00112       k--;
00113    } while (k > 0U);
00114 }
00115 
00116 /* Prepares data for inverse cfft */
00117 void merge_rfft_f32(
00118 arm_rfft_fast_instance_f32 * S,
00119 float32_t * p, float32_t * pOut)
00120 {
00121    uint32_t  k;                             /* Loop Counter                     */
00122    float32_t twR, twI;                      /* RFFT Twiddle coefficients        */
00123    float32_t *pCoeff = S->pTwiddleRFFT;     /* Points to RFFT Twiddle factors   */
00124    float32_t *pA = p;                       /* increasing pointer               */
00125    float32_t *pB = p;                       /* decreasing pointer               */
00126    float32_t xAR, xAI, xBR, xBI;            /* temporary variables              */
00127    float32_t t1a, t1b, r, s, t, u;          /* temporary variables              */
00128 
00129    k = (S->Sint).fftLen - 1;
00130 
00131    xAR = pA[0];
00132    xAI = pA[1];
00133 
00134    pCoeff += 2 ;
00135 
00136    *pOut++ = 0.5f * ( xAR + xAI );
00137    *pOut++ = 0.5f * ( xAR - xAI );
00138 
00139    pB  =  p + 2*k ;
00140    pA +=  2    ;
00141 
00142    while (k > 0U)
00143    {
00144       /* G is half of the frequency complex spectrum */
00145       //for k = 2:N
00146       //    Xk(k) = 1/2 * (G(k) + conj(G(N-k+2)) + Tw(k)*( G(k) - conj(G(N-k+2))));
00147       xBI =   pB[1]    ;
00148       xBR =   pB[0]    ;
00149       xAR =  pA[0];
00150       xAI =  pA[1];
00151 
00152       twR = *pCoeff++;
00153       twI = *pCoeff++;
00154 
00155       t1a = xAR - xBR ;
00156       t1b = xAI + xBI ;
00157 
00158       r = twR * t1a;
00159       s = twI * t1b;
00160       t = twI * t1a;
00161       u = twR * t1b;
00162 
00163       // real(tw * (xA - xB)) = twR * (xAR - xBR) - twI * (xAI - xBI);
00164       // imag(tw * (xA - xB)) = twI * (xAR - xBR) + twR * (xAI - xBI);
00165       *pOut++ = 0.5f * (xAR + xBR - r - s ); //xAR
00166       *pOut++ = 0.5f * (xAI - xBI + t - u ); //xAI
00167 
00168       pA += 2;
00169       pB -= 2;
00170       k--;
00171    }
00172 
00173 }
00174 
00175 /**
00176 * @ingroup groupTransforms
00177 */
00178 
00179 /**
00180  * @defgroup RealFFT Real FFT Functions
00181  *
00182  * \par
00183  * The CMSIS DSP library includes specialized algorithms for computing the
00184  * FFT of real data sequences.  The FFT is defined over complex data but
00185  * in many applications the input is real.  Real FFT algorithms take advantage
00186  * of the symmetry properties of the FFT and have a speed advantage over complex
00187  * algorithms of the same length.
00188  * \par
00189  * The Fast RFFT algorith relays on the mixed radix CFFT that save processor usage.
00190  * \par
00191  * The real length N forward FFT of a sequence is computed using the steps shown below.
00192  * \par
00193  * \image html RFFT.gif "Real Fast Fourier Transform"
00194  * \par
00195  * The real sequence is initially treated as if it were complex to perform a CFFT.
00196  * Later, a processing stage reshapes the data to obtain half of the frequency spectrum
00197  * in complex format. Except the first complex number that contains the two real numbers
00198  * X[0] and X[N/2] all the data is complex. In other words, the first complex sample
00199  * contains two real values packed.
00200  * \par
00201  * The input for the inverse RFFT should keep the same format as the output of the
00202  * forward RFFT. A first processing stage pre-process the data to later perform an
00203  * inverse CFFT.
00204  * \par
00205  * \image html RIFFT.gif "Real Inverse Fast Fourier Transform"
00206  * \par
00207  * The algorithms for floating-point, Q15, and Q31 data are slightly different
00208  * and we describe each algorithm in turn.
00209  * \par Floating-point
00210  * The main functions are arm_rfft_fast_f32() and arm_rfft_fast_init_f32().
00211  * The older functions arm_rfft_f32() and arm_rfft_init_f32() have been
00212  * deprecated but are still documented.
00213  * \par
00214  * The FFT of a real N-point sequence has even symmetry in the frequency
00215  * domain. The second half of the data equals the conjugate of the first
00216  * half flipped in frequency. Looking at the data, we see that we can
00217  * uniquely represent the FFT using only N/2 complex numbers. These are
00218  * packed into the output array in alternating real and imaginary
00219  * components:
00220  * \par
00221  * X = { real[0], imag[0], real[1], imag[1], real[2], imag[2] ...
00222  * real[(N/2)-1], imag[(N/2)-1 }
00223  * \par
00224  * It happens that the first complex number (real[0], imag[0]) is actually
00225  * all real. real[0] represents the DC offset, and imag[0] should be 0.
00226  * (real[1], imag[1]) is the fundamental frequency, (real[2], imag[2]) is
00227  * the first harmonic and so on.
00228  * \par
00229  * The real FFT functions pack the frequency domain data in this fashion.
00230  * The forward transform outputs the data in this form and the inverse
00231  * transform expects input data in this form. The function always performs
00232  * the needed bitreversal so that the input and output data is always in
00233  * normal order. The functions support lengths of [32, 64, 128, ..., 4096]
00234  * samples.
00235  * \par Q15 and Q31
00236  * The real algorithms are defined in a similar manner and utilize N/2 complex
00237  * transforms behind the scenes.
00238  * \par
00239  * The complex transforms used internally include scaling to prevent fixed-point
00240  * overflows.  The overall scaling equals 1/(fftLen/2).
00241  * \par
00242  * A separate instance structure must be defined for each transform used but
00243  * twiddle factor and bit reversal tables can be reused.
00244  * \par
00245  * There is also an associated initialization function for each data type.
00246  * The initialization function performs the following operations:
00247  * - Sets the values of the internal structure fields.
00248  * - Initializes twiddle factor table and bit reversal table pointers.
00249  * - Initializes the internal complex FFT data structure.
00250  * \par
00251  * Use of the initialization function is optional.
00252  * However, if the initialization function is used, then the instance structure
00253  * cannot be placed into a const data section. To place an instance structure
00254  * into a const data section, the instance structure should be manually
00255  * initialized as follows:
00256  * <pre>
00257  *arm_rfft_instance_q31 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
00258  *arm_rfft_instance_q15 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
00259  * </pre>
00260  * where <code>fftLenReal</code> is the length of the real transform;
00261  * <code>fftLenBy2</code> length of  the internal complex transform.
00262  * <code>ifftFlagR</code> Selects forward (=0) or inverse (=1) transform.
00263  * <code>bitReverseFlagR</code> Selects bit reversed output (=0) or normal order
00264  * output (=1).
00265  * <code>twidCoefRModifier</code> stride modifier for the twiddle factor table.
00266  * The value is based on the FFT length;
00267  * <code>pTwiddleAReal</code>points to the A array of twiddle coefficients;
00268  * <code>pTwiddleBReal</code>points to the B array of twiddle coefficients;
00269  * <code>pCfft</code> points to the CFFT Instance structure. The CFFT structure
00270  * must also be initialized.  Refer to arm_cfft_radix4_f32() for details regarding
00271  * static initialization of the complex FFT instance structure.
00272  */
00273 
00274 /**
00275 * @addtogroup RealFFT
00276 * @{
00277 */
00278 
00279 /**
00280 * @brief Processing function for the floating-point real FFT.
00281 * @param[in]  *S              points to an arm_rfft_fast_instance_f32 structure.
00282 * @param[in]  *p              points to the input buffer.
00283 * @param[in]  *pOut           points to the output buffer.
00284 * @param[in]  ifftFlag        RFFT if flag is 0, RIFFT if flag is 1
00285 * @return none.
00286 */
00287 
00288 void arm_rfft_fast_f32(
00289 arm_rfft_fast_instance_f32 * S,
00290 float32_t * p, float32_t * pOut,
00291 uint8_t ifftFlag)
00292 {
00293    arm_cfft_instance_f32 * Sint = &(S->Sint);
00294    Sint->fftLen = S->fftLenRFFT / 2;
00295 
00296    /* Calculation of Real FFT */
00297    if (ifftFlag)
00298    {
00299       /*  Real FFT compression */
00300       merge_rfft_f32(S, p, pOut);
00301 
00302       /* Complex radix-4 IFFT process */
00303       arm_cfft_f32( Sint, pOut, ifftFlag, 1);
00304    }
00305    else
00306    {
00307       /* Calculation of RFFT of input */
00308       arm_cfft_f32( Sint, p, ifftFlag, 1);
00309 
00310       /*  Real FFT extraction */
00311       stage_rfft_f32(S, p, pOut);
00312    }
00313 }
00314 
00315 /**
00316 * @} end of RealFFT group
00317 */
00318