查看 | 历史 | 提交 | 提交对比
QuakeGod
2024-07-27 842bb64195f958b050867c50db66fc0aa413dafb
提交 | 用户 | age
483170 1 /* ----------------------------------------------------------------------
Q 2 * Copyright (C) 2010-2014 ARM Limited. All rights reserved.
3 *
4 * $Date:        19. March 2015
5 * $Revision:     V.1.4.5
6 *
7 * Project:         CMSIS DSP Library
8 * Title:        arm_rfft_f32.c
9 *
10 * Description:    RFFT & RIFFT Floating point process function
11 *
12 * Target Processor: Cortex-M4/Cortex-M3/Cortex-M0
13 *
14 * Redistribution and use in source and binary forms, with or without
15 * modification, are permitted provided that the following conditions
16 * are met:
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
22 *     distribution.
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.
26 *
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 * -------------------------------------------------------------------- */
40
41 #include "arm_math.h"
42
43 void stage_rfft_f32(
44   arm_rfft_fast_instance_f32 * S,
45   float32_t * p, float32_t * pOut)
46 {
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              */
55
56
57    k = (S->Sint).fftLen - 1;                    
58
59    /* Pack first and last sample of the frequency domain together */
60
61    xBR = pB[0];
62    xBI = pB[1];
63    xAR = pA[0];
64    xAI = pA[1];
65
66    twR = *pCoeff++ ;
67    twI = *pCoeff++ ;
68    
69    // U1 = XA(1) + XB(1); % It is real
70    t1a = xBR + xAR  ;
71    
72    // U2 = XB(1) - XA(1); % It is imaginary
73    t1b = xBI + xAI  ;
74
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 );
79
80    // XA(1) = 1/2*( U1 - imag(U2) +  i*( U1 +imag(U2) ));
81    pB  = p + 2*k;
82    pA += 2;
83
84    do
85    {
86       /*
87          function X = my_split_rfft(X, ifftFlag)
88          % X is a series of real numbers
89          L  = length(X);
90          XC = X(1:2:end) +i*X(2:2:end);
91          XA = fft(XC);
92          XB = conj(XA([1 end:-1:2]));
93          TW = i*exp(-2*pi*i*[0:L/2-1]/L).';
94          for l = 2:L/2
95             XA(l) = 1/2 * (XA(l) + XB(l) + TW(l) * (XB(l) - XA(l)));
96          end
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))));
98          X = XA;
99       */
100
101       xBI = pB[1];
102       xBR = pB[0];
103       xAR = pA[0];
104       xAI = pA[1];
105
106       twR = *pCoeff++;
107       twI = *pCoeff++;
108
109       t1a = xBR - xAR ;
110       t1b = xBI + xAI ;
111
112       // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
113       // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
114       p0 = twR * t1a;
115       p1 = twI * t1a;
116       p2 = twR * t1b;
117       p3 = twI * t1b;
118
119       *pOut++ = 0.5f * (xAR + xBR + p0 + p3 ); //xAR
120       *pOut++ = 0.5f * (xAI - xBI + p1 - p2 ); //xAI
121
122       pA += 2;
123       pB -= 2;
124       k--;
125    } while(k > 0u);
126 }
127
128 /* Prepares data for inverse cfft */
129 void merge_rfft_f32(
130 arm_rfft_fast_instance_f32 * S,
131 float32_t * p, float32_t * pOut)
132 {
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              */
140
141    k = (S->Sint).fftLen - 1;                    
142
143    xAR = pA[0];
144    xAI = pA[1];
145
146    pCoeff += 2 ;
147
148    *pOut++ = 0.5f * ( xAR + xAI );
149    *pOut++ = 0.5f * ( xAR - xAI );
150
151    pB  =  p + 2*k ;
152    pA +=  2       ;
153
154    while(k > 0u)
155    {
156       /* G is half of the frequency complex spectrum */
157       //for k = 2:N
158       //    Xk(k) = 1/2 * (G(k) + conj(G(N-k+2)) + Tw(k)*( G(k) - conj(G(N-k+2))));
159       xBI =   pB[1]    ;
160       xBR =   pB[0]    ;
161       xAR =  pA[0];
162       xAI =  pA[1];
163
164       twR = *pCoeff++;
165       twI = *pCoeff++;
166
167       t1a = xAR - xBR ;
168       t1b = xAI + xBI ;
169
170       r = twR * t1a;
171       s = twI * t1b;
172       t = twI * t1a;
173       u = twR * t1b;
174
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
179
180       pA += 2;
181       pB -= 2;
182       k--;
183    }
184
185 }
186
187 /**
188 * @ingroup groupTransforms
189 */
190
191 /**
192  * @defgroup Fast Real FFT Functions
193  *
194  * \par
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.
200  * \par
201  * The Fast RFFT algorith relays on the mixed radix CFFT that save processor usage.
202  * \par
203  * The real length N forward FFT of a sequence is computed using the steps shown below.
204  * \par
205  * \image html RFFT.gif "Real Fast Fourier Transform"
206  * \par
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.
212  * \par
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
215  * inverse CFFT.
216  * \par    
217  * \image html RIFFT.gif "Real Inverse Fast Fourier Transform"    
218  * \par    
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.
226  * \par
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:
230  * <pre>
231  *X[0] - real data
232  *X[1] - complex data
233  *X[2] - complex data
234  *... 
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]
239  *... 
240  *X[fftLen-1] - conjugate of X[1]
241  * </pre>
242  * Looking at the data, we see that we can uniquely represent the FFT using only
243  * <pre>
244  *N/2+1 samples:
245  *X[0] - real data
246  *X[1] - complex data
247  *X[2] - complex data
248  *... 
249  *X[fftLen/2-1] - complex data
250  *X[fftLen/2] - real data
251  * </pre>
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:
255  * <pre>
256  *X[0],X[N/2] - packed real data: X[0] + jX[N/2]
257  *X[1] - complex data
258  *X[2] - complex data
259  *... 
260  *X[fftLen/2-1] - complex data
261  * </pre>
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.
267  * \par
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
270  * transform.
271  * \par Q15 and Q31
272  * The real algorithms are defined in a similar manner and utilize N/2 complex
273  * transforms behind the scenes.  
274  * \par
275  * The complex transforms used internally include scaling to prevent fixed-point
276  * overflows.  The overall scaling equals 1/(fftLen/2).
277  * \par
278  * A separate instance structure must be defined for each transform used but 
279  * twiddle factor and bit reversal tables can be reused.
280  * \par
281  * There is also an associated initialization function for each data type. 
282  * The initialization function performs the following operations:
283  * - Sets the values of the internal structure fields.   
284  * - Initializes twiddle factor table and bit reversal table pointers.
285  * - Initializes the internal complex FFT data structure.
286  * \par   
287  * Use of the initialization function is optional.   
288  * However, if the initialization function is used, then the instance structure 
289  * cannot be placed into a const data section. To place an instance structure 
290  * into a const data section, the instance structure should be manually 
291  * initialized as follows:
292  * <pre>
293  *arm_rfft_instance_q31 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};    
294  *arm_rfft_instance_q15 S = {fftLenReal, fftLenBy2, ifftFlagR, bitReverseFlagR, twidCoefRModifier, pTwiddleAReal, pTwiddleBReal, pCfft};    
295  * </pre>
296  * where <code>fftLenReal</code> is the length of the real transform;
297  * <code>fftLenBy2</code> length of  the internal complex transform.
298  * <code>ifftFlagR</code> Selects forward (=0) or inverse (=1) transform.
299  * <code>bitReverseFlagR</code> Selects bit reversed output (=0) or normal order
300  * output (=1).
301  * <code>twidCoefRModifier</code> stride modifier for the twiddle factor table.
302  * The value is based on the FFT length;
303  * <code>pTwiddleAReal</code>points to the A array of twiddle coefficients; 
304  * <code>pTwiddleBReal</code>points to the B array of twiddle coefficients;    
305  * <code>pCfft</code> points to the CFFT Instance structure. The CFFT structure
306  * must also be initialized.  Refer to arm_cfft_radix4_f32() for details regarding    
307  * static initialization of the complex FFT instance structure.    
308  */
309
310 /**
311 * @addtogroup RealFFT
312 * @{
313 */
314
315 /**
316 * @brief Processing function for the floating-point real FFT.
317 * @param[in]  *S              points to an arm_rfft_fast_instance_f32 structure.
318 * @param[in]  *p              points to the input buffer.
319 * @param[in]  *pOut           points to the output buffer.
320 * @param[in]  ifftFlag        RFFT if flag is 0, RIFFT if flag is 1
321 * @return none.
322 */
323
324 void arm_rfft_fast_f32(
325 arm_rfft_fast_instance_f32 * S,
326 float32_t * p, float32_t * pOut,
327 uint8_t ifftFlag)
328 {
329    arm_cfft_instance_f32 * Sint = &(S->Sint);
330    Sint->fftLen = S->fftLenRFFT / 2;
331
332    /* Calculation of Real FFT */
333    if(ifftFlag)
334    {
335       /*  Real FFT compression */
336       merge_rfft_f32(S, p, pOut);
337
338       /* Complex radix-4 IFFT process */
339       arm_cfft_f32( Sint, pOut, ifftFlag, 1);
340    }
341    else
342    {
343       /* Calculation of RFFT of input */
344       arm_cfft_f32( Sint, p, ifftFlag, 1);
345    
346       /*  Real FFT extraction */
347       stage_rfft_f32(S, p, pOut);
348    }
349 }
350
351 /**    
352 * @} end of RealFFT group    
353 */