• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * Copyright (C) 2008-2009 SVOX AG, Baslerstr. 30, 8048 Zuerich, Switzerland
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *     http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 /**
17  * @file picofftsg.c
18  *
19  * FFT/DCT related data types, constants and functions in Pico
20  *
21  * Copyright (C) 2008-2009 SVOX AG, Baslerstr. 30, 8048 Zuerich, Switzerland
22  * All rights reserved.
23  *
24  * History:
25  * - 2009-04-20 -- initial version
26  *
27 */
28 
29 #include "picoos.h"
30 #include "picofftsg.h"
31 #include "picodbg.h"
32 
33 #ifdef __cplusplus
34 extern "C" {
35 #endif
36 
37 /**
38  * @addtogroup picofft
39  * ---------------------------------------------------\n
40  * <b> Fast Fourier/Cosine/Sine Transform </b>\n
41  * Adapted from http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html (Copyright Takuya OOURA, 1996-2001)\n
42  * ---------------------------------------------------\n
43 
44   Overall features
45   - dimension   :one
46   - data length :power of 2
47   - decimation  :frequency
48   - radix       :split-radix
49   - data        :inplace
50   - table       :not use
51 
52   functions
53   - cdft: Complex Discrete Fourier Transform
54   - rdft: Real Discrete Fourier Transform
55   - ddct: Discrete Cosine Transform
56   - ddst: Discrete Sine Transform
57   - dfct: Cosine Transform of RDFT (Real Symmetric DFT)
58   - dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
59 
60   function prototypes
61   - void cdft(picoos_int32, picoos_int32, PICOFFTSG_FFTTYPE *);
62   - void rdft(picoos_int32, picoos_int32, PICOFFTSG_FFTTYPE *);
63   - void ddct(picoos_int32, picoos_int32, PICOFFTSG_FFTTYPE *);
64   - void ddst(picoos_int32, picoos_int32, PICOFFTSG_FFTTYPE *);
65   - void dfct(picoos_int32, PICOFFTSG_FFTTYPE *);
66   - void dfst(picoos_int32, PICOFFTSG_FFTTYPE *);
67 
68   <b>Complex DFT (Discrete Fourier Transform)</b>
69 
70   [definition]
71   - <case1>
72       - X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
73   - <case2>
74       - X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
75   - (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
76 
77   [usage]
78   - <case1>
79       - cdft(2*n, 1, a);
80   - <case2>
81       - cdft(2*n, -1, a);
82 
83   [parameters]
84   - 2*n            :data length (picoos_int32)
85   - n >= 1, n = power of 2
86   - a[0...2*n-1]   :input/output data (PICOFFTSG_FFTTYPE *)
87   - input data
88       - a[2*j] = Re(x[j]),
89       - a[2*j+1] = Im(x[j]), 0<=j<n
90   - output data
91       - a[2*k] = Re(X[k]),
92       - a[2*k+1] = Im(X[k]), 0<=k<n
93 
94   [remark]
95   - Inverse of cdft(2*n, -1, a); is
96       -cdft(2*n, 1, a);
97           - for (j = 0; j <= 2 * n - 1; j++) {
98               - a[j] *= 1.0 / n;
99           - }
100 
101 
102   <b> Real DFT / Inverse of Real DFT </b>
103 
104   [definition]
105   - <case1> RDFT
106       - R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
107       - I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
108   -    <case2> IRDFT (excluding scale)
109       - a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
110       - sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
111       - sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
112 
113   [usage]
114   - <case1>
115       - rdft(n, 1, a);
116   -    <case2>
117       - rdft(n, -1, a);
118 
119   [parameters]
120   - n              :data length (picoos_int32)
121   - n >= 2, n = power of 2
122   - a[0...n-1]     :input/output data (PICOFFTSG_FFTTYPE *)
123   - <case1>
124       - output data
125           - a[2*k] = R[k], 0<=k<n/2
126           - a[2*k+1] = I[k], 0<k<n/2
127           - a[1] = R[n/2]
128   - <case2>
129       - input data
130           - a[2*j] = R[j], 0<=j<n/2
131           - a[2*j+1] = I[j], 0<j<n/2
132           - a[1] = R[n/2]
133 
134   [remark]
135   - Inverse of rdft(n, 1, a); is
136   - rdft(n, -1, a);
137       - for (j = 0; j <= n - 1; j++) {
138           - a[j] *= 2.0 / n;
139      -}
140 
141 
142   <b> DCT (Discrete Cosine Transform) / Inverse of DCT</b>
143 
144   [definition]
145   - <case1> IDCT (excluding scale)
146       - C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
147   - <case2> DCT
148       - C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
149 
150   [usage]
151   - <case1>
152     - ddct(n, 1, a);
153   - <case2>
154     - ddct(n, -1, a);
155 
156   [parameters]
157   - n              :data length (picoos_int32)
158   - n >= 2, n = power of 2
159   - a[0...n-1]     :input/output data (PICOFFTSG_FFTTYPE *)
160   - output data
161       - a[k] = C[k], 0<=k<n
162 
163   [remark]
164   - Inverse of ddct(n, -1, a); is
165   - a[0] *= 0.5;
166   - ddct(n, 1, a);
167   - for (j = 0; j <= n - 1; j++) {
168       - a[j] *= 2.0 / n;
169   - }
170 
171   <b> DST (Discrete Sine Transform) / Inverse of DST</b>
172 
173   [definition]
174   - <case1> IDST (excluding scale)
175       - S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
176   - <case2> DST
177       - S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
178 
179   [usage]
180   - <case1>
181       - ddst(n, 1, a);
182   - <case2>
183       - ddst(n, -1, a);
184 
185   [parameters]
186   - n              :data length (picoos_int32)
187   - n >= 2, n = power of 2
188   - a[0...n-1]     :input/output data (PICOFFTSG_FFTTYPE *)
189   - <case1>
190       - input data
191           - a[j] = A[j], 0<j<n
192           - a[0] = A[n]
193       - output data
194           - a[k] = S[k], 0<=k<n
195   - <case2>
196     - output data
197         - a[k] = S[k], 0<k<n
198         - a[0] = S[n]
199 
200   [remark]
201   - Inverse of ddst(n, -1, a); is
202   - a[0] *= 0.5;
203   - ddst(n, 1, a);
204   - for (j = 0; j <= n - 1; j++) {
205       - a[j] *= 2.0 / n;
206   - }
207 
208   <b> Cosine Transform of RDFT (Real Symmetric DFT)</b>
209 
210   [definition]
211   - C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
212 
213   [usage]
214   - dfct(n, a);
215 
216   [parameters]
217   - n              :data length - 1 (picoos_int32)
218   - n >= 2, n = power of 2
219   - a[0...n]       :input/output data (PICOFFTSG_FFTTYPE *)
220 
221   - output data
222     - a[k] = C[k], 0<=k<=n
223 
224   [remark]
225   - Inverse of a[0] *= 0.5; a[n] *= 0.5; dfct(n, a); is
226   - a[0] *= 0.5;
227   - a[n] *= 0.5;
228   - dfct(n, a);
229   - for (j = 0; j <= n; j++) {
230       - a[j] *= 2.0 / n;
231   - }
232 
233   <b> Sine Transform of RDFT (Real Anti-symmetric DFT)</b>
234 
235   [definition]
236   - S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
237 
238   [usage]
239   - dfst(n, a);
240 
241   [parameters]
242   - n              :data length + 1 (picoos_int32)
243   - n >= 2, n = power of 2
244   - a[0...n-1]     :input/output data (PICOFFTSG_FFTTYPE *)
245   - output data
246       - a[k] = S[k], 0<k<n
247       - (a[0] is used for work area)
248 
249   [remark]
250   - Inverse of dfst(n, a); is
251   - dfst(n, a);
252   - for (j = 1; j <= n - 1; j++) {
253     - a[j] *= 2.0 / n;
254   - }
255 
256 */
257 
258 /* fixed point multiplier for weights */
259 #define PICODSP_WGT_SHIFT  (0x20000000)  /* 2^29 */
260 #define PICOFFTSG_WGT_SHIFT2 (0x10000000)  /* PICODSP_WGT_SHIFT/2 */
261 #define PICOFFTSG_WGT_N_SHIFT (29)  /* 2^29 */
262 /* fixed point known constants */
263 #ifndef WR5000  /* cos(M_PI_2*0.5000) */
264 #define WR5000      (PICOFFTSG_FFTTYPE)(0.707106781186547524400844362104849039284835937688*PICODSP_WGT_SHIFT)
265 #ifndef WR2500  /* cos(M_PI_2*0.2500) */
266 #define WR2500      (PICOFFTSG_FFTTYPE)(0.923879532511286756128183189396788286822416625863*PICODSP_WGT_SHIFT)
267 #endif
268 #ifndef WI2500  /* sin(M_PI_2*0.2500) */
269 #define WI2500      (PICOFFTSG_FFTTYPE)(0.382683432365089771728459984030398866761344562485*PICODSP_WGT_SHIFT)
270 #endif
271 #ifndef WR1250  /* cos(M_PI_2*0.1250) */
272 #define WR1250      (PICOFFTSG_FFTTYPE)(0.980785280403230449126182236134239036973933730893*PICODSP_WGT_SHIFT)
273 #endif
274 #ifndef WI1250  /* sin(M_PI_2*0.1250) */
275 #define WI1250      (PICOFFTSG_FFTTYPE)(0.195090322016128267848284868477022240927691617751*PICODSP_WGT_SHIFT)
276 #endif
277 #ifndef WR3750  /* cos(M_PI_2*0.3750) */
278 #define WR3750      (PICOFFTSG_FFTTYPE)(0.831469612302545237078788377617905756738560811987*PICODSP_WGT_SHIFT)
279 #endif
280 #ifndef WI3750  /* sin(M_PI_2*0.3750) */
281 #define WI3750      (PICOFFTSG_FFTTYPE)(0.555570233019602224742830813948532874374937190754*PICODSP_WGT_SHIFT)
282 #endif
283 
284 #else
285 
286  /*#ifndef M_PI_2
287    #define M_PI_2      (PICOFFTSG_FFTTYPE)1.570796326794896619231321691639751442098584699687
288    #endif*/
289 #ifndef WR5000  /* cos(M_PI_2*0.5000) */
290 #define WR5000      (PICOFFTSG_FFTTYPE)0.707106781186547524400844362104849039284835937688
291 #endif
292 #ifndef WR2500  /* cos(M_PI_2*0.2500) */
293 #define WR2500      (PICOFFTSG_FFTTYPE)0.923879532511286756128183189396788286822416625863
294 #endif
295 #ifndef WI2500  /* sin(M_PI_2*0.2500) */
296 #define WI2500      (PICOFFTSG_FFTTYPE)0.382683432365089771728459984030398866761344562485
297 #endif
298 #ifndef WR1250  /* cos(M_PI_2*0.1250) */
299 #define WR1250      (PICOFFTSG_FFTTYPE)0.980785280403230449126182236134239036973933730893
300 #endif
301 #ifndef WI1250  /* sin(M_PI_2*0.1250) */
302 #define WI1250      (PICOFFTSG_FFTTYPE)0.195090322016128267848284868477022240927691617751
303 #endif
304 #ifndef WR3750  /* cos(M_PI_2*0.3750) */
305 #define WR3750      (PICOFFTSG_FFTTYPE)0.831469612302545237078788377617905756738560811987
306 #endif
307 #ifndef WI3750  /* sin(M_PI_2*0.3750) */
308 #define WI3750      (PICOFFTSG_FFTTYPE)0.555570233019602224742830813948532874374937190754
309 #endif
310 
311 #endif
312 
313 #ifndef CDFT_LOOP_DIV  /* control of the CDFT's speed & tolerance */
314 #define CDFT_LOOP_DIV 32
315 #define CDFT_LOOP_DIV_4 128
316 #endif
317 
318 #ifndef RDFT_LOOP_DIV  /* control of the RDFT's speed & tolerance */
319 #define RDFT_LOOP_DIV 64
320 #define RDFT_LOOP_DIV4 256
321 #define RDFT_LOOP_DIV_4 256
322 #endif
323 
324 #ifndef DCST_LOOP_DIV  /* control of the DCT,DST's speed & tolerance */
325 #define DCST_LOOP_DIV 64
326 #define DCST_LOOP_DIV2 128
327 #endif
328 
329 
330 #define POW1 (0x1)
331 #define POW2 (0x2)
332 #define POW3 (0x4)
333 #define POW4 (0x8)
334 #define POW5 (0x10)
335 #define POW6 (0x20)
336 #define POW7 (0x40)
337 #define POW8 (0x80)
338 #define POW9 (0x100)
339 #define POW10 (0x200)
340 #define POW11 (0x400)
341 #define POW12 (0x800)
342 #define POW13 (0x1000)
343 #define POW14 (0x2000)
344 #define POW15 (0x4000)
345 #define POW16 (0x8000)
346 #define POW17 (0x10000)
347 #define POW18 (0x20000)
348 #define POW19 (0x40000)
349 #define POW20 (0x80000)
350 #define POW21 (0x100000)
351 #define POW22 (0x200000)
352 #define POW23 (0x400000)
353 #define POW24 (0x800000)
354 #define POW25 (0x1000000)
355 #define POW26 (0x2000000)
356 #define POW27 (0x4000000)
357 #define POW28 (0x8000000)
358 #define POW29 (0x10000000)
359 #define POW30 (0x20000000)
360 #define POW31 (0x40000000)
361 /**************
362  * useful macros
363  ************** */
364 #define picofftsg_highestBitPos(x) (x==0?0:(x>=POW17?(x>=POW25?(x>=POW29?(x>=POW31?31:(x>=POW30?30:29)):(x>=POW27?(x>=POW28?28:27):(x>=POW26?26:25))):(x>=POW21?(x>=POW23?(x>=POW24?24:23):(x>=POW22?22:21)):(x>=POW19?(x>=POW20?20:19):(x>=POW18?18:17)))):(x>=POW9?(x>=POW13?(x>=POW15?(x>=POW16?16:15):(x>=POW14?14:13)):(x>=POW11?(x>=POW12?12:11):(x>=POW10?10:9))):(x>=POW5?(x>=POW7?(x>=POW8?8:7):(x>=POW6?6:5)):(x>=POW3?(x>=POW4?4:3):(x>=POW2?2:1))))))
365 #define picofftsg_highestBit(x) (x==0?0:(x<0?(zz=-x,(zz>=POW17?(zz>=POW25?(zz>=POW29?(zz>=POW31?31:(zz>=POW30?30:29)):(zz>=POW27?(zz>=POW28?28:27):(zz>=POW26?26:25))):(zz>=POW21?(zz>=POW23?(zz>=POW24?24:23):(zz>=POW22?22:21)):(zz>=POW19?(zz>=POW20?20:19):(zz>=POW18?18:17)))):(zz>=POW9?(zz>=POW13?(zz>=POW15?(zz>=POW16?16:15):(zz>=POW14?14:13)):(zz>=POW11?(zz>=POW12?12:11):(zz>=POW10?10:9))):(zz>=POW5?(zz>=POW7?(zz>=POW8?8:7):(zz>=POW6?6:5)):(zz>=POW3?(zz>=POW4?4:3):(zz>=POW2?2:1)))))):(x>=POW17?(x>=POW25?(x>=POW29?(x>=POW31?31:(x>=POW30?30:29)):(x>=POW27?(x>=POW28?28:27):(x>=POW26?26:25))):(x>=POW21?(x>=POW23?(x>=POW24?24:23):(x>=POW22?22:21)):(x>=POW19?(x>=POW20?20:19):(x>=POW18?18:17)))):(x>=POW9?(x>=POW13?(x>=POW15?(x>=POW16?16:15):(x>=POW14?14:13)):(x>=POW11?(x>=POW12?12:11):(x>=POW10?10:9))):(x>=POW5?(x>=POW7?(x>=POW8?8:7):(x>=POW6?6:5)):(x>=POW3?(x>=POW4?4:3):(x>=POW2?2:1)))))))
366 #define Mult_W_W picofftsg_mult_w_w
367 
368 
369 /* ***********************************************************************************************/
370 /* forward declarations */
371 /* ***********************************************************************************************/
372 static PICOFFTSG_FFTTYPE picofftsg_mult_w_w(PICOFFTSG_FFTTYPE x1, PICOFFTSG_FFTTYPE y1);
373 static PICOFFTSG_FFTTYPE picofftsg_mult_w_a(PICOFFTSG_FFTTYPE x1, PICOFFTSG_FFTTYPE y1);
374 
375 
376 static void cftfsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
377 static void cftbsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
378 static void rftfsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
379 static void rftbsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
380 
381 static void cftfsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
382 static void cftbsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
383 static void rftfsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
384 static void rftbsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
385 static void dctsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
386 static void dctsub4(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
387 
388 static void ddct(picoos_int32 n, picoos_int32 isgn, PICOFFTSG_FFTTYPE *a);
389 static void bitrv1(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
390 
391 static void bitrv2(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
392 static void bitrv216(PICOFFTSG_FFTTYPE *a);
393 static void bitrv208(PICOFFTSG_FFTTYPE *a);
394 static void cftmdl1(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
395 static void cftrec4(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
396 static void cftleaf(picoos_int32 n, picoos_int32 isplt, PICOFFTSG_FFTTYPE *a);
397 static void cftfx41(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
398 static void cftf161(PICOFFTSG_FFTTYPE *a);
399 static void cftf081(PICOFFTSG_FFTTYPE *a);
400 static void cftf040(PICOFFTSG_FFTTYPE *a);
401 static void cftx020(PICOFFTSG_FFTTYPE *a);
402 
403 void bitrv2conj(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
404 void bitrv216neg(PICOFFTSG_FFTTYPE *a);
405 void bitrv208neg(PICOFFTSG_FFTTYPE *a);
406 void cftb1st(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
407 void cftrec4(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
408 void cftleaf(picoos_int32 n, picoos_int32 isplt, PICOFFTSG_FFTTYPE *a);
409 void cftfx41(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
410 void cftf161(PICOFFTSG_FFTTYPE *a);
411 void cftf081(PICOFFTSG_FFTTYPE *a);
412 void cftb040(PICOFFTSG_FFTTYPE *a);
413 void cftx020(PICOFFTSG_FFTTYPE *a);
414 
415 static picoos_int32 cfttree(picoos_int32 n, picoos_int32 j, picoos_int32 k, PICOFFTSG_FFTTYPE *a);
416 static void cftleaf(picoos_int32 n, picoos_int32 isplt, PICOFFTSG_FFTTYPE *a);
417 static void cftmdl1(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
418 
419 static void cftmdl1(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
420 static void cftmdl2(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
421 
422 static void cftmdl1(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
423 static void cftmdl2(picoos_int32 n, PICOFFTSG_FFTTYPE *a);
424 static void cftf161(PICOFFTSG_FFTTYPE *a);
425 static void cftf162(PICOFFTSG_FFTTYPE *a);
426 static void cftf081(PICOFFTSG_FFTTYPE *a);
427 static void cftf082(PICOFFTSG_FFTTYPE *a);
428 
429 static void cftf161(PICOFFTSG_FFTTYPE *a);
430 static void cftf162(PICOFFTSG_FFTTYPE *a);
431 static void cftf081(PICOFFTSG_FFTTYPE *a);
432 static void cftf082(PICOFFTSG_FFTTYPE *a);
433 
434 /* ***********************************************************************************************/
435 /* Exported functions */
436 /* ***********************************************************************************************/
rdft(picoos_int32 n,picoos_int32 isgn,PICOFFTSG_FFTTYPE * a)437 void rdft(picoos_int32 n, picoos_int32 isgn, PICOFFTSG_FFTTYPE *a)
438 {
439     PICOFFTSG_FFTTYPE xi;
440 
441     if (isgn >= 0) {
442         if (n > 4) {
443             cftfsub(n, a);
444             rftfsub(n, a);
445         } else if (n == 4) {
446             cftfsub(n, a);
447         }
448         xi = a[0] - a[1];
449         a[0] += a[1];
450         a[1] = xi;
451     } else {
452         a[1] =  (a[0] - a[1]) / 2;
453         a[0] -= a[1];
454         if (n > 4) {
455             rftbsub(n, a);
456             cftbsub(n, a);
457         } else if (n == 4) {
458             cftbsub(n, a);
459         }
460     }
461 
462 }
463 
464 
norm_result(picoos_int32 m2,PICOFFTSG_FFTTYPE * tmpX,PICOFFTSG_FFTTYPE * norm_window)465 picoos_single norm_result(picoos_int32 m2, PICOFFTSG_FFTTYPE *tmpX, PICOFFTSG_FFTTYPE *norm_window)
466 {
467     picoos_int16 nI;
468     PICOFFTSG_FFTTYPE a,b, E;
469 
470     E = (picoos_int32)0;
471     for (nI=0; nI<m2; nI++) {
472         a = (norm_window[nI]>>18) * ((tmpX[nI]>0) ? tmpX[nI]>>11 : -((-tmpX[nI])>>11));
473         tmpX[nI] = a;
474         b = (a>=0?a:-a)  >> 18;
475         E += (b*b);
476     }
477 
478     if (E>0) {
479         return (picoos_single)sqrt((double)E/16.0)/m2;
480     }
481     else {
482         return 0.0;
483     }
484 }
485 
ddct(picoos_int32 n,picoos_int32 isgn,PICOFFTSG_FFTTYPE * a)486 void ddct(picoos_int32 n, picoos_int32 isgn, PICOFFTSG_FFTTYPE *a)
487 {
488     picoos_int32 j;
489     PICOFFTSG_FFTTYPE xr;
490 
491     if (isgn < 0) {
492         xr = a[n - 1];
493         for (j = n - 2; j >= 2; j -= 2) {
494             a[j + 1] = a[j] - a[j - 1];
495             a[j] += a[j - 1];
496         }
497         a[1] = a[0] - xr;
498         a[0] += xr;
499         if (n > 4) {
500             rftbsub(n, a);
501             cftbsub(n, a);
502         } else if (n == 4) {
503             cftbsub(n, a);
504         }
505     }
506     if (n > 4) {
507         dctsub(n, a);
508     } else {
509         dctsub4(n, a);
510     }
511     if (isgn >= 0) {
512         if (n > 4) {
513             cftfsub(n, a);
514             rftfsub(n, a);
515         } else if (n == 4) {
516             cftfsub(n, a);
517         }
518 
519 
520         xr = a[0] - a[1];
521         a[0] += a[1];
522         for (j = 2; j < n; j += 2) {
523             a[j - 1] = a[j] - a[j + 1];
524             a[j] += a[j + 1];
525         }
526         a[n - 1] = xr;
527     }
528 }
529 
dfct_nmf(picoos_int32 n,picoos_int32 * a)530 void dfct_nmf(picoos_int32 n, picoos_int32 *a)
531 {
532     picoos_int32 j, k, m, mh;
533     PICOFFTSG_FFTTYPE xr, xi, yr, yi, an;
534     PICOFFTSG_FFTTYPE *aj, *ak, *amj, *amk;
535 
536     m = n >> 1;
537     for (j = 0; j < m; j++) {
538         k = n - j;
539         xr = a[j] + a[k];
540         a[j] -= a[k];
541         a[k] = xr;
542     }
543     an = a[n];
544     while (m >= 2) {
545         ddct(m, 1, a);
546         if (m > 2) {
547             bitrv1(m, a);
548         }
549         mh = m >> 1;
550         xi = a[m];
551         a[m] = a[0];
552         a[0] = an - xi;
553         an += xi;
554         k = m-1;
555         aj = a + 1; ak = a + k; amj = aj + m; amk = ak + m;
556         for (j = 1; j < mh; j++, aj++, ak--, amj++, amk--) {
557             xr = *amk;
558             xi = *amj;
559             yr = *aj;
560             yi = *ak;
561             *amj = yr;
562             *amk = yi;
563             *aj = xr - xi;
564             *ak = xr + xi;
565         }
566         xr = *aj;
567         *aj = *amj;
568         *amj = xr;
569 
570         m = mh;
571     }
572 
573     xi = a[1];
574     a[1] = a[0];
575     a[0] = an + xi;
576     a[n] = an - xi;
577     if (n > 2) {
578         bitrv1(n, a);
579     }
580 
581 }
582 
583 /* ***********************************************************************************************/
584 /* internal routines */
585 /* ***********************************************************************************************/
586 /*
587   mult two numbers which are guaranteed to be in the range -1..1
588   shift right as little as possible before mult, and the rest after the mult
589   Also, shift bigger number more - lose less accuracy
590  */
picofftsg_mult_w_w(PICOFFTSG_FFTTYPE x1,PICOFFTSG_FFTTYPE y1)591 static PICOFFTSG_FFTTYPE picofftsg_mult_w_w(PICOFFTSG_FFTTYPE x1, PICOFFTSG_FFTTYPE y1)
592 {
593     PICOFFTSG_FFTTYPE x, y;
594     x = x1>=0 ? x1>>15 : -((-x1)>>15);
595     y = y1>=0 ? y1>>14 : -((-y1)>>14);
596     return  x * y;
597 }
598 
picofftsg_mult_w_a(PICOFFTSG_FFTTYPE x1,PICOFFTSG_FFTTYPE y1)599 static PICOFFTSG_FFTTYPE picofftsg_mult_w_a(PICOFFTSG_FFTTYPE x1, PICOFFTSG_FFTTYPE y1)
600 {
601     PICOFFTSG_FFTTYPE x, y;
602 
603 
604     x = x1>=0 ? x1>>15 : -((-x1)>>15);
605     y = y1>=0 ? y1>>14 : -((-y1)>>14);
606     return x * y;
607 }
608 
cftfsub(picoos_int32 n,PICOFFTSG_FFTTYPE * a)609 static void cftfsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
610 {
611 
612     if (n > 8) {
613         if (n > 32) {
614             cftmdl1(n, a);
615                 if (n > 512) {
616                     cftrec4(n, a);
617                 } else if (n > 128) {
618                     cftleaf(n, 1, a);
619                 } else {
620                     cftfx41(n, a);
621                 }
622             bitrv2(n, a);
623         } else if (n == 32) {
624             cftf161(a);
625             bitrv216(a);
626         } else {
627             cftf081(a);
628             bitrv208(a);
629         }
630     } else if (n == 8) {
631         cftf040(a);
632     } else if (n == 4) {
633         cftx020(a);
634     }
635 }
636 
637 
cftbsub(picoos_int32 n,PICOFFTSG_FFTTYPE * a)638 void cftbsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
639 {
640     if (n > 8) {
641         if (n > 32) {
642             cftb1st(n, a);
643                 if (n > 512) {
644                     cftrec4(n, a);
645                 } else if (n > 128) {
646                     cftleaf(n, 1, a);
647                 } else {
648                     cftfx41(n, a);
649                 }
650             bitrv2conj(n, a);
651         } else if (n == 32) {
652             cftf161(a);
653             bitrv216neg(a);
654         } else {
655             cftf081(a);
656             bitrv208neg(a);
657         }
658     } else if (n == 8) {
659         cftb040(a);
660     } else if (n == 4) {
661         cftx020(a);
662     }
663 }
664 
665 /* **************************************************************************************************/
666 
667 /* **************************************************************************************************/
bitrv2(picoos_int32 n,PICOFFTSG_FFTTYPE * a)668 void bitrv2(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
669 {
670     picoos_int32 j0, k0, j1, k1, l, m, i, j, k, nh, m2;
671     PICOFFTSG_FFTTYPE xr, xi, yr, yi;
672 
673     m = 4;
674     for (l = n >> 2; l > 8; l >>= 2) {
675         m <<= 1;
676     }
677     m2 = m + m;
678     nh = n >> 1;
679     if (l == 8) {
680         j0 = 0;
681         for (k0 = 0; k0 < m; k0 += 4) {
682             k = k0;
683             for (j = j0; j < j0 + k0; j += 4) {
684                 xr = a[j];
685                 xi = a[j + 1];
686                 yr = a[k];
687                 yi = a[k + 1];
688                 a[j] = yr;
689                 a[j + 1] = yi;
690                 a[k] = xr;
691                 a[k + 1] = xi;
692                 j1 = j + m;
693                 k1 = k + m2;
694                 xr = a[j1];
695                 xi = a[j1 + 1];
696                 yr = a[k1];
697                 yi = a[k1 + 1];
698                 a[j1] = yr;
699                 a[j1 + 1] = yi;
700                 a[k1] = xr;
701                 a[k1 + 1] = xi;
702                 j1 += m;
703                 k1 -= m;
704                 xr = a[j1];
705                 xi = a[j1 + 1];
706                 yr = a[k1];
707                 yi = a[k1 + 1];
708                 a[j1] = yr;
709                 a[j1 + 1] = yi;
710                 a[k1] = xr;
711                 a[k1 + 1] = xi;
712                 j1 += m;
713                 k1 += m2;
714                 xr = a[j1];
715                 xi = a[j1 + 1];
716                 yr = a[k1];
717                 yi = a[k1 + 1];
718                 a[j1] = yr;
719                 a[j1 + 1] = yi;
720                 a[k1] = xr;
721                 a[k1 + 1] = xi;
722                 j1 += nh;
723                 k1 += 2;
724                 xr = a[j1];
725                 xi = a[j1 + 1];
726                 yr = a[k1];
727                 yi = a[k1 + 1];
728                 a[j1] = yr;
729                 a[j1 + 1] = yi;
730                 a[k1] = xr;
731                 a[k1 + 1] = xi;
732                 j1 -= m;
733                 k1 -= m2;
734                 xr = a[j1];
735                 xi = a[j1 + 1];
736                 yr = a[k1];
737                 yi = a[k1 + 1];
738                 a[j1] = yr;
739                 a[j1 + 1] = yi;
740                 a[k1] = xr;
741                 a[k1 + 1] = xi;
742                 j1 -= m;
743                 k1 += m;
744                 xr = a[j1];
745                 xi = a[j1 + 1];
746                 yr = a[k1];
747                 yi = a[k1 + 1];
748                 a[j1] = yr;
749                 a[j1 + 1] = yi;
750                 a[k1] = xr;
751                 a[k1 + 1] = xi;
752                 j1 -= m;
753                 k1 -= m2;
754                 xr = a[j1];
755                 xi = a[j1 + 1];
756                 yr = a[k1];
757                 yi = a[k1 + 1];
758                 a[j1] = yr;
759                 a[j1 + 1] = yi;
760                 a[k1] = xr;
761                 a[k1 + 1] = xi;
762                 j1 += 2;
763                 k1 += nh;
764                 xr = a[j1];
765                 xi = a[j1 + 1];
766                 yr = a[k1];
767                 yi = a[k1 + 1];
768                 a[j1] = yr;
769                 a[j1 + 1] = yi;
770                 a[k1] = xr;
771                 a[k1 + 1] = xi;
772                 j1 += m;
773                 k1 += m2;
774                 xr = a[j1];
775                 xi = a[j1 + 1];
776                 yr = a[k1];
777                 yi = a[k1 + 1];
778                 a[j1] = yr;
779                 a[j1 + 1] = yi;
780                 a[k1] = xr;
781                 a[k1 + 1] = xi;
782                 j1 += m;
783                 k1 -= m;
784                 xr = a[j1];
785                 xi = a[j1 + 1];
786                 yr = a[k1];
787                 yi = a[k1 + 1];
788                 a[j1] = yr;
789                 a[j1 + 1] = yi;
790                 a[k1] = xr;
791                 a[k1 + 1] = xi;
792                 j1 += m;
793                 k1 += m2;
794                 xr = a[j1];
795                 xi = a[j1 + 1];
796                 yr = a[k1];
797                 yi = a[k1 + 1];
798                 a[j1] = yr;
799                 a[j1 + 1] = yi;
800                 a[k1] = xr;
801                 a[k1 + 1] = xi;
802                 j1 -= nh;
803                 k1 -= 2;
804                 xr = a[j1];
805                 xi = a[j1 + 1];
806                 yr = a[k1];
807                 yi = a[k1 + 1];
808                 a[j1] = yr;
809                 a[j1 + 1] = yi;
810                 a[k1] = xr;
811                 a[k1 + 1] = xi;
812                 j1 -= m;
813                 k1 -= m2;
814                 xr = a[j1];
815                 xi = a[j1 + 1];
816                 yr = a[k1];
817                 yi = a[k1 + 1];
818                 a[j1] = yr;
819                 a[j1 + 1] = yi;
820                 a[k1] = xr;
821                 a[k1 + 1] = xi;
822                 j1 -= m;
823                 k1 += m;
824                 xr = a[j1];
825                 xi = a[j1 + 1];
826                 yr = a[k1];
827                 yi = a[k1 + 1];
828                 a[j1] = yr;
829                 a[j1 + 1] = yi;
830                 a[k1] = xr;
831                 a[k1 + 1] = xi;
832                 j1 -= m;
833                 k1 -= m2;
834                 xr = a[j1];
835                 xi = a[j1 + 1];
836                 yr = a[k1];
837                 yi = a[k1 + 1];
838                 a[j1] = yr;
839                 a[j1 + 1] = yi;
840                 a[k1] = xr;
841                 a[k1 + 1] = xi;
842                 for (i = nh >> 1; i > (k ^= i); i >>= 1) {
843                     /* Avoid warning*/
844                 };
845             }
846             k1 = j0 + k0;
847             j1 = k1 + 2;
848             k1 += nh;
849             xr = a[j1];
850             xi = a[j1 + 1];
851             yr = a[k1];
852             yi = a[k1 + 1];
853             a[j1] = yr;
854             a[j1 + 1] = yi;
855             a[k1] = xr;
856             a[k1 + 1] = xi;
857             j1 += m;
858             k1 += m2;
859             xr = a[j1];
860             xi = a[j1 + 1];
861             yr = a[k1];
862             yi = a[k1 + 1];
863             a[j1] = yr;
864             a[j1 + 1] = yi;
865             a[k1] = xr;
866             a[k1 + 1] = xi;
867             j1 += m;
868             k1 -= m;
869             xr = a[j1];
870             xi = a[j1 + 1];
871             yr = a[k1];
872             yi = a[k1 + 1];
873             a[j1] = yr;
874             a[j1 + 1] = yi;
875             a[k1] = xr;
876             a[k1 + 1] = xi;
877             j1 -= 2;
878             k1 -= nh;
879             xr = a[j1];
880             xi = a[j1 + 1];
881             yr = a[k1];
882             yi = a[k1 + 1];
883             a[j1] = yr;
884             a[j1 + 1] = yi;
885             a[k1] = xr;
886             a[k1 + 1] = xi;
887             j1 += nh + 2;
888             k1 += nh + 2;
889             xr = a[j1];
890             xi = a[j1 + 1];
891             yr = a[k1];
892             yi = a[k1 + 1];
893             a[j1] = yr;
894             a[j1 + 1] = yi;
895             a[k1] = xr;
896             a[k1 + 1] = xi;
897             j1 -= nh - m;
898             k1 += m2 - 2;
899             xr = a[j1];
900             xi = a[j1 + 1];
901             yr = a[k1];
902             yi = a[k1 + 1];
903             a[j1] = yr;
904             a[j1 + 1] = yi;
905             a[k1] = xr;
906             a[k1 + 1] = xi;
907             for (i = nh >> 1; i > (j0 ^= i); i >>= 1) {
908                 /* Avoid warning */
909             }
910         }
911     } else {
912         j0 = 0;
913         for (k0 = 0; k0 < m; k0 += 4) {
914             k = k0;
915             for (j = j0; j < j0 + k0; j += 4) {
916                 xr = a[j];
917                 xi = a[j + 1];
918                 yr = a[k];
919                 yi = a[k + 1];
920                 a[j] = yr;
921                 a[j + 1] = yi;
922                 a[k] = xr;
923                 a[k + 1] = xi;
924                 j1 = j + m;
925                 k1 = k + m;
926                 xr = a[j1];
927                 xi = a[j1 + 1];
928                 yr = a[k1];
929                 yi = a[k1 + 1];
930                 a[j1] = yr;
931                 a[j1 + 1] = yi;
932                 a[k1] = xr;
933                 a[k1 + 1] = xi;
934                 j1 += nh;
935                 k1 += 2;
936                 xr = a[j1];
937                 xi = a[j1 + 1];
938                 yr = a[k1];
939                 yi = a[k1 + 1];
940                 a[j1] = yr;
941                 a[j1 + 1] = yi;
942                 a[k1] = xr;
943                 a[k1 + 1] = xi;
944                 j1 -= m;
945                 k1 -= m;
946                 xr = a[j1];
947                 xi = a[j1 + 1];
948                 yr = a[k1];
949                 yi = a[k1 + 1];
950                 a[j1] = yr;
951                 a[j1 + 1] = yi;
952                 a[k1] = xr;
953                 a[k1 + 1] = xi;
954                 j1 += 2;
955                 k1 += nh;
956                 xr = a[j1];
957                 xi = a[j1 + 1];
958                 yr = a[k1];
959                 yi = a[k1 + 1];
960                 a[j1] = yr;
961                 a[j1 + 1] = yi;
962                 a[k1] = xr;
963                 a[k1 + 1] = xi;
964                 j1 += m;
965                 k1 += m;
966                 xr = a[j1];
967                 xi = a[j1 + 1];
968                 yr = a[k1];
969                 yi = a[k1 + 1];
970                 a[j1] = yr;
971                 a[j1 + 1] = yi;
972                 a[k1] = xr;
973                 a[k1 + 1] = xi;
974                 j1 -= nh;
975                 k1 -= 2;
976                 xr = a[j1];
977                 xi = a[j1 + 1];
978                 yr = a[k1];
979                 yi = a[k1 + 1];
980                 a[j1] = yr;
981                 a[j1 + 1] = yi;
982                 a[k1] = xr;
983                 a[k1 + 1] = xi;
984                 j1 -= m;
985                 k1 -= m;
986                 xr = a[j1];
987                 xi = a[j1 + 1];
988                 yr = a[k1];
989                 yi = a[k1 + 1];
990                 a[j1] = yr;
991                 a[j1 + 1] = yi;
992                 a[k1] = xr;
993                 a[k1 + 1] = xi;
994                 for (i = nh >> 1; i > (k ^= i); i >>= 1){
995                     /* Avoid warning */
996                 }
997             }
998             k1 = j0 + k0;
999             j1 = k1 + 2;
1000             k1 += nh;
1001             xr = a[j1];
1002             xi = a[j1 + 1];
1003             yr = a[k1];
1004             yi = a[k1 + 1];
1005             a[j1] = yr;
1006             a[j1 + 1] = yi;
1007             a[k1] = xr;
1008             a[k1 + 1] = xi;
1009             j1 += m;
1010             k1 += m;
1011             xr = a[j1];
1012             xi = a[j1 + 1];
1013             yr = a[k1];
1014             yi = a[k1 + 1];
1015             a[j1] = yr;
1016             a[j1 + 1] = yi;
1017             a[k1] = xr;
1018             a[k1 + 1] = xi;
1019             for (i = nh >> 1; i > (j0 ^= i); i >>= 1){
1020                 /* Avoid  warning */
1021             }
1022         }
1023     }
1024 }
1025 
1026 
bitrv2conj(picoos_int32 n,PICOFFTSG_FFTTYPE * a)1027 void bitrv2conj(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
1028 {
1029     picoos_int32 j0, k0, j1, k1, l, m, i, j, k, nh, m2;
1030     PICOFFTSG_FFTTYPE xr, xi, yr, yi;
1031 
1032 
1033     m = 4;
1034     for (l = n >> 2; l > 8; l >>= 2) {
1035         m <<= 1;
1036     }
1037     m2 = m + m;
1038     nh = n >> 1;
1039     if (l == 8) {
1040         j0 = 0;
1041         for (k0 = 0; k0 < m; k0 += 4) {
1042             k = k0;
1043             for (j = j0; j < j0 + k0; j += 4) {
1044                 xr = a[j];
1045                 xi = -a[j + 1];
1046                 yr = a[k];
1047                 yi = -a[k + 1];
1048                 a[j] = yr;
1049                 a[j + 1] = yi;
1050                 a[k] = xr;
1051                 a[k + 1] = xi;
1052                 j1 = j + m;
1053                 k1 = k + m2;
1054                 xr = a[j1];
1055                 xi = -a[j1 + 1];
1056                 yr = a[k1];
1057                 yi = -a[k1 + 1];
1058                 a[j1] = yr;
1059                 a[j1 + 1] = yi;
1060                 a[k1] = xr;
1061                 a[k1 + 1] = xi;
1062                 j1 += m;
1063                 k1 -= m;
1064                 xr = a[j1];
1065                 xi = -a[j1 + 1];
1066                 yr = a[k1];
1067                 yi = -a[k1 + 1];
1068                 a[j1] = yr;
1069                 a[j1 + 1] = yi;
1070                 a[k1] = xr;
1071                 a[k1 + 1] = xi;
1072                 j1 += m;
1073                 k1 += m2;
1074                 xr = a[j1];
1075                 xi = -a[j1 + 1];
1076                 yr = a[k1];
1077                 yi = -a[k1 + 1];
1078                 a[j1] = yr;
1079                 a[j1 + 1] = yi;
1080                 a[k1] = xr;
1081                 a[k1 + 1] = xi;
1082                 j1 += nh;
1083                 k1 += 2;
1084                 xr = a[j1];
1085                 xi = -a[j1 + 1];
1086                 yr = a[k1];
1087                 yi = -a[k1 + 1];
1088                 a[j1] = yr;
1089                 a[j1 + 1] = yi;
1090                 a[k1] = xr;
1091                 a[k1 + 1] = xi;
1092                 j1 -= m;
1093                 k1 -= m2;
1094                 xr = a[j1];
1095                 xi = -a[j1 + 1];
1096                 yr = a[k1];
1097                 yi = -a[k1 + 1];
1098                 a[j1] = yr;
1099                 a[j1 + 1] = yi;
1100                 a[k1] = xr;
1101                 a[k1 + 1] = xi;
1102                 j1 -= m;
1103                 k1 += m;
1104                 xr = a[j1];
1105                 xi = -a[j1 + 1];
1106                 yr = a[k1];
1107                 yi = -a[k1 + 1];
1108                 a[j1] = yr;
1109                 a[j1 + 1] = yi;
1110                 a[k1] = xr;
1111                 a[k1 + 1] = xi;
1112                 j1 -= m;
1113                 k1 -= m2;
1114                 xr = a[j1];
1115                 xi = -a[j1 + 1];
1116                 yr = a[k1];
1117                 yi = -a[k1 + 1];
1118                 a[j1] = yr;
1119                 a[j1 + 1] = yi;
1120                 a[k1] = xr;
1121                 a[k1 + 1] = xi;
1122                 j1 += 2;
1123                 k1 += nh;
1124                 xr = a[j1];
1125                 xi = -a[j1 + 1];
1126                 yr = a[k1];
1127                 yi = -a[k1 + 1];
1128                 a[j1] = yr;
1129                 a[j1 + 1] = yi;
1130                 a[k1] = xr;
1131                 a[k1 + 1] = xi;
1132                 j1 += m;
1133                 k1 += m2;
1134                 xr = a[j1];
1135                 xi = -a[j1 + 1];
1136                 yr = a[k1];
1137                 yi = -a[k1 + 1];
1138                 a[j1] = yr;
1139                 a[j1 + 1] = yi;
1140                 a[k1] = xr;
1141                 a[k1 + 1] = xi;
1142                 j1 += m;
1143                 k1 -= m;
1144                 xr = a[j1];
1145                 xi = -a[j1 + 1];
1146                 yr = a[k1];
1147                 yi = -a[k1 + 1];
1148                 a[j1] = yr;
1149                 a[j1 + 1] = yi;
1150                 a[k1] = xr;
1151                 a[k1 + 1] = xi;
1152                 j1 += m;
1153                 k1 += m2;
1154                 xr = a[j1];
1155                 xi = -a[j1 + 1];
1156                 yr = a[k1];
1157                 yi = -a[k1 + 1];
1158                 a[j1] = yr;
1159                 a[j1 + 1] = yi;
1160                 a[k1] = xr;
1161                 a[k1 + 1] = xi;
1162                 j1 -= nh;
1163                 k1 -= 2;
1164                 xr = a[j1];
1165                 xi = -a[j1 + 1];
1166                 yr = a[k1];
1167                 yi = -a[k1 + 1];
1168                 a[j1] = yr;
1169                 a[j1 + 1] = yi;
1170                 a[k1] = xr;
1171                 a[k1 + 1] = xi;
1172                 j1 -= m;
1173                 k1 -= m2;
1174                 xr = a[j1];
1175                 xi = -a[j1 + 1];
1176                 yr = a[k1];
1177                 yi = -a[k1 + 1];
1178                 a[j1] = yr;
1179                 a[j1 + 1] = yi;
1180                 a[k1] = xr;
1181                 a[k1 + 1] = xi;
1182                 j1 -= m;
1183                 k1 += m;
1184                 xr = a[j1];
1185                 xi = -a[j1 + 1];
1186                 yr = a[k1];
1187                 yi = -a[k1 + 1];
1188                 a[j1] = yr;
1189                 a[j1 + 1] = yi;
1190                 a[k1] = xr;
1191                 a[k1 + 1] = xi;
1192                 j1 -= m;
1193                 k1 -= m2;
1194                 xr = a[j1];
1195                 xi = -a[j1 + 1];
1196                 yr = a[k1];
1197                 yi = -a[k1 + 1];
1198                 a[j1] = yr;
1199                 a[j1 + 1] = yi;
1200                 a[k1] = xr;
1201                 a[k1 + 1] = xi;
1202                 for (i = nh >> 1; i > (k ^= i); i >>= 1) {
1203                     /* Avoid warning */
1204                 }
1205             }
1206             k1 = j0 + k0;
1207             j1 = k1 + 2;
1208             k1 += nh;
1209             a[j1 - 1] = -a[j1 - 1];
1210             xr = a[j1];
1211             xi = -a[j1 + 1];
1212             yr = a[k1];
1213             yi = -a[k1 + 1];
1214             a[j1] = yr;
1215             a[j1 + 1] = yi;
1216             a[k1] = xr;
1217             a[k1 + 1] = xi;
1218             a[k1 + 3] = -a[k1 + 3];
1219             j1 += m;
1220             k1 += m2;
1221             xr = a[j1];
1222             xi = -a[j1 + 1];
1223             yr = a[k1];
1224             yi = -a[k1 + 1];
1225             a[j1] = yr;
1226             a[j1 + 1] = yi;
1227             a[k1] = xr;
1228             a[k1 + 1] = xi;
1229             j1 += m;
1230             k1 -= m;
1231             xr = a[j1];
1232             xi = -a[j1 + 1];
1233             yr = a[k1];
1234             yi = -a[k1 + 1];
1235             a[j1] = yr;
1236             a[j1 + 1] = yi;
1237             a[k1] = xr;
1238             a[k1 + 1] = xi;
1239             j1 -= 2;
1240             k1 -= nh;
1241             xr = a[j1];
1242             xi = -a[j1 + 1];
1243             yr = a[k1];
1244             yi = -a[k1 + 1];
1245             a[j1] = yr;
1246             a[j1 + 1] = yi;
1247             a[k1] = xr;
1248             a[k1 + 1] = xi;
1249             j1 += nh + 2;
1250             k1 += nh + 2;
1251             xr = a[j1];
1252             xi = -a[j1 + 1];
1253             yr = a[k1];
1254             yi = -a[k1 + 1];
1255             a[j1] = yr;
1256             a[j1 + 1] = yi;
1257             a[k1] = xr;
1258             a[k1 + 1] = xi;
1259             j1 -= nh - m;
1260             k1 += m2 - 2;
1261             a[j1 - 1] = -a[j1 - 1];
1262             xr = a[j1];
1263             xi = -a[j1 + 1];
1264             yr = a[k1];
1265             yi = -a[k1 + 1];
1266             a[j1] = yr;
1267             a[j1 + 1] = yi;
1268             a[k1] = xr;
1269             a[k1 + 1] = xi;
1270             a[k1 + 3] = -a[k1 + 3];
1271             for (i = nh >> 1; i > (j0 ^= i); i >>= 1) { /* Avoid warning*/
1272 
1273             }
1274         }
1275     } else {
1276         j0 = 0;
1277         for (k0 = 0; k0 < m; k0 += 4) {
1278             k = k0;
1279             for (j = j0; j < j0 + k0; j += 4) {
1280                 xr = a[j];
1281                 xi = -a[j + 1];
1282                 yr = a[k];
1283                 yi = -a[k + 1];
1284                 a[j] = yr;
1285                 a[j + 1] = yi;
1286                 a[k] = xr;
1287                 a[k + 1] = xi;
1288                 j1 = j + m;
1289                 k1 = k + m;
1290                 xr = a[j1];
1291                 xi = -a[j1 + 1];
1292                 yr = a[k1];
1293                 yi = -a[k1 + 1];
1294                 a[j1] = yr;
1295                 a[j1 + 1] = yi;
1296                 a[k1] = xr;
1297                 a[k1 + 1] = xi;
1298                 j1 += nh;
1299                 k1 += 2;
1300                 xr = a[j1];
1301                 xi = -a[j1 + 1];
1302                 yr = a[k1];
1303                 yi = -a[k1 + 1];
1304                 a[j1] = yr;
1305                 a[j1 + 1] = yi;
1306                 a[k1] = xr;
1307                 a[k1 + 1] = xi;
1308                 j1 -= m;
1309                 k1 -= m;
1310                 xr = a[j1];
1311                 xi = -a[j1 + 1];
1312                 yr = a[k1];
1313                 yi = -a[k1 + 1];
1314                 a[j1] = yr;
1315                 a[j1 + 1] = yi;
1316                 a[k1] = xr;
1317                 a[k1 + 1] = xi;
1318                 j1 += 2;
1319                 k1 += nh;
1320                 xr = a[j1];
1321                 xi = -a[j1 + 1];
1322                 yr = a[k1];
1323                 yi = -a[k1 + 1];
1324                 a[j1] = yr;
1325                 a[j1 + 1] = yi;
1326                 a[k1] = xr;
1327                 a[k1 + 1] = xi;
1328                 j1 += m;
1329                 k1 += m;
1330                 xr = a[j1];
1331                 xi = -a[j1 + 1];
1332                 yr = a[k1];
1333                 yi = -a[k1 + 1];
1334                 a[j1] = yr;
1335                 a[j1 + 1] = yi;
1336                 a[k1] = xr;
1337                 a[k1 + 1] = xi;
1338                 j1 -= nh;
1339                 k1 -= 2;
1340                 xr = a[j1];
1341                 xi = -a[j1 + 1];
1342                 yr = a[k1];
1343                 yi = -a[k1 + 1];
1344                 a[j1] = yr;
1345                 a[j1 + 1] = yi;
1346                 a[k1] = xr;
1347                 a[k1 + 1] = xi;
1348                 j1 -= m;
1349                 k1 -= m;
1350                 xr = a[j1];
1351                 xi = -a[j1 + 1];
1352                 yr = a[k1];
1353                 yi = -a[k1 + 1];
1354                 a[j1] = yr;
1355                 a[j1 + 1] = yi;
1356                 a[k1] = xr;
1357                 a[k1 + 1] = xi;
1358                 for (i = nh >> 1; i > (k ^= i); i >>= 1) {
1359                     /* Avoid warning*/
1360                 }
1361             }
1362             k1 = j0 + k0;
1363             j1 = k1 + 2;
1364             k1 += nh;
1365             a[j1 - 1] = -a[j1 - 1];
1366             xr = a[j1];
1367             xi = -a[j1 + 1];
1368             yr = a[k1];
1369             yi = -a[k1 + 1];
1370             a[j1] = yr;
1371             a[j1 + 1] = yi;
1372             a[k1] = xr;
1373             a[k1 + 1] = xi;
1374             a[k1 + 3] = -a[k1 + 3];
1375             j1 += m;
1376             k1 += m;
1377             a[j1 - 1] = -a[j1 - 1];
1378             xr = a[j1];
1379             xi = -a[j1 + 1];
1380             yr = a[k1];
1381             yi = -a[k1 + 1];
1382             a[j1] = yr;
1383             a[j1 + 1] = yi;
1384             a[k1] = xr;
1385             a[k1 + 1] = xi;
1386             a[k1 + 3] = -a[k1 + 3];
1387             for (i = nh >> 1; i > (j0 ^= i); i >>= 1) {
1388                 /* Avoid warning*/
1389             }
1390         }
1391     }
1392 }
1393 
1394 
bitrv216(PICOFFTSG_FFTTYPE * a)1395 void bitrv216(PICOFFTSG_FFTTYPE *a)
1396 {
1397     PICOFFTSG_FFTTYPE x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1398         x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i,
1399         x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
1400 
1401     x1r = a[2];
1402     x1i = a[3];
1403     x2r = a[4];
1404     x2i = a[5];
1405     x3r = a[6];
1406     x3i = a[7];
1407     x4r = a[8];
1408     x4i = a[9];
1409     x5r = a[10];
1410     x5i = a[11];
1411     x7r = a[14];
1412     x7i = a[15];
1413     x8r = a[16];
1414     x8i = a[17];
1415     x10r = a[20];
1416     x10i = a[21];
1417     x11r = a[22];
1418     x11i = a[23];
1419     x12r = a[24];
1420     x12i = a[25];
1421     x13r = a[26];
1422     x13i = a[27];
1423     x14r = a[28];
1424     x14i = a[29];
1425     a[2] = x8r;
1426     a[3] = x8i;
1427     a[4] = x4r;
1428     a[5] = x4i;
1429     a[6] = x12r;
1430     a[7] = x12i;
1431     a[8] = x2r;
1432     a[9] = x2i;
1433     a[10] = x10r;
1434     a[11] = x10i;
1435     a[14] = x14r;
1436     a[15] = x14i;
1437     a[16] = x1r;
1438     a[17] = x1i;
1439     a[20] = x5r;
1440     a[21] = x5i;
1441     a[22] = x13r;
1442     a[23] = x13i;
1443     a[24] = x3r;
1444     a[25] = x3i;
1445     a[26] = x11r;
1446     a[27] = x11i;
1447     a[28] = x7r;
1448     a[29] = x7i;
1449 }
1450 
1451 
bitrv216neg(PICOFFTSG_FFTTYPE * a)1452 void bitrv216neg(PICOFFTSG_FFTTYPE *a)
1453 {
1454     PICOFFTSG_FFTTYPE x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1455         x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i,
1456         x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i,
1457         x13r, x13i, x14r, x14i, x15r, x15i;
1458 
1459     x1r = a[2];
1460     x1i = a[3];
1461     x2r = a[4];
1462     x2i = a[5];
1463     x3r = a[6];
1464     x3i = a[7];
1465     x4r = a[8];
1466     x4i = a[9];
1467     x5r = a[10];
1468     x5i = a[11];
1469     x6r = a[12];
1470     x6i = a[13];
1471     x7r = a[14];
1472     x7i = a[15];
1473     x8r = a[16];
1474     x8i = a[17];
1475     x9r = a[18];
1476     x9i = a[19];
1477     x10r = a[20];
1478     x10i = a[21];
1479     x11r = a[22];
1480     x11i = a[23];
1481     x12r = a[24];
1482     x12i = a[25];
1483     x13r = a[26];
1484     x13i = a[27];
1485     x14r = a[28];
1486     x14i = a[29];
1487     x15r = a[30];
1488     x15i = a[31];
1489     a[2] = x15r;
1490     a[3] = x15i;
1491     a[4] = x7r;
1492     a[5] = x7i;
1493     a[6] = x11r;
1494     a[7] = x11i;
1495     a[8] = x3r;
1496     a[9] = x3i;
1497     a[10] = x13r;
1498     a[11] = x13i;
1499     a[12] = x5r;
1500     a[13] = x5i;
1501     a[14] = x9r;
1502     a[15] = x9i;
1503     a[16] = x1r;
1504     a[17] = x1i;
1505     a[18] = x14r;
1506     a[19] = x14i;
1507     a[20] = x6r;
1508     a[21] = x6i;
1509     a[22] = x10r;
1510     a[23] = x10i;
1511     a[24] = x2r;
1512     a[25] = x2i;
1513     a[26] = x12r;
1514     a[27] = x12i;
1515     a[28] = x4r;
1516     a[29] = x4i;
1517     a[30] = x8r;
1518     a[31] = x8i;
1519 }
1520 
1521 
bitrv208(PICOFFTSG_FFTTYPE * a)1522 void bitrv208(PICOFFTSG_FFTTYPE *a)
1523 {
1524     PICOFFTSG_FFTTYPE x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
1525 
1526     x1r = a[2];
1527     x1i = a[3];
1528     x3r = a[6];
1529     x3i = a[7];
1530     x4r = a[8];
1531     x4i = a[9];
1532     x6r = a[12];
1533     x6i = a[13];
1534     a[2] = x4r;
1535     a[3] = x4i;
1536     a[6] = x6r;
1537     a[7] = x6i;
1538     a[8] = x1r;
1539     a[9] = x1i;
1540     a[12] = x3r;
1541     a[13] = x3i;
1542 }
1543 
1544 
bitrv208neg(PICOFFTSG_FFTTYPE * a)1545 void bitrv208neg(PICOFFTSG_FFTTYPE *a)
1546 {
1547     PICOFFTSG_FFTTYPE x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1548         x5r, x5i, x6r, x6i, x7r, x7i;
1549 
1550     x1r = a[2];
1551     x1i = a[3];
1552     x2r = a[4];
1553     x2i = a[5];
1554     x3r = a[6];
1555     x3i = a[7];
1556     x4r = a[8];
1557     x4i = a[9];
1558     x5r = a[10];
1559     x5i = a[11];
1560     x6r = a[12];
1561     x6i = a[13];
1562     x7r = a[14];
1563     x7i = a[15];
1564     a[2] = x7r;
1565     a[3] = x7i;
1566     a[4] = x3r;
1567     a[5] = x3i;
1568     a[6] = x5r;
1569     a[7] = x5i;
1570     a[8] = x1r;
1571     a[9] = x1i;
1572     a[10] = x6r;
1573     a[11] = x6i;
1574     a[12] = x2r;
1575     a[13] = x2i;
1576     a[14] = x4r;
1577     a[15] = x4i;
1578 }
1579 
1580 
bitrv1(picoos_int32 n,PICOFFTSG_FFTTYPE * a)1581 void bitrv1(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
1582 {
1583     picoos_int32 j0, k0, j1, k1, l, m, i, j, k, nh;
1584     PICOFFTSG_FFTTYPE x;
1585     nh = n >> 1;
1586     x = a[1];
1587     a[1] = a[nh];
1588     a[nh] = x;
1589     m = 2;
1590     for (l = n >> 2; l > 2; l >>= 2) {
1591         m <<= 1;
1592     }
1593     if (l == 2) {
1594         j1 = m + 1;
1595         k1 = m + nh;
1596         x = a[j1];
1597         a[j1] = a[k1];
1598         a[k1] = x;
1599         j0 = 0;
1600         for (k0 = 2; k0 < m; k0 += 2) {
1601             for (i = nh >> 1; i > (j0 ^= i); i >>= 1) {
1602                 /* Avoid warning*/
1603             }
1604             k = k0;
1605             for (j = j0; j < j0 + k0; j += 2) {
1606                 x = a[j];
1607                 a[j] = a[k];
1608                 a[k] = x;
1609                 j1 = j + m;
1610                 k1 = k + m;
1611                 x = a[j1];
1612                 a[j1] = a[k1];
1613                 a[k1] = x;
1614                 j1 += nh;
1615                 k1++;
1616                 x = a[j1];
1617                 a[j1] = a[k1];
1618                 a[k1] = x;
1619                 j1 -= m;
1620                 k1 -= m;
1621                 x = a[j1];
1622                 a[j1] = a[k1];
1623                 a[k1] = x;
1624                 j1++;
1625                 k1 += nh;
1626                 x = a[j1];
1627                 a[j1] = a[k1];
1628                 a[k1] = x;
1629                 j1 += m;
1630                 k1 += m;
1631                 x = a[j1];
1632                 a[j1] = a[k1];
1633                 a[k1] = x;
1634                 j1 -= nh;
1635                 k1--;
1636                 x = a[j1];
1637                 a[j1] = a[k1];
1638                 a[k1] = x;
1639                 j1 -= m;
1640                 k1 -= m;
1641                 x = a[j1];
1642                 a[j1] = a[k1];
1643                 a[k1] = x;
1644                 for (i = nh >> 1; i > (k ^= i); i >>= 1) {
1645                     /* Avoid warning*/
1646                 }
1647             }
1648             k1 = j0 + k0;
1649             j1 = k1 + 1;
1650             k1 += nh;
1651             x = a[j1];
1652             a[j1] = a[k1];
1653             a[k1] = x;
1654             j1 += m;
1655             k1 += m;
1656             x = a[j1];
1657             a[j1] = a[k1];
1658             a[k1] = x;
1659         }
1660     } else {
1661         j0 = 0;
1662         for (k0 = 2; k0 < m; k0 += 2) {
1663             for (i = nh >> 1; i > (j0 ^= i); i >>= 1) {
1664                 /* Avoid warning*/
1665             }
1666             k = k0;
1667             for (j = j0; j < j0 + k0; j += 2) {
1668                 x = a[j];
1669                 a[j] = a[k];
1670                 a[k] = x;
1671                 j1 = j + nh;
1672                 k1 = k + 1;
1673                 x = a[j1];
1674                 a[j1] = a[k1];
1675                 a[k1] = x;
1676                 j1++;
1677                 k1 += nh;
1678                 x = a[j1];
1679                 a[j1] = a[k1];
1680                 a[k1] = x;
1681                 j1 -= nh;
1682                 k1--;
1683                 x = a[j1];
1684                 a[j1] = a[k1];
1685                 a[k1] = x;
1686                 for (i = nh >> 1; i > (k ^= i); i >>= 1) {
1687                     /* Avoid warning*/
1688                 }
1689             }
1690             k1 = j0 + k0;
1691             j1 = k1 + 1;
1692             k1 += nh;
1693             x = a[j1];
1694             a[j1] = a[k1];
1695             a[k1] = x;
1696         }
1697     }
1698 }
1699 
1700 
1701 /* **************************************************************************************************/
1702 
1703 /* **************************************************************************************************/
1704 
cftb1st(picoos_int32 n,PICOFFTSG_FFTTYPE * a)1705 void cftb1st(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
1706 {
1707     picoos_int32 i, i0, j, j0, j1, j2, j3, m, mh;
1708     PICOFFTSG_FFTTYPE wk1r, wk1i, wk3r, wk3i,
1709         wd1r, wd1i, wd3r, wd3i, ss1, ss3;
1710     PICOFFTSG_FFTTYPE x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
1711 
1712     mh = n >> 3;
1713     m = 2 * mh;
1714     j1 = m;
1715     j2 = j1 + m;
1716     j3 = j2 + m;
1717     x0r = a[0] + a[j2];
1718     x0i = -a[1] - a[j2 + 1];
1719     x1r = a[0] - a[j2];
1720     x1i = -a[1] + a[j2 + 1];
1721     x2r = a[j1] + a[j3];
1722     x2i = a[j1 + 1] + a[j3 + 1];
1723     x3r = a[j1] - a[j3];
1724     x3i = a[j1 + 1] - a[j3 + 1];
1725     a[0] = x0r + x2r;
1726     a[1] = x0i - x2i;
1727     a[j1] = x0r - x2r;
1728     a[j1 + 1] = x0i + x2i;
1729     a[j2] = x1r + x3i;
1730     a[j2 + 1] = x1i + x3r;
1731     a[j3] = x1r - x3i;
1732     a[j3 + 1] = x1i - x3r;
1733     wd1r = PICODSP_WGT_SHIFT;
1734     wd1i = 0;
1735     wd3r = PICODSP_WGT_SHIFT;
1736     wd3i = 0;
1737 
1738     wk1r  = (PICOFFTSG_FFTTYPE) (0.998795449734  *PICODSP_WGT_SHIFT);
1739     wk1i  = (PICOFFTSG_FFTTYPE) (0.049067676067  *PICODSP_WGT_SHIFT);
1740     ss1   = (PICOFFTSG_FFTTYPE) (0.098135352135  *PICODSP_WGT_SHIFT);
1741     wk3i  = (PICOFFTSG_FFTTYPE) (-0.146730467677 *PICODSP_WGT_SHIFT);
1742     wk3r  = (PICOFFTSG_FFTTYPE) (0.989176511765  *PICODSP_WGT_SHIFT);
1743     ss3   = (PICOFFTSG_FFTTYPE) (-0.293460935354 *PICODSP_WGT_SHIFT);
1744 
1745     i = 0;
1746     for (;;) {
1747         i0 = i + CDFT_LOOP_DIV_4;
1748         if (i0 > mh - 4) {
1749             i0 = mh - 4;
1750         }
1751         for (j = i + 2; j < i0; j += 4) {
1752 
1753             wd1r -= Mult_W_W(ss1, wk1i);
1754             wd1i += Mult_W_W(ss1, wk1r);
1755             wd3r -= Mult_W_W(ss3, wk3i);
1756             wd3i += Mult_W_W(ss3, wk3r);
1757 
1758             j1 = j + m;
1759             j2 = j1 + m;
1760             j3 = j2 + m;
1761             x0r = a[j] + a[j2];
1762             x0i = -a[j + 1] - a[j2 + 1];
1763             x1r = a[j] - a[j2];
1764             x1i = -a[j + 1] + a[j2 + 1];
1765             x2r = a[j1] + a[j3];
1766             x2i = a[j1 + 1] + a[j3 + 1];
1767             x3r = a[j1] - a[j3];
1768             x3i = a[j1 + 1] - a[j3 + 1];
1769             a[j] = x0r + x2r;
1770             a[j + 1] = x0i - x2i;
1771             a[j1] = x0r - x2r;
1772             a[j1 + 1] = x0i + x2i;
1773             x0r = x1r + x3i;
1774             x0i = x1i + x3r;
1775             a[j2] = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
1776             a[j2 + 1] = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
1777             x0r = x1r - x3i;
1778             x0i = x1i - x3r;
1779             a[j3] = Mult_W_W(wk3r, x0r) + Mult_W_W(wk3i, x0i);
1780             a[j3 + 1] = Mult_W_W(wk3r, x0i) - Mult_W_W(wk3i, x0r);
1781             x0r = a[j + 2] + a[j2 + 2];
1782             x0i = -a[j + 3] - a[j2 + 3];
1783             x1r = a[j + 2] - a[j2 + 2];
1784             x1i = -a[j + 3] + a[j2 + 3];
1785             x2r = a[j1 + 2] + a[j3 + 2];
1786             x2i = a[j1 + 3] + a[j3 + 3];
1787             x3r = a[j1 + 2] - a[j3 + 2];
1788             x3i = a[j1 + 3] - a[j3 + 3];
1789             a[j + 2] = x0r + x2r;
1790             a[j + 3] = x0i - x2i;
1791             a[j1 + 2] = x0r - x2r;
1792             a[j1 + 3] = x0i + x2i;
1793             x0r = x1r + x3i;
1794             x0i = x1i + x3r;
1795             a[j2 + 2] = Mult_W_W(wd1r, x0r) - Mult_W_W(wd1i, x0i);
1796             a[j2 + 3] = Mult_W_W(wd1r, x0i) + Mult_W_W(wd1i, x0r);
1797             x0r = x1r - x3i;
1798             x0i = x1i - x3r;
1799             a[j3 + 2] = Mult_W_W(wd3r, x0r) + Mult_W_W(wd3i, x0i);
1800             a[j3 + 3] = Mult_W_W(wd3r, x0i) - Mult_W_W(wd3i, x0r);
1801             j0 = m - j;
1802             j1 = j0 + m;
1803             j2 = j1 + m;
1804             j3 = j2 + m;
1805             x0r = a[j0] + a[j2];
1806             x0i = -a[j0 + 1] - a[j2 + 1];
1807             x1r = a[j0] - a[j2];
1808             x1i = -a[j0 + 1] + a[j2 + 1];
1809             x2r = a[j1] + a[j3];
1810             x2i = a[j1 + 1] + a[j3 + 1];
1811             x3r = a[j1] - a[j3];
1812             x3i = a[j1 + 1] - a[j3 + 1];
1813             a[j0] = x0r + x2r;
1814             a[j0 + 1] = x0i - x2i;
1815             a[j1] = x0r - x2r;
1816             a[j1 + 1] = x0i + x2i;
1817             x0r = x1r + x3i;
1818             x0i = x1i + x3r;
1819             a[j2] = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
1820             a[j2 + 1] = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
1821             x0r = x1r - x3i;
1822             x0i = x1i - x3r;
1823             a[j3] = Mult_W_W(wk3i, x0r) + Mult_W_W(wk3r, x0i);
1824             a[j3 + 1] = Mult_W_W(wk3i, x0i) - Mult_W_W(wk3r, x0r);
1825             x0r = a[j0 - 2] + a[j2 - 2];
1826             x0i = -a[j0 - 1] - a[j2 - 1];
1827             x1r = a[j0 - 2] - a[j2 - 2];
1828             x1i = -a[j0 - 1] + a[j2 - 1];
1829             x2r = a[j1 - 2] + a[j3 - 2];
1830             x2i = a[j1 - 1] + a[j3 - 1];
1831             x3r = a[j1 - 2] - a[j3 - 2];
1832             x3i = a[j1 - 1] - a[j3 - 1];
1833             a[j0 - 2] = x0r + x2r;
1834             a[j0 - 1] = x0i - x2i;
1835             a[j1 - 2] = x0r - x2r;
1836             a[j1 - 1] = x0i + x2i;
1837             x0r = x1r + x3i;
1838             x0i = x1i + x3r;
1839             a[j2 - 2] = Mult_W_W(wd1i, x0r) - Mult_W_W(wd1r, x0i);
1840             a[j2 - 1] = Mult_W_W(wd1i, x0i) + Mult_W_W(wd1r, x0r);
1841             x0r = x1r - x3i;
1842             x0i = x1i - x3r;
1843             a[j3 - 2] = Mult_W_W(wd3i, x0r) + Mult_W_W(wd3r, x0i);
1844             a[j3 - 1] = Mult_W_W(wd3i, x0i) - Mult_W_W(wd3r, x0r);
1845             wk1r -= Mult_W_W(ss1, wd1i);
1846             wk1i += Mult_W_W(ss1, wd1r);
1847             wk3r -= Mult_W_W(ss3, wd3i);
1848             wk3i += Mult_W_W(ss3, wd3r);
1849         }
1850         if (i0 == mh - 4) {
1851             break;
1852         }
1853     }
1854     wd1r = WR5000;
1855     j0 = mh;
1856     j1 = j0 + m;
1857     j2 = j1 + m;
1858     j3 = j2 + m;
1859     x0r = a[j0 - 2] + a[j2 - 2];
1860     x0i = -a[j0 - 1] - a[j2 - 1];
1861     x1r = a[j0 - 2] - a[j2 - 2];
1862     x1i = -a[j0 - 1] + a[j2 - 1];
1863     x2r = a[j1 - 2] + a[j3 - 2];
1864     x2i = a[j1 - 1] + a[j3 - 1];
1865     x3r = a[j1 - 2] - a[j3 - 2];
1866     x3i = a[j1 - 1] - a[j3 - 1];
1867     a[j0 - 2] = x0r + x2r;
1868     a[j0 - 1] = x0i - x2i;
1869     a[j1 - 2] = x0r - x2r;
1870     a[j1 - 1] = x0i + x2i;
1871     x0r = x1r + x3i;
1872     x0i = x1i + x3r;
1873     a[j2 - 2] = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
1874     a[j2 - 1] = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
1875     x0r = x1r - x3i;
1876     x0i = x1i - x3r;
1877     a[j3 - 2] = Mult_W_W(wk3r, x0r) + Mult_W_W(wk3i, x0i);
1878     a[j3 - 1] = Mult_W_W(wk3r, x0i) - Mult_W_W(wk3i, x0r);
1879     x0r = a[j0] + a[j2];
1880     x0i = -a[j0 + 1] - a[j2 + 1];
1881     x1r = a[j0] - a[j2];
1882     x1i = -a[j0 + 1] + a[j2 + 1];
1883     x2r = a[j1] + a[j3];
1884     x2i = a[j1 + 1] + a[j3 + 1];
1885     x3r = a[j1] - a[j3];
1886     x3i = a[j1 + 1] - a[j3 + 1];
1887     a[j0] = x0r + x2r;
1888     a[j0 + 1] = x0i - x2i;
1889     a[j1] = x0r - x2r;
1890     a[j1 + 1] = x0i + x2i;
1891     x0r = x1r + x3i;
1892     x0i = x1i + x3r;
1893     a[j2] = picofftsg_mult_w_a(wd1r, (x0r - x0i));
1894     a[j2 + 1] = picofftsg_mult_w_a(wd1r, (x0i + x0r));
1895     x0r = x1r - x3i;
1896     x0i = x1i - x3r;
1897     a[j3] = -picofftsg_mult_w_a(wd1r, (x0r + x0i));
1898     a[j3 + 1] = -picofftsg_mult_w_a(wd1r, (x0i - x0r));
1899     x0r = a[j0 + 2] + a[j2 + 2];
1900     x0i = -a[j0 + 3] - a[j2 + 3];
1901     x1r = a[j0 + 2] - a[j2 + 2];
1902     x1i = -a[j0 + 3] + a[j2 + 3];
1903     x2r = a[j1 + 2] + a[j3 + 2];
1904     x2i = a[j1 + 3] + a[j3 + 3];
1905     x3r = a[j1 + 2] - a[j3 + 2];
1906     x3i = a[j1 + 3] - a[j3 + 3];
1907     a[j0 + 2] = x0r + x2r;
1908     a[j0 + 3] = x0i - x2i;
1909     a[j1 + 2] = x0r - x2r;
1910     a[j1 + 3] = x0i + x2i;
1911     x0r = x1r + x3i;
1912     x0i = x1i + x3r;
1913     a[j2 + 2] = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
1914     a[j2 + 3] = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
1915     x0r = x1r - x3i;
1916     x0i = x1i - x3r;
1917     a[j3 + 2] = Mult_W_W(wk3i, x0r) + Mult_W_W(wk3r, x0i);
1918     a[j3 + 3] = Mult_W_W(wk3i, x0i) - Mult_W_W(wk3r, x0r);
1919 }
1920 
cftrec4(picoos_int32 n,PICOFFTSG_FFTTYPE * a)1921 void cftrec4(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
1922 {
1923     picoos_int32 isplt, j, k, m;
1924 
1925     m = n;
1926     while (m > 512) {
1927         m >>= 2;
1928         cftmdl1(m, &a[n - m]);
1929     }
1930     cftleaf(m, 1, &a[n - m]);
1931     k = 0;
1932     for (j = n - m; j > 0; j -= m) {
1933         k++;
1934         isplt = cfttree(m, j, k, a);
1935         cftleaf(m, isplt, &a[j - m]);
1936     }
1937 }
1938 
1939 
cfttree(picoos_int32 n,picoos_int32 j,picoos_int32 k,PICOFFTSG_FFTTYPE * a)1940 picoos_int32 cfttree(picoos_int32 n, picoos_int32 j, picoos_int32 k, PICOFFTSG_FFTTYPE *a)
1941 {
1942     picoos_int32 i, isplt, m;
1943 
1944     if ((k & 3) != 0) {
1945         isplt = k & 1;
1946         if (isplt != 0) {
1947             cftmdl1(n, &a[j - n]);
1948         } else {
1949             cftmdl2(n, &a[j - n]);
1950         }
1951     } else {
1952         m = n;
1953         for (i = k; (i & 3) == 0; i >>= 2) {
1954             m <<= 2;
1955         }
1956         isplt = i & 1;
1957         if (isplt != 0) {
1958             while (m > 128) {
1959                 cftmdl1(m, &a[j - m]);
1960                 m >>= 2;
1961             }
1962         } else {
1963             while (m > 128) {
1964                 cftmdl2(m, &a[j - m]);
1965                 m >>= 2;
1966             }
1967         }
1968     }
1969     return isplt;
1970 }
1971 
1972 
cftleaf(picoos_int32 n,picoos_int32 isplt,PICOFFTSG_FFTTYPE * a)1973 void cftleaf(picoos_int32 n, picoos_int32 isplt, PICOFFTSG_FFTTYPE *a)
1974 {
1975 
1976     if (n == 512) {
1977         cftmdl1(128, a);
1978         cftf161(a);
1979         cftf162(&a[32]);
1980         cftf161(&a[64]);
1981         cftf161(&a[96]);
1982         cftmdl2(128, &a[128]);
1983         cftf161(&a[128]);
1984         cftf162(&a[160]);
1985         cftf161(&a[192]);
1986         cftf162(&a[224]);
1987         cftmdl1(128, &a[256]);
1988         cftf161(&a[256]);
1989         cftf162(&a[288]);
1990         cftf161(&a[320]);
1991         cftf161(&a[352]);
1992         if (isplt != 0) {
1993             cftmdl1(128, &a[384]);
1994             cftf161(&a[480]);
1995         } else {
1996             cftmdl2(128, &a[384]);
1997             cftf162(&a[480]);
1998         }
1999         cftf161(&a[384]);
2000         cftf162(&a[416]);
2001         cftf161(&a[448]);
2002     } else {
2003         cftmdl1(64, a);
2004         cftf081(a);
2005         cftf082(&a[16]);
2006         cftf081(&a[32]);
2007         cftf081(&a[48]);
2008         cftmdl2(64, &a[64]);
2009         cftf081(&a[64]);
2010         cftf082(&a[80]);
2011         cftf081(&a[96]);
2012         cftf082(&a[112]);
2013         cftmdl1(64, &a[128]);
2014         cftf081(&a[128]);
2015         cftf082(&a[144]);
2016         cftf081(&a[160]);
2017         cftf081(&a[176]);
2018         if (isplt != 0) {
2019             cftmdl1(64, &a[192]);
2020             cftf081(&a[240]);
2021         } else {
2022             cftmdl2(64, &a[192]);
2023             cftf082(&a[240]);
2024         }
2025         cftf081(&a[192]);
2026         cftf082(&a[208]);
2027         cftf081(&a[224]);
2028     }
2029 }
2030 
2031 
cftmdl1(picoos_int32 n,PICOFFTSG_FFTTYPE * a)2032 void cftmdl1(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
2033 {
2034     picoos_int32 i, i0, j, j0, j1, j2, j3, m, mh;
2035     PICOFFTSG_FFTTYPE wk1r, wk1i, wk3r, wk3i,
2036         wd1r, wd1i, wd3r, wd3i, ss1, ss3;
2037     PICOFFTSG_FFTTYPE x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2038 
2039     mh = n >> 3;
2040     m = 2 * mh;
2041     j1 = m;
2042     j2 = j1 + m;
2043     j3 = j2 + m;
2044     x0r = a[0] + a[j2];
2045     x0i = a[1] + a[j2 + 1];
2046     x1r = a[0] - a[j2];
2047     x1i = a[1] - a[j2 + 1];
2048     x2r = a[j1] + a[j3];
2049     x2i = a[j1 + 1] + a[j3 + 1];
2050     x3r = a[j1] - a[j3];
2051     x3i = a[j1 + 1] - a[j3 + 1];
2052     a[0] = x0r + x2r;
2053     a[1] = x0i + x2i;
2054     a[j1] = x0r - x2r;
2055     a[j1 + 1] = x0i - x2i;
2056     a[j2] = x1r - x3i;
2057     a[j2 + 1] = x1i + x3r;
2058     a[j3] = x1r + x3i;
2059     a[j3 + 1] = x1i - x3r;
2060     wd1r = (PICOFFTSG_FFTTYPE)(PICODSP_WGT_SHIFT);
2061     wd1i = 0;
2062     wd3r = (PICOFFTSG_FFTTYPE)(PICODSP_WGT_SHIFT);
2063     wd3i = 0;
2064     wk1r  =  (PICOFFTSG_FFTTYPE) (0.980785250664  *PICODSP_WGT_SHIFT);
2065     wk1i  =  (PICOFFTSG_FFTTYPE) (0.195090323687  *PICODSP_WGT_SHIFT);
2066     ss1   =  (PICOFFTSG_FFTTYPE) (0.390180647373  *PICODSP_WGT_SHIFT);
2067     wk3i  =  (PICOFFTSG_FFTTYPE) (-0.555570185184 *PICODSP_WGT_SHIFT);
2068     wk3r  =  (PICOFFTSG_FFTTYPE) (0.831469595432  *PICODSP_WGT_SHIFT);
2069     ss3   =  (PICOFFTSG_FFTTYPE) (-1.111140370369 *PICODSP_WGT_SHIFT);
2070 
2071     i = 0;
2072     for (;;) {
2073         i0 = i + CDFT_LOOP_DIV_4;
2074         if (i0 > mh - 4) {
2075             i0 = mh - 4;
2076         }
2077         for (j = i + 2; j < i0; j += 4) {
2078             wd1r -= Mult_W_W(ss1, wk1i);
2079             wd1i += Mult_W_W(ss1, wk1r);
2080             wd3r -= Mult_W_W(ss3, wk3i);
2081             wd3i += Mult_W_W(ss3, wk3r);
2082             j1 = j + m;
2083             j2 = j1 + m;
2084             j3 = j2 + m;
2085             x0r = a[j] + a[j2];
2086             x0i = a[j + 1] + a[j2 + 1];
2087             x1r = a[j] - a[j2];
2088             x1i = a[j + 1] - a[j2 + 1];
2089             x2r = a[j1] + a[j3];
2090             x2i = a[j1 + 1] + a[j3 + 1];
2091             x3r = a[j1] - a[j3];
2092             x3i = a[j1 + 1] - a[j3 + 1];
2093             a[j] = x0r + x2r;
2094             a[j + 1] = x0i + x2i;
2095             a[j1] = x0r - x2r;
2096             a[j1 + 1] = x0i - x2i;
2097             x0r = x1r - x3i;
2098             x0i = x1i + x3r;
2099             a[j2] = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2100             a[j2 + 1] = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2101             x0r = x1r + x3i;
2102             x0i = x1i - x3r;
2103             a[j3] = Mult_W_W(wk3r, x0r) + Mult_W_W(wk3i, x0i);
2104             a[j3 + 1] = Mult_W_W(wk3r, x0i) - Mult_W_W(wk3i, x0r);
2105             x0r = a[j + 2] + a[j2 + 2];
2106             x0i = a[j + 3] + a[j2 + 3];
2107             x1r = a[j + 2] - a[j2 + 2];
2108             x1i = a[j + 3] - a[j2 + 3];
2109             x2r = a[j1 + 2] + a[j3 + 2];
2110             x2i = a[j1 + 3] + a[j3 + 3];
2111             x3r = a[j1 + 2] - a[j3 + 2];
2112             x3i = a[j1 + 3] - a[j3 + 3];
2113             a[j + 2] = x0r + x2r;
2114             a[j + 3] = x0i + x2i;
2115             a[j1 + 2] = x0r - x2r;
2116             a[j1 + 3] = x0i - x2i;
2117             x0r = x1r - x3i;
2118             x0i = x1i + x3r;
2119             a[j2 + 2] = Mult_W_W(wd1r, x0r) - Mult_W_W(wd1i, x0i);
2120             a[j2 + 3] = Mult_W_W(wd1r, x0i) + Mult_W_W(wd1i, x0r);
2121             x0r = x1r + x3i;
2122             x0i = x1i - x3r;
2123             a[j3 + 2] = Mult_W_W(wd3r, x0r) + Mult_W_W(wd3i, x0i);
2124             a[j3 + 3] = Mult_W_W(wd3r, x0i) - Mult_W_W(wd3i, x0r);
2125             j0 = m - j;
2126             j1 = j0 + m;
2127             j2 = j1 + m;
2128             j3 = j2 + m;
2129             x0r = a[j0] + a[j2];
2130             x0i = a[j0 + 1] + a[j2 + 1];
2131             x1r = a[j0] - a[j2];
2132             x1i = a[j0 + 1] - a[j2 + 1];
2133             x2r = a[j1] + a[j3];
2134             x2i = a[j1 + 1] + a[j3 + 1];
2135             x3r = a[j1] - a[j3];
2136             x3i = a[j1 + 1] - a[j3 + 1];
2137             a[j0] = x0r + x2r;
2138             a[j0 + 1] = x0i + x2i;
2139             a[j1] = x0r - x2r;
2140             a[j1 + 1] = x0i - x2i;
2141             x0r = x1r - x3i;
2142             x0i = x1i + x3r;
2143             a[j2] = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
2144             a[j2 + 1] = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
2145             x0r = x1r + x3i;
2146             x0i = x1i - x3r;
2147             a[j3] = Mult_W_W(wk3i, x0r) + Mult_W_W(wk3r, x0i);
2148             a[j3 + 1] = Mult_W_W(wk3i, x0i) - Mult_W_W(wk3r, x0r);
2149             x0r = a[j0 - 2] + a[j2 - 2];
2150             x0i = a[j0 - 1] + a[j2 - 1];
2151             x1r = a[j0 - 2] - a[j2 - 2];
2152             x1i = a[j0 - 1] - a[j2 - 1];
2153             x2r = a[j1 - 2] + a[j3 - 2];
2154             x2i = a[j1 - 1] + a[j3 - 1];
2155             x3r = a[j1 - 2] - a[j3 - 2];
2156             x3i = a[j1 - 1] - a[j3 - 1];
2157             a[j0 - 2] = x0r + x2r;
2158             a[j0 - 1] = x0i + x2i;
2159             a[j1 - 2] = x0r - x2r;
2160             a[j1 - 1] = x0i - x2i;
2161             x0r = x1r - x3i;
2162             x0i = x1i + x3r;
2163             a[j2 - 2] = Mult_W_W(wd1i, x0r) - Mult_W_W(wd1r, x0i);
2164             a[j2 - 1] = Mult_W_W(wd1i, x0i) + Mult_W_W(wd1r, x0r);
2165             x0r = x1r + x3i;
2166             x0i = x1i - x3r;
2167             a[j3 - 2] = Mult_W_W(wd3i, x0r) + Mult_W_W(wd3r, x0i);
2168             a[j3 - 1] = Mult_W_W(wd3i, x0i) - Mult_W_W(wd3r, x0r);
2169             wk1r -= Mult_W_W(ss1, wd1i);
2170             wk1i += Mult_W_W(ss1, wd1r);
2171             wk3r -= Mult_W_W(ss3, wd3i);
2172             wk3i += Mult_W_W(ss3, wd3r);
2173         }
2174         if (i0 == mh - 4) {
2175             break;
2176         }
2177     }
2178     wd1r = WR5000;
2179     j0 = mh;
2180     j1 = j0 + m;
2181     j2 = j1 + m;
2182     j3 = j2 + m;
2183     x0r = a[j0 - 2] + a[j2 - 2];
2184     x0i = a[j0 - 1] + a[j2 - 1];
2185     x1r = a[j0 - 2] - a[j2 - 2];
2186     x1i = a[j0 - 1] - a[j2 - 1];
2187     x2r = a[j1 - 2] + a[j3 - 2];
2188     x2i = a[j1 - 1] + a[j3 - 1];
2189     x3r = a[j1 - 2] - a[j3 - 2];
2190     x3i = a[j1 - 1] - a[j3 - 1];
2191     a[j0 - 2] = x0r + x2r;
2192     a[j0 - 1] = x0i + x2i;
2193     a[j1 - 2] = x0r - x2r;
2194     a[j1 - 1] = x0i - x2i;
2195     x0r = x1r - x3i;
2196     x0i = x1i + x3r;
2197     a[j2 - 2] = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2198     a[j2 - 1] = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2199     x0r = x1r + x3i;
2200     x0i = x1i - x3r;
2201     a[j3 - 2] = Mult_W_W(wk3r, x0r) + Mult_W_W(wk3i, x0i);
2202     a[j3 - 1] = Mult_W_W(wk3r, x0i) - Mult_W_W(wk3i, x0r);
2203     x0r = a[j0] + a[j2];
2204     x0i = a[j0 + 1] + a[j2 + 1];
2205     x1r = a[j0] - a[j2];
2206     x1i = a[j0 + 1] - a[j2 + 1];
2207     x2r = a[j1] + a[j3];
2208     x2i = a[j1 + 1] + a[j3 + 1];
2209     x3r = a[j1] - a[j3];
2210     x3i = a[j1 + 1] - a[j3 + 1];
2211     a[j0] = x0r + x2r;
2212     a[j0 + 1] = x0i + x2i;
2213     a[j1] = x0r - x2r;
2214     a[j1 + 1] = x0i - x2i;
2215     x0r = x1r - x3i;
2216     x0i = x1i + x3r;
2217     a[j2] = picofftsg_mult_w_a(wd1r, (x0r - x0i));
2218     a[j2 + 1] = picofftsg_mult_w_a(wd1r, (x0i + x0r));
2219     x0r = x1r + x3i;
2220     x0i = x1i - x3r;
2221     a[j3] = -picofftsg_mult_w_a(wd1r, (x0r + x0i));
2222     a[j3 + 1] = -picofftsg_mult_w_a(wd1r, (x0i - x0r));
2223     x0r = a[j0 + 2] + a[j2 + 2];
2224     x0i = a[j0 + 3] + a[j2 + 3];
2225     x1r = a[j0 + 2] - a[j2 + 2];
2226     x1i = a[j0 + 3] - a[j2 + 3];
2227     x2r = a[j1 + 2] + a[j3 + 2];
2228     x2i = a[j1 + 3] + a[j3 + 3];
2229     x3r = a[j1 + 2] - a[j3 + 2];
2230     x3i = a[j1 + 3] - a[j3 + 3];
2231     a[j0 + 2] = x0r + x2r;
2232     a[j0 + 3] = x0i + x2i;
2233     a[j1 + 2] = x0r - x2r;
2234     a[j1 + 3] = x0i - x2i;
2235     x0r = x1r - x3i;
2236     x0i = x1i + x3r;
2237     a[j2 + 2] = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
2238     a[j2 + 3] = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
2239     x0r = x1r + x3i;
2240     x0i = x1i - x3r;
2241     a[j3 + 2] = Mult_W_W(wk3i, x0r) + Mult_W_W(wk3r, x0i);
2242     a[j3 + 3] = Mult_W_W(wk3i, x0i) - Mult_W_W(wk3r, x0r);
2243 }
2244 
2245 
cftmdl2(picoos_int32 n,PICOFFTSG_FFTTYPE * a)2246 void cftmdl2(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
2247 {
2248     picoos_int32 i, i0, j, j0, j1, j2, j3, m, mh;
2249     PICOFFTSG_FFTTYPE wn4r, wk1r, wk1i, wk3r, wk3i,
2250         wl1r, wl1i, wl3r, wl3i, wd1r, wd1i, wd3r, wd3i,
2251         we1r, we1i, we3r, we3i, ss1, ss3;
2252     PICOFFTSG_FFTTYPE x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
2253 
2254     mh = n >> 3;
2255     m = 2 * mh;
2256     wn4r = WR5000;
2257     j1 = m;
2258     j2 = j1 + m;
2259     j3 = j2 + m;
2260     x0r = a[0] - a[j2 + 1];
2261     x0i = a[1] + a[j2];
2262     x1r = a[0] + a[j2 + 1];
2263     x1i = a[1] - a[j2];
2264     x2r = a[j1] - a[j3 + 1];
2265     x2i = a[j1 + 1] + a[j3];
2266     x3r = a[j1] + a[j3 + 1];
2267     x3i = a[j1 + 1] - a[j3];
2268     y0r = picofftsg_mult_w_a(wn4r, (x2r - x2i));
2269     y0i = picofftsg_mult_w_a(wn4r, (x2i + x2r));
2270     a[0] = x0r + y0r;
2271     a[1] = x0i + y0i;
2272     a[j1] = x0r - y0r;
2273     a[j1 + 1] = x0i - y0i;
2274     y0r = picofftsg_mult_w_a(wn4r, (x3r - x3i));
2275     y0i = picofftsg_mult_w_a(wn4r, (x3i + x3r));
2276     a[j2] = x1r - y0i;
2277     a[j2 + 1] = x1i + y0r;
2278     a[j3] = x1r + y0i;
2279     a[j3 + 1] = x1i - y0r;
2280     wl1r = PICODSP_WGT_SHIFT;
2281     wl1i = 0;
2282     wl3r = PICODSP_WGT_SHIFT;
2283     wl3i = 0;
2284     we1r = wn4r;
2285     we1i = wn4r;
2286     we3r = -wn4r;
2287     we3i = -wn4r;
2288 
2289     wk1r  =  (PICOFFTSG_FFTTYPE)(0.995184719563  *PICODSP_WGT_SHIFT);
2290     wk1i  =  (PICOFFTSG_FFTTYPE)(0.098017141223  *PICODSP_WGT_SHIFT);
2291     wd1r  =  (PICOFFTSG_FFTTYPE)(0.634393274784  *PICODSP_WGT_SHIFT);
2292     wd1i  =  (PICOFFTSG_FFTTYPE)(0.773010432720  *PICODSP_WGT_SHIFT);
2293     ss1   =  (PICOFFTSG_FFTTYPE)(0.196034282446  *PICODSP_WGT_SHIFT);
2294     wk3i  =  (PICOFFTSG_FFTTYPE)(-0.290284663439 *PICODSP_WGT_SHIFT);
2295     wk3r  =  (PICOFFTSG_FFTTYPE)(0.956940352917  *PICODSP_WGT_SHIFT);
2296     ss3   =  (PICOFFTSG_FFTTYPE)(-0.580569326878 *PICODSP_WGT_SHIFT);
2297     wd3r  =  (PICOFFTSG_FFTTYPE)(-0.881921231747 *PICODSP_WGT_SHIFT);
2298     wd3i  =  (PICOFFTSG_FFTTYPE)(-0.471396744251 *PICODSP_WGT_SHIFT);
2299 
2300     i = 0;
2301     for (;;) {
2302         i0 = i + 4 * CDFT_LOOP_DIV;
2303         if (i0 > mh - 4) {
2304             i0 = mh - 4;
2305         }
2306         for (j = i + 2; j < i0; j += 4) {
2307             wl1r -= Mult_W_W(ss1, wk1i);
2308             wl1i += Mult_W_W(ss1, wk1r);
2309             wl3r -= Mult_W_W(ss3, wk3i);
2310             wl3i += Mult_W_W(ss3, wk3r);
2311             we1r -= Mult_W_W(ss1, wd1i);
2312             we1i += Mult_W_W(ss1, wd1r);
2313             we3r -= Mult_W_W(ss3, wd3i);
2314             we3i += Mult_W_W(ss3, wd3r);
2315             j1 = j + m;
2316             j2 = j1 + m;
2317             j3 = j2 + m;
2318             x0r = a[j] - a[j2 + 1];
2319             x0i = a[j + 1] + a[j2];
2320             x1r = a[j] + a[j2 + 1];
2321             x1i = a[j + 1] - a[j2];
2322             x2r = a[j1] - a[j3 + 1];
2323             x2i = a[j1 + 1] + a[j3];
2324             x3r = a[j1] + a[j3 + 1];
2325             x3i = a[j1 + 1] - a[j3];
2326             y0r = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2327             y0i = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2328             y2r = Mult_W_W(wd1r, x2r) - Mult_W_W(wd1i, x2i);
2329             y2i = Mult_W_W(wd1r, x2i) + Mult_W_W(wd1i, x2r);
2330             a[j] = y0r + y2r;
2331             a[j + 1] = y0i + y2i;
2332             a[j1] = y0r - y2r;
2333             a[j1 + 1] = y0i - y2i;
2334             y0r = Mult_W_W(wk3r, x1r) + Mult_W_W(wk3i, x1i);
2335             y0i = Mult_W_W(wk3r, x1i) - Mult_W_W(wk3i, x1r);
2336             y2r = Mult_W_W(wd3r, x3r) + Mult_W_W(wd3i, x3i);
2337             y2i = Mult_W_W(wd3r, x3i) - Mult_W_W(wd3i, x3r);
2338             a[j2] = y0r + y2r;
2339             a[j2 + 1] = y0i + y2i;
2340             a[j3] = y0r - y2r;
2341             a[j3 + 1] = y0i - y2i;
2342             x0r = a[j + 2] - a[j2 + 3];
2343             x0i = a[j + 3] + a[j2 + 2];
2344             x1r = a[j + 2] + a[j2 + 3];
2345             x1i = a[j + 3] - a[j2 + 2];
2346             x2r = a[j1 + 2] - a[j3 + 3];
2347             x2i = a[j1 + 3] + a[j3 + 2];
2348             x3r = a[j1 + 2] + a[j3 + 3];
2349             x3i = a[j1 + 3] - a[j3 + 2];
2350             y0r = Mult_W_W(wl1r, x0r) - Mult_W_W(wl1i, x0i);
2351             y0i = Mult_W_W(wl1r, x0i) + Mult_W_W(wl1i, x0r);
2352             y2r = Mult_W_W(we1r, x2r) - Mult_W_W(we1i, x2i);
2353             y2i = Mult_W_W(we1r, x2i) + Mult_W_W(we1i, x2r);
2354             a[j + 2] = y0r + y2r;
2355             a[j + 3] = y0i + y2i;
2356             a[j1 + 2] = y0r - y2r;
2357             a[j1 + 3] = y0i - y2i;
2358             y0r = Mult_W_W(wl3r, x1r) + Mult_W_W(wl3i, x1i);
2359             y0i = Mult_W_W(wl3r, x1i) - Mult_W_W(wl3i, x1r);
2360             y2r = Mult_W_W(we3r, x3r) + Mult_W_W(we3i, x3i);
2361             y2i = Mult_W_W(we3r, x3i) - Mult_W_W(we3i, x3r);
2362             a[j2 + 2] = y0r + y2r;
2363             a[j2 + 3] = y0i + y2i;
2364             a[j3 + 2] = y0r - y2r;
2365             a[j3 + 3] = y0i - y2i;
2366             j0 = m - j;
2367             j1 = j0 + m;
2368             j2 = j1 + m;
2369             j3 = j2 + m;
2370             x0r = a[j0] - a[j2 + 1];
2371             x0i = a[j0 + 1] + a[j2];
2372             x1r = a[j0] + a[j2 + 1];
2373             x1i = a[j0 + 1] - a[j2];
2374             x2r = a[j1] - a[j3 + 1];
2375             x2i = a[j1 + 1] + a[j3];
2376             x3r = a[j1] + a[j3 + 1];
2377             x3i = a[j1 + 1] - a[j3];
2378             y0r = Mult_W_W(wd1i, x0r) - Mult_W_W(wd1r, x0i);
2379             y0i = Mult_W_W(wd1i, x0i) + Mult_W_W(wd1r, x0r);
2380             y2r = Mult_W_W(wk1i, x2r) - Mult_W_W(wk1r, x2i);
2381             y2i = Mult_W_W(wk1i, x2i) + Mult_W_W(wk1r, x2r);
2382             a[j0] = y0r + y2r;
2383             a[j0 + 1] = y0i + y2i;
2384             a[j1] = y0r - y2r;
2385             a[j1 + 1] = y0i - y2i;
2386             y0r = Mult_W_W(wd3i, x1r) + Mult_W_W(wd3r, x1i);
2387             y0i = Mult_W_W(wd3i, x1i) - Mult_W_W(wd3r, x1r);
2388             y2r = Mult_W_W(wk3i, x3r) + Mult_W_W(wk3r, x3i);
2389             y2i = Mult_W_W(wk3i, x3i) - Mult_W_W(wk3r, x3r);
2390             a[j2] = y0r + y2r;
2391             a[j2 + 1] = y0i + y2i;
2392             a[j3] = y0r - y2r;
2393             a[j3 + 1] = y0i - y2i;
2394             x0r = a[j0 - 2] - a[j2 - 1];
2395             x0i = a[j0 - 1] + a[j2 - 2];
2396             x1r = a[j0 - 2] + a[j2 - 1];
2397             x1i = a[j0 - 1] - a[j2 - 2];
2398             x2r = a[j1 - 2] - a[j3 - 1];
2399             x2i = a[j1 - 1] + a[j3 - 2];
2400             x3r = a[j1 - 2] + a[j3 - 1];
2401             x3i = a[j1 - 1] - a[j3 - 2];
2402             y0r = Mult_W_W(we1i, x0r) - Mult_W_W(we1r, x0i);
2403             y0i = Mult_W_W(we1i, x0i) + Mult_W_W(we1r, x0r);
2404             y2r = Mult_W_W(wl1i, x2r) - Mult_W_W(wl1r, x2i);
2405             y2i = Mult_W_W(wl1i, x2i) + Mult_W_W(wl1r, x2r);
2406             a[j0 - 2] = y0r + y2r;
2407             a[j0 - 1] = y0i + y2i;
2408             a[j1 - 2] = y0r - y2r;
2409             a[j1 - 1] = y0i - y2i;
2410             y0r = Mult_W_W(we3i, x1r) + Mult_W_W(we3r, x1i);
2411             y0i = Mult_W_W(we3i, x1i) - Mult_W_W(we3r, x1r);
2412             y2r = Mult_W_W(wl3i, x3r) + Mult_W_W(wl3r, x3i);
2413             y2i = Mult_W_W(wl3i, x3i) - Mult_W_W(wl3r, x3r);
2414             a[j2 - 2] = y0r + y2r;
2415             a[j2 - 1] = y0i + y2i;
2416             a[j3 - 2] = y0r - y2r;
2417             a[j3 - 1] = y0i - y2i;
2418             wk1r -= Mult_W_W(ss1, wl1i);
2419             wk1i += Mult_W_W(ss1, wl1r);
2420             wk3r -= Mult_W_W(ss3, wl3i);
2421             wk3i += Mult_W_W(ss3, wl3r);
2422             wd1r -= Mult_W_W(ss1, we1i);
2423             wd1i += Mult_W_W(ss1, we1r);
2424             wd3r -= Mult_W_W(ss3, we3i);
2425             wd3i += Mult_W_W(ss3, we3r);
2426         }
2427         if (i0 == mh - 4) {
2428             break;
2429         }
2430     }
2431     wl1r = WR2500;
2432     wl1i = WI2500;
2433     j0 = mh;
2434     j1 = j0 + m;
2435     j2 = j1 + m;
2436     j3 = j2 + m;
2437     x0r = a[j0 - 2] - a[j2 - 1];
2438     x0i = a[j0 - 1] + a[j2 - 2];
2439     x1r = a[j0 - 2] + a[j2 - 1];
2440     x1i = a[j0 - 1] - a[j2 - 2];
2441     x2r = a[j1 - 2] - a[j3 - 1];
2442     x2i = a[j1 - 1] + a[j3 - 2];
2443     x3r = a[j1 - 2] + a[j3 - 1];
2444     x3i = a[j1 - 1] - a[j3 - 2];
2445     y0r = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2446     y0i = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2447     y2r = Mult_W_W(wd1r, x2r) - Mult_W_W(wd1i, x2i);
2448     y2i = Mult_W_W(wd1r, x2i) + Mult_W_W(wd1i, x2r);
2449     a[j0 - 2] = y0r + y2r;
2450     a[j0 - 1] = y0i + y2i;
2451     a[j1 - 2] = y0r - y2r;
2452     a[j1 - 1] = y0i - y2i;
2453     y0r = Mult_W_W(wk3r, x1r) + Mult_W_W(wk3i, x1i);
2454     y0i = Mult_W_W(wk3r, x1i) - Mult_W_W(wk3i, x1r);
2455     y2r = Mult_W_W(wd3r, x3r) + Mult_W_W(wd3i, x3i);
2456     y2i = Mult_W_W(wd3r, x3i) - Mult_W_W(wd3i, x3r);
2457     a[j2 - 2] = y0r + y2r;
2458     a[j2 - 1] = y0i + y2i;
2459     a[j3 - 2] = y0r - y2r;
2460     a[j3 - 1] = y0i - y2i;
2461     x0r = a[j0] - a[j2 + 1];
2462     x0i = a[j0 + 1] + a[j2];
2463     x1r = a[j0] + a[j2 + 1];
2464     x1i = a[j0 + 1] - a[j2];
2465     x2r = a[j1] - a[j3 + 1];
2466     x2i = a[j1 + 1] + a[j3];
2467     x3r = a[j1] + a[j3 + 1];
2468     x3i = a[j1 + 1] - a[j3];
2469     y0r = Mult_W_W(wl1r, x0r) - Mult_W_W(wl1i, x0i);
2470     y0i = Mult_W_W(wl1r, x0i) + Mult_W_W(wl1i, x0r);
2471     y2r = Mult_W_W(wl1i, x2r) - Mult_W_W(wl1r, x2i);
2472     y2i = Mult_W_W(wl1i, x2i) + Mult_W_W(wl1r, x2r);
2473     a[j0] = y0r + y2r;
2474     a[j0 + 1] = y0i + y2i;
2475     a[j1] = y0r - y2r;
2476     a[j1 + 1] = y0i - y2i;
2477     y0r = Mult_W_W(wl1i, x1r) - Mult_W_W(wl1r, x1i);
2478     y0i = Mult_W_W(wl1i, x1i) + Mult_W_W(wl1r, x1r);
2479     y2r = Mult_W_W(wl1r, x3r) - Mult_W_W(wl1i, x3i);
2480     y2i = Mult_W_W(wl1r, x3i) + Mult_W_W(wl1i, x3r);
2481     a[j2] = y0r - y2r;
2482     a[j2 + 1] = y0i - y2i;
2483     a[j3] = y0r + y2r;
2484     a[j3 + 1] = y0i + y2i;
2485     x0r = a[j0 + 2] - a[j2 + 3];
2486     x0i = a[j0 + 3] + a[j2 + 2];
2487     x1r = a[j0 + 2] + a[j2 + 3];
2488     x1i = a[j0 + 3] - a[j2 + 2];
2489     x2r = a[j1 + 2] - a[j3 + 3];
2490     x2i = a[j1 + 3] + a[j3 + 2];
2491     x3r = a[j1 + 2] + a[j3 + 3];
2492     x3i = a[j1 + 3] - a[j3 + 2];
2493     y0r = Mult_W_W(wd1i, x0r) - Mult_W_W(wd1r, x0i);
2494     y0i = Mult_W_W(wd1i, x0i) + Mult_W_W(wd1r, x0r);
2495     y2r = Mult_W_W(wk1i, x2r) - Mult_W_W(wk1r, x2i);
2496     y2i = Mult_W_W(wk1i, x2i) + Mult_W_W(wk1r, x2r);
2497     a[j0 + 2] = y0r + y2r;
2498     a[j0 + 3] = y0i + y2i;
2499     a[j1 + 2] = y0r - y2r;
2500     a[j1 + 3] = y0i - y2i;
2501     y0r = Mult_W_W(wd3i, x1r) + Mult_W_W(wd3r, x1i);
2502     y0i = Mult_W_W(wd3i, x1i) - Mult_W_W(wd3r, x1r);
2503     y2r = Mult_W_W(wk3i, x3r) + Mult_W_W(wk3r, x3i);
2504     y2i = Mult_W_W(wk3i, x3i) - Mult_W_W(wk3r, x3r);
2505     a[j2 + 2] = y0r + y2r;
2506     a[j2 + 3] = y0i + y2i;
2507     a[j3 + 2] = y0r - y2r;
2508     a[j3 + 3] = y0i - y2i;
2509 }
2510 
2511 
cftfx41(picoos_int32 n,PICOFFTSG_FFTTYPE * a)2512 void cftfx41(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
2513 {
2514 
2515     if (n == 128) {
2516         cftf161(a);
2517         cftf162(&a[32]);
2518         cftf161(&a[64]);
2519         cftf161(&a[96]);
2520     } else {
2521         cftf081(a);
2522         cftf082(&a[16]);
2523         cftf081(&a[32]);
2524         cftf081(&a[48]);
2525     }
2526 }
2527 
2528 
cftf161(PICOFFTSG_FFTTYPE * a)2529 void cftf161(PICOFFTSG_FFTTYPE *a)
2530 {
2531     PICOFFTSG_FFTTYPE wn4r, wk1r, wk1i,
2532         x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
2533         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2534         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2535         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2536         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2537 
2538     wn4r = WR5000;
2539     wk1r = WR2500;
2540     wk1i = WI2500;
2541     x0r = a[0] + a[16];
2542     x0i = a[1] + a[17];
2543     x1r = a[0] - a[16];
2544     x1i = a[1] - a[17];
2545     x2r = a[8] + a[24];
2546     x2i = a[9] + a[25];
2547     x3r = a[8] - a[24];
2548     x3i = a[9] - a[25];
2549     y0r = x0r + x2r;
2550     y0i = x0i + x2i;
2551     y4r = x0r - x2r;
2552     y4i = x0i - x2i;
2553     y8r = x1r - x3i;
2554     y8i = x1i + x3r;
2555     y12r = x1r + x3i;
2556     y12i = x1i - x3r;
2557     x0r = a[2] + a[18];
2558     x0i = a[3] + a[19];
2559     x1r = a[2] - a[18];
2560     x1i = a[3] - a[19];
2561     x2r = a[10] + a[26];
2562     x2i = a[11] + a[27];
2563     x3r = a[10] - a[26];
2564     x3i = a[11] - a[27];
2565     y1r = x0r + x2r;
2566     y1i = x0i + x2i;
2567     y5r = x0r - x2r;
2568     y5i = x0i - x2i;
2569     x0r = x1r - x3i;
2570     x0i = x1i + x3r;
2571     y9r = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2572     y9i = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2573     x0r = x1r + x3i;
2574     x0i = x1i - x3r;
2575     y13r = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
2576     y13i = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
2577     x0r = a[4] + a[20];
2578     x0i = a[5] + a[21];
2579     x1r = a[4] - a[20];
2580     x1i = a[5] - a[21];
2581     x2r = a[12] + a[28];
2582     x2i = a[13] + a[29];
2583     x3r = a[12] - a[28];
2584     x3i = a[13] - a[29];
2585     y2r = x0r + x2r;
2586     y2i = x0i + x2i;
2587     y6r = x0r - x2r;
2588     y6i = x0i - x2i;
2589     x0r = x1r - x3i;
2590     x0i = x1i + x3r;
2591     y10r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2592     y10i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2593     x0r = x1r + x3i;
2594     x0i = x1i - x3r;
2595     y14r = picofftsg_mult_w_a(wn4r, (x0r + x0i));
2596     y14i = picofftsg_mult_w_a(wn4r, (x0i - x0r));
2597     x0r = a[6] + a[22];
2598     x0i = a[7] + a[23];
2599     x1r = a[6] - a[22];
2600     x1i = a[7] - a[23];
2601     x2r = a[14] + a[30];
2602     x2i = a[15] + a[31];
2603     x3r = a[14] - a[30];
2604     x3i = a[15] - a[31];
2605     y3r = x0r + x2r;
2606     y3i = x0i + x2i;
2607     y7r = x0r - x2r;
2608     y7i = x0i - x2i;
2609     x0r = x1r - x3i;
2610     x0i = x1i + x3r;
2611     y11r = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
2612     y11i = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
2613     x0r = x1r + x3i;
2614     x0i = x1i - x3r;
2615     y15r = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2616     y15i = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2617     x0r = y12r - y14r;
2618     x0i = y12i - y14i;
2619     x1r = y12r + y14r;
2620     x1i = y12i + y14i;
2621     x2r = y13r - y15r;
2622     x2i = y13i - y15i;
2623     x3r = y13r + y15r;
2624     x3i = y13i + y15i;
2625     a[24] = x0r + x2r;
2626     a[25] = x0i + x2i;
2627     a[26] = x0r - x2r;
2628     a[27] = x0i - x2i;
2629     a[28] = x1r - x3i;
2630     a[29] = x1i + x3r;
2631     a[30] = x1r + x3i;
2632     a[31] = x1i - x3r;
2633     x0r = y8r + y10r;
2634     x0i = y8i + y10i;
2635     x1r = y8r - y10r;
2636     x1i = y8i - y10i;
2637     x2r = y9r + y11r;
2638     x2i = y9i + y11i;
2639     x3r = y9r - y11r;
2640     x3i = y9i - y11i;
2641     a[16] = x0r + x2r;
2642     a[17] = x0i + x2i;
2643     a[18] = x0r - x2r;
2644     a[19] = x0i - x2i;
2645     a[20] = x1r - x3i;
2646     a[21] = x1i + x3r;
2647     a[22] = x1r + x3i;
2648     a[23] = x1i - x3r;
2649     x0r = y5r - y7i;
2650     x0i = y5i + y7r;
2651     x2r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2652     x2i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2653     x0r = y5r + y7i;
2654     x0i = y5i - y7r;
2655     x3r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2656     x3i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2657     x0r = y4r - y6i;
2658     x0i = y4i + y6r;
2659     x1r = y4r + y6i;
2660     x1i = y4i - y6r;
2661     a[8] = x0r + x2r;
2662     a[9] = x0i + x2i;
2663     a[10] = x0r - x2r;
2664     a[11] = x0i - x2i;
2665     a[12] = x1r - x3i;
2666     a[13] = x1i + x3r;
2667     a[14] = x1r + x3i;
2668     a[15] = x1i - x3r;
2669     x0r = y0r + y2r;
2670     x0i = y0i + y2i;
2671     x1r = y0r - y2r;
2672     x1i = y0i - y2i;
2673     x2r = y1r + y3r;
2674     x2i = y1i + y3i;
2675     x3r = y1r - y3r;
2676     x3i = y1i - y3i;
2677     a[0] = x0r + x2r;
2678     a[1] = x0i + x2i;
2679     a[2] = x0r - x2r;
2680     a[3] = x0i - x2i;
2681     a[4] = x1r - x3i;
2682     a[5] = x1i + x3r;
2683     a[6] = x1r + x3i;
2684     a[7] = x1i - x3r;
2685 }
2686 
2687 
cftf162(PICOFFTSG_FFTTYPE * a)2688 void cftf162(PICOFFTSG_FFTTYPE *a)
2689 {
2690     PICOFFTSG_FFTTYPE wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
2691         x0r, x0i, x1r, x1i, x2r, x2i,
2692         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2693         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2694         y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2695         y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2696 
2697     wn4r = WR5000;
2698     wk1r = WR1250;
2699     wk1i = WI1250;
2700     wk2r = WR2500;
2701     wk2i = WI2500;
2702     wk3r = WR3750;
2703     wk3i = WI3750;
2704     x1r = a[0] - a[17];
2705     x1i = a[1] + a[16];
2706     x0r = a[8] - a[25];
2707     x0i = a[9] + a[24];
2708     x2r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2709     x2i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2710     y0r = x1r + x2r;
2711     y0i = x1i + x2i;
2712     y4r = x1r - x2r;
2713     y4i = x1i - x2i;
2714     x1r = a[0] + a[17];
2715     x1i = a[1] - a[16];
2716     x0r = a[8] + a[25];
2717     x0i = a[9] - a[24];
2718     x2r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2719     x2i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2720     y8r = x1r - x2i;
2721     y8i = x1i + x2r;
2722     y12r = x1r + x2i;
2723     y12i = x1i - x2r;
2724     x0r = a[2] - a[19];
2725     x0i = a[3] + a[18];
2726     x1r = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2727     x1i = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2728     x0r = a[10] - a[27];
2729     x0i = a[11] + a[26];
2730     x2r = Mult_W_W(wk3i, x0r) - Mult_W_W(wk3r, x0i);
2731     x2i = Mult_W_W(wk3i, x0i) + Mult_W_W(wk3r, x0r);
2732     y1r = x1r + x2r;
2733     y1i = x1i + x2i;
2734     y5r = x1r - x2r;
2735     y5i = x1i - x2i;
2736     x0r = a[2] + a[19];
2737     x0i = a[3] - a[18];
2738     x1r = Mult_W_W(wk3r, x0r) - Mult_W_W(wk3i, x0i);
2739     x1i = Mult_W_W(wk3r, x0i) + Mult_W_W(wk3i, x0r);
2740     x0r = a[10] + a[27];
2741     x0i = a[11] - a[26];
2742     x2r = Mult_W_W(wk1r, x0r) + Mult_W_W(wk1i, x0i);
2743     x2i = Mult_W_W(wk1r, x0i) - Mult_W_W(wk1i, x0r);
2744     y9r = x1r - x2r;
2745     y9i = x1i - x2i;
2746     y13r = x1r + x2r;
2747     y13i = x1i + x2i;
2748     x0r = a[4] - a[21];
2749     x0i = a[5] + a[20];
2750     x1r = Mult_W_W(wk2r, x0r) - Mult_W_W(wk2i, x0i);
2751     x1i = Mult_W_W(wk2r, x0i) + Mult_W_W(wk2i, x0r);
2752     x0r = a[12] - a[29];
2753     x0i = a[13] + a[28];
2754     x2r = Mult_W_W(wk2i, x0r) - Mult_W_W(wk2r, x0i);
2755     x2i = Mult_W_W(wk2i, x0i) + Mult_W_W(wk2r, x0r);
2756     y2r = x1r + x2r;
2757     y2i = x1i + x2i;
2758     y6r = x1r - x2r;
2759     y6i = x1i - x2i;
2760     x0r = a[4] + a[21];
2761     x0i = a[5] - a[20];
2762     x1r = Mult_W_W(wk2i, x0r) - Mult_W_W(wk2r, x0i);
2763     x1i = Mult_W_W(wk2i, x0i) + Mult_W_W(wk2r, x0r);
2764     x0r = a[12] + a[29];
2765     x0i = a[13] - a[28];
2766     x2r = Mult_W_W(wk2r, x0r) - Mult_W_W(wk2i, x0i);
2767     x2i = Mult_W_W(wk2r, x0i) + Mult_W_W(wk2i, x0r);
2768     y10r = x1r - x2r;
2769     y10i = x1i - x2i;
2770     y14r = x1r + x2r;
2771     y14i = x1i + x2i;
2772     x0r = a[6] - a[23];
2773     x0i = a[7] + a[22];
2774     x1r = Mult_W_W(wk3r, x0r) - Mult_W_W(wk3i, x0i);
2775     x1i = Mult_W_W(wk3r, x0i) + Mult_W_W(wk3i, x0r);
2776     x0r = a[14] - a[31];
2777     x0i = a[15] + a[30];
2778     x2r = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
2779     x2i = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
2780     y3r = x1r + x2r;
2781     y3i = x1i + x2i;
2782     y7r = x1r - x2r;
2783     y7i = x1i - x2i;
2784     x0r = a[6] + a[23];
2785     x0i = a[7] - a[22];
2786     x1r = Mult_W_W(wk1i, x0r) + Mult_W_W(wk1r, x0i);
2787     x1i = Mult_W_W(wk1i, x0i) - Mult_W_W(wk1r, x0r);
2788     x0r = a[14] + a[31];
2789     x0i = a[15] - a[30];
2790     x2r = Mult_W_W(wk3i, x0r) - Mult_W_W(wk3r, x0i);
2791     x2i = Mult_W_W(wk3i, x0i) + Mult_W_W(wk3r, x0r);
2792     y11r = x1r + x2r;
2793     y11i = x1i + x2i;
2794     y15r = x1r - x2r;
2795     y15i = x1i - x2i;
2796     x1r = y0r + y2r;
2797     x1i = y0i + y2i;
2798     x2r = y1r + y3r;
2799     x2i = y1i + y3i;
2800     a[0] = x1r + x2r;
2801     a[1] = x1i + x2i;
2802     a[2] = x1r - x2r;
2803     a[3] = x1i - x2i;
2804     x1r = y0r - y2r;
2805     x1i = y0i - y2i;
2806     x2r = y1r - y3r;
2807     x2i = y1i - y3i;
2808     a[4] = x1r - x2i;
2809     a[5] = x1i + x2r;
2810     a[6] = x1r + x2i;
2811     a[7] = x1i - x2r;
2812     x1r = y4r - y6i;
2813     x1i = y4i + y6r;
2814     x0r = y5r - y7i;
2815     x0i = y5i + y7r;
2816     x2r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2817     x2i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2818     a[8] = x1r + x2r;
2819     a[9] = x1i + x2i;
2820     a[10] = x1r - x2r;
2821     a[11] = x1i - x2i;
2822     x1r = y4r + y6i;
2823     x1i = y4i - y6r;
2824     x0r = y5r + y7i;
2825     x0i = y5i - y7r;
2826     x2r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2827     x2i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2828     a[12] = x1r - x2i;
2829     a[13] = x1i + x2r;
2830     a[14] = x1r + x2i;
2831     a[15] = x1i - x2r;
2832     x1r = y8r + y10r;
2833     x1i = y8i + y10i;
2834     x2r = y9r - y11r;
2835     x2i = y9i - y11i;
2836     a[16] = x1r + x2r;
2837     a[17] = x1i + x2i;
2838     a[18] = x1r - x2r;
2839     a[19] = x1i - x2i;
2840     x1r = y8r - y10r;
2841     x1i = y8i - y10i;
2842     x2r = y9r + y11r;
2843     x2i = y9i + y11i;
2844     a[20] = x1r - x2i;
2845     a[21] = x1i + x2r;
2846     a[22] = x1r + x2i;
2847     a[23] = x1i - x2r;
2848     x1r = y12r - y14i;
2849     x1i = y12i + y14r;
2850     x0r = y13r + y15i;
2851     x0i = y13i - y15r;
2852     x2r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2853     x2i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2854     a[24] = x1r + x2r;
2855     a[25] = x1i + x2i;
2856     a[26] = x1r - x2r;
2857     a[27] = x1i - x2i;
2858     x1r = y12r + y14i;
2859     x1i = y12i - y14r;
2860     x0r = y13r - y15i;
2861     x0i = y13i + y15r;
2862     x2r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2863     x2i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2864     a[28] = x1r - x2i;
2865     a[29] = x1i + x2r;
2866     a[30] = x1r + x2i;
2867     a[31] = x1i - x2r;
2868 }
2869 
2870 
cftf081(PICOFFTSG_FFTTYPE * a)2871 void cftf081(PICOFFTSG_FFTTYPE *a)
2872 {
2873     PICOFFTSG_FFTTYPE wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
2874         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2875         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
2876 
2877     wn4r = WR5000;
2878     x0r = a[0] + a[8];
2879     x0i = a[1] + a[9];
2880     x1r = a[0] - a[8];
2881     x1i = a[1] - a[9];
2882     x2r = a[4] + a[12];
2883     x2i = a[5] + a[13];
2884     x3r = a[4] - a[12];
2885     x3i = a[5] - a[13];
2886     y0r = x0r + x2r;
2887     y0i = x0i + x2i;
2888     y2r = x0r - x2r;
2889     y2i = x0i - x2i;
2890     y1r = x1r - x3i;
2891     y1i = x1i + x3r;
2892     y3r = x1r + x3i;
2893     y3i = x1i - x3r;
2894     x0r = a[2] + a[10];
2895     x0i = a[3] + a[11];
2896     x1r = a[2] - a[10];
2897     x1i = a[3] - a[11];
2898     x2r = a[6] + a[14];
2899     x2i = a[7] + a[15];
2900     x3r = a[6] - a[14];
2901     x3i = a[7] - a[15];
2902     y4r = x0r + x2r;
2903     y4i = x0i + x2i;
2904     y6r = x0r - x2r;
2905     y6i = x0i - x2i;
2906     x0r = x1r - x3i;
2907     x0i = x1i + x3r;
2908     x2r = x1r + x3i;
2909     x2i = x1i - x3r;
2910     y5r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2911     y5i = picofftsg_mult_w_a(wn4r, (x0r + x0i));
2912     y7r = picofftsg_mult_w_a(wn4r, (x2r - x2i));
2913     y7i = picofftsg_mult_w_a(wn4r, (x2r + x2i));
2914     a[8] = y1r + y5r;
2915     a[9] = y1i + y5i;
2916     a[10] = y1r - y5r;
2917     a[11] = y1i - y5i;
2918     a[12] = y3r - y7i;
2919     a[13] = y3i + y7r;
2920     a[14] = y3r + y7i;
2921     a[15] = y3i - y7r;
2922     a[0] = y0r + y4r;
2923     a[1] = y0i + y4i;
2924     a[2] = y0r - y4r;
2925     a[3] = y0i - y4i;
2926     a[4] = y2r - y6i;
2927     a[5] = y2i + y6r;
2928     a[6] = y2r + y6i;
2929     a[7] = y2i - y6r;
2930 }
2931 
2932 
cftf082(PICOFFTSG_FFTTYPE * a)2933 void cftf082(PICOFFTSG_FFTTYPE *a)
2934 {
2935     PICOFFTSG_FFTTYPE wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i,
2936         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2937         y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
2938 
2939     wn4r = WR5000;
2940     wk1r = WR2500;
2941     wk1i = WI2500;
2942     y0r = a[0] - a[9];
2943     y0i = a[1] + a[8];
2944     y1r = a[0] + a[9];
2945     y1i = a[1] - a[8];
2946     x0r = a[4] - a[13];
2947     x0i = a[5] + a[12];
2948     y2r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2949     y2i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2950     x0r = a[4] + a[13];
2951     x0i = a[5] - a[12];
2952     y3r = picofftsg_mult_w_a(wn4r, (x0r - x0i));
2953     y3i = picofftsg_mult_w_a(wn4r, (x0i + x0r));
2954     x0r = a[2] - a[11];
2955     x0i = a[3] + a[10];
2956     y4r = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2957     y4i = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2958     x0r = a[2] + a[11];
2959     x0i = a[3] - a[10];
2960     y5r = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
2961     y5i = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
2962     x0r = a[6] - a[15];
2963     x0i = a[7] + a[14];
2964     y6r = Mult_W_W(wk1i, x0r) - Mult_W_W(wk1r, x0i);
2965     y6i = Mult_W_W(wk1i, x0i) + Mult_W_W(wk1r, x0r);
2966     x0r = a[6] + a[15];
2967     x0i = a[7] - a[14];
2968     y7r = Mult_W_W(wk1r, x0r) - Mult_W_W(wk1i, x0i);
2969     y7i = Mult_W_W(wk1r, x0i) + Mult_W_W(wk1i, x0r);
2970     x0r = y0r + y2r;
2971     x0i = y0i + y2i;
2972     x1r = y4r + y6r;
2973     x1i = y4i + y6i;
2974     a[0] = x0r + x1r;
2975     a[1] = x0i + x1i;
2976     a[2] = x0r - x1r;
2977     a[3] = x0i - x1i;
2978     x0r = y0r - y2r;
2979     x0i = y0i - y2i;
2980     x1r = y4r - y6r;
2981     x1i = y4i - y6i;
2982     a[4] = x0r - x1i;
2983     a[5] = x0i + x1r;
2984     a[6] = x0r + x1i;
2985     a[7] = x0i - x1r;
2986     x0r = y1r - y3i;
2987     x0i = y1i + y3r;
2988     x1r = y5r - y7r;
2989     x1i = y5i - y7i;
2990     a[8] = x0r + x1r;
2991     a[9] = x0i + x1i;
2992     a[10] = x0r - x1r;
2993     a[11] = x0i - x1i;
2994     x0r = y1r + y3i;
2995     x0i = y1i - y3r;
2996     x1r = y5r + y7r;
2997     x1i = y5i + y7i;
2998     a[12] = x0r - x1i;
2999     a[13] = x0i + x1r;
3000     a[14] = x0r + x1i;
3001     a[15] = x0i - x1r;
3002 }
3003 
3004 
cftf040(PICOFFTSG_FFTTYPE * a)3005 void cftf040(PICOFFTSG_FFTTYPE *a)
3006 {
3007     PICOFFTSG_FFTTYPE x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3008 
3009     x0r = a[0] + a[4];
3010     x0i = a[1] + a[5];
3011     x1r = a[0] - a[4];
3012     x1i = a[1] - a[5];
3013     x2r = a[2] + a[6];
3014     x2i = a[3] + a[7];
3015     x3r = a[2] - a[6];
3016     x3i = a[3] - a[7];
3017     a[0] = x0r + x2r;
3018     a[1] = x0i + x2i;
3019     a[2] = x1r - x3i;
3020     a[3] = x1i + x3r;
3021     a[4] = x0r - x2r;
3022     a[5] = x0i - x2i;
3023     a[6] = x1r + x3i;
3024     a[7] = x1i - x3r;
3025 }
3026 
3027 
cftb040(PICOFFTSG_FFTTYPE * a)3028 void cftb040(PICOFFTSG_FFTTYPE *a)
3029 {
3030     PICOFFTSG_FFTTYPE x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3031 
3032     x0r = a[0] + a[4];
3033     x0i = a[1] + a[5];
3034     x1r = a[0] - a[4];
3035     x1i = a[1] - a[5];
3036     x2r = a[2] + a[6];
3037     x2i = a[3] + a[7];
3038     x3r = a[2] - a[6];
3039     x3i = a[3] - a[7];
3040     a[0] = x0r + x2r;
3041     a[1] = x0i + x2i;
3042     a[2] = x1r + x3i;
3043     a[3] = x1i - x3r;
3044     a[4] = x0r - x2r;
3045     a[5] = x0i - x2i;
3046     a[6] = x1r - x3i;
3047     a[7] = x1i + x3r;
3048 }
3049 
3050 
cftx020(PICOFFTSG_FFTTYPE * a)3051 void cftx020(PICOFFTSG_FFTTYPE *a)
3052 {
3053     PICOFFTSG_FFTTYPE x0r, x0i;
3054 
3055     x0r = a[0] - a[2];
3056     x0i = a[1] - a[3];
3057     a[0] += a[2];
3058     a[1] += a[3];
3059     a[2] = x0r;
3060     a[3] = x0i;
3061 }
3062 
3063 
rftfsub(picoos_int32 n,PICOFFTSG_FFTTYPE * a)3064 void rftfsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
3065 {
3066     picoos_int32 i, i0, j, k;
3067     PICOFFTSG_FFTTYPE w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
3068 
3069     wkr = 0;
3070     wki = 0;
3071 
3072     switch (n) {
3073         case 8 :
3074             wdi=(PICOFFTSG_FFTTYPE)(0.353553414345*PICODSP_WGT_SHIFT); wdr=(PICOFFTSG_FFTTYPE)(0.146446630359*PICODSP_WGT_SHIFT);
3075             w1r=(PICOFFTSG_FFTTYPE)(0.707106709480*PICODSP_WGT_SHIFT); w1i=(PICOFFTSG_FFTTYPE)(0.707106828690*PICODSP_WGT_SHIFT);
3076             ss =(PICOFFTSG_FFTTYPE)(1.414213657379*PICODSP_WGT_SHIFT); break;
3077         case 16 :
3078             wdi=(PICOFFTSG_FFTTYPE)(0.191341713071*PICODSP_WGT_SHIFT); wdr=(PICOFFTSG_FFTTYPE)(0.038060232997*PICODSP_WGT_SHIFT);
3079             w1r=(PICOFFTSG_FFTTYPE)(0.923879504204*PICODSP_WGT_SHIFT); w1i=(PICOFFTSG_FFTTYPE)(0.382683426142*PICODSP_WGT_SHIFT);
3080             ss =(PICOFFTSG_FFTTYPE)(0.765366852283*PICODSP_WGT_SHIFT); break;
3081         case 32 :
3082             wdi=(PICOFFTSG_FFTTYPE)(0.097545161843*PICODSP_WGT_SHIFT); wdr=(PICOFFTSG_FFTTYPE)(0.009607359767*PICODSP_WGT_SHIFT);
3083             w1r=(PICOFFTSG_FFTTYPE)(0.980785250664*PICODSP_WGT_SHIFT); w1i=(PICOFFTSG_FFTTYPE)(0.195090323687*PICODSP_WGT_SHIFT);
3084             ss =(PICOFFTSG_FFTTYPE)(0.390180647373*PICODSP_WGT_SHIFT); break;
3085         case 64 :
3086             wdi=(PICOFFTSG_FFTTYPE)(0.049008570611*PICODSP_WGT_SHIFT); wdr=(PICOFFTSG_FFTTYPE)(0.002407636726*PICODSP_WGT_SHIFT);
3087             w1r=(PICOFFTSG_FFTTYPE)(0.995184719563*PICODSP_WGT_SHIFT); w1i=(PICOFFTSG_FFTTYPE)(0.098017141223*PICODSP_WGT_SHIFT);
3088             ss =(PICOFFTSG_FFTTYPE)(0.196034282446*PICODSP_WGT_SHIFT); break;
3089         default :
3090             wdr = 0; wdi = 0; ss = 0;
3091             break;
3092     }
3093 
3094     i = n >> 1;
3095     for (;;) {
3096         i0 = i - RDFT_LOOP_DIV_4;
3097         if (i0 < 4) {
3098             i0 = 4;
3099         }
3100         for (j = i - 4; j >= i0; j -= 4) {
3101             k = n - j;
3102             xr = a[j + 2] - a[k - 2];
3103             xi = a[j + 3] + a[k - 1];
3104             yr = Mult_W_W(wdr, xr) - Mult_W_W(wdi, xi);
3105             yi = Mult_W_W(wdr, xi) + Mult_W_W(wdi, xr);
3106             a[j + 2] -= yr;
3107             a[j + 3] -= yi;
3108             a[k - 2] += yr;
3109             a[k - 1] -= yi;
3110             wkr += Mult_W_W(ss, wdi);
3111             wki += picofftsg_mult_w_w(ss, (PICOFFTSG_WGT_SHIFT2 - wdr));
3112             xr = a[j] - a[k];
3113             xi = a[j + 1] + a[k + 1];
3114             yr = Mult_W_W(wkr, xr) - Mult_W_W(wki, xi);
3115             yi = Mult_W_W(wkr, xi) + Mult_W_W(wki, xr);
3116             a[j] -= yr;
3117             a[j + 1] -= yi;
3118             a[k] += yr;
3119             a[k + 1] -= yi;
3120             wdr += Mult_W_W(ss, wki);
3121             wdi += picofftsg_mult_w_w(ss, (PICOFFTSG_WGT_SHIFT2 - wkr));
3122         }
3123         if (i0 == 4) {
3124             break;
3125         }
3126     }
3127 
3128     xr = a[2] - a[n - 2];
3129     xi = a[3] + a[n - 1];
3130     yr = Mult_W_W(wdr, xr) - Mult_W_W(wdi, xi);
3131     yi = Mult_W_W(wdr, xi) + Mult_W_W(wdi, xr);
3132 
3133     a[2] -= yr;
3134     a[3] -= yi;
3135     a[n - 2] += yr;
3136     a[n - 1] -= yi;
3137 
3138 
3139 }
3140 
3141 
rftbsub(picoos_int32 n,PICOFFTSG_FFTTYPE * a)3142 void rftbsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
3143 {
3144     picoos_int32 i, i0, j, k;
3145     PICOFFTSG_FFTTYPE w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
3146     wkr = 0;
3147     wki = 0;
3148     wdi=(PICOFFTSG_FFTTYPE)(0.012270614505*PICODSP_WGT_SHIFT);
3149     wdr=(PICOFFTSG_FFTTYPE)(0.000150590655*PICODSP_WGT_SHIFT);
3150     w1r=(PICOFFTSG_FFTTYPE)(0.999698817730*PICODSP_WGT_SHIFT);
3151     w1i=(PICOFFTSG_FFTTYPE)(0.024541229010*PICODSP_WGT_SHIFT);
3152     ss=(PICOFFTSG_FFTTYPE)(0.049082458019*PICODSP_WGT_SHIFT);
3153 
3154     i = n >> 1;
3155     for (;;) {
3156         i0 = i - RDFT_LOOP_DIV4;
3157         if (i0 < 4) {
3158             i0 = 4;
3159         }
3160         for (j = i - 4; j >= i0; j -= 4) {
3161             k = n - j;
3162             xr = a[j + 2] - a[k - 2];
3163             xi = a[j + 3] + a[k - 1];
3164             yr = Mult_W_W(wdr, xr) + Mult_W_W(wdi, xi);
3165             yi = Mult_W_W(wdr, xi) - Mult_W_W(wdi, xr);
3166             a[j + 2] -= yr;
3167             a[j + 3] -= yi;
3168             a[k - 2] += yr;
3169             a[k - 1] -= yi;
3170             wkr += Mult_W_W(ss, wdi);
3171             wki += picofftsg_mult_w_w(ss, (PICOFFTSG_WGT_SHIFT2 - wdr));
3172             xr = a[j] - a[k];
3173             xi = a[j + 1] + a[k + 1];
3174             yr = Mult_W_W(wkr, xr) + Mult_W_W(wki, xi);
3175             yi = Mult_W_W(wkr, xi) - Mult_W_W(wki, xr);
3176             a[j] -= yr;
3177             a[j + 1] -= yi;
3178             a[k] += yr;
3179             a[k + 1] -= yi;
3180             wdr += Mult_W_W(ss, wki);
3181             wdi += picofftsg_mult_w_w(ss, (PICOFFTSG_WGT_SHIFT2 - wkr));
3182         }
3183         if (i0 == 4) {
3184             break;
3185         }
3186     }
3187     xr = a[2] - a[n - 2];
3188     xi = a[3] + a[n - 1];
3189     yr = Mult_W_W(wdr, xr) + Mult_W_W(wdi, xi);
3190     yi = Mult_W_W(wdr, xi) - Mult_W_W(wdi, xr);
3191     a[2] -= yr;
3192     a[3] -= yi;
3193     a[n - 2] += yr;
3194     a[n - 1] -= yi;
3195 }
3196 
3197 
dctsub(picoos_int32 n,PICOFFTSG_FFTTYPE * a)3198 void dctsub(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
3199 {
3200     picoos_int32 i, i0, j, k, m;
3201     PICOFFTSG_FFTTYPE w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
3202     wkr = (PICOFFTSG_FFTTYPE)(0.5*PICODSP_WGT_SHIFT);
3203     wki = (PICOFFTSG_FFTTYPE)(0.5*PICODSP_WGT_SHIFT);
3204 
3205     switch (n) {
3206         case 8 :  wdi=(PICOFFTSG_FFTTYPE)(0.587937772274*PICODSP_WGT_SHIFT);
3207             wdr=(PICOFFTSG_FFTTYPE)(0.392847478390*PICODSP_WGT_SHIFT); w1r=(PICOFFTSG_FFTTYPE)(0.980785250664*PICODSP_WGT_SHIFT);
3208             w1i=(PICOFFTSG_FFTTYPE)(0.195090323687*PICODSP_WGT_SHIFT); ss =(PICOFFTSG_FFTTYPE)(0.390180647373*PICODSP_WGT_SHIFT); break;
3209         case 16 : wdi=(PICOFFTSG_FFTTYPE)(0.546600937843*PICODSP_WGT_SHIFT);
3210             wdr=(PICOFFTSG_FFTTYPE)(0.448583781719*PICODSP_WGT_SHIFT); w1r=(PICOFFTSG_FFTTYPE)(0.995184719563*PICODSP_WGT_SHIFT);
3211             w1i=(PICOFFTSG_FFTTYPE)(0.098017141223*PICODSP_WGT_SHIFT); ss =(PICOFFTSG_FFTTYPE)(0.196034282446*PICODSP_WGT_SHIFT); break;
3212         case 32 : wdi=(PICOFFTSG_FFTTYPE)(0.523931562901*PICODSP_WGT_SHIFT);
3213             wdr=(PICOFFTSG_FFTTYPE)(0.474863886833*PICODSP_WGT_SHIFT); w1r=(PICOFFTSG_FFTTYPE)(0.998795449734*PICODSP_WGT_SHIFT);
3214             w1i=(PICOFFTSG_FFTTYPE)(0.049067676067*PICODSP_WGT_SHIFT); ss =(PICOFFTSG_FFTTYPE)(0.098135352135*PICODSP_WGT_SHIFT); break;
3215         case 64 : wdi=(PICOFFTSG_FFTTYPE)(0.512120008469*PICODSP_WGT_SHIFT);
3216             wdr=(PICOFFTSG_FFTTYPE)(0.487578809261*PICODSP_WGT_SHIFT); w1r=(PICOFFTSG_FFTTYPE)(0.999698817730*PICODSP_WGT_SHIFT);
3217             w1i=(PICOFFTSG_FFTTYPE)(0.024541229010*PICODSP_WGT_SHIFT); ss =(PICOFFTSG_FFTTYPE)(0.049082458019*PICODSP_WGT_SHIFT); break;
3218         default:
3219             wdr = 0; wdi = 0; ss = 0; break;
3220     }
3221 
3222     m = n >> 1;
3223     i = 0;
3224     for (;;) {
3225         i0 = i + DCST_LOOP_DIV2;
3226         if (i0 > m - 2) {
3227             i0 = m - 2;
3228         }
3229         for (j = i + 2; j <= i0; j += 2) {
3230             k = n - j;
3231             xr = picofftsg_mult_w_a(wdi, a[j - 1]) - picofftsg_mult_w_a(wdr, a[k + 1]);
3232             xi = picofftsg_mult_w_a(wdr, a[j - 1]) + picofftsg_mult_w_a(wdi, a[k + 1]);
3233             wkr -= Mult_W_W(ss, wdi);
3234             wki += Mult_W_W(ss, wdr);
3235             yr = Mult_W_W(wki, a[j]) - Mult_W_W(wkr, a[k]);
3236             yi = Mult_W_W(wkr, a[j]) + Mult_W_W(wki, a[k]);
3237             wdr -= Mult_W_W(ss, wki);
3238             wdi += Mult_W_W(ss, wkr);
3239             a[k + 1] = xr;
3240             a[k] = yr;
3241             a[j - 1] = xi;
3242             a[j] = yi;
3243         }
3244         if (i0 == m - 2) {
3245             break;
3246         }
3247     }
3248     xr = picofftsg_mult_w_a(wdi, a[m - 1]) - picofftsg_mult_w_a(wdr, a[m + 1]);
3249     a[m - 1] = picofftsg_mult_w_a(wdr, a[m - 1]) + picofftsg_mult_w_a(wdi, a[m + 1]);
3250     a[m + 1] = xr;
3251     a[m] = Mult_W_W(WR5000, a[m]);
3252 }
3253 
3254 
dctsub4(picoos_int32 n,PICOFFTSG_FFTTYPE * a)3255 void dctsub4(picoos_int32 n, PICOFFTSG_FFTTYPE *a)
3256 {
3257     picoos_int32 m;
3258     PICOFFTSG_FFTTYPE wki, wdr, wdi, xr;
3259 
3260     wki = WR5000;
3261     m = n >> 1;
3262     if (m == 2) {
3263         wdr = Mult_W_W(wki, WI2500);
3264         wdi = Mult_W_W(wki, WR2500);
3265         xr = Mult_W_W(wdi, a[1]) - Mult_W_W(wdr, a[3]);
3266         a[1] = Mult_W_W(wdr, a[1]) + Mult_W_W(wdi, a[3]);
3267         a[3] = xr;
3268     }
3269     a[m] = Mult_W_W(wki, a[m]);
3270 }
3271 
3272 #ifdef __cplusplus
3273 }
3274 #endif
3275