1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_rfft_fast_f64.c
4 * Description: RFFT & RIFFT Double precision Floating point process function
5 *
6 * $Date: 23 April 2021
7 * $Revision: V1.9.0
8 *
9 * Target Processor: Cortex-M and Cortex-A cores
10 * -------------------------------------------------------------------- */
11 /*
12 * Copyright (C) 2010-2021 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 "dsp/transform_functions.h"
30
stage_rfft_f64(const arm_rfft_fast_instance_f64 * S,float64_t * p,float64_t * pOut)31 void stage_rfft_f64(
32 const arm_rfft_fast_instance_f64 * S,
33 float64_t * p,
34 float64_t * pOut)
35 {
36 uint32_t k; /* Loop Counter */
37 float64_t twR, twI; /* RFFT Twiddle coefficients */
38 const float64_t * pCoeff = S->pTwiddleRFFT; /* Points to RFFT Twiddle factors */
39 float64_t *pA = p; /* increasing pointer */
40 float64_t *pB = p; /* decreasing pointer */
41 float64_t xAR, xAI, xBR, xBI; /* temporary variables */
42 float64_t t1a, t1b; /* temporary variables */
43 float64_t p0, p1, p2, p3; /* temporary variables */
44
45
46 k = (S->Sint).fftLen - 1;
47
48 /* Pack first and last sample of the frequency domain together */
49
50 xBR = pB[0];
51 xBI = pB[1];
52 xAR = pA[0];
53 xAI = pA[1];
54
55 twR = *pCoeff++ ;
56 twI = *pCoeff++ ;
57
58 // U1 = XA(1) + XB(1); % It is real
59 t1a = xBR + xAR ;
60
61 // U2 = XB(1) - XA(1); % It is imaginary
62 t1b = xBI + xAI ;
63
64 // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
65 // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
66 *pOut++ = 0.5 * ( t1a + t1b );
67 *pOut++ = 0.5 * ( t1a - t1b );
68
69 // XA(1) = 1/2*( U1 - imag(U2) + i*( U1 +imag(U2) ));
70 pB = p + 2*k;
71 pA += 2;
72
73 do
74 {
75 /*
76 function X = my_split_rfft(X, ifftFlag)
77 % X is a series of real numbers
78 L = length(X);
79 XC = X(1:2:end) +i*X(2:2:end);
80 XA = fft(XC);
81 XB = conj(XA([1 end:-1:2]));
82 TW = i*exp(-2*pi*i*[0:L/2-1]/L).';
83 for l = 2:L/2
84 XA(l) = 1/2 * (XA(l) + XB(l) + TW(l) * (XB(l) - XA(l)));
85 end
86 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))));
87 X = XA;
88 */
89
90 xBI = pB[1];
91 xBR = pB[0];
92 xAR = pA[0];
93 xAI = pA[1];
94
95 twR = *pCoeff++;
96 twI = *pCoeff++;
97
98 t1a = xBR - xAR ;
99 t1b = xBI + xAI ;
100
101 // real(tw * (xB - xA)) = twR * (xBR - xAR) - twI * (xBI - xAI);
102 // imag(tw * (xB - xA)) = twI * (xBR - xAR) + twR * (xBI - xAI);
103 p0 = twR * t1a;
104 p1 = twI * t1a;
105 p2 = twR * t1b;
106 p3 = twI * t1b;
107
108 *pOut++ = 0.5 * (xAR + xBR + p0 + p3 ); //xAR
109 *pOut++ = 0.5 * (xAI - xBI + p1 - p2 ); //xAI
110
111 pA += 2;
112 pB -= 2;
113 k--;
114 } while (k > 0U);
115 }
116
117 /* Prepares data for inverse cfft */
merge_rfft_f64(const arm_rfft_fast_instance_f64 * S,float64_t * p,float64_t * pOut)118 void merge_rfft_f64(
119 const arm_rfft_fast_instance_f64 * S,
120 float64_t * p,
121 float64_t * pOut)
122 {
123 uint32_t k; /* Loop Counter */
124 float64_t twR, twI; /* RFFT Twiddle coefficients */
125 const float64_t *pCoeff = S->pTwiddleRFFT; /* Points to RFFT Twiddle factors */
126 float64_t *pA = p; /* increasing pointer */
127 float64_t *pB = p; /* decreasing pointer */
128 float64_t xAR, xAI, xBR, xBI; /* temporary variables */
129 float64_t t1a, t1b, r, s, t, u; /* temporary variables */
130
131 k = (S->Sint).fftLen - 1;
132
133 xAR = pA[0];
134 xAI = pA[1];
135
136 pCoeff += 2 ;
137
138 *pOut++ = 0.5 * ( xAR + xAI );
139 *pOut++ = 0.5 * ( xAR - xAI );
140
141 pB = p + 2*k ;
142 pA += 2 ;
143
144 while (k > 0U)
145 {
146 /* G is half of the frequency complex spectrum */
147 //for k = 2:N
148 // Xk(k) = 1/2 * (G(k) + conj(G(N-k+2)) + Tw(k)*( G(k) - conj(G(N-k+2))));
149 xBI = pB[1] ;
150 xBR = pB[0] ;
151 xAR = pA[0];
152 xAI = pA[1];
153
154 twR = *pCoeff++;
155 twI = *pCoeff++;
156
157 t1a = xAR - xBR ;
158 t1b = xAI + xBI ;
159
160 r = twR * t1a;
161 s = twI * t1b;
162 t = twI * t1a;
163 u = twR * t1b;
164
165 // real(tw * (xA - xB)) = twR * (xAR - xBR) - twI * (xAI - xBI);
166 // imag(tw * (xA - xB)) = twI * (xAR - xBR) + twR * (xAI - xBI);
167 *pOut++ = 0.5 * (xAR + xBR - r - s ); //xAR
168 *pOut++ = 0.5 * (xAI - xBI + t - u ); //xAI
169
170 pA += 2;
171 pB -= 2;
172 k--;
173 }
174
175 }
176
177 /**
178 @ingroup groupTransforms
179 */
180
181
182 /**
183 @addtogroup RealFFT
184 @{
185 */
186
187 /**
188 @brief Processing function for the Double Precision floating-point real FFT.
189 @param[in] S points to an arm_rfft_fast_instance_f64 structure
190 @param[in] p points to input buffer (Source buffer is modified by this function.)
191 @param[in] pOut points to output buffer
192 @param[in] ifftFlag
193 - value = 0: RFFT
194 - value = 1: RIFFT
195 @return none
196 */
197
arm_rfft_fast_f64(arm_rfft_fast_instance_f64 * S,float64_t * p,float64_t * pOut,uint8_t ifftFlag)198 void arm_rfft_fast_f64(
199 arm_rfft_fast_instance_f64 * S,
200 float64_t * p,
201 float64_t * pOut,
202 uint8_t ifftFlag)
203 {
204 arm_cfft_instance_f64 * Sint = &(S->Sint);
205 Sint->fftLen = S->fftLenRFFT / 2;
206
207 /* Calculation of Real FFT */
208 if (ifftFlag)
209 {
210 /* Real FFT compression */
211 merge_rfft_f64(S, p, pOut);
212
213 /* Complex radix-4 IFFT process */
214 arm_cfft_f64( Sint, pOut, ifftFlag, 1);
215 }
216 else
217 {
218 /* Calculation of RFFT of input */
219 arm_cfft_f64( Sint, p, ifftFlag, 1);
220
221 /* Real FFT extraction */
222 stage_rfft_f64(S, p, pOut);
223 }
224 }
225
226 /**
227 * @} end of RealFFT group
228 */
229
