1 /*
2  *  Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3  *
4  *  Use of this source code is governed by a BSD-style license
5  *  that can be found in the LICENSE file in the root of the source
6  *  tree. An additional intellectual property rights grant can be found
7  *  in the file PATENTS.  All contributing project authors may
8  *  be found in the AUTHORS file in the root of the source tree.
9  */
10 
11 
12 /*
13  * This file contains the function WebRtcSpl_ComplexFFT().
14  * The description header can be found in signal_processing_library.h
15  *
16  */
17 
18 #include "signal_processing_library.h"
19 
20 #define CFFTSFT 14
21 #define CFFTRND 1
22 #define CFFTRND2 16384
23 
24 #define CIFFTSFT 14
25 #define CIFFTRND 1
26 
27 static const WebRtc_Word16 kSinTable1024[] = {
28       0,    201,    402,    603,    804,   1005,   1206,   1406,
29    1607,   1808,   2009,   2209,   2410,   2610,   2811,   3011,
30    3211,   3411,   3611,   3811,   4011,   4210,   4409,   4608,
31    4807,   5006,   5205,   5403,   5601,   5799,   5997,   6195,
32    6392,   6589,   6786,   6982,   7179,   7375,   7571,   7766,
33    7961,   8156,   8351,   8545,   8739,   8932,   9126,   9319,
34    9511,   9703,   9895,  10087,  10278,  10469,  10659,  10849,
35   11038,  11227,  11416,  11604,  11792,  11980,  12166,  12353,
36   12539,  12724,  12909,  13094,  13278,  13462,  13645,  13827,
37   14009,  14191,  14372,  14552,  14732,  14911,  15090,  15268,
38   15446,  15623,  15799,  15975,  16150,  16325,  16499,  16672,
39   16845,  17017,  17189,  17360,  17530,  17699,  17868,  18036,
40   18204,  18371,  18537,  18702,  18867,  19031,  19194,  19357,
41   19519,  19680,  19840,  20000,  20159,  20317,  20474,  20631,
42   20787,  20942,  21096,  21249,  21402,  21554,  21705,  21855,
43   22004,  22153,  22301,  22448,  22594,  22739,  22883,  23027,
44   23169,  23311,  23452,  23592,  23731,  23869,  24006,  24143,
45   24278,  24413,  24546,  24679,  24811,  24942,  25072,  25201,
46   25329,  25456,  25582,  25707,  25831,  25954,  26077,  26198,
47   26318,  26437,  26556,  26673,  26789,  26905,  27019,  27132,
48   27244,  27355,  27466,  27575,  27683,  27790,  27896,  28001,
49   28105,  28208,  28309,  28410,  28510,  28608,  28706,  28802,
50   28897,  28992,  29085,  29177,  29268,  29358,  29446,  29534,
51   29621,  29706,  29790,  29873,  29955,  30036,  30116,  30195,
52   30272,  30349,  30424,  30498,  30571,  30643,  30713,  30783,
53   30851,  30918,  30984,  31049,
54   31113,  31175,  31236,  31297,
55   31356,  31413,  31470,  31525,  31580,  31633,  31684,  31735,
56   31785,  31833,  31880,  31926,  31970,  32014,  32056,  32097,
57   32137,  32176,  32213,  32249,  32284,  32318,  32350,  32382,
58   32412,  32441,  32468,  32495,  32520,  32544,  32567,  32588,
59   32609,  32628,  32646,  32662,  32678,  32692,  32705,  32717,
60   32727,  32736,  32744,  32751,  32757,  32761,  32764,  32766,
61   32767,  32766,  32764,  32761,  32757,  32751,  32744,  32736,
62   32727,  32717,  32705,  32692,  32678,  32662,  32646,  32628,
63   32609,  32588,  32567,  32544,  32520,  32495,  32468,  32441,
64   32412,  32382,  32350,  32318,  32284,  32249,  32213,  32176,
65   32137,  32097,  32056,  32014,  31970,  31926,  31880,  31833,
66   31785,  31735,  31684,  31633,  31580,  31525,  31470,  31413,
67   31356,  31297,  31236,  31175,  31113,  31049,  30984,  30918,
68   30851,  30783,  30713,  30643,  30571,  30498,  30424,  30349,
69   30272,  30195,  30116,  30036,  29955,  29873,  29790,  29706,
70   29621,  29534,  29446,  29358,  29268,  29177,  29085,  28992,
71   28897,  28802,  28706,  28608,  28510,  28410,  28309,  28208,
72   28105,  28001,  27896,  27790,  27683,  27575,  27466,  27355,
73   27244,  27132,  27019,  26905,  26789,  26673,  26556,  26437,
74   26318,  26198,  26077,  25954,  25831,  25707,  25582,  25456,
75   25329,  25201,  25072,  24942,  24811,  24679,  24546,  24413,
76   24278,  24143,  24006,  23869,  23731,  23592,  23452,  23311,
77   23169,  23027,  22883,  22739,  22594,  22448,  22301,  22153,
78   22004,  21855,  21705,  21554,  21402,  21249,  21096,  20942,
79   20787,  20631,  20474,  20317,  20159,  20000,  19840,  19680,
80   19519,  19357,  19194,  19031,  18867,  18702,  18537,  18371,
81   18204,  18036,  17868,  17699,  17530,  17360,  17189,  17017,
82   16845,  16672,  16499,  16325,  16150,  15975,  15799,  15623,
83   15446,  15268,  15090,  14911,  14732,  14552,  14372,  14191,
84   14009,  13827,  13645,  13462,  13278,  13094,  12909,  12724,
85   12539,  12353,  12166,  11980,  11792,  11604,  11416,  11227,
86   11038,  10849,  10659,  10469,  10278,  10087,   9895,   9703,
87    9511,   9319,   9126,   8932,   8739,   8545,   8351,   8156,
88    7961,   7766,   7571,   7375,   7179,   6982,   6786,   6589,
89    6392,   6195,   5997,   5799,   5601,   5403,   5205,   5006,
90    4807,   4608,   4409,   4210,   4011,   3811,   3611,   3411,
91    3211,   3011,   2811,   2610,   2410,   2209,   2009,   1808,
92    1607,   1406,   1206,   1005,    804,    603,    402,    201,
93       0,   -201,   -402,   -603,   -804,  -1005,  -1206,  -1406,
94   -1607,  -1808,  -2009,  -2209,  -2410,  -2610,  -2811,  -3011,
95   -3211,  -3411,  -3611,  -3811,  -4011,  -4210,  -4409,  -4608,
96   -4807,  -5006,  -5205,  -5403,  -5601,  -5799,  -5997,  -6195,
97   -6392,  -6589,  -6786,  -6982,  -7179,  -7375,  -7571,  -7766,
98   -7961,  -8156,  -8351,  -8545,  -8739,  -8932,  -9126,  -9319,
99   -9511,  -9703,  -9895, -10087, -10278, -10469, -10659, -10849,
100  -11038, -11227, -11416, -11604, -11792, -11980, -12166, -12353,
101  -12539, -12724, -12909, -13094, -13278, -13462, -13645, -13827,
102  -14009, -14191, -14372, -14552, -14732, -14911, -15090, -15268,
103  -15446, -15623, -15799, -15975, -16150, -16325, -16499, -16672,
104  -16845, -17017, -17189, -17360, -17530, -17699, -17868, -18036,
105  -18204, -18371, -18537, -18702, -18867, -19031, -19194, -19357,
106  -19519, -19680, -19840, -20000, -20159, -20317, -20474, -20631,
107  -20787, -20942, -21096, -21249, -21402, -21554, -21705, -21855,
108  -22004, -22153, -22301, -22448, -22594, -22739, -22883, -23027,
109  -23169, -23311, -23452, -23592, -23731, -23869, -24006, -24143,
110  -24278, -24413, -24546, -24679, -24811, -24942, -25072, -25201,
111  -25329, -25456, -25582, -25707, -25831, -25954, -26077, -26198,
112  -26318, -26437, -26556, -26673, -26789, -26905, -27019, -27132,
113  -27244, -27355, -27466, -27575, -27683, -27790, -27896, -28001,
114  -28105, -28208, -28309, -28410, -28510, -28608, -28706, -28802,
115  -28897, -28992, -29085, -29177, -29268, -29358, -29446, -29534,
116  -29621, -29706, -29790, -29873, -29955, -30036, -30116, -30195,
117  -30272, -30349, -30424, -30498, -30571, -30643, -30713, -30783,
118  -30851, -30918, -30984, -31049, -31113, -31175, -31236, -31297,
119  -31356, -31413, -31470, -31525, -31580, -31633, -31684, -31735,
120  -31785, -31833, -31880, -31926, -31970, -32014, -32056, -32097,
121  -32137, -32176, -32213, -32249, -32284, -32318, -32350, -32382,
122  -32412, -32441, -32468, -32495, -32520, -32544, -32567, -32588,
123  -32609, -32628, -32646, -32662, -32678, -32692, -32705, -32717,
124  -32727, -32736, -32744, -32751, -32757, -32761, -32764, -32766,
125  -32767, -32766, -32764, -32761, -32757, -32751, -32744, -32736,
126  -32727, -32717, -32705, -32692, -32678, -32662, -32646, -32628,
127  -32609, -32588, -32567, -32544, -32520, -32495, -32468, -32441,
128  -32412, -32382, -32350, -32318, -32284, -32249, -32213, -32176,
129  -32137, -32097, -32056, -32014, -31970, -31926, -31880, -31833,
130  -31785, -31735, -31684, -31633, -31580, -31525, -31470, -31413,
131  -31356, -31297, -31236, -31175, -31113, -31049, -30984, -30918,
132  -30851, -30783, -30713, -30643, -30571, -30498, -30424, -30349,
133  -30272, -30195, -30116, -30036, -29955, -29873, -29790, -29706,
134  -29621, -29534, -29446, -29358, -29268, -29177, -29085, -28992,
135  -28897, -28802, -28706, -28608, -28510, -28410, -28309, -28208,
136  -28105, -28001, -27896, -27790, -27683, -27575, -27466, -27355,
137  -27244, -27132, -27019, -26905, -26789, -26673, -26556, -26437,
138  -26318, -26198, -26077, -25954, -25831, -25707, -25582, -25456,
139  -25329, -25201, -25072, -24942, -24811, -24679, -24546, -24413,
140  -24278, -24143, -24006, -23869, -23731, -23592, -23452, -23311,
141  -23169, -23027, -22883, -22739, -22594, -22448, -22301, -22153,
142  -22004, -21855, -21705, -21554, -21402, -21249, -21096, -20942,
143  -20787, -20631, -20474, -20317, -20159, -20000, -19840, -19680,
144  -19519, -19357, -19194, -19031, -18867, -18702, -18537, -18371,
145  -18204, -18036, -17868, -17699, -17530, -17360, -17189, -17017,
146  -16845, -16672, -16499, -16325, -16150, -15975, -15799, -15623,
147  -15446, -15268, -15090, -14911, -14732, -14552, -14372, -14191,
148  -14009, -13827, -13645, -13462, -13278, -13094, -12909, -12724,
149  -12539, -12353, -12166, -11980, -11792, -11604, -11416, -11227,
150  -11038, -10849, -10659, -10469, -10278, -10087,  -9895,  -9703,
151   -9511,  -9319,  -9126,  -8932,  -8739,  -8545,  -8351,  -8156,
152   -7961,  -7766,  -7571,  -7375,  -7179,  -6982,  -6786,  -6589,
153   -6392,  -6195,  -5997,  -5799,  -5601,  -5403,  -5205,  -5006,
154   -4807,  -4608,  -4409,  -4210,  -4011,  -3811,  -3611,  -3411,
155   -3211,  -3011,  -2811,  -2610,  -2410,  -2209,  -2009,  -1808,
156   -1607,  -1406,  -1206,  -1005,   -804,   -603,   -402,   -201
157 };
158 
WebRtcSpl_ComplexFFT(WebRtc_Word16 frfi[],int stages,int mode)159 int WebRtcSpl_ComplexFFT(WebRtc_Word16 frfi[], int stages, int mode)
160 {
161     int i, j, l, k, istep, n, m;
162     WebRtc_Word16 wr, wi;
163     WebRtc_Word32 tr32, ti32, qr32, qi32;
164 
165     /* The 1024-value is a constant given from the size of kSinTable1024[],
166      * and should not be changed depending on the input parameter 'stages'
167      */
168     n = 1 << stages;
169     if (n > 1024)
170         return -1;
171 
172     l = 1;
173     k = 10 - 1; /* Constant for given kSinTable1024[]. Do not change
174          depending on the input parameter 'stages' */
175 
176     if (mode == 0)
177     {
178         // mode==0: Low-complexity and Low-accuracy mode
179         while (l < n)
180         {
181             istep = l << 1;
182 
183             for (m = 0; m < l; ++m)
184             {
185                 j = m << k;
186 
187                 /* The 256-value is a constant given as 1/4 of the size of
188                  * kSinTable1024[], and should not be changed depending on the input
189                  * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
190                  */
191                 wr = kSinTable1024[j + 256];
192                 wi = -kSinTable1024[j];
193 
194                 for (i = m; i < n; i += istep)
195                 {
196                     j = i + l;
197 
198                     tr32 = WEBRTC_SPL_RSHIFT_W32((WEBRTC_SPL_MUL_16_16(wr, frfi[2 * j])
199                             - WEBRTC_SPL_MUL_16_16(wi, frfi[2 * j + 1])), 15);
200 
201                     ti32 = WEBRTC_SPL_RSHIFT_W32((WEBRTC_SPL_MUL_16_16(wr, frfi[2 * j + 1])
202                             + WEBRTC_SPL_MUL_16_16(wi, frfi[2 * j])), 15);
203 
204                     qr32 = (WebRtc_Word32)frfi[2 * i];
205                     qi32 = (WebRtc_Word32)frfi[2 * i + 1];
206                     frfi[2 * j] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(qr32 - tr32, 1);
207                     frfi[2 * j + 1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(qi32 - ti32, 1);
208                     frfi[2 * i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(qr32 + tr32, 1);
209                     frfi[2 * i + 1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(qi32 + ti32, 1);
210                 }
211             }
212 
213             --k;
214             l = istep;
215 
216         }
217 
218     } else
219     {
220         // mode==1: High-complexity and High-accuracy mode
221         while (l < n)
222         {
223             istep = l << 1;
224 
225             for (m = 0; m < l; ++m)
226             {
227                 j = m << k;
228 
229                 /* The 256-value is a constant given as 1/4 of the size of
230                  * kSinTable1024[], and should not be changed depending on the input
231                  * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
232                  */
233                 wr = kSinTable1024[j + 256];
234                 wi = -kSinTable1024[j];
235 
236 #ifdef WEBRTC_ARCH_ARM_V7A
237                 WebRtc_Word32 wri;
238                 WebRtc_Word32 frfi_r;
239                 __asm__("pkhbt %0, %1, %2, lsl #16" : "=r"(wri) :
240                     "r"((WebRtc_Word32)wr), "r"((WebRtc_Word32)wi));
241 #endif
242 
243                 for (i = m; i < n; i += istep)
244                 {
245                     j = i + l;
246 
247 #ifdef WEBRTC_ARCH_ARM_V7A
248                     __asm__("pkhbt %0, %1, %2, lsl #16" : "=r"(frfi_r) :
249                         "r"((WebRtc_Word32)frfi[2*j]), "r"((WebRtc_Word32)frfi[2*j +1]));
250                     __asm__("smlsd %0, %1, %2, %3" : "=r"(tr32) :
251                         "r"(wri), "r"(frfi_r), "r"(CFFTRND));
252                     __asm__("smladx %0, %1, %2, %3" : "=r"(ti32) :
253                         "r"(wri), "r"(frfi_r), "r"(CFFTRND));
254 
255 #else
256                     tr32 = WEBRTC_SPL_MUL_16_16(wr, frfi[2 * j])
257                             - WEBRTC_SPL_MUL_16_16(wi, frfi[2 * j + 1]) + CFFTRND;
258 
259                     ti32 = WEBRTC_SPL_MUL_16_16(wr, frfi[2 * j + 1])
260                             + WEBRTC_SPL_MUL_16_16(wi, frfi[2 * j]) + CFFTRND;
261 #endif
262 
263                     tr32 = WEBRTC_SPL_RSHIFT_W32(tr32, 15 - CFFTSFT);
264                     ti32 = WEBRTC_SPL_RSHIFT_W32(ti32, 15 - CFFTSFT);
265 
266                     qr32 = ((WebRtc_Word32)frfi[2 * i]) << CFFTSFT;
267                     qi32 = ((WebRtc_Word32)frfi[2 * i + 1]) << CFFTSFT;
268 
269                     frfi[2 * j] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(
270                             (qr32 - tr32 + CFFTRND2), 1 + CFFTSFT);
271                     frfi[2 * j + 1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(
272                             (qi32 - ti32 + CFFTRND2), 1 + CFFTSFT);
273                     frfi[2 * i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(
274                             (qr32 + tr32 + CFFTRND2), 1 + CFFTSFT);
275                     frfi[2 * i + 1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(
276                             (qi32 + ti32 + CFFTRND2), 1 + CFFTSFT);
277                 }
278             }
279 
280             --k;
281             l = istep;
282         }
283     }
284     return 0;
285 }
286 
WebRtcSpl_ComplexIFFT(WebRtc_Word16 frfi[],int stages,int mode)287 int WebRtcSpl_ComplexIFFT(WebRtc_Word16 frfi[], int stages, int mode)
288 {
289     int i, j, l, k, istep, n, m, scale, shift;
290     WebRtc_Word16 wr, wi;
291     WebRtc_Word32 tr32, ti32, qr32, qi32;
292     WebRtc_Word32 tmp32, round2;
293 
294     /* The 1024-value is a constant given from the size of kSinTable1024[],
295      * and should not be changed depending on the input parameter 'stages'
296      */
297     n = 1 << stages;
298     if (n > 1024)
299         return -1;
300 
301     scale = 0;
302 
303     l = 1;
304     k = 10 - 1; /* Constant for given kSinTable1024[]. Do not change
305          depending on the input parameter 'stages' */
306 
307     while (l < n)
308     {
309         // variable scaling, depending upon data
310         shift = 0;
311         round2 = 8192;
312 
313         tmp32 = (WebRtc_Word32)WebRtcSpl_MaxAbsValueW16(frfi, 2 * n);
314         if (tmp32 > 13573)
315         {
316             shift++;
317             scale++;
318             round2 <<= 1;
319         }
320         if (tmp32 > 27146)
321         {
322             shift++;
323             scale++;
324             round2 <<= 1;
325         }
326 
327         istep = l << 1;
328 
329         if (mode == 0)
330         {
331             // mode==0: Low-complexity and Low-accuracy mode
332             for (m = 0; m < l; ++m)
333             {
334                 j = m << k;
335 
336                 /* The 256-value is a constant given as 1/4 of the size of
337                  * kSinTable1024[], and should not be changed depending on the input
338                  * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
339                  */
340                 wr = kSinTable1024[j + 256];
341                 wi = kSinTable1024[j];
342 
343                 for (i = m; i < n; i += istep)
344                 {
345                     j = i + l;
346 
347                     tr32 = WEBRTC_SPL_RSHIFT_W32((WEBRTC_SPL_MUL_16_16_RSFT(wr, frfi[2 * j], 0)
348                             - WEBRTC_SPL_MUL_16_16_RSFT(wi, frfi[2 * j + 1], 0)), 15);
349 
350                     ti32 = WEBRTC_SPL_RSHIFT_W32(
351                             (WEBRTC_SPL_MUL_16_16_RSFT(wr, frfi[2 * j + 1], 0)
352                                     + WEBRTC_SPL_MUL_16_16_RSFT(wi,frfi[2*j],0)), 15);
353 
354                     qr32 = (WebRtc_Word32)frfi[2 * i];
355                     qi32 = (WebRtc_Word32)frfi[2 * i + 1];
356                     frfi[2 * j] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(qr32 - tr32, shift);
357                     frfi[2 * j + 1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(qi32 - ti32, shift);
358                     frfi[2 * i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(qr32 + tr32, shift);
359                     frfi[2 * i + 1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(qi32 + ti32, shift);
360                 }
361             }
362         } else
363         {
364             // mode==1: High-complexity and High-accuracy mode
365 
366             for (m = 0; m < l; ++m)
367             {
368                 j = m << k;
369 
370                 /* The 256-value is a constant given as 1/4 of the size of
371                  * kSinTable1024[], and should not be changed depending on the input
372                  * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
373                  */
374                 wr = kSinTable1024[j + 256];
375                 wi = kSinTable1024[j];
376 
377 #ifdef WEBRTC_ARCH_ARM_V7A
378                 WebRtc_Word32 wri;
379                 WebRtc_Word32 frfi_r;
380                 __asm__("pkhbt %0, %1, %2, lsl #16" : "=r"(wri) :
381                     "r"((WebRtc_Word32)wr), "r"((WebRtc_Word32)wi));
382 #endif
383 
384                 for (i = m; i < n; i += istep)
385                 {
386                     j = i + l;
387 
388 #ifdef WEBRTC_ARCH_ARM_V7A
389                     __asm__("pkhbt %0, %1, %2, lsl #16" : "=r"(frfi_r) :
390                         "r"((WebRtc_Word32)frfi[2*j]), "r"((WebRtc_Word32)frfi[2*j +1]));
391                     __asm__("smlsd %0, %1, %2, %3" : "=r"(tr32) :
392                         "r"(wri), "r"(frfi_r), "r"(CIFFTRND));
393                     __asm__("smladx %0, %1, %2, %3" : "=r"(ti32) :
394                         "r"(wri), "r"(frfi_r), "r"(CIFFTRND));
395 #else
396 
397                     tr32 = WEBRTC_SPL_MUL_16_16(wr, frfi[2 * j])
398                             - WEBRTC_SPL_MUL_16_16(wi, frfi[2 * j + 1]) + CIFFTRND;
399 
400                     ti32 = WEBRTC_SPL_MUL_16_16(wr, frfi[2 * j + 1])
401                             + WEBRTC_SPL_MUL_16_16(wi, frfi[2 * j]) + CIFFTRND;
402 #endif
403                     tr32 = WEBRTC_SPL_RSHIFT_W32(tr32, 15 - CIFFTSFT);
404                     ti32 = WEBRTC_SPL_RSHIFT_W32(ti32, 15 - CIFFTSFT);
405 
406                     qr32 = ((WebRtc_Word32)frfi[2 * i]) << CIFFTSFT;
407                     qi32 = ((WebRtc_Word32)frfi[2 * i + 1]) << CIFFTSFT;
408 
409                     frfi[2 * j] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((qr32 - tr32+round2),
410                                                                        shift+CIFFTSFT);
411                     frfi[2 * j + 1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(
412                             (qi32 - ti32 + round2), shift + CIFFTSFT);
413                     frfi[2 * i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((qr32 + tr32 + round2),
414                                                                        shift + CIFFTSFT);
415                     frfi[2 * i + 1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(
416                             (qi32 + ti32 + round2), shift + CIFFTSFT);
417                 }
418             }
419 
420         }
421         --k;
422         l = istep;
423     }
424     return scale;
425 }
426