1 /* ----------------------------------------------------------------------
2 * Copyright (C) 2010-2013 ARM Limited. All rights reserved.
4 * $Date: 17. January 2013
7 * Project: CMSIS DSP Library
8 * Title: arm_rfft_f32.c
10 * Description: RFFT & RIFFT Floating point process 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 * -------------------------------------------------------------------- */
44 arm_rfft_fast_instance_f32
* S
,
45 float32_t
* p
, float32_t
* pOut
)
47 uint32_t k
; /* Loop Counter */
48 float32_t twR
, twI
; /* RFFT Twiddle coefficients */
49 float32_t
* pCoeff
= S
->pTwiddleRFFT
; /* Points to RFFT Twiddle factors */
50 float32_t
*pA
= p
; /* increasing pointer */
51 float32_t
*pB
= p
; /* decreasing pointer */
52 float32_t xAR
, xAI
, xBR
, xBI
; /* temporary variables */
53 float32_t t1a
, t1b
; /* temporary variables */
54 float32_t p0
, p1
, p2
, p3
; /* temporary variables */
57 k
= (S
->Sint
).fftLen
- 1;
59 /* Pack first and last sample of the frequency domain together */
69 // U1 = XA(1) + XB(1); % It is real
72 // U2 = XB(1) - XA(1); % It is imaginary
75 // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
76 // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
77 *pOut
++ = 0.5f
* ( t1a
+ t1b
);
78 *pOut
++ = 0.5f
* ( t1a
- t1b
);
80 // XA(1) = 1/2*( U1 - imag(U2) + i*( U1 +imag(U2) ));
87 function X = my_split_rfft(X, ifftFlag)
88 % X is a series of real numbers
90 XC = X(1:2:end) +i*X(2:2:end);
92 XB = conj(XA([1 end:-1:2]));
93 TW = i*exp(-2*pi*i*[0:L/2-1]/L).';
95 XA(l) = 1/2 * (XA(l) + XB(l) + TW(l) * (XB(l) - XA(l)));
97 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))));
112 // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
113 // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
119 *pOut
++ = 0.5f
* (xAR
+ xBR
+ p0
+ p3
); //xAR
120 *pOut
++ = 0.5f
* (xAI
- xBI
+ p1
- p2
); //xAI
128 /* Prepares data for inverse cfft */
130 arm_rfft_fast_instance_f32
* S
,
131 float32_t
* p
, float32_t
* pOut
)
133 uint32_t k
; /* Loop Counter */
134 float32_t twR
, twI
; /* RFFT Twiddle coefficients */
135 float32_t
*pCoeff
= S
->pTwiddleRFFT
; /* Points to RFFT Twiddle factors */
136 float32_t
*pA
= p
; /* increasing pointer */
137 float32_t
*pB
= p
; /* decreasing pointer */
138 float32_t xAR
, xAI
, xBR
, xBI
; /* temporary variables */
139 float32_t t1a
, t1b
, r
, s
, t
, u
; /* temporary variables */
141 k
= (S
->Sint
).fftLen
- 1;
148 *pOut
++ = 0.5f
* ( xAR
+ xAI
);
149 *pOut
++ = 0.5f
* ( xAR
- xAI
);
156 /* G is half of the frequency complex spectrum */
158 // Xk(k) = 1/2 * (G(k) + conj(G(N-k+2)) + Tw(k)*( G(k) - conj(G(N-k+2))));
175 // real(tw * (xA - xB)) = twR * (xAR - xBR) - twI * (xAI - xBI);
176 // imag(tw * (xA - xB)) = twI * (xAR - xBR) + twR * (xAI - xBI);
177 *pOut
++ = 0.5f
* (xAR
+ xBR
- r
- s
); //xAR
178 *pOut
++ = 0.5f
* (xAI
- xBI
+ t
- u
); //xAI
188 * @ingroup groupTransforms
192 * @defgroup Fast Real FFT Functions
195 * The CMSIS DSP library includes specialized algorithms for computing the
196 * FFT of real data sequences. The FFT is defined over complex data but
197 * in many applications the input is real. Real FFT algorithms take advantage
198 * of the symmetry properties of the FFT and have a speed advantage over complex
199 * algorithms of the same length.
201 * The Fast RFFT algorith relays on the mixed radix CFFT that save processor usage.
203 * The real length N forward FFT of a sequence is computed using the steps shown below.
205 * \image html RFFT.gif "Real Fast Fourier Transform"
207 * The real sequence is initially treated as if it were complex to perform a CFFT.
208 * Later, a processing stage reshapes the data to obtain half of the frequency spectrum
209 * in complex format. Except the first complex number that contains the two real numbers
210 * X[0] and X[N/2] all the data is complex. In other words, the first complex sample
211 * contains two real values packed.
213 * The input for the inverse RFFT should keep the same format as the output of the
214 * forward RFFT. A first processing stage pre-process the data to later perform an
217 * \image html RIFFT.gif "Real Inverse Fast Fourier Transform"
219 * The algorithms for floating-point, Q15, and Q31 data are slightly different
220 * and we describe each algorithm in turn.
221 * \par Floating-point
222 * The main functions are <code>arm_rfft_fast_f32()</code>
223 * and <code>arm_rfft_fast_init_f32()</code>. The older functions
224 * <code>arm_rfft_f32()</code> and <code>arm_rfft_init_f32()</code> have been
225 * deprecated but are still documented.
227 * The FFT of a real N-point sequence has even symmetry in the frequency
228 * domain. The second half of the data equals the conjugate of the first half
229 * flipped in frequency:
235 *X[fftLen/2-1] - complex data
236 *X[fftLen/2] - real data
237 *X[fftLen/2+1] - conjugate of X[fftLen/2-1]
238 *X[fftLen/2+2] - conjugate of X[fftLen/2-2]
240 *X[fftLen-1] - conjugate of X[1]
242 * Looking at the data, we see that we can uniquely represent the FFT using only
249 *X[fftLen/2-1] - complex data
250 *X[fftLen/2] - real data
252 * Looking more closely we see that the first and last samples are real valued.
253 * They can be packed together and we can thus represent the FFT of an N-point
254 * real sequence by N/2 complex values:
256 *X[0],X[N/2] - packed real data: X[0] + jX[N/2]
260 *X[fftLen/2-1] - complex data
262 * The real FFT functions pack the frequency domain data in this fashion. The
263 * forward transform outputs the data in this form and the inverse transform
264 * expects input data in this form. The function always performs the needed
265 * bitreversal so that the input and output data is always in normal order. The
266 * functions support lengths of [32, 64, 128, ..., 4096] samples.
268 * The forward and inverse real FFT functions apply the standard FFT scaling; no
269 * scaling on the forward transform and 1/fftLen scaling on the inverse
272 * The real algorithms are defined in a similar manner and utilize N/2 complex
273 * transforms behind the scenes. In the case of fixed-point data, a radix-4
274 * complex transform is performed and this limits the allows sequence lengths to
275 * 128, 512, and 2048 samples.
277 * TBD. We need to document input and output order of data.
279 * The complex transforms used internally include scaling to prevent fixed-point
280 * overflows. The overall scaling equals 1/(fftLen/2).
282 * A separate instance structure must be defined for each transform used but
283 * twiddle factor and bit reversal tables can be reused.
285 * There is also an associated initialization function for each data type.
286 * The initialization function performs the following operations:
287 * - Sets the values of the internal structure fields.
288 * - Initializes twiddle factor table and bit reversal table pointers.
289 * - Initializes the internal complex FFT data structure.
291 * Use of the initialization function is optional.
292 * However, if the initialization function is used, then the instance structure
293 * cannot be placed into a const data section. To place an instance structure
294 * into a const data section, the instance structure should be manually
295 * initialized as follows:
297 *arm_rfft_instance_q31 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
298 *arm_rfft_instance_q15 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
300 * where <code>fftLenReal</code> is the length of the real transform;
301 * <code>fftLenBy2</code> length of the internal complex transform.
302 * <code>ifftFlagR</code> Selects forward (=0) or inverse (=1) transform.
303 * <code>bitReverseFlagR</code> Selects bit reversed output (=0) or normal order
305 * <code>twidCoefRModifier</code> stride modifier for the twiddle factor table.
306 * The value is based on the FFT length;
307 * <code>pTwiddleAReal</code>points to the A array of twiddle coefficients;
308 * <code>pTwiddleBReal</code>points to the B array of twiddle coefficients;
309 * <code>pCfft</code> points to the CFFT Instance structure. The CFFT structure
310 * must also be initialized. Refer to arm_cfft_radix4_f32() for details regarding
311 * static initialization of the complex FFT instance structure.
315 * @addtogroup RealFFT
320 * @brief Processing function for the floating-point real FFT.
321 * @param[in] *S points to an arm_rfft_fast_instance_f32 structure.
322 * @param[in] *p points to the input buffer.
323 * @param[in] *pOut points to an arm_rfft_fast_instance_f32 structure.
324 * @param[in] ifftFlag RFFT if flag is 0, RIFFT if flag is 1
328 void arm_rfft_fast_f32(
329 arm_rfft_fast_instance_f32
* S
,
330 float32_t
* p
, float32_t
* pOut
,
333 arm_cfft_instance_f32
* Sint
= &(S
->Sint
);
334 Sint
->fftLen
= S
->fftLenRFFT
/ 2;
336 /* Calculation of Real FFT */
339 /* Real FFT comression */
340 merge_rfft_f32(S
, p
, pOut
);
342 /* Complex radix-4 IFFT process */
343 arm_cfft_f32( Sint
, pOut
, ifftFlag
, 1);
347 /* Calculation of RFFT of input */
348 arm_cfft_f32( Sint
, p
, ifftFlag
, 1);
350 /* Real FFT extraction */
351 stage_rfft_f32(S
, p
, pOut
);