Diff: TransformFunctions/arm_cfft_radix8_f32.c

Revision:

0:3d9c67d97d6f

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/TransformFunctions/arm_cfft_radix8_f32.c	Mon Jul 28 15:03:15 2014 +0000
@@ -0,0 +1,384 @@
+/* ----------------------------------------------------------------------    
+* Copyright (C) 2010-2014 ARM Limited. All rights reserved.    
+*    
+* $Date:        12. March 2014  
+* $Revision: 	V1.4.3  
+*    
+* Project: 	    CMSIS DSP Library    
+* Title:	    arm_cfft_radix8_f32.c    
+*    
+* Description:	Radix-8 Decimation in Frequency CFFT & CIFFT Floating point processing 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"
+
+/**    
+* @ingroup groupTransforms    
+*/
+
+/**    
+* @defgroup Radix8_CFFT_CIFFT Radix-8 Complex FFT Functions    
+*    
+* \par    
+* 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).    
+* Computational complexity of CFFT reduces drastically when compared to DFT.    
+* \par    
+* This set of functions implements CFFT/CIFFT    
+* for floating-point data types.  The functions operates on in-place buffer which uses same buffer for input and output.    
+* Complex input is stored in input buffer in an interleaved fashion.    
+*    
+* \par    
+* The functions operate on blocks of input and output data and each call to the function processes    
+* <code>2*fftLen</code> samples through the transform.  <code>pSrc</code>  points to In-place arrays containing <code>2*fftLen</code> values.    
+* \par   
+* 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.    
+* <pre> {real[0], imag[0], real[1], imag[1],..} </pre>    
+*    
+* \par Lengths supported by the transform:   
+* \par    
+* Internally, the function utilize a Radix-8 decimation in frequency(DIF) algorithm    
+* and the size of the FFT supported are of the lengths [ 64, 512, 4096].   
+*     
+*    
+* \par Algorithm:    
+*    
+* <b>Complex Fast Fourier Transform:</b>    
+* \par     
+* Input real and imaginary data:    
+* <pre>    
+* x(n) = xa + j * ya    
+* x(n+N/4 ) = xb + j * yb    
+* x(n+N/2 ) = xc + j * yc    
+* x(n+3N 4) = xd + j * yd    
+* </pre>    
+* where N is length of FFT    
+* \par    
+* Output real and imaginary data:    
+* <pre>    
+* X(4r) = xa'+ j * ya'    
+* X(4r+1) = xb'+ j * yb'    
+* X(4r+2) = xc'+ j * yc'    
+* X(4r+3) = xd'+ j * yd'    
+* </pre>    
+* \par    
+* Twiddle factors for Radix-8 FFT:    
+* <pre>    
+* Wn = co1 + j * (- si1)    
+* W2n = co2 + j * (- si2)    
+* W3n = co3 + j * (- si3)    
+* </pre>    
+*    
+* \par    
+* \image html CFFT.gif "Radix-8 Decimation-in Frequency Complex Fast Fourier Transform"    
+*    
+* \par    
+* Output from Radix-8 CFFT Results in Digit reversal order. Interchange middle two branches of every butterfly results in Bit reversed output.    
+* \par    
+* <b> Butterfly CFFT equations:</b>    
+* <pre>    
+* xa' = xa + xb + xc + xd    
+* ya' = ya + yb + yc + yd    
+* xc' = (xa+yb-xc-yd)* co1 + (ya-xb-yc+xd)* (si1)    
+* yc' = (ya-xb-yc+xd)* co1 - (xa+yb-xc-yd)* (si1)    
+* xb' = (xa-xb+xc-xd)* co2 + (ya-yb+yc-yd)* (si2)    
+* yb' = (ya-yb+yc-yd)* co2 - (xa-xb+xc-xd)* (si2)    
+* xd' = (xa-yb-xc+yd)* co3 + (ya+xb-yc-xd)* (si3)    
+* yd' = (ya+xb-yc-xd)* co3 - (xa-yb-xc+yd)* (si3)    
+* </pre>    
+*    
+* \par    
+* 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);    
+* <code>bitReverseFlag</code> Flag for selection of output order(Set bitReverseFlag to output in normal order otherwise output in bit reversed order);     
+* <code>pTwiddle</code>points to array of twiddle coefficients; <code>pBitRevTable</code> points to the array of bit reversal table.    
+* <code>twidCoefModifier</code> modifier for twiddle factor table which supports all FFT lengths with same table;     
+* <code>pBitRevTable</code> modifier for bit reversal table which supports all FFT lengths with same table.    
+* <code>onebyfftLen</code> value of 1/fftLen to calculate CIFFT;    
+*   
+* \par Fixed-Point Behavior    
+* Care must be taken when using the fixed-point versions of the CFFT/CIFFT function.    
+* Refer to the function specific documentation below for usage guidelines.    
+*/
+
+
+/*    
+* @brief  Core function for the floating-point CFFT butterfly process.   
+* @param[in, out] *pSrc            points to the in-place buffer of floating-point data type.   
+* @param[in]      fftLen           length of the FFT.   
+* @param[in]      *pCoef           points to the twiddle coefficient buffer.   
+* @param[in]      twidCoefModifier twiddle coefficient modifier that supports different size FFTs with the same twiddle factor table.   
+* @return none.   
+*/
+
+void arm_radix8_butterfly_f32(
+float32_t * pSrc,
+uint16_t fftLen,
+const float32_t * pCoef,
+uint16_t twidCoefModifier)
+{
+   uint32_t ia1, ia2, ia3, ia4, ia5, ia6, ia7;
+   uint32_t i1, i2, i3, i4, i5, i6, i7, i8;
+   uint32_t id;
+   uint32_t n1, n2, j;
+   
+   float32_t r1, r2, r3, r4, r5, r6, r7, r8;
+   float32_t t1, t2;
+   float32_t s1, s2, s3, s4, s5, s6, s7, s8;
+   float32_t p1, p2, p3, p4;
+   float32_t co2, co3, co4, co5, co6, co7, co8;
+   float32_t si2, si3, si4, si5, si6, si7, si8;
+   const float32_t C81 = 0.70710678118f;
+
+   n2 = fftLen;
+   
+   do 
+   {
+      n1 = n2;
+      n2 = n2 >> 3;
+      i1 = 0;
+      
+      do
+      {
+         i2 = i1 + n2;
+         i3 = i2 + n2;
+         i4 = i3 + n2;
+         i5 = i4 + n2;
+         i6 = i5 + n2;
+         i7 = i6 + n2;
+         i8 = i7 + n2;
+         r1 = pSrc[2 * i1] + pSrc[2 * i5];
+         r5 = pSrc[2 * i1] - pSrc[2 * i5];
+         r2 = pSrc[2 * i2] + pSrc[2 * i6];
+         r6 = pSrc[2 * i2] - pSrc[2 * i6];
+         r3 = pSrc[2 * i3] + pSrc[2 * i7];
+         r7 = pSrc[2 * i3] - pSrc[2 * i7];
+         r4 = pSrc[2 * i4] + pSrc[2 * i8];
+         r8 = pSrc[2 * i4] - pSrc[2 * i8];
+         t1 = r1 - r3;
+         r1 = r1 + r3;
+         r3 = r2 - r4;
+         r2 = r2 + r4;
+         pSrc[2 * i1] = r1 + r2;   
+         pSrc[2 * i5] = r1 - r2;
+         r1 = pSrc[2 * i1 + 1] + pSrc[2 * i5 + 1];
+         s5 = pSrc[2 * i1 + 1] - pSrc[2 * i5 + 1];
+         r2 = pSrc[2 * i2 + 1] + pSrc[2 * i6 + 1];
+         s6 = pSrc[2 * i2 + 1] - pSrc[2 * i6 + 1];
+         s3 = pSrc[2 * i3 + 1] + pSrc[2 * i7 + 1];
+         s7 = pSrc[2 * i3 + 1] - pSrc[2 * i7 + 1];
+         r4 = pSrc[2 * i4 + 1] + pSrc[2 * i8 + 1];
+         s8 = pSrc[2 * i4 + 1] - pSrc[2 * i8 + 1];
+         t2 = r1 - s3;
+         r1 = r1 + s3;
+         s3 = r2 - r4;
+         r2 = r2 + r4;
+         pSrc[2 * i1 + 1] = r1 + r2;
+         pSrc[2 * i5 + 1] = r1 - r2;
+         pSrc[2 * i3]     = t1 + s3;
+         pSrc[2 * i7]     = t1 - s3;
+         pSrc[2 * i3 + 1] = t2 - r3;
+         pSrc[2 * i7 + 1] = t2 + r3;
+         r1 = (r6 - r8) * C81;
+         r6 = (r6 + r8) * C81;
+         r2 = (s6 - s8) * C81;
+         s6 = (s6 + s8) * C81;
+         t1 = r5 - r1;
+         r5 = r5 + r1;
+         r8 = r7 - r6;
+         r7 = r7 + r6;
+         t2 = s5 - r2;
+         s5 = s5 + r2;
+         s8 = s7 - s6;
+         s7 = s7 + s6;
+         pSrc[2 * i2]     = r5 + s7;
+         pSrc[2 * i8]     = r5 - s7;
+         pSrc[2 * i6]     = t1 + s8;
+         pSrc[2 * i4]     = t1 - s8;
+         pSrc[2 * i2 + 1] = s5 - r7;
+         pSrc[2 * i8 + 1] = s5 + r7;
+         pSrc[2 * i6 + 1] = t2 - r8;
+         pSrc[2 * i4 + 1] = t2 + r8;
+         
+         i1 += n1;
+      } while(i1 < fftLen);
+      
+      if(n2 < 8)
+         break;
+      
+      ia1 = 0;
+      j = 1;
+      
+      do
+      {      
+         /*  index calculation for the coefficients */
+         id  = ia1 + twidCoefModifier;
+         ia1 = id;
+         ia2 = ia1 + id;
+         ia3 = ia2 + id;
+         ia4 = ia3 + id;
+         ia5 = ia4 + id;
+         ia6 = ia5 + id;
+         ia7 = ia6 + id;
+                  
+         co2 = pCoef[2 * ia1];
+         co3 = pCoef[2 * ia2];
+         co4 = pCoef[2 * ia3];
+         co5 = pCoef[2 * ia4];
+         co6 = pCoef[2 * ia5];
+         co7 = pCoef[2 * ia6];
+         co8 = pCoef[2 * ia7];
+         si2 = pCoef[2 * ia1 + 1];
+         si3 = pCoef[2 * ia2 + 1];
+         si4 = pCoef[2 * ia3 + 1];
+         si5 = pCoef[2 * ia4 + 1];
+         si6 = pCoef[2 * ia5 + 1];
+         si7 = pCoef[2 * ia6 + 1];
+         si8 = pCoef[2 * ia7 + 1];         
+         
+         i1 = j;
+         
+         do
+         {
+            /*  index calculation for the input */
+            i2 = i1 + n2;
+            i3 = i2 + n2;
+            i4 = i3 + n2;
+            i5 = i4 + n2;
+            i6 = i5 + n2;
+            i7 = i6 + n2;
+            i8 = i7 + n2;
+            r1 = pSrc[2 * i1] + pSrc[2 * i5];
+            r5 = pSrc[2 * i1] - pSrc[2 * i5];
+            r2 = pSrc[2 * i2] + pSrc[2 * i6];
+            r6 = pSrc[2 * i2] - pSrc[2 * i6];
+            r3 = pSrc[2 * i3] + pSrc[2 * i7];
+            r7 = pSrc[2 * i3] - pSrc[2 * i7];
+            r4 = pSrc[2 * i4] + pSrc[2 * i8];
+            r8 = pSrc[2 * i4] - pSrc[2 * i8];
+            t1 = r1 - r3;
+            r1 = r1 + r3;
+            r3 = r2 - r4;
+            r2 = r2 + r4;
+            pSrc[2 * i1] = r1 + r2;
+            r2 = r1 - r2;
+            s1 = pSrc[2 * i1 + 1] + pSrc[2 * i5 + 1];
+            s5 = pSrc[2 * i1 + 1] - pSrc[2 * i5 + 1];
+            s2 = pSrc[2 * i2 + 1] + pSrc[2 * i6 + 1];
+            s6 = pSrc[2 * i2 + 1] - pSrc[2 * i6 + 1];
+            s3 = pSrc[2 * i3 + 1] + pSrc[2 * i7 + 1];
+            s7 = pSrc[2 * i3 + 1] - pSrc[2 * i7 + 1];
+            s4 = pSrc[2 * i4 + 1] + pSrc[2 * i8 + 1];
+            s8 = pSrc[2 * i4 + 1] - pSrc[2 * i8 + 1];
+            t2 = s1 - s3;
+            s1 = s1 + s3;
+            s3 = s2 - s4;
+            s2 = s2 + s4;
+            r1 = t1 + s3;
+            t1 = t1 - s3;
+            pSrc[2 * i1 + 1] = s1 + s2;
+            s2 = s1 - s2;
+            s1 = t2 - r3;
+            t2 = t2 + r3;
+            p1 = co5 * r2;
+            p2 = si5 * s2;
+            p3 = co5 * s2;
+            p4 = si5 * r2;
+            pSrc[2 * i5]     = p1 + p2;
+            pSrc[2 * i5 + 1] = p3 - p4;
+            p1 = co3 * r1;
+            p2 = si3 * s1;
+            p3 = co3 * s1;
+            p4 = si3 * r1;
+            pSrc[2 * i3]     = p1 + p2;
+            pSrc[2 * i3 + 1] = p3 - p4;
+            p1 = co7 * t1;
+            p2 = si7 * t2;
+            p3 = co7 * t2;
+            p4 = si7 * t1;
+            pSrc[2 * i7]     = p1 + p2;
+            pSrc[2 * i7 + 1] = p3 - p4;
+            r1 = (r6 - r8) * C81;
+            r6 = (r6 + r8) * C81;
+            s1 = (s6 - s8) * C81;
+            s6 = (s6 + s8) * C81;
+            t1 = r5 - r1;
+            r5 = r5 + r1;
+            r8 = r7 - r6;
+            r7 = r7 + r6;
+            t2 = s5 - s1;
+            s5 = s5 + s1;
+            s8 = s7 - s6;
+            s7 = s7 + s6;
+            r1 = r5 + s7;
+            r5 = r5 - s7;
+            r6 = t1 + s8;
+            t1 = t1 - s8;
+            s1 = s5 - r7;
+            s5 = s5 + r7;
+            s6 = t2 - r8;
+            t2 = t2 + r8;
+            p1 = co2 * r1;
+            p2 = si2 * s1;
+            p3 = co2 * s1;
+            p4 = si2 * r1;
+            pSrc[2 * i2]     = p1 + p2;
+            pSrc[2 * i2 + 1] = p3 - p4;
+            p1 = co8 * r5;
+            p2 = si8 * s5;
+            p3 = co8 * s5;
+            p4 = si8 * r5;
+            pSrc[2 * i8]     = p1 + p2;
+            pSrc[2 * i8 + 1] = p3 - p4;
+            p1 = co6 * r6;
+            p2 = si6 * s6;
+            p3 = co6 * s6;
+            p4 = si6 * r6;
+            pSrc[2 * i6]     = p1 + p2;
+            pSrc[2 * i6 + 1] = p3 - p4;
+            p1 = co4 * t1;
+            p2 = si4 * t2;
+            p3 = co4 * t2;
+            p4 = si4 * t1;
+            pSrc[2 * i4]     = p1 + p2;
+            pSrc[2 * i4 + 1] = p3 - p4;
+            
+            i1 += n1;
+         } while(i1 < fftLen);
+         
+         j++;
+      } while(j < n2);
+      
+      twidCoefModifier <<= 3;
+   } while(n2 > 7);   
+}
+
+/**    
+* @} end of Radix8_CFFT_CIFFT group    
+*/