arm_cfft_radix8_f32.c

00001 /* ----------------------------------------------------------------------    
00002 * Copyright (C) 2010-2014 ARM Limited. All rights reserved.    
00003 *    
00004 * $Date:        19. March 2015 
00005 * $Revision:    V.1.4.5  
00006 *    
00007 * Project:      CMSIS DSP Library    
00008 * Title:        arm_cfft_radix8_f32.c    
00009 *    
00010 * Description:  Radix-8 Decimation in Frequency CFFT & CIFFT Floating point processing function        
00011 *    
00012 * Target Processor: Cortex-M4/Cortex-M3/Cortex-M0
00013 *  
00014 * Redistribution and use in source and binary forms, with or without 
00015 * modification, are permitted provided that the following conditions
00016 * are met:
00017 *   - Redistributions of source code must retain the above copyright
00018 *     notice, this list of conditions and the following disclaimer.
00019 *   - Redistributions in binary form must reproduce the above copyright
00020 *     notice, this list of conditions and the following disclaimer in
00021 *     the documentation and/or other materials provided with the 
00022 *     distribution.
00023 *   - Neither the name of ARM LIMITED nor the names of its contributors
00024 *     may be used to endorse or promote products derived from this
00025 *     software without specific prior written permission.
00026 *
00027 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00028 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00029 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
00030 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 
00031 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
00032 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
00033 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
00034 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
00035 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00036 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
00037 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
00038 * POSSIBILITY OF SUCH DAMAGE.      
00039 * -------------------------------------------------------------------- */
00040 
00041 #include "arm_math.h"
00042 
00043 /**    
00044 * @ingroup groupTransforms    
00045 */
00046 
00047 /**    
00048 * @defgroup Radix8_CFFT_CIFFT Radix-8 Complex FFT Functions    
00049 *    
00050 * \par    
00051 * Complex Fast Fourier Transform(CFFT) and Complex Inverse Fast Fourier Transform(CIFFT) is an efficient algorithm to compute Discrete Fourier Transform(DFT) and Inverse Discrete Fourier Transform(IDFT).    
00052 * Computational complexity of CFFT reduces drastically when compared to DFT.    
00053 * \par    
00054 * This set of functions implements CFFT/CIFFT    
00055 * for floating-point data types.  The functions operates on in-place buffer which uses same buffer for input and output.    
00056 * Complex input is stored in input buffer in an interleaved fashion.    
00057 *    
00058 * \par    
00059 * The functions operate on blocks of input and output data and each call to the function processes    
00060 * <code>2*fftLen</code> samples through the transform.  <code>pSrc</code>  points to In-place arrays containing <code>2*fftLen</code> values.    
00061 * \par   
00062 * The <code>pSrc</code> points to the array of in-place buffer of size <code>2*fftLen</code> and inputs and outputs are stored in an interleaved fashion as shown below.    
00063 * <pre> {real[0], imag[0], real[1], imag[1],..} </pre>    
00064 *    
00065 * \par Lengths supported by the transform:   
00066 * \par    
00067 * Internally, the function utilize a Radix-8 decimation in frequency(DIF) algorithm    
00068 * and the size of the FFT supported are of the lengths [ 64, 512, 4096].   
00069 *     
00070 *    
00071 * \par Algorithm:    
00072 *    
00073 * <b>Complex Fast Fourier Transform:</b>    
00074 * \par     
00075 * Input real and imaginary data:    
00076 * <pre>    
00077 * x(n) = xa + j * ya    
00078 * x(n+N/4 ) = xb + j * yb    
00079 * x(n+N/2 ) = xc + j * yc    
00080 * x(n+3N 4) = xd + j * yd    
00081 * </pre>    
00082 * where N is length of FFT    
00083 * \par    
00084 * Output real and imaginary data:    
00085 * <pre>    
00086 * X(4r) = xa'+ j * ya'    
00087 * X(4r+1) = xb'+ j * yb'    
00088 * X(4r+2) = xc'+ j * yc'    
00089 * X(4r+3) = xd'+ j * yd'    
00090 * </pre>    
00091 * \par    
00092 * Twiddle factors for Radix-8 FFT:    
00093 * <pre>    
00094 * Wn = co1 + j * (- si1)    
00095 * W2n = co2 + j * (- si2)    
00096 * W3n = co3 + j * (- si3)    
00097 * </pre>    
00098 *    
00099 * \par    
00100 * \image html CFFT.gif "Radix-8 Decimation-in Frequency Complex Fast Fourier Transform"    
00101 *    
00102 * \par    
00103 * Output from Radix-8 CFFT Results in Digit reversal order. Interchange middle two branches of every butterfly results in Bit reversed output.    
00104 * \par    
00105 * <b> Butterfly CFFT equations:</b>    
00106 * <pre>    
00107 * xa' = xa + xb + xc + xd    
00108 * ya' = ya + yb + yc + yd    
00109 * xc' = (xa+yb-xc-yd)* co1 + (ya-xb-yc+xd)* (si1)    
00110 * yc' = (ya-xb-yc+xd)* co1 - (xa+yb-xc-yd)* (si1)    
00111 * xb' = (xa-xb+xc-xd)* co2 + (ya-yb+yc-yd)* (si2)    
00112 * yb' = (ya-yb+yc-yd)* co2 - (xa-xb+xc-xd)* (si2)    
00113 * xd' = (xa-yb-xc+yd)* co3 + (ya+xb-yc-xd)* (si3)    
00114 * yd' = (ya+xb-yc-xd)* co3 - (xa-yb-xc+yd)* (si3)    
00115 * </pre>    
00116 *    
00117 * \par    
00118 * where <code>fftLen</code> length of CFFT/CIFFT; <code>ifftFlag</code> Flag for selection of CFFT or CIFFT(Set ifftFlag to calculate CIFFT otherwise calculates CFFT);    
00119 * <code>bitReverseFlag</code> Flag for selection of output order(Set bitReverseFlag to output in normal order otherwise output in bit reversed order);     
00120 * <code>pTwiddle</code>points to array of twiddle coefficients; <code>pBitRevTable</code> points to the array of bit reversal table.    
00121 * <code>twidCoefModifier</code> modifier for twiddle factor table which supports all FFT lengths with same table;     
00122 * <code>pBitRevTable</code> modifier for bit reversal table which supports all FFT lengths with same table.    
00123 * <code>onebyfftLen</code> value of 1/fftLen to calculate CIFFT;    
00124 *   
00125 * \par Fixed-Point Behavior    
00126 * Care must be taken when using the fixed-point versions of the CFFT/CIFFT function.    
00127 * Refer to the function specific documentation below for usage guidelines.    
00128 */
00129 
00130 
00131 /*    
00132 * @brief  Core function for the floating-point CFFT butterfly process.   
00133 * @param[in, out] *pSrc            points to the in-place buffer of floating-point data type.   
00134 * @param[in]      fftLen           length of the FFT.   
00135 * @param[in]      *pCoef           points to the twiddle coefficient buffer.   
00136 * @param[in]      twidCoefModifier twiddle coefficient modifier that supports different size FFTs with the same twiddle factor table.   
00137 * @return none.   
00138 */
00139 
00140 void arm_radix8_butterfly_f32(
00141 float32_t * pSrc,
00142 uint16_t fftLen,
00143 const float32_t * pCoef,
00144 uint16_t twidCoefModifier)
00145 {
00146    uint32_t ia1, ia2, ia3, ia4, ia5, ia6, ia7;
00147    uint32_t i1, i2, i3, i4, i5, i6, i7, i8;
00148    uint32_t id;
00149    uint32_t n1, n2, j;
00150    
00151    float32_t r1, r2, r3, r4, r5, r6, r7, r8;
00152    float32_t t1, t2;
00153    float32_t s1, s2, s3, s4, s5, s6, s7, s8;
00154    float32_t p1, p2, p3, p4;
00155    float32_t co2, co3, co4, co5, co6, co7, co8;
00156    float32_t si2, si3, si4, si5, si6, si7, si8;
00157    const float32_t C81 = 0.70710678118f;
00158 
00159    n2 = fftLen;
00160    
00161    do 
00162    {
00163       n1 = n2;
00164       n2 = n2 >> 3;
00165       i1 = 0;
00166       
00167       do
00168       {
00169          i2 = i1 + n2;
00170          i3 = i2 + n2;
00171          i4 = i3 + n2;
00172          i5 = i4 + n2;
00173          i6 = i5 + n2;
00174          i7 = i6 + n2;
00175          i8 = i7 + n2;
00176          r1 = pSrc[2 * i1] + pSrc[2 * i5];
00177          r5 = pSrc[2 * i1] - pSrc[2 * i5];
00178          r2 = pSrc[2 * i2] + pSrc[2 * i6];
00179          r6 = pSrc[2 * i2] - pSrc[2 * i6];
00180          r3 = pSrc[2 * i3] + pSrc[2 * i7];
00181          r7 = pSrc[2 * i3] - pSrc[2 * i7];
00182          r4 = pSrc[2 * i4] + pSrc[2 * i8];
00183          r8 = pSrc[2 * i4] - pSrc[2 * i8];
00184          t1 = r1 - r3;
00185          r1 = r1 + r3;
00186          r3 = r2 - r4;
00187          r2 = r2 + r4;
00188          pSrc[2 * i1] = r1 + r2;   
00189          pSrc[2 * i5] = r1 - r2;
00190          r1 = pSrc[2 * i1 + 1] + pSrc[2 * i5 + 1];
00191          s5 = pSrc[2 * i1 + 1] - pSrc[2 * i5 + 1];
00192          r2 = pSrc[2 * i2 + 1] + pSrc[2 * i6 + 1];
00193          s6 = pSrc[2 * i2 + 1] - pSrc[2 * i6 + 1];
00194          s3 = pSrc[2 * i3 + 1] + pSrc[2 * i7 + 1];
00195          s7 = pSrc[2 * i3 + 1] - pSrc[2 * i7 + 1];
00196          r4 = pSrc[2 * i4 + 1] + pSrc[2 * i8 + 1];
00197          s8 = pSrc[2 * i4 + 1] - pSrc[2 * i8 + 1];
00198          t2 = r1 - s3;
00199          r1 = r1 + s3;
00200          s3 = r2 - r4;
00201          r2 = r2 + r4;
00202          pSrc[2 * i1 + 1] = r1 + r2;
00203          pSrc[2 * i5 + 1] = r1 - r2;
00204          pSrc[2 * i3]     = t1 + s3;
00205          pSrc[2 * i7]     = t1 - s3;
00206          pSrc[2 * i3 + 1] = t2 - r3;
00207          pSrc[2 * i7 + 1] = t2 + r3;
00208          r1 = (r6 - r8) * C81;
00209          r6 = (r6 + r8) * C81;
00210          r2 = (s6 - s8) * C81;
00211          s6 = (s6 + s8) * C81;
00212          t1 = r5 - r1;
00213          r5 = r5 + r1;
00214          r8 = r7 - r6;
00215          r7 = r7 + r6;
00216          t2 = s5 - r2;
00217          s5 = s5 + r2;
00218          s8 = s7 - s6;
00219          s7 = s7 + s6;
00220          pSrc[2 * i2]     = r5 + s7;
00221          pSrc[2 * i8]     = r5 - s7;
00222          pSrc[2 * i6]     = t1 + s8;
00223          pSrc[2 * i4]     = t1 - s8;
00224          pSrc[2 * i2 + 1] = s5 - r7;
00225          pSrc[2 * i8 + 1] = s5 + r7;
00226          pSrc[2 * i6 + 1] = t2 - r8;
00227          pSrc[2 * i4 + 1] = t2 + r8;
00228          
00229          i1 += n1;
00230       } while(i1 < fftLen);
00231       
00232       if(n2 < 8)
00233          break;
00234       
00235       ia1 = 0;
00236       j = 1;
00237       
00238       do
00239       {      
00240          /*  index calculation for the coefficients */
00241          id  = ia1 + twidCoefModifier;
00242          ia1 = id;
00243          ia2 = ia1 + id;
00244          ia3 = ia2 + id;
00245          ia4 = ia3 + id;
00246          ia5 = ia4 + id;
00247          ia6 = ia5 + id;
00248          ia7 = ia6 + id;
00249                   
00250          co2 = pCoef[2 * ia1];
00251          co3 = pCoef[2 * ia2];
00252          co4 = pCoef[2 * ia3];
00253          co5 = pCoef[2 * ia4];
00254          co6 = pCoef[2 * ia5];
00255          co7 = pCoef[2 * ia6];
00256          co8 = pCoef[2 * ia7];
00257          si2 = pCoef[2 * ia1 + 1];
00258          si3 = pCoef[2 * ia2 + 1];
00259          si4 = pCoef[2 * ia3 + 1];
00260          si5 = pCoef[2 * ia4 + 1];
00261          si6 = pCoef[2 * ia5 + 1];
00262          si7 = pCoef[2 * ia6 + 1];
00263          si8 = pCoef[2 * ia7 + 1];         
00264          
00265          i1 = j;
00266          
00267          do
00268          {
00269             /*  index calculation for the input */
00270             i2 = i1 + n2;
00271             i3 = i2 + n2;
00272             i4 = i3 + n2;
00273             i5 = i4 + n2;
00274             i6 = i5 + n2;
00275             i7 = i6 + n2;
00276             i8 = i7 + n2;
00277             r1 = pSrc[2 * i1] + pSrc[2 * i5];
00278             r5 = pSrc[2 * i1] - pSrc[2 * i5];
00279             r2 = pSrc[2 * i2] + pSrc[2 * i6];
00280             r6 = pSrc[2 * i2] - pSrc[2 * i6];
00281             r3 = pSrc[2 * i3] + pSrc[2 * i7];
00282             r7 = pSrc[2 * i3] - pSrc[2 * i7];
00283             r4 = pSrc[2 * i4] + pSrc[2 * i8];
00284             r8 = pSrc[2 * i4] - pSrc[2 * i8];
00285             t1 = r1 - r3;
00286             r1 = r1 + r3;
00287             r3 = r2 - r4;
00288             r2 = r2 + r4;
00289             pSrc[2 * i1] = r1 + r2;
00290             r2 = r1 - r2;
00291             s1 = pSrc[2 * i1 + 1] + pSrc[2 * i5 + 1];
00292             s5 = pSrc[2 * i1 + 1] - pSrc[2 * i5 + 1];
00293             s2 = pSrc[2 * i2 + 1] + pSrc[2 * i6 + 1];
00294             s6 = pSrc[2 * i2 + 1] - pSrc[2 * i6 + 1];
00295             s3 = pSrc[2 * i3 + 1] + pSrc[2 * i7 + 1];
00296             s7 = pSrc[2 * i3 + 1] - pSrc[2 * i7 + 1];
00297             s4 = pSrc[2 * i4 + 1] + pSrc[2 * i8 + 1];
00298             s8 = pSrc[2 * i4 + 1] - pSrc[2 * i8 + 1];
00299             t2 = s1 - s3;
00300             s1 = s1 + s3;
00301             s3 = s2 - s4;
00302             s2 = s2 + s4;
00303             r1 = t1 + s3;
00304             t1 = t1 - s3;
00305             pSrc[2 * i1 + 1] = s1 + s2;
00306             s2 = s1 - s2;
00307             s1 = t2 - r3;
00308             t2 = t2 + r3;
00309             p1 = co5 * r2;
00310             p2 = si5 * s2;
00311             p3 = co5 * s2;
00312             p4 = si5 * r2;
00313             pSrc[2 * i5]     = p1 + p2;
00314             pSrc[2 * i5 + 1] = p3 - p4;
00315             p1 = co3 * r1;
00316             p2 = si3 * s1;
00317             p3 = co3 * s1;
00318             p4 = si3 * r1;
00319             pSrc[2 * i3]     = p1 + p2;
00320             pSrc[2 * i3 + 1] = p3 - p4;
00321             p1 = co7 * t1;
00322             p2 = si7 * t2;
00323             p3 = co7 * t2;
00324             p4 = si7 * t1;
00325             pSrc[2 * i7]     = p1 + p2;
00326             pSrc[2 * i7 + 1] = p3 - p4;
00327             r1 = (r6 - r8) * C81;
00328             r6 = (r6 + r8) * C81;
00329             s1 = (s6 - s8) * C81;
00330             s6 = (s6 + s8) * C81;
00331             t1 = r5 - r1;
00332             r5 = r5 + r1;
00333             r8 = r7 - r6;
00334             r7 = r7 + r6;
00335             t2 = s5 - s1;
00336             s5 = s5 + s1;
00337             s8 = s7 - s6;
00338             s7 = s7 + s6;
00339             r1 = r5 + s7;
00340             r5 = r5 - s7;
00341             r6 = t1 + s8;
00342             t1 = t1 - s8;
00343             s1 = s5 - r7;
00344             s5 = s5 + r7;
00345             s6 = t2 - r8;
00346             t2 = t2 + r8;
00347             p1 = co2 * r1;
00348             p2 = si2 * s1;
00349             p3 = co2 * s1;
00350             p4 = si2 * r1;
00351             pSrc[2 * i2]     = p1 + p2;
00352             pSrc[2 * i2 + 1] = p3 - p4;
00353             p1 = co8 * r5;
00354             p2 = si8 * s5;
00355             p3 = co8 * s5;
00356             p4 = si8 * r5;
00357             pSrc[2 * i8]     = p1 + p2;
00358             pSrc[2 * i8 + 1] = p3 - p4;
00359             p1 = co6 * r6;
00360             p2 = si6 * s6;
00361             p3 = co6 * s6;
00362             p4 = si6 * r6;
00363             pSrc[2 * i6]     = p1 + p2;
00364             pSrc[2 * i6 + 1] = p3 - p4;
00365             p1 = co4 * t1;
00366             p2 = si4 * t2;
00367             p3 = co4 * t2;
00368             p4 = si4 * t1;
00369             pSrc[2 * i4]     = p1 + p2;
00370             pSrc[2 * i4 + 1] = p3 - p4;
00371             
00372             i1 += n1;
00373          } while(i1 < fftLen);
00374          
00375          j++;
00376       } while(j < n2);
00377       
00378       twidCoefModifier <<= 3;
00379    } while(n2 > 7);   
00380 }
00381 
00382 /**    
00383 * @} end of Radix8_CFFT_CIFFT group    
00384 */