]>
git.gir.st - tmk_keyboard.git/blob - tmk_core/tool/mbed/mbed-sdk/libraries/dsp/cmsis_dsp/TransformFunctions/arm_cfft_f32.c
1 /* ----------------------------------------------------------------------
2 * Copyright (C) 2010-2013 ARM Limited. All rights reserved.
4 * $Date: 17. January 2013
7 * Project: CMSIS DSP Library
8 * Title: arm_cfft_f32.c
10 * Description: Combined Radix Decimation in Frequency CFFT Floating point processing function
12 * Target Processor: Cortex-M4/Cortex-M3/Cortex-M0
14 * Redistribution and use in source and binary forms, with or without
15 * modification, are permitted provided that the following conditions
17 * - Redistributions of source code must retain the above copyright
18 * notice, this list of conditions and the following disclaimer.
19 * - Redistributions in binary form must reproduce the above copyright
20 * notice, this list of conditions and the following disclaimer in
21 * the documentation and/or other materials provided with the
23 * - Neither the name of ARM LIMITED nor the names of its contributors
24 * may be used to endorse or promote products derived from this
25 * software without specific prior written permission.
27 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
28 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
29 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
30 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
31 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
32 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
33 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
34 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
35 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
36 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
37 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
38 * POSSIBILITY OF SUCH DAMAGE.
39 * -------------------------------------------------------------------- */
43 #include "arm_common_tables.h"
45 extern void arm_radix8_butterfly_f32 (
48 const float32_t
* pCoef
,
49 uint16_t twidCoefModifier
);
51 extern void arm_bitreversal_32 (
53 const uint16_t bitRevLen
,
54 const uint16_t * pBitRevTable
);
57 * @ingroup groupTransforms
61 * @defgroup ComplexFFT Complex FFT Functions
64 * The Fast Fourier Transform (FFT) is an efficient algorithm for computing the
65 * Discrete Fourier Transform (DFT). The FFT can be orders of magnitude faster
66 * than the DFT, especially for long lengths.
67 * The algorithms described in this section
68 * operate on complex data. A separate set of functions is devoted to handling
71 * There are separate algorithms for handling floating-point, Q15, and Q31 data
72 * types. The algorithms available for each data type are described next.
74 * The FFT functions operate in-place. That is, the array holding the input data
75 * will also be used to hold the corresponding result. The input data is complex
76 * and contains <code>2*fftLen</code> interleaved values as shown below.
77 * <pre> {real[0], imag[0], real[1], imag[1],..} </pre>
78 * The FFT result will be contained in the same array and the frequency domain
79 * values will have the same interleaving.
82 * The floating-point complex FFT uses a mixed-radix algorithm. Multiple radix-8
83 * stages are performed along with a single radix-2 or radix-4 stage, as needed.
84 * The algorithm supports lengths of [16, 32, 64, ..., 4096] and each length uses
85 * a different twiddle factor table.
87 * The function uses the standard FFT definition and output values may grow by a
88 * factor of <code>fftLen</code> when computing the forward transform. The
89 * inverse transform includes a scale of <code>1/fftLen</code> as part of the
90 * calculation and this matches the textbook definition of the inverse FFT.
92 * Preinitialized data structures containing twiddle factors and bit reversal
93 * tables are provided and defined in <code>arm_const_structs.h</code>. Include
94 * this header in your function and then pass one of the constant structures as
95 * an argument to arm_cfft_f32. For example:
97 * <code>arm_cfft_f32(arm_cfft_sR_f32_len64, pSrc, 1, 1)</code>
99 * computes a 64-point inverse complex FFT including bit reversal.
100 * The data structures are treated as constant data and not modified during the
101 * calculation. The same data structure can be reused for multiple transforms
102 * including mixing forward and inverse transforms.
104 * Earlier releases of the library provided separate radix-2 and radix-4
105 * algorithms that operated on floating-point data. These functions are still
106 * provided but are deprecated. The older functions are slower and less general
107 * than the new functions.
109 * An example of initialization of the constants for the arm_cfft_f32 function follows:
111 * const static arm_cfft_instance_f32 *S;
115 * S = & arm_cfft_sR_f32_len16;
118 * S = & arm_cfft_sR_f32_len32;
121 * S = & arm_cfft_sR_f32_len64;
124 * S = & arm_cfft_sR_f32_len128;
127 * S = & arm_cfft_sR_f32_len256;
130 * S = & arm_cfft_sR_f32_len512;
133 * S = & arm_cfft_sR_f32_len1024;
136 * S = & arm_cfft_sR_f32_len2048;
139 * S = & arm_cfft_sR_f32_len4096;
143 * The library provides radix-2 and radix-4 FFT algorithms for fixed-point data. The
144 * radix-2 algorithm supports lengths of [16, 32, 64, ..., 4096]. The radix-4
145 * algorithm supports lengths of [16, 64, 256, ..., 4096]. When possible, you
146 * should use the radix-4 algorithm since it is faster than the radix-2 of the
149 * The forward FFTs include scaling in order to prevent results from overflowing.
150 * Intermediate results are scaled down during each butterfly stage. In the
151 * radix-2 algorithm, a scale of 0.5 is applied during each butterfly. In the
152 * radix-4 algorithm, a scale of 0.25 is applied. The scaling applies to both
153 * the forward and the inverse FFTs. Thus the forward FFT contains an additional
154 * scale factor of <code>1/fftLen</code> as compared to the standard textbook
155 * definition of the FFT. The inverse FFT also scales down during each butterfly
156 * stage and this corresponds to the standard textbook definition.
158 * A separate instance structure must be defined for each transform used but
159 * twiddle factor and bit reversal tables can be reused.
161 * There is also an associated initialization function for each data type.
162 * The initialization function performs the following operations:
163 * - Sets the values of the internal structure fields.
164 * - Initializes twiddle factor table and bit reversal table pointers.
166 * Use of the initialization function is optional.
167 * However, if the initialization function is used, then the instance structure
168 * cannot be placed into a const data section. To place an instance structure
169 * into a const data section, the instance structure should be manually
170 * initialized as follows:
172 *arm_cfft_radix2_instance_q31 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};
173 *arm_cfft_radix2_instance_q15 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};
174 *arm_cfft_radix4_instance_q31 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};
175 *arm_cfft_radix4_instance_q15 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};
176 *arm_cfft_instance_f32 S = {fftLen, pTwiddle, pBitRevTable, bitRevLength};
179 * where <code>fftLen</code> length of CFFT/CIFFT; <code>ifftFlag</code> Flag for
180 * selection of forward or inverse transform. When ifftFlag is set the inverse
181 * transform is calculated.
182 * <code>bitReverseFlag</code> Flag for selection of output order (Set bitReverseFlag to output in normal order otherwise output in bit reversed order);
183 * <code>pTwiddle</code>points to array of twiddle coefficients; <code>pBitRevTable</code> points to the bit reversal table.
184 * <code>twidCoefModifier</code> modifier for twiddle factor table which supports all FFT lengths with same table;
185 * <code>pBitRevTable</code> modifier for bit reversal table which supports all FFT lengths with same table.
186 * <code>onebyfftLen</code> value of 1/fftLen to calculate CIFFT;
188 * The Q15 and Q31 FFT functions use a large bit reversal and twiddle factor
189 * table. The tables are defined for the maximum length transform and a subset
190 * of the coefficients are used in shorter transforms.
194 void arm_cfft_radix8by2_f32 ( arm_cfft_instance_f32
* S
, float32_t
* p1
)
196 uint32_t L
= S
-> fftLen
;
197 float32_t
* pCol1
, * pCol2
, * pMid1
, * pMid2
;
198 float32_t
* p2
= p1
+ L
;
199 const float32_t
* tw
= ( float32_t
*) S
-> pTwiddle
;
200 float32_t t1
[ 4 ], t2
[ 4 ], t3
[ 4 ], t4
[ 4 ], twR
, twI
;
201 float32_t m0
, m1
, m2
, m3
;
209 // Initialize mid pointers
213 // do two dot Fourier transform
214 for ( l
= L
>> 2 ; l
> 0 ; l
-- )
236 * p1
++ = t1
[ 0 ] + t2
[ 0 ];
237 * p1
++ = t1
[ 1 ] + t2
[ 1 ];
238 * p1
++ = t1
[ 2 ] + t2
[ 2 ];
239 * p1
++ = t1
[ 3 ] + t2
[ 3 ]; // col 1
241 t2
[ 0 ] = t1
[ 0 ] - t2
[ 0 ];
242 t2
[ 1 ] = t1
[ 1 ] - t2
[ 1 ];
243 t2
[ 2 ] = t1
[ 2 ] - t2
[ 2 ];
244 t2
[ 3 ] = t1
[ 3 ] - t2
[ 3 ]; // for col 2
246 * pMid1
++ = t3
[ 0 ] + t4
[ 0 ];
247 * pMid1
++ = t3
[ 1 ] + t4
[ 1 ];
248 * pMid1
++ = t3
[ 2 ] + t4
[ 2 ];
249 * pMid1
++ = t3
[ 3 ] + t4
[ 3 ]; // col 1
251 t4
[ 0 ] = t4
[ 0 ] - t3
[ 0 ];
252 t4
[ 1 ] = t4
[ 1 ] - t3
[ 1 ];
253 t4
[ 2 ] = t4
[ 2 ] - t3
[ 2 ];
254 t4
[ 3 ] = t4
[ 3 ] - t3
[ 3 ]; // for col 2
259 // multiply by twiddle factors
265 // R = R * Tr - I * Ti
267 // I = I * Tr + R * Ti
270 // use vertical symmetry
271 // 0.9988 - 0.0491i <==> -0.0491 - 0.9988i
301 arm_radix8_butterfly_f32 ( pCol1
, L
, ( float32_t
*) S
-> pTwiddle
, 2u );
303 arm_radix8_butterfly_f32 ( pCol2
, L
, ( float32_t
*) S
-> pTwiddle
, 2u );
307 void arm_cfft_radix8by4_f32 ( arm_cfft_instance_f32
* S
, float32_t
* p1
)
309 uint32_t L
= S
-> fftLen
>> 1 ;
310 float32_t
* pCol1
, * pCol2
, * pCol3
, * pCol4
, * pEnd1
, * pEnd2
, * pEnd3
, * pEnd4
;
311 const float32_t
* tw2
, * tw3
, * tw4
;
312 float32_t
* p2
= p1
+ L
;
313 float32_t
* p3
= p2
+ L
;
314 float32_t
* p4
= p3
+ L
;
315 float32_t t2
[ 4 ], t3
[ 4 ], t4
[ 4 ], twR
, twI
;
316 float32_t p1ap3_0
, p1sp3_0
, p1ap3_1
, p1sp3_1
;
317 float32_t m0
, m1
, m2
, m3
;
318 uint32_t l
, twMod2
, twMod3
, twMod4
;
320 pCol1
= p1
; // points to real values by default
324 pEnd1
= p2
- 1 ; // points to imaginary values by default
329 tw2
= tw3
= tw4
= ( float32_t
*) S
-> pTwiddle
;
333 // do four dot Fourier transform
340 p1ap3_0
= p1
[ 0 ] + p3
[ 0 ];
341 p1sp3_0
= p1
[ 0 ] - p3
[ 0 ];
342 p1ap3_1
= p1
[ 1 ] + p3
[ 1 ];
343 p1sp3_1
= p1
[ 1 ] - p3
[ 1 ];
346 t2
[ 0 ] = p1sp3_0
+ p2
[ 1 ] - p4
[ 1 ];
347 t2
[ 1 ] = p1sp3_1
- p2
[ 0 ] + p4
[ 0 ];
349 t3
[ 0 ] = p1ap3_0
- p2
[ 0 ] - p4
[ 0 ];
350 t3
[ 1 ] = p1ap3_1
- p2
[ 1 ] - p4
[ 1 ];
352 t4
[ 0 ] = p1sp3_0
- p2
[ 1 ] + p4
[ 1 ];
353 t4
[ 1 ] = p1sp3_1
+ p2
[ 0 ] - p4
[ 0 ];
355 * p1
++ = p1ap3_0
+ p2
[ 0 ] + p4
[ 0 ];
356 * p1
++ = p1ap3_1
+ p2
[ 1 ] + p4
[ 1 ];
358 // Twiddle factors are ones
370 for ( l
= ( L
- 2 ) >> 1 ; l
> 0 ; l
-- )
374 p1ap3_0
= p1
[ 0 ] + p3
[ 0 ];
375 p1sp3_0
= p1
[ 0 ] - p3
[ 0 ];
376 p1ap3_1
= p1
[ 1 ] + p3
[ 1 ];
377 p1sp3_1
= p1
[ 1 ] - p3
[ 1 ];
379 t2
[ 0 ] = p1sp3_0
+ p2
[ 1 ] - p4
[ 1 ];
380 t2
[ 1 ] = p1sp3_1
- p2
[ 0 ] + p4
[ 0 ];
382 t3
[ 0 ] = p1ap3_0
- p2
[ 0 ] - p4
[ 0 ];
383 t3
[ 1 ] = p1ap3_1
- p2
[ 1 ] - p4
[ 1 ];
385 t4
[ 0 ] = p1sp3_0
- p2
[ 1 ] + p4
[ 1 ];
386 t4
[ 1 ] = p1sp3_1
+ p2
[ 0 ] - p4
[ 0 ];
388 * p1
++ = p1ap3_0
+ p2
[ 0 ] + p4
[ 0 ];
389 * p1
++ = p1ap3_1
+ p2
[ 1 ] + p4
[ 1 ];
392 p1ap3_1
= pEnd1
[- 1 ] + pEnd3
[- 1 ];
393 p1sp3_1
= pEnd1
[- 1 ] - pEnd3
[- 1 ];
394 p1ap3_0
= pEnd1
[ 0 ] + pEnd3
[ 0 ];
395 p1sp3_0
= pEnd1
[ 0 ] - pEnd3
[ 0 ];
397 t2
[ 2 ] = pEnd2
[ 0 ] - pEnd4
[ 0 ] + p1sp3_1
;
398 t2
[ 3 ] = pEnd1
[ 0 ] - pEnd3
[ 0 ] - pEnd2
[- 1 ] + pEnd4
[- 1 ];
400 t3
[ 2 ] = p1ap3_1
- pEnd2
[- 1 ] - pEnd4
[- 1 ];
401 t3
[ 3 ] = p1ap3_0
- pEnd2
[ 0 ] - pEnd4
[ 0 ];
403 t4
[ 2 ] = pEnd2
[ 0 ] - pEnd4
[ 0 ] - p1sp3_1
;
404 t4
[ 3 ] = pEnd4
[- 1 ] - pEnd2
[- 1 ] - p1sp3_0
;
406 * pEnd1
-- = p1ap3_0
+ pEnd2
[ 0 ] + pEnd4
[ 0 ];
407 * pEnd1
-- = p1ap3_1
+ pEnd2
[- 1 ] + pEnd4
[- 1 ];
410 // read twiddle factors
413 // multiply by twiddle factors
414 // let Z1 = a + i(b), Z2 = c + i(d)
415 // => Z1 * Z2 = (a*c - b*d) + i(b*c + a*d)
424 // use vertical symmetry col 2
425 // 0.9997 - 0.0245i <==> 0.0245 - 0.9997i
447 // use vertical symmetry col 3
448 // 0.9988 - 0.0491i <==> -0.9988 - 0.0491i
470 // use vertical symmetry col 4
471 // 0.9973 - 0.0736i <==> -0.0736 + 0.9973i
483 // Twiddle factors are
484 // 1.0000 0.7071-0.7071i -1.0000i -0.7071-0.7071i
485 p1ap3_0
= p1
[ 0 ] + p3
[ 0 ];
486 p1sp3_0
= p1
[ 0 ] - p3
[ 0 ];
487 p1ap3_1
= p1
[ 1 ] + p3
[ 1 ];
488 p1sp3_1
= p1
[ 1 ] - p3
[ 1 ];
491 t2
[ 0 ] = p1sp3_0
+ p2
[ 1 ] - p4
[ 1 ];
492 t2
[ 1 ] = p1sp3_1
- p2
[ 0 ] + p4
[ 0 ];
494 t3
[ 0 ] = p1ap3_0
- p2
[ 0 ] - p4
[ 0 ];
495 t3
[ 1 ] = p1ap3_1
- p2
[ 1 ] - p4
[ 1 ];
497 t4
[ 0 ] = p1sp3_0
- p2
[ 1 ] + p4
[ 1 ];
498 t4
[ 1 ] = p1sp3_1
+ p2
[ 0 ] - p4
[ 0 ];
500 * p1
++ = p1ap3_0
+ p2
[ 0 ] + p4
[ 0 ];
501 * p1
++ = p1ap3_1
+ p2
[ 1 ] + p4
[ 1 ];
538 arm_radix8_butterfly_f32 ( pCol1
, L
, ( float32_t
*) S
-> pTwiddle
, 4u );
540 arm_radix8_butterfly_f32 ( pCol2
, L
, ( float32_t
*) S
-> pTwiddle
, 4u );
542 arm_radix8_butterfly_f32 ( pCol3
, L
, ( float32_t
*) S
-> pTwiddle
, 4u );
544 arm_radix8_butterfly_f32 ( pCol4
, L
, ( float32_t
*) S
-> pTwiddle
, 4u );
549 * @addtogroup ComplexFFT
555 * @brief Processing function for the floating-point complex FFT.
556 * @param[in] *S points to an instance of the floating-point CFFT structure.
557 * @param[in, out] *p1 points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.
558 * @param[in] ifftFlag flag that selects forward (ifftFlag=0) or inverse (ifftFlag=1) transform.
559 * @param[in] bitReverseFlag flag that enables (bitReverseFlag=1) or disables (bitReverseFlag=0) bit reversal of output.
564 const arm_cfft_instance_f32
* S
,
567 uint8_t bitReverseFlag
)
570 uint32_t L
= S
-> fftLen
, l
;
571 float32_t invL
, * pSrc
;
575 /* Conjugate input data */
587 arm_cfft_radix8by2_f32 ( ( arm_cfft_instance_f32
*) S
, p1
);
592 arm_cfft_radix8by4_f32 ( ( arm_cfft_instance_f32
*) S
, p1
);
597 arm_radix8_butterfly_f32 ( p1
, L
, ( float32_t
*) S
-> pTwiddle
, 1 );
602 arm_bitreversal_32 (( uint32_t *) p1
, S
-> bitRevLength
, S
-> pBitRevTable
);
606 invL
= 1.0 f
/( float32_t
) L
;
607 /* Conjugate and scale output data */
611 * pSrc
= -(* pSrc
) * invL
;