Joep D / CMSIS_DSP_5

The CMSIS DSP 5 library

Dependents:   Nucleo-Heart-Rate ejercicioVrms2 PROYECTOFINAL ejercicioVrms ... more

functions/TransformFunctions/arm_cfft_f32.c@3:4098b9d3d571, 2018-06-21 (annotated)

Committer:
xorjoep
Date:
Thu Jun 21 11:56:27 2018 +0000
Revision:
3:4098b9d3d571
Parent:
1:24714b45cd1b 
headers is a folder not a library

Who changed what in which revision?

UserRevisionLine numberNew contents of line
xorjoep 1:24714b45cd1b 1 /* ----------------------------------------------------------------------
xorjoep 1:24714b45cd1b 2 * Project: CMSIS DSP Library
xorjoep 1:24714b45cd1b 3 * Title: arm_cfft_f32.c
xorjoep 1:24714b45cd1b 4 * Description: Combined Radix Decimation in Frequency CFFT Floating point processing function
xorjoep 1:24714b45cd1b 5 *
xorjoep 1:24714b45cd1b 6 * $Date: 27. January 2017
xorjoep 1:24714b45cd1b 7 * $Revision: V.1.5.1
xorjoep 1:24714b45cd1b 8 *
xorjoep 1:24714b45cd1b 9 * Target Processor: Cortex-M cores
xorjoep 1:24714b45cd1b 10 * -------------------------------------------------------------------- */
xorjoep 1:24714b45cd1b 11 /*
xorjoep 1:24714b45cd1b 12 * Copyright (C) 2010-2017 ARM Limited or its affiliates. All rights reserved.
xorjoep 1:24714b45cd1b 13 *
xorjoep 1:24714b45cd1b 14 * SPDX-License-Identifier: Apache-2.0
xorjoep 1:24714b45cd1b 15 *
xorjoep 1:24714b45cd1b 16 * Licensed under the Apache License, Version 2.0 (the License); you may
xorjoep 1:24714b45cd1b 17 * not use this file except in compliance with the License.
xorjoep 1:24714b45cd1b 18 * You may obtain a copy of the License at
xorjoep 1:24714b45cd1b 19 *
xorjoep 1:24714b45cd1b 20 * www.apache.org/licenses/LICENSE-2.0
xorjoep 1:24714b45cd1b 21 *
xorjoep 1:24714b45cd1b 22 * Unless required by applicable law or agreed to in writing, software
xorjoep 1:24714b45cd1b 23 * distributed under the License is distributed on an AS IS BASIS, WITHOUT
xorjoep 1:24714b45cd1b 24 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
xorjoep 1:24714b45cd1b 25 * See the License for the specific language governing permissions and
xorjoep 1:24714b45cd1b 26 * limitations under the License.
xorjoep 1:24714b45cd1b 27 */
xorjoep 1:24714b45cd1b 28
xorjoep 1:24714b45cd1b 29 #include "arm_math.h"
xorjoep 1:24714b45cd1b 30 #include "arm_common_tables.h"
xorjoep 1:24714b45cd1b 31
xorjoep 1:24714b45cd1b 32 extern void arm_radix8_butterfly_f32(
xorjoep 1:24714b45cd1b 33 float32_t * pSrc,
xorjoep 1:24714b45cd1b 34 uint16_t fftLen,
xorjoep 1:24714b45cd1b 35 const float32_t * pCoef,
xorjoep 1:24714b45cd1b 36 uint16_t twidCoefModifier);
xorjoep 1:24714b45cd1b 37
xorjoep 1:24714b45cd1b 38 extern void arm_bitreversal_32(
xorjoep 1:24714b45cd1b 39 uint32_t * pSrc,
xorjoep 1:24714b45cd1b 40 const uint16_t bitRevLen,
xorjoep 1:24714b45cd1b 41 const uint16_t * pBitRevTable);
xorjoep 1:24714b45cd1b 42
xorjoep 1:24714b45cd1b 43 /**
xorjoep 1:24714b45cd1b 44 * @ingroup groupTransforms
xorjoep 1:24714b45cd1b 45 */
xorjoep 1:24714b45cd1b 46
xorjoep 1:24714b45cd1b 47 /**
xorjoep 1:24714b45cd1b 48 * @defgroup ComplexFFT Complex FFT Functions
xorjoep 1:24714b45cd1b 49 *
xorjoep 1:24714b45cd1b 50 * \par
xorjoep 1:24714b45cd1b 51 * The Fast Fourier Transform (FFT) is an efficient algorithm for computing the
xorjoep 1:24714b45cd1b 52 * Discrete Fourier Transform (DFT). The FFT can be orders of magnitude faster
xorjoep 1:24714b45cd1b 53 * than the DFT, especially for long lengths.
xorjoep 1:24714b45cd1b 54 * The algorithms described in this section
xorjoep 1:24714b45cd1b 55 * operate on complex data. A separate set of functions is devoted to handling
xorjoep 1:24714b45cd1b 56 * of real sequences.
xorjoep 1:24714b45cd1b 57 * \par
xorjoep 1:24714b45cd1b 58 * There are separate algorithms for handling floating-point, Q15, and Q31 data
xorjoep 1:24714b45cd1b 59 * types. The algorithms available for each data type are described next.
xorjoep 1:24714b45cd1b 60 * \par
xorjoep 1:24714b45cd1b 61 * The FFT functions operate in-place. That is, the array holding the input data
xorjoep 1:24714b45cd1b 62 * will also be used to hold the corresponding result. The input data is complex
xorjoep 1:24714b45cd1b 63 * and contains <code>2*fftLen</code> interleaved values as shown below.
xorjoep 1:24714b45cd1b 64 * <pre> {real[0], imag[0], real[1], imag[1],..} </pre>
xorjoep 1:24714b45cd1b 65 * The FFT result will be contained in the same array and the frequency domain
xorjoep 1:24714b45cd1b 66 * values will have the same interleaving.
xorjoep 1:24714b45cd1b 67 *
xorjoep 1:24714b45cd1b 68 * \par Floating-point
xorjoep 1:24714b45cd1b 69 * The floating-point complex FFT uses a mixed-radix algorithm. Multiple radix-8
xorjoep 1:24714b45cd1b 70 * stages are performed along with a single radix-2 or radix-4 stage, as needed.
xorjoep 1:24714b45cd1b 71 * The algorithm supports lengths of [16, 32, 64, ..., 4096] and each length uses
xorjoep 1:24714b45cd1b 72 * a different twiddle factor table.
xorjoep 1:24714b45cd1b 73 * \par
xorjoep 1:24714b45cd1b 74 * The function uses the standard FFT definition and output values may grow by a
xorjoep 1:24714b45cd1b 75 * factor of <code>fftLen</code> when computing the forward transform. The
xorjoep 1:24714b45cd1b 76 * inverse transform includes a scale of <code>1/fftLen</code> as part of the
xorjoep 1:24714b45cd1b 77 * calculation and this matches the textbook definition of the inverse FFT.
xorjoep 1:24714b45cd1b 78 * \par
xorjoep 1:24714b45cd1b 79 * Pre-initialized data structures containing twiddle factors and bit reversal
xorjoep 1:24714b45cd1b 80 * tables are provided and defined in <code>arm_const_structs.h</code>. Include
xorjoep 1:24714b45cd1b 81 * this header in your function and then pass one of the constant structures as
xorjoep 1:24714b45cd1b 82 * an argument to arm_cfft_f32. For example:
xorjoep 1:24714b45cd1b 83 * \par
xorjoep 1:24714b45cd1b 84 * <code>arm_cfft_f32(arm_cfft_sR_f32_len64, pSrc, 1, 1)</code>
xorjoep 1:24714b45cd1b 85 * \par
xorjoep 1:24714b45cd1b 86 * computes a 64-point inverse complex FFT including bit reversal.
xorjoep 1:24714b45cd1b 87 * The data structures are treated as constant data and not modified during the
xorjoep 1:24714b45cd1b 88 * calculation. The same data structure can be reused for multiple transforms
xorjoep 1:24714b45cd1b 89 * including mixing forward and inverse transforms.
xorjoep 1:24714b45cd1b 90 * \par
xorjoep 1:24714b45cd1b 91 * Earlier releases of the library provided separate radix-2 and radix-4
xorjoep 1:24714b45cd1b 92 * algorithms that operated on floating-point data. These functions are still
xorjoep 1:24714b45cd1b 93 * provided but are deprecated. The older functions are slower and less general
xorjoep 1:24714b45cd1b 94 * than the new functions.
xorjoep 1:24714b45cd1b 95 * \par
xorjoep 1:24714b45cd1b 96 * An example of initialization of the constants for the arm_cfft_f32 function follows:
xorjoep 1:24714b45cd1b 97 * \code
xorjoep 1:24714b45cd1b 98 * const static arm_cfft_instance_f32 *S;
xorjoep 1:24714b45cd1b 99 * ...
xorjoep 1:24714b45cd1b 100 * switch (length) {
xorjoep 1:24714b45cd1b 101 * case 16:
xorjoep 1:24714b45cd1b 102 * S = &arm_cfft_sR_f32_len16;
xorjoep 1:24714b45cd1b 103 * break;
xorjoep 1:24714b45cd1b 104 * case 32:
xorjoep 1:24714b45cd1b 105 * S = &arm_cfft_sR_f32_len32;
xorjoep 1:24714b45cd1b 106 * break;
xorjoep 1:24714b45cd1b 107 * case 64:
xorjoep 1:24714b45cd1b 108 * S = &arm_cfft_sR_f32_len64;
xorjoep 1:24714b45cd1b 109 * break;
xorjoep 1:24714b45cd1b 110 * case 128:
xorjoep 1:24714b45cd1b 111 * S = &arm_cfft_sR_f32_len128;
xorjoep 1:24714b45cd1b 112 * break;
xorjoep 1:24714b45cd1b 113 * case 256:
xorjoep 1:24714b45cd1b 114 * S = &arm_cfft_sR_f32_len256;
xorjoep 1:24714b45cd1b 115 * break;
xorjoep 1:24714b45cd1b 116 * case 512:
xorjoep 1:24714b45cd1b 117 * S = &arm_cfft_sR_f32_len512;
xorjoep 1:24714b45cd1b 118 * break;
xorjoep 1:24714b45cd1b 119 * case 1024:
xorjoep 1:24714b45cd1b 120 * S = &arm_cfft_sR_f32_len1024;
xorjoep 1:24714b45cd1b 121 * break;
xorjoep 1:24714b45cd1b 122 * case 2048:
xorjoep 1:24714b45cd1b 123 * S = &arm_cfft_sR_f32_len2048;
xorjoep 1:24714b45cd1b 124 * break;
xorjoep 1:24714b45cd1b 125 * case 4096:
xorjoep 1:24714b45cd1b 126 * S = &arm_cfft_sR_f32_len4096;
xorjoep 1:24714b45cd1b 127 * break;
xorjoep 1:24714b45cd1b 128 * }
xorjoep 1:24714b45cd1b 129 * \endcode
xorjoep 1:24714b45cd1b 130 * \par Q15 and Q31
xorjoep 1:24714b45cd1b 131 * The floating-point complex FFT uses a mixed-radix algorithm. Multiple radix-4
xorjoep 1:24714b45cd1b 132 * stages are performed along with a single radix-2 stage, as needed.
xorjoep 1:24714b45cd1b 133 * The algorithm supports lengths of [16, 32, 64, ..., 4096] and each length uses
xorjoep 1:24714b45cd1b 134 * a different twiddle factor table.
xorjoep 1:24714b45cd1b 135 * \par
xorjoep 1:24714b45cd1b 136 * The function uses the standard FFT definition and output values may grow by a
xorjoep 1:24714b45cd1b 137 * factor of <code>fftLen</code> when computing the forward transform. The
xorjoep 1:24714b45cd1b 138 * inverse transform includes a scale of <code>1/fftLen</code> as part of the
xorjoep 1:24714b45cd1b 139 * calculation and this matches the textbook definition of the inverse FFT.
xorjoep 1:24714b45cd1b 140 * \par
xorjoep 1:24714b45cd1b 141 * Pre-initialized data structures containing twiddle factors and bit reversal
xorjoep 1:24714b45cd1b 142 * tables are provided and defined in <code>arm_const_structs.h</code>. Include
xorjoep 1:24714b45cd1b 143 * this header in your function and then pass one of the constant structures as
xorjoep 1:24714b45cd1b 144 * an argument to arm_cfft_q31. For example:
xorjoep 1:24714b45cd1b 145 * \par
xorjoep 1:24714b45cd1b 146 * <code>arm_cfft_q31(arm_cfft_sR_q31_len64, pSrc, 1, 1)</code>
xorjoep 1:24714b45cd1b 147 * \par
xorjoep 1:24714b45cd1b 148 * computes a 64-point inverse complex FFT including bit reversal.
xorjoep 1:24714b45cd1b 149 * The data structures are treated as constant data and not modified during the
xorjoep 1:24714b45cd1b 150 * calculation. The same data structure can be reused for multiple transforms
xorjoep 1:24714b45cd1b 151 * including mixing forward and inverse transforms.
xorjoep 1:24714b45cd1b 152 * \par
xorjoep 1:24714b45cd1b 153 * Earlier releases of the library provided separate radix-2 and radix-4
xorjoep 1:24714b45cd1b 154 * algorithms that operated on floating-point data. These functions are still
xorjoep 1:24714b45cd1b 155 * provided but are deprecated. The older functions are slower and less general
xorjoep 1:24714b45cd1b 156 * than the new functions.
xorjoep 1:24714b45cd1b 157 * \par
xorjoep 1:24714b45cd1b 158 * An example of initialization of the constants for the arm_cfft_q31 function follows:
xorjoep 1:24714b45cd1b 159 * \code
xorjoep 1:24714b45cd1b 160 * const static arm_cfft_instance_q31 *S;
xorjoep 1:24714b45cd1b 161 * ...
xorjoep 1:24714b45cd1b 162 * switch (length) {
xorjoep 1:24714b45cd1b 163 * case 16:
xorjoep 1:24714b45cd1b 164 * S = &arm_cfft_sR_q31_len16;
xorjoep 1:24714b45cd1b 165 * break;
xorjoep 1:24714b45cd1b 166 * case 32:
xorjoep 1:24714b45cd1b 167 * S = &arm_cfft_sR_q31_len32;
xorjoep 1:24714b45cd1b 168 * break;
xorjoep 1:24714b45cd1b 169 * case 64:
xorjoep 1:24714b45cd1b 170 * S = &arm_cfft_sR_q31_len64;
xorjoep 1:24714b45cd1b 171 * break;
xorjoep 1:24714b45cd1b 172 * case 128:
xorjoep 1:24714b45cd1b 173 * S = &arm_cfft_sR_q31_len128;
xorjoep 1:24714b45cd1b 174 * break;
xorjoep 1:24714b45cd1b 175 * case 256:
xorjoep 1:24714b45cd1b 176 * S = &arm_cfft_sR_q31_len256;
xorjoep 1:24714b45cd1b 177 * break;
xorjoep 1:24714b45cd1b 178 * case 512:
xorjoep 1:24714b45cd1b 179 * S = &arm_cfft_sR_q31_len512;
xorjoep 1:24714b45cd1b 180 * break;
xorjoep 1:24714b45cd1b 181 * case 1024:
xorjoep 1:24714b45cd1b 182 * S = &arm_cfft_sR_q31_len1024;
xorjoep 1:24714b45cd1b 183 * break;
xorjoep 1:24714b45cd1b 184 * case 2048:
xorjoep 1:24714b45cd1b 185 * S = &arm_cfft_sR_q31_len2048;
xorjoep 1:24714b45cd1b 186 * break;
xorjoep 1:24714b45cd1b 187 * case 4096:
xorjoep 1:24714b45cd1b 188 * S = &arm_cfft_sR_q31_len4096;
xorjoep 1:24714b45cd1b 189 * break;
xorjoep 1:24714b45cd1b 190 * }
xorjoep 1:24714b45cd1b 191 * \endcode
xorjoep 1:24714b45cd1b 192 *
xorjoep 1:24714b45cd1b 193 */
xorjoep 1:24714b45cd1b 194
xorjoep 1:24714b45cd1b 195 void arm_cfft_radix8by2_f32( arm_cfft_instance_f32 * S, float32_t * p1)
xorjoep 1:24714b45cd1b 196 {
xorjoep 1:24714b45cd1b 197 uint32_t L = S->fftLen;
xorjoep 1:24714b45cd1b 198 float32_t * pCol1, * pCol2, * pMid1, * pMid2;
xorjoep 1:24714b45cd1b 199 float32_t * p2 = p1 + L;
xorjoep 1:24714b45cd1b 200 const float32_t * tw = (float32_t *) S->pTwiddle;
xorjoep 1:24714b45cd1b 201 float32_t t1[4], t2[4], t3[4], t4[4], twR, twI;
xorjoep 1:24714b45cd1b 202 float32_t m0, m1, m2, m3;
xorjoep 1:24714b45cd1b 203 uint32_t l;
xorjoep 1:24714b45cd1b 204
xorjoep 1:24714b45cd1b 205 pCol1 = p1;
xorjoep 1:24714b45cd1b 206 pCol2 = p2;
xorjoep 1:24714b45cd1b 207
xorjoep 1:24714b45cd1b 208 // Define new length
xorjoep 1:24714b45cd1b 209 L >>= 1;
xorjoep 1:24714b45cd1b 210 // Initialize mid pointers
xorjoep 1:24714b45cd1b 211 pMid1 = p1 + L;
xorjoep 1:24714b45cd1b 212 pMid2 = p2 + L;
xorjoep 1:24714b45cd1b 213
xorjoep 1:24714b45cd1b 214 // do two dot Fourier transform
xorjoep 1:24714b45cd1b 215 for ( l = L >> 2; l > 0; l-- )
xorjoep 1:24714b45cd1b 216 {
xorjoep 1:24714b45cd1b 217 t1[0] = p1[0];
xorjoep 1:24714b45cd1b 218 t1[1] = p1[1];
xorjoep 1:24714b45cd1b 219 t1[2] = p1[2];
xorjoep 1:24714b45cd1b 220 t1[3] = p1[3];
xorjoep 1:24714b45cd1b 221
xorjoep 1:24714b45cd1b 222 t2[0] = p2[0];
xorjoep 1:24714b45cd1b 223 t2[1] = p2[1];
xorjoep 1:24714b45cd1b 224 t2[2] = p2[2];
xorjoep 1:24714b45cd1b 225 t2[3] = p2[3];
xorjoep 1:24714b45cd1b 226
xorjoep 1:24714b45cd1b 227 t3[0] = pMid1[0];
xorjoep 1:24714b45cd1b 228 t3[1] = pMid1[1];
xorjoep 1:24714b45cd1b 229 t3[2] = pMid1[2];
xorjoep 1:24714b45cd1b 230 t3[3] = pMid1[3];
xorjoep 1:24714b45cd1b 231
xorjoep 1:24714b45cd1b 232 t4[0] = pMid2[0];
xorjoep 1:24714b45cd1b 233 t4[1] = pMid2[1];
xorjoep 1:24714b45cd1b 234 t4[2] = pMid2[2];
xorjoep 1:24714b45cd1b 235 t4[3] = pMid2[3];
xorjoep 1:24714b45cd1b 236
xorjoep 1:24714b45cd1b 237 *p1++ = t1[0] + t2[0];
xorjoep 1:24714b45cd1b 238 *p1++ = t1[1] + t2[1];
xorjoep 1:24714b45cd1b 239 *p1++ = t1[2] + t2[2];
xorjoep 1:24714b45cd1b 240 *p1++ = t1[3] + t2[3]; // col 1
xorjoep 1:24714b45cd1b 241
xorjoep 1:24714b45cd1b 242 t2[0] = t1[0] - t2[0];
xorjoep 1:24714b45cd1b 243 t2[1] = t1[1] - t2[1];
xorjoep 1:24714b45cd1b 244 t2[2] = t1[2] - t2[2];
xorjoep 1:24714b45cd1b 245 t2[3] = t1[3] - t2[3]; // for col 2
xorjoep 1:24714b45cd1b 246
xorjoep 1:24714b45cd1b 247 *pMid1++ = t3[0] + t4[0];
xorjoep 1:24714b45cd1b 248 *pMid1++ = t3[1] + t4[1];
xorjoep 1:24714b45cd1b 249 *pMid1++ = t3[2] + t4[2];
xorjoep 1:24714b45cd1b 250 *pMid1++ = t3[3] + t4[3]; // col 1
xorjoep 1:24714b45cd1b 251
xorjoep 1:24714b45cd1b 252 t4[0] = t4[0] - t3[0];
xorjoep 1:24714b45cd1b 253 t4[1] = t4[1] - t3[1];
xorjoep 1:24714b45cd1b 254 t4[2] = t4[2] - t3[2];
xorjoep 1:24714b45cd1b 255 t4[3] = t4[3] - t3[3]; // for col 2
xorjoep 1:24714b45cd1b 256
xorjoep 1:24714b45cd1b 257 twR = *tw++;
xorjoep 1:24714b45cd1b 258 twI = *tw++;
xorjoep 1:24714b45cd1b 259
xorjoep 1:24714b45cd1b 260 // multiply by twiddle factors
xorjoep 1:24714b45cd1b 261 m0 = t2[0] * twR;
xorjoep 1:24714b45cd1b 262 m1 = t2[1] * twI;
xorjoep 1:24714b45cd1b 263 m2 = t2[1] * twR;
xorjoep 1:24714b45cd1b 264 m3 = t2[0] * twI;
xorjoep 1:24714b45cd1b 265
xorjoep 1:24714b45cd1b 266 // R = R * Tr - I * Ti
xorjoep 1:24714b45cd1b 267 *p2++ = m0 + m1;
xorjoep 1:24714b45cd1b 268 // I = I * Tr + R * Ti
xorjoep 1:24714b45cd1b 269 *p2++ = m2 - m3;
xorjoep 1:24714b45cd1b 270
xorjoep 1:24714b45cd1b 271 // use vertical symmetry
xorjoep 1:24714b45cd1b 272 // 0.9988 - 0.0491i <==> -0.0491 - 0.9988i
xorjoep 1:24714b45cd1b 273 m0 = t4[0] * twI;
xorjoep 1:24714b45cd1b 274 m1 = t4[1] * twR;
xorjoep 1:24714b45cd1b 275 m2 = t4[1] * twI;
xorjoep 1:24714b45cd1b 276 m3 = t4[0] * twR;
xorjoep 1:24714b45cd1b 277
xorjoep 1:24714b45cd1b 278 *pMid2++ = m0 - m1;
xorjoep 1:24714b45cd1b 279 *pMid2++ = m2 + m3;
xorjoep 1:24714b45cd1b 280
xorjoep 1:24714b45cd1b 281 twR = *tw++;
xorjoep 1:24714b45cd1b 282 twI = *tw++;
xorjoep 1:24714b45cd1b 283
xorjoep 1:24714b45cd1b 284 m0 = t2[2] * twR;
xorjoep 1:24714b45cd1b 285 m1 = t2[3] * twI;
xorjoep 1:24714b45cd1b 286 m2 = t2[3] * twR;
xorjoep 1:24714b45cd1b 287 m3 = t2[2] * twI;
xorjoep 1:24714b45cd1b 288
xorjoep 1:24714b45cd1b 289 *p2++ = m0 + m1;
xorjoep 1:24714b45cd1b 290 *p2++ = m2 - m3;
xorjoep 1:24714b45cd1b 291
xorjoep 1:24714b45cd1b 292 m0 = t4[2] * twI;
xorjoep 1:24714b45cd1b 293 m1 = t4[3] * twR;
xorjoep 1:24714b45cd1b 294 m2 = t4[3] * twI;
xorjoep 1:24714b45cd1b 295 m3 = t4[2] * twR;
xorjoep 1:24714b45cd1b 296
xorjoep 1:24714b45cd1b 297 *pMid2++ = m0 - m1;
xorjoep 1:24714b45cd1b 298 *pMid2++ = m2 + m3;
xorjoep 1:24714b45cd1b 299 }
xorjoep 1:24714b45cd1b 300
xorjoep 1:24714b45cd1b 301 // first col
xorjoep 1:24714b45cd1b 302 arm_radix8_butterfly_f32( pCol1, L, (float32_t *) S->pTwiddle, 2U);
xorjoep 1:24714b45cd1b 303 // second col
xorjoep 1:24714b45cd1b 304 arm_radix8_butterfly_f32( pCol2, L, (float32_t *) S->pTwiddle, 2U);
xorjoep 1:24714b45cd1b 305 }
xorjoep 1:24714b45cd1b 306
xorjoep 1:24714b45cd1b 307 void arm_cfft_radix8by4_f32( arm_cfft_instance_f32 * S, float32_t * p1)
xorjoep 1:24714b45cd1b 308 {
xorjoep 1:24714b45cd1b 309 uint32_t L = S->fftLen >> 1;
xorjoep 1:24714b45cd1b 310 float32_t * pCol1, *pCol2, *pCol3, *pCol4, *pEnd1, *pEnd2, *pEnd3, *pEnd4;
xorjoep 1:24714b45cd1b 311 const float32_t *tw2, *tw3, *tw4;
xorjoep 1:24714b45cd1b 312 float32_t * p2 = p1 + L;
xorjoep 1:24714b45cd1b 313 float32_t * p3 = p2 + L;
xorjoep 1:24714b45cd1b 314 float32_t * p4 = p3 + L;
xorjoep 1:24714b45cd1b 315 float32_t t2[4], t3[4], t4[4], twR, twI;
xorjoep 1:24714b45cd1b 316 float32_t p1ap3_0, p1sp3_0, p1ap3_1, p1sp3_1;
xorjoep 1:24714b45cd1b 317 float32_t m0, m1, m2, m3;
xorjoep 1:24714b45cd1b 318 uint32_t l, twMod2, twMod3, twMod4;
xorjoep 1:24714b45cd1b 319
xorjoep 1:24714b45cd1b 320 pCol1 = p1; // points to real values by default
xorjoep 1:24714b45cd1b 321 pCol2 = p2;
xorjoep 1:24714b45cd1b 322 pCol3 = p3;
xorjoep 1:24714b45cd1b 323 pCol4 = p4;
xorjoep 1:24714b45cd1b 324 pEnd1 = p2 - 1; // points to imaginary values by default
xorjoep 1:24714b45cd1b 325 pEnd2 = p3 - 1;
xorjoep 1:24714b45cd1b 326 pEnd3 = p4 - 1;
xorjoep 1:24714b45cd1b 327 pEnd4 = pEnd3 + L;
xorjoep 1:24714b45cd1b 328
xorjoep 1:24714b45cd1b 329 tw2 = tw3 = tw4 = (float32_t *) S->pTwiddle;
xorjoep 1:24714b45cd1b 330
xorjoep 1:24714b45cd1b 331 L >>= 1;
xorjoep 1:24714b45cd1b 332
xorjoep 1:24714b45cd1b 333 // do four dot Fourier transform
xorjoep 1:24714b45cd1b 334
xorjoep 1:24714b45cd1b 335 twMod2 = 2;
xorjoep 1:24714b45cd1b 336 twMod3 = 4;
xorjoep 1:24714b45cd1b 337 twMod4 = 6;
xorjoep 1:24714b45cd1b 338
xorjoep 1:24714b45cd1b 339 // TOP
xorjoep 1:24714b45cd1b 340 p1ap3_0 = p1[0] + p3[0];
xorjoep 1:24714b45cd1b 341 p1sp3_0 = p1[0] - p3[0];
xorjoep 1:24714b45cd1b 342 p1ap3_1 = p1[1] + p3[1];
xorjoep 1:24714b45cd1b 343 p1sp3_1 = p1[1] - p3[1];
xorjoep 1:24714b45cd1b 344
xorjoep 1:24714b45cd1b 345 // col 2
xorjoep 1:24714b45cd1b 346 t2[0] = p1sp3_0 + p2[1] - p4[1];
xorjoep 1:24714b45cd1b 347 t2[1] = p1sp3_1 - p2[0] + p4[0];
xorjoep 1:24714b45cd1b 348 // col 3
xorjoep 1:24714b45cd1b 349 t3[0] = p1ap3_0 - p2[0] - p4[0];
xorjoep 1:24714b45cd1b 350 t3[1] = p1ap3_1 - p2[1] - p4[1];
xorjoep 1:24714b45cd1b 351 // col 4
xorjoep 1:24714b45cd1b 352 t4[0] = p1sp3_0 - p2[1] + p4[1];
xorjoep 1:24714b45cd1b 353 t4[1] = p1sp3_1 + p2[0] - p4[0];
xorjoep 1:24714b45cd1b 354 // col 1
xorjoep 1:24714b45cd1b 355 *p1++ = p1ap3_0 + p2[0] + p4[0];
xorjoep 1:24714b45cd1b 356 *p1++ = p1ap3_1 + p2[1] + p4[1];
xorjoep 1:24714b45cd1b 357
xorjoep 1:24714b45cd1b 358 // Twiddle factors are ones
xorjoep 1:24714b45cd1b 359 *p2++ = t2[0];
xorjoep 1:24714b45cd1b 360 *p2++ = t2[1];
xorjoep 1:24714b45cd1b 361 *p3++ = t3[0];
xorjoep 1:24714b45cd1b 362 *p3++ = t3[1];
xorjoep 1:24714b45cd1b 363 *p4++ = t4[0];
xorjoep 1:24714b45cd1b 364 *p4++ = t4[1];
xorjoep 1:24714b45cd1b 365
xorjoep 1:24714b45cd1b 366 tw2 += twMod2;
xorjoep 1:24714b45cd1b 367 tw3 += twMod3;
xorjoep 1:24714b45cd1b 368 tw4 += twMod4;
xorjoep 1:24714b45cd1b 369
xorjoep 1:24714b45cd1b 370 for (l = (L - 2) >> 1; l > 0; l-- )
xorjoep 1:24714b45cd1b 371 {
xorjoep 1:24714b45cd1b 372 // TOP
xorjoep 1:24714b45cd1b 373 p1ap3_0 = p1[0] + p3[0];
xorjoep 1:24714b45cd1b 374 p1sp3_0 = p1[0] - p3[0];
xorjoep 1:24714b45cd1b 375 p1ap3_1 = p1[1] + p3[1];
xorjoep 1:24714b45cd1b 376 p1sp3_1 = p1[1] - p3[1];
xorjoep 1:24714b45cd1b 377 // col 2
xorjoep 1:24714b45cd1b 378 t2[0] = p1sp3_0 + p2[1] - p4[1];
xorjoep 1:24714b45cd1b 379 t2[1] = p1sp3_1 - p2[0] + p4[0];
xorjoep 1:24714b45cd1b 380 // col 3
xorjoep 1:24714b45cd1b 381 t3[0] = p1ap3_0 - p2[0] - p4[0];
xorjoep 1:24714b45cd1b 382 t3[1] = p1ap3_1 - p2[1] - p4[1];
xorjoep 1:24714b45cd1b 383 // col 4
xorjoep 1:24714b45cd1b 384 t4[0] = p1sp3_0 - p2[1] + p4[1];
xorjoep 1:24714b45cd1b 385 t4[1] = p1sp3_1 + p2[0] - p4[0];
xorjoep 1:24714b45cd1b 386 // col 1 - top
xorjoep 1:24714b45cd1b 387 *p1++ = p1ap3_0 + p2[0] + p4[0];
xorjoep 1:24714b45cd1b 388 *p1++ = p1ap3_1 + p2[1] + p4[1];
xorjoep 1:24714b45cd1b 389
xorjoep 1:24714b45cd1b 390 // BOTTOM
xorjoep 1:24714b45cd1b 391 p1ap3_1 = pEnd1[-1] + pEnd3[-1];
xorjoep 1:24714b45cd1b 392 p1sp3_1 = pEnd1[-1] - pEnd3[-1];
xorjoep 1:24714b45cd1b 393 p1ap3_0 = pEnd1[0] + pEnd3[0];
xorjoep 1:24714b45cd1b 394 p1sp3_0 = pEnd1[0] - pEnd3[0];
xorjoep 1:24714b45cd1b 395 // col 2
xorjoep 1:24714b45cd1b 396 t2[2] = pEnd2[0] - pEnd4[0] + p1sp3_1;
xorjoep 1:24714b45cd1b 397 t2[3] = pEnd1[0] - pEnd3[0] - pEnd2[-1] + pEnd4[-1];
xorjoep 1:24714b45cd1b 398 // col 3
xorjoep 1:24714b45cd1b 399 t3[2] = p1ap3_1 - pEnd2[-1] - pEnd4[-1];
xorjoep 1:24714b45cd1b 400 t3[3] = p1ap3_0 - pEnd2[0] - pEnd4[0];
xorjoep 1:24714b45cd1b 401 // col 4
xorjoep 1:24714b45cd1b 402 t4[2] = pEnd2[0] - pEnd4[0] - p1sp3_1;
xorjoep 1:24714b45cd1b 403 t4[3] = pEnd4[-1] - pEnd2[-1] - p1sp3_0;
xorjoep 1:24714b45cd1b 404 // col 1 - Bottom
xorjoep 1:24714b45cd1b 405 *pEnd1-- = p1ap3_0 + pEnd2[0] + pEnd4[0];
xorjoep 1:24714b45cd1b 406 *pEnd1-- = p1ap3_1 + pEnd2[-1] + pEnd4[-1];
xorjoep 1:24714b45cd1b 407
xorjoep 1:24714b45cd1b 408 // COL 2
xorjoep 1:24714b45cd1b 409 // read twiddle factors
xorjoep 1:24714b45cd1b 410 twR = *tw2++;
xorjoep 1:24714b45cd1b 411 twI = *tw2++;
xorjoep 1:24714b45cd1b 412 // multiply by twiddle factors
xorjoep 1:24714b45cd1b 413 // let Z1 = a + i(b), Z2 = c + i(d)
xorjoep 1:24714b45cd1b 414 // => Z1 * Z2 = (a*c - b*d) + i(b*c + a*d)
xorjoep 1:24714b45cd1b 415
xorjoep 1:24714b45cd1b 416 // Top
xorjoep 1:24714b45cd1b 417 m0 = t2[0] * twR;
xorjoep 1:24714b45cd1b 418 m1 = t2[1] * twI;
xorjoep 1:24714b45cd1b 419 m2 = t2[1] * twR;
xorjoep 1:24714b45cd1b 420 m3 = t2[0] * twI;
xorjoep 1:24714b45cd1b 421
xorjoep 1:24714b45cd1b 422 *p2++ = m0 + m1;
xorjoep 1:24714b45cd1b 423 *p2++ = m2 - m3;
xorjoep 1:24714b45cd1b 424 // use vertical symmetry col 2
xorjoep 1:24714b45cd1b 425 // 0.9997 - 0.0245i <==> 0.0245 - 0.9997i
xorjoep 1:24714b45cd1b 426 // Bottom
xorjoep 1:24714b45cd1b 427 m0 = t2[3] * twI;
xorjoep 1:24714b45cd1b 428 m1 = t2[2] * twR;
xorjoep 1:24714b45cd1b 429 m2 = t2[2] * twI;
xorjoep 1:24714b45cd1b 430 m3 = t2[3] * twR;
xorjoep 1:24714b45cd1b 431
xorjoep 1:24714b45cd1b 432 *pEnd2-- = m0 - m1;
xorjoep 1:24714b45cd1b 433 *pEnd2-- = m2 + m3;
xorjoep 1:24714b45cd1b 434
xorjoep 1:24714b45cd1b 435 // COL 3
xorjoep 1:24714b45cd1b 436 twR = tw3[0];
xorjoep 1:24714b45cd1b 437 twI = tw3[1];
xorjoep 1:24714b45cd1b 438 tw3 += twMod3;
xorjoep 1:24714b45cd1b 439 // Top
xorjoep 1:24714b45cd1b 440 m0 = t3[0] * twR;
xorjoep 1:24714b45cd1b 441 m1 = t3[1] * twI;
xorjoep 1:24714b45cd1b 442 m2 = t3[1] * twR;
xorjoep 1:24714b45cd1b 443 m3 = t3[0] * twI;
xorjoep 1:24714b45cd1b 444
xorjoep 1:24714b45cd1b 445 *p3++ = m0 + m1;
xorjoep 1:24714b45cd1b 446 *p3++ = m2 - m3;
xorjoep 1:24714b45cd1b 447 // use vertical symmetry col 3
xorjoep 1:24714b45cd1b 448 // 0.9988 - 0.0491i <==> -0.9988 - 0.0491i
xorjoep 1:24714b45cd1b 449 // Bottom
xorjoep 1:24714b45cd1b 450 m0 = -t3[3] * twR;
xorjoep 1:24714b45cd1b 451 m1 = t3[2] * twI;
xorjoep 1:24714b45cd1b 452 m2 = t3[2] * twR;
xorjoep 1:24714b45cd1b 453 m3 = t3[3] * twI;
xorjoep 1:24714b45cd1b 454
xorjoep 1:24714b45cd1b 455 *pEnd3-- = m0 - m1;
xorjoep 1:24714b45cd1b 456 *pEnd3-- = m3 - m2;
xorjoep 1:24714b45cd1b 457
xorjoep 1:24714b45cd1b 458 // COL 4
xorjoep 1:24714b45cd1b 459 twR = tw4[0];
xorjoep 1:24714b45cd1b 460 twI = tw4[1];
xorjoep 1:24714b45cd1b 461 tw4 += twMod4;
xorjoep 1:24714b45cd1b 462 // Top
xorjoep 1:24714b45cd1b 463 m0 = t4[0] * twR;
xorjoep 1:24714b45cd1b 464 m1 = t4[1] * twI;
xorjoep 1:24714b45cd1b 465 m2 = t4[1] * twR;
xorjoep 1:24714b45cd1b 466 m3 = t4[0] * twI;
xorjoep 1:24714b45cd1b 467
xorjoep 1:24714b45cd1b 468 *p4++ = m0 + m1;
xorjoep 1:24714b45cd1b 469 *p4++ = m2 - m3;
xorjoep 1:24714b45cd1b 470 // use vertical symmetry col 4
xorjoep 1:24714b45cd1b 471 // 0.9973 - 0.0736i <==> -0.0736 + 0.9973i
xorjoep 1:24714b45cd1b 472 // Bottom
xorjoep 1:24714b45cd1b 473 m0 = t4[3] * twI;
xorjoep 1:24714b45cd1b 474 m1 = t4[2] * twR;
xorjoep 1:24714b45cd1b 475 m2 = t4[2] * twI;
xorjoep 1:24714b45cd1b 476 m3 = t4[3] * twR;
xorjoep 1:24714b45cd1b 477
xorjoep 1:24714b45cd1b 478 *pEnd4-- = m0 - m1;
xorjoep 1:24714b45cd1b 479 *pEnd4-- = m2 + m3;
xorjoep 1:24714b45cd1b 480 }
xorjoep 1:24714b45cd1b 481
xorjoep 1:24714b45cd1b 482 //MIDDLE
xorjoep 1:24714b45cd1b 483 // Twiddle factors are
xorjoep 1:24714b45cd1b 484 // 1.0000 0.7071-0.7071i -1.0000i -0.7071-0.7071i
xorjoep 1:24714b45cd1b 485 p1ap3_0 = p1[0] + p3[0];
xorjoep 1:24714b45cd1b 486 p1sp3_0 = p1[0] - p3[0];
xorjoep 1:24714b45cd1b 487 p1ap3_1 = p1[1] + p3[1];
xorjoep 1:24714b45cd1b 488 p1sp3_1 = p1[1] - p3[1];
xorjoep 1:24714b45cd1b 489
xorjoep 1:24714b45cd1b 490 // col 2
xorjoep 1:24714b45cd1b 491 t2[0] = p1sp3_0 + p2[1] - p4[1];
xorjoep 1:24714b45cd1b 492 t2[1] = p1sp3_1 - p2[0] + p4[0];
xorjoep 1:24714b45cd1b 493 // col 3
xorjoep 1:24714b45cd1b 494 t3[0] = p1ap3_0 - p2[0] - p4[0];
xorjoep 1:24714b45cd1b 495 t3[1] = p1ap3_1 - p2[1] - p4[1];
xorjoep 1:24714b45cd1b 496 // col 4
xorjoep 1:24714b45cd1b 497 t4[0] = p1sp3_0 - p2[1] + p4[1];
xorjoep 1:24714b45cd1b 498 t4[1] = p1sp3_1 + p2[0] - p4[0];
xorjoep 1:24714b45cd1b 499 // col 1 - Top
xorjoep 1:24714b45cd1b 500 *p1++ = p1ap3_0 + p2[0] + p4[0];
xorjoep 1:24714b45cd1b 501 *p1++ = p1ap3_1 + p2[1] + p4[1];
xorjoep 1:24714b45cd1b 502
xorjoep 1:24714b45cd1b 503 // COL 2
xorjoep 1:24714b45cd1b 504 twR = tw2[0];
xorjoep 1:24714b45cd1b 505 twI = tw2[1];
xorjoep 1:24714b45cd1b 506
xorjoep 1:24714b45cd1b 507 m0 = t2[0] * twR;
xorjoep 1:24714b45cd1b 508 m1 = t2[1] * twI;
xorjoep 1:24714b45cd1b 509 m2 = t2[1] * twR;
xorjoep 1:24714b45cd1b 510 m3 = t2[0] * twI;
xorjoep 1:24714b45cd1b 511
xorjoep 1:24714b45cd1b 512 *p2++ = m0 + m1;
xorjoep 1:24714b45cd1b 513 *p2++ = m2 - m3;
xorjoep 1:24714b45cd1b 514 // COL 3
xorjoep 1:24714b45cd1b 515 twR = tw3[0];
xorjoep 1:24714b45cd1b 516 twI = tw3[1];
xorjoep 1:24714b45cd1b 517
xorjoep 1:24714b45cd1b 518 m0 = t3[0] * twR;
xorjoep 1:24714b45cd1b 519 m1 = t3[1] * twI;
xorjoep 1:24714b45cd1b 520 m2 = t3[1] * twR;
xorjoep 1:24714b45cd1b 521 m3 = t3[0] * twI;
xorjoep 1:24714b45cd1b 522
xorjoep 1:24714b45cd1b 523 *p3++ = m0 + m1;
xorjoep 1:24714b45cd1b 524 *p3++ = m2 - m3;
xorjoep 1:24714b45cd1b 525 // COL 4
xorjoep 1:24714b45cd1b 526 twR = tw4[0];
xorjoep 1:24714b45cd1b 527 twI = tw4[1];
xorjoep 1:24714b45cd1b 528
xorjoep 1:24714b45cd1b 529 m0 = t4[0] * twR;
xorjoep 1:24714b45cd1b 530 m1 = t4[1] * twI;
xorjoep 1:24714b45cd1b 531 m2 = t4[1] * twR;
xorjoep 1:24714b45cd1b 532 m3 = t4[0] * twI;
xorjoep 1:24714b45cd1b 533
xorjoep 1:24714b45cd1b 534 *p4++ = m0 + m1;
xorjoep 1:24714b45cd1b 535 *p4++ = m2 - m3;
xorjoep 1:24714b45cd1b 536
xorjoep 1:24714b45cd1b 537 // first col
xorjoep 1:24714b45cd1b 538 arm_radix8_butterfly_f32( pCol1, L, (float32_t *) S->pTwiddle, 4U);
xorjoep 1:24714b45cd1b 539 // second col
xorjoep 1:24714b45cd1b 540 arm_radix8_butterfly_f32( pCol2, L, (float32_t *) S->pTwiddle, 4U);
xorjoep 1:24714b45cd1b 541 // third col
xorjoep 1:24714b45cd1b 542 arm_radix8_butterfly_f32( pCol3, L, (float32_t *) S->pTwiddle, 4U);
xorjoep 1:24714b45cd1b 543 // fourth col
xorjoep 1:24714b45cd1b 544 arm_radix8_butterfly_f32( pCol4, L, (float32_t *) S->pTwiddle, 4U);
xorjoep 1:24714b45cd1b 545 }
xorjoep 1:24714b45cd1b 546
xorjoep 1:24714b45cd1b 547 /**
xorjoep 1:24714b45cd1b 548 * @addtogroup ComplexFFT
xorjoep 1:24714b45cd1b 549 * @{
xorjoep 1:24714b45cd1b 550 */
xorjoep 1:24714b45cd1b 551
xorjoep 1:24714b45cd1b 552 /**
xorjoep 1:24714b45cd1b 553 * @details
xorjoep 1:24714b45cd1b 554 * @brief Processing function for the floating-point complex FFT.
xorjoep 1:24714b45cd1b 555 * @param[in] *S points to an instance of the floating-point CFFT structure.
xorjoep 1:24714b45cd1b 556 * @param[in, out] *p1 points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.
xorjoep 1:24714b45cd1b 557 * @param[in] ifftFlag flag that selects forward (ifftFlag=0) or inverse (ifftFlag=1) transform.
xorjoep 1:24714b45cd1b 558 * @param[in] bitReverseFlag flag that enables (bitReverseFlag=1) or disables (bitReverseFlag=0) bit reversal of output.
xorjoep 1:24714b45cd1b 559 * @return none.
xorjoep 1:24714b45cd1b 560 */
xorjoep 1:24714b45cd1b 561
xorjoep 1:24714b45cd1b 562 void arm_cfft_f32(
xorjoep 1:24714b45cd1b 563 const arm_cfft_instance_f32 * S,
xorjoep 1:24714b45cd1b 564 float32_t * p1,
xorjoep 1:24714b45cd1b 565 uint8_t ifftFlag,
xorjoep 1:24714b45cd1b 566 uint8_t bitReverseFlag)
xorjoep 1:24714b45cd1b 567 {
xorjoep 1:24714b45cd1b 568 uint32_t L = S->fftLen, l;
xorjoep 1:24714b45cd1b 569 float32_t invL, * pSrc;
xorjoep 1:24714b45cd1b 570
xorjoep 1:24714b45cd1b 571 if (ifftFlag == 1U)
xorjoep 1:24714b45cd1b 572 {
xorjoep 1:24714b45cd1b 573 /* Conjugate input data */
xorjoep 1:24714b45cd1b 574 pSrc = p1 + 1;
xorjoep 1:24714b45cd1b 575 for(l=0; l<L; l++)
xorjoep 1:24714b45cd1b 576 {
xorjoep 1:24714b45cd1b 577 *pSrc = -*pSrc;
xorjoep 1:24714b45cd1b 578 pSrc += 2;
xorjoep 1:24714b45cd1b 579 }
xorjoep 1:24714b45cd1b 580 }
xorjoep 1:24714b45cd1b 581
xorjoep 1:24714b45cd1b 582 switch (L)
xorjoep 1:24714b45cd1b 583 {
xorjoep 1:24714b45cd1b 584 case 16:
xorjoep 1:24714b45cd1b 585 case 128:
xorjoep 1:24714b45cd1b 586 case 1024:
xorjoep 1:24714b45cd1b 587 arm_cfft_radix8by2_f32 ( (arm_cfft_instance_f32 *) S, p1);
xorjoep 1:24714b45cd1b 588 break;
xorjoep 1:24714b45cd1b 589 case 32:
xorjoep 1:24714b45cd1b 590 case 256:
xorjoep 1:24714b45cd1b 591 case 2048:
xorjoep 1:24714b45cd1b 592 arm_cfft_radix8by4_f32 ( (arm_cfft_instance_f32 *) S, p1);
xorjoep 1:24714b45cd1b 593 break;
xorjoep 1:24714b45cd1b 594 case 64:
xorjoep 1:24714b45cd1b 595 case 512:
xorjoep 1:24714b45cd1b 596 case 4096:
xorjoep 1:24714b45cd1b 597 arm_radix8_butterfly_f32( p1, L, (float32_t *) S->pTwiddle, 1);
xorjoep 1:24714b45cd1b 598 break;
xorjoep 1:24714b45cd1b 599 }
xorjoep 1:24714b45cd1b 600
xorjoep 1:24714b45cd1b 601 if ( bitReverseFlag )
xorjoep 1:24714b45cd1b 602 arm_bitreversal_32((uint32_t*)p1,S->bitRevLength,S->pBitRevTable);
xorjoep 1:24714b45cd1b 603
xorjoep 1:24714b45cd1b 604 if (ifftFlag == 1U)
xorjoep 1:24714b45cd1b 605 {
xorjoep 1:24714b45cd1b 606 invL = 1.0f/(float32_t)L;
xorjoep 1:24714b45cd1b 607 /* Conjugate and scale output data */
xorjoep 1:24714b45cd1b 608 pSrc = p1;
xorjoep 1:24714b45cd1b 609 for(l=0; l<L; l++)
xorjoep 1:24714b45cd1b 610 {
xorjoep 1:24714b45cd1b 611 *pSrc++ *= invL ;
xorjoep 1:24714b45cd1b 612 *pSrc = -(*pSrc) * invL;
xorjoep 1:24714b45cd1b 613 pSrc++;
xorjoep 1:24714b45cd1b 614 }
xorjoep 1:24714b45cd1b 615 }
xorjoep 1:24714b45cd1b 616 }
xorjoep 1:24714b45cd1b 617
xorjoep 1:24714b45cd1b 618 /**
xorjoep 1:24714b45cd1b 619 * @} end of ComplexFFT group
xorjoep 1:24714b45cd1b 620 */