1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_rfft_f32.c
4 * Description: RFFT & RIFFT Floating point process function
5 *
6 * $Date: 27. January 2017
7 * $Revision: V.1.5.1
8 *
9 * Target Processor: Cortex-M cores
10 * -------------------------------------------------------------------- */
11 /*
12 * Copyright (C) 2010-2017 ARM Limited or its affiliates. All rights reserved.
13 *
14 * SPDX-License-Identifier: Apache-2.0
15 *
16 * Licensed under the Apache License, Version 2.0 (the License); you may
17 * not use this file except in compliance with the License.
18 * You may obtain a copy of the License at
19 *
20 * www.apache.org/licenses/LICENSE-2.0
21 *
22 * Unless required by applicable law or agreed to in writing, software
23 * distributed under the License is distributed on an AS IS BASIS, WITHOUT
24 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
25 * See the License for the specific language governing permissions and
26 * limitations under the License.
27 */
28
29 #include "arm_math.h"
30
stage_rfft_f32(arm_rfft_fast_instance_f32 * S,float32_t * p,float32_t * pOut)31 void stage_rfft_f32(
32 arm_rfft_fast_instance_f32 * S,
33 float32_t * p, float32_t * pOut)
34 {
35 uint32_t k; /* Loop Counter */
36 float32_t twR, twI; /* RFFT Twiddle coefficients */
37 float32_t * pCoeff = S->pTwiddleRFFT; /* Points to RFFT Twiddle factors */
38 float32_t *pA = p; /* increasing pointer */
39 float32_t *pB = p; /* decreasing pointer */
40 float32_t xAR, xAI, xBR, xBI; /* temporary variables */
41 float32_t t1a, t1b; /* temporary variables */
42 float32_t p0, p1, p2, p3; /* temporary variables */
43
44
45 k = (S->Sint).fftLen - 1;
46
47 /* Pack first and last sample of the frequency domain together */
48
49 xBR = pB[0];
50 xBI = pB[1];
51 xAR = pA[0];
52 xAI = pA[1];
53
54 twR = *pCoeff++ ;
55 twI = *pCoeff++ ;
56
57 // U1 = XA(1) + XB(1); % It is real
58 t1a = xBR + xAR ;
59
60 // U2 = XB(1) - XA(1); % It is imaginary
61 t1b = xBI + xAI ;
62
63 // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
64 // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
65 *pOut++ = 0.5f * ( t1a + t1b );
66 *pOut++ = 0.5f * ( t1a - t1b );
67
68 // XA(1) = 1/2*( U1 - imag(U2) + i*( U1 +imag(U2) ));
69 pB = p + 2*k;
70 pA += 2;
71
72 do
73 {
74 /*
75 function X = my_split_rfft(X, ifftFlag)
76 % X is a series of real numbers
77 L = length(X);
78 XC = X(1:2:end) +i*X(2:2:end);
79 XA = fft(XC);
80 XB = conj(XA([1 end:-1:2]));
81 TW = i*exp(-2*pi*i*[0:L/2-1]/L).';
82 for l = 2:L/2
83 XA(l) = 1/2 * (XA(l) + XB(l) + TW(l) * (XB(l) - XA(l)));
84 end
85 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))));
86 X = XA;
87 */
88
89 xBI = pB[1];
90 xBR = pB[0];
91 xAR = pA[0];
92 xAI = pA[1];
93
94 twR = *pCoeff++;
95 twI = *pCoeff++;
96
97 t1a = xBR - xAR ;
98 t1b = xBI + xAI ;
99
100 // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
101 // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
102 p0 = twR * t1a;
103 p1 = twI * t1a;
104 p2 = twR * t1b;
105 p3 = twI * t1b;
106
107 *pOut++ = 0.5f * (xAR + xBR + p0 + p3 ); //xAR
108 *pOut++ = 0.5f * (xAI - xBI + p1 - p2 ); //xAI
109
110 pA += 2;
111 pB -= 2;
112 k--;
113 } while (k > 0u);
114 }
115
116 /* Prepares data for inverse cfft */
merge_rfft_f32(arm_rfft_fast_instance_f32 * S,float32_t * p,float32_t * pOut)117 void merge_rfft_f32(
118 arm_rfft_fast_instance_f32 * S,
119 float32_t * p, float32_t * pOut)
120 {
121 uint32_t k; /* Loop Counter */
122 float32_t twR, twI; /* RFFT Twiddle coefficients */
123 float32_t *pCoeff = S->pTwiddleRFFT; /* Points to RFFT Twiddle factors */
124 float32_t *pA = p; /* increasing pointer */
125 float32_t *pB = p; /* decreasing pointer */
126 float32_t xAR, xAI, xBR, xBI; /* temporary variables */
127 float32_t t1a, t1b, r, s, t, u; /* temporary variables */
128
129 k = (S->Sint).fftLen - 1;
130
131 xAR = pA[0];
132 xAI = pA[1];
133
134 pCoeff += 2 ;
135
136 *pOut++ = 0.5f * ( xAR + xAI );
137 *pOut++ = 0.5f * ( xAR - xAI );
138
139 pB = p + 2*k ;
140 pA += 2 ;
141
142 while (k > 0u)
143 {
144 /* G is half of the frequency complex spectrum */
145 //for k = 2:N
146 // Xk(k) = 1/2 * (G(k) + conj(G(N-k+2)) + Tw(k)*( G(k) - conj(G(N-k+2))));
147 xBI = pB[1] ;
148 xBR = pB[0] ;
149 xAR = pA[0];
150 xAI = pA[1];
151
152 twR = *pCoeff++;
153 twI = *pCoeff++;
154
155 t1a = xAR - xBR ;
156 t1b = xAI + xBI ;
157
158 r = twR * t1a;
159 s = twI * t1b;
160 t = twI * t1a;
161 u = twR * t1b;
162
163 // real(tw * (xA - xB)) = twR * (xAR - xBR) - twI * (xAI - xBI);
164 // imag(tw * (xA - xB)) = twI * (xAR - xBR) + twR * (xAI - xBI);
165 *pOut++ = 0.5f * (xAR + xBR - r - s ); //xAR
166 *pOut++ = 0.5f * (xAI - xBI + t - u ); //xAI
167
168 pA += 2;
169 pB -= 2;
170 k--;
171 }
172
173 }
174
175 /**
176 * @ingroup groupTransforms
177 */
178
179 /**
180 * @defgroup RealFFT Real FFT Functions
181 *
182 * \par
183 * The CMSIS DSP library includes specialized algorithms for computing the
184 * FFT of real data sequences. The FFT is defined over complex data but
185 * in many applications the input is real. Real FFT algorithms take advantage
186 * of the symmetry properties of the FFT and have a speed advantage over complex
187 * algorithms of the same length.
188 * \par
189 * The Fast RFFT algorith relays on the mixed radix CFFT that save processor usage.
190 * \par
191 * The real length N forward FFT of a sequence is computed using the steps shown below.
192 * \par
193 * \image html RFFT.gif "Real Fast Fourier Transform"
194 * \par
195 * The real sequence is initially treated as if it were complex to perform a CFFT.
196 * Later, a processing stage reshapes the data to obtain half of the frequency spectrum
197 * in complex format. Except the first complex number that contains the two real numbers
198 * X[0] and X[N/2] all the data is complex. In other words, the first complex sample
199 * contains two real values packed.
200 * \par
201 * The input for the inverse RFFT should keep the same format as the output of the
202 * forward RFFT. A first processing stage pre-process the data to later perform an
203 * inverse CFFT.
204 * \par
205 * \image html RIFFT.gif "Real Inverse Fast Fourier Transform"
206 * \par
207 * The algorithms for floating-point, Q15, and Q31 data are slightly different
208 * and we describe each algorithm in turn.
209 * \par Floating-point
210 * The main functions are arm_rfft_fast_f32() and arm_rfft_fast_init_f32().
211 * The older functions arm_rfft_f32() and arm_rfft_init_f32() have been
212 * deprecated but are still documented.
213 * \par
214 * The FFT of a real N-point sequence has even symmetry in the frequency
215 * domain. The second half of the data equals the conjugate of the first
216 * half flipped in frequency. Looking at the data, we see that we can
217 * uniquely represent the FFT using only N/2 complex numbers. These are
218 * packed into the output array in alternating real and imaginary
219 * components:
220 * \par
221 * X = { real[0], imag[0], real[1], imag[1], real[2], imag[2] ...
222 * real[(N/2)-1], imag[(N/2)-1 }
223 * \par
224 * It happens that the first complex number (real[0], imag[0]) is actually
225 * all real. real[0] represents the DC offset, and imag[0] should be 0.
226 * (real[1], imag[1]) is the fundamental frequency, (real[2], imag[2]) is
227 * the first harmonic and so on.
228 * \par
229 * The real FFT functions pack the frequency domain data in this fashion.
230 * The forward transform outputs the data in this form and the inverse
231 * transform expects input data in this form. The function always performs
232 * the needed bitreversal so that the input and output data is always in
233 * normal order. The functions support lengths of [32, 64, 128, ..., 4096]
234 * samples.
235 * \par Q15 and Q31
236 * The real algorithms are defined in a similar manner and utilize N/2 complex
237 * transforms behind the scenes.
238 * \par
239 * The complex transforms used internally include scaling to prevent fixed-point
240 * overflows. The overall scaling equals 1/(fftLen/2).
241 * \par
242 * A separate instance structure must be defined for each transform used but
243 * twiddle factor and bit reversal tables can be reused.
244 * \par
245 * There is also an associated initialization function for each data type.
246 * The initialization function performs the following operations:
247 * - Sets the values of the internal structure fields.
248 * - Initializes twiddle factor table and bit reversal table pointers.
249 * - Initializes the internal complex FFT data structure.
250 * \par
251 * Use of the initialization function is optional.
252 * However, if the initialization function is used, then the instance structure
253 * cannot be placed into a const data section. To place an instance structure
254 * into a const data section, the instance structure should be manually
255 * initialized as follows:
256 * <pre>
257 *arm_rfft_instance_q31 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
258 *arm_rfft_instance_q15 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};
259 * </pre>
260 * where <code>fftLenReal</code> is the length of the real transform;
261 * <code>fftLenBy2</code> length of the internal complex transform.
262 * <code>ifftFlagR</code> Selects forward (=0) or inverse (=1) transform.
263 * <code>bitReverseFlagR</code> Selects bit reversed output (=0) or normal order
264 * output (=1).
265 * <code>twidCoefRModifier</code> stride modifier for the twiddle factor table.
266 * The value is based on the FFT length;
267 * <code>pTwiddleAReal</code>points to the A array of twiddle coefficients;
268 * <code>pTwiddleBReal</code>points to the B array of twiddle coefficients;
269 * <code>pCfft</code> points to the CFFT Instance structure. The CFFT structure
270 * must also be initialized. Refer to arm_cfft_radix4_f32() for details regarding
271 * static initialization of the complex FFT instance structure.
272 */
273
274 /**
275 * @addtogroup RealFFT
276 * @{
277 */
278
279 /**
280 * @brief Processing function for the floating-point real FFT.
281 * @param[in] *S points to an arm_rfft_fast_instance_f32 structure.
282 * @param[in] *p points to the input buffer.
283 * @param[in] *pOut points to the output buffer.
284 * @param[in] ifftFlag RFFT if flag is 0, RIFFT if flag is 1
285 * @return none.
286 */
287
arm_rfft_fast_f32(arm_rfft_fast_instance_f32 * S,float32_t * p,float32_t * pOut,uint8_t ifftFlag)288 void arm_rfft_fast_f32(
289 arm_rfft_fast_instance_f32 * S,
290 float32_t * p, float32_t * pOut,
291 uint8_t ifftFlag)
292 {
293 arm_cfft_instance_f32 * Sint = &(S->Sint);
294 Sint->fftLen = S->fftLenRFFT / 2;
295
296 /* Calculation of Real FFT */
297 if (ifftFlag)
298 {
299 /* Real FFT compression */
300 merge_rfft_f32(S, p, pOut);
301
302 /* Complex radix-4 IFFT process */
303 arm_cfft_f32( Sint, pOut, ifftFlag, 1);
304 }
305 else
306 {
307 /* Calculation of RFFT of input */
308 arm_cfft_f32( Sint, p, ifftFlag, 1);
309
310 /* Real FFT extraction */
311 stage_rfft_f32(S, p, pOut);
312 }
313 }
314
315 /**
316 * @} end of RealFFT group
317 */
318