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 "webrtc/common_audio/signal_processing/complex_fft_tables.h"
19#include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
20
21#define CFFTSFT 14
22#define CFFTRND 1
23#define CFFTRND2 16384
24
25#define CIFFTSFT 14
26#define CIFFTRND 1
27
28
29int WebRtcSpl_ComplexFFT(int16_t frfi[], int stages, int mode)
30{
31    int i, j, l, k, istep, n, m;
32    int16_t wr, wi;
33    int32_t tr32, ti32, qr32, qi32;
34
35    /* The 1024-value is a constant given from the size of kSinTable1024[],
36     * and should not be changed depending on the input parameter 'stages'
37     */
38    n = 1 << stages;
39    if (n > 1024)
40        return -1;
41
42    l = 1;
43    k = 10 - 1; /* Constant for given kSinTable1024[]. Do not change
44         depending on the input parameter 'stages' */
45
46    if (mode == 0)
47    {
48        // mode==0: Low-complexity and Low-accuracy mode
49        while (l < n)
50        {
51            istep = l << 1;
52
53            for (m = 0; m < l; ++m)
54            {
55                j = m << k;
56
57                /* The 256-value is a constant given as 1/4 of the size of
58                 * kSinTable1024[], and should not be changed depending on the input
59                 * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
60                 */
61                wr = kSinTable1024[j + 256];
62                wi = -kSinTable1024[j];
63
64                for (i = m; i < n; i += istep)
65                {
66                    j = i + l;
67
68                    tr32 = (wr * frfi[2 * j] - wi * frfi[2 * j + 1]) >> 15;
69
70                    ti32 = (wr * frfi[2 * j + 1] + wi * frfi[2 * j]) >> 15;
71
72                    qr32 = (int32_t)frfi[2 * i];
73                    qi32 = (int32_t)frfi[2 * i + 1];
74                    frfi[2 * j] = (int16_t)((qr32 - tr32) >> 1);
75                    frfi[2 * j + 1] = (int16_t)((qi32 - ti32) >> 1);
76                    frfi[2 * i] = (int16_t)((qr32 + tr32) >> 1);
77                    frfi[2 * i + 1] = (int16_t)((qi32 + ti32) >> 1);
78                }
79            }
80
81            --k;
82            l = istep;
83
84        }
85
86    } else
87    {
88        // mode==1: High-complexity and High-accuracy mode
89        while (l < n)
90        {
91            istep = l << 1;
92
93            for (m = 0; m < l; ++m)
94            {
95                j = m << k;
96
97                /* The 256-value is a constant given as 1/4 of the size of
98                 * kSinTable1024[], and should not be changed depending on the input
99                 * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
100                 */
101                wr = kSinTable1024[j + 256];
102                wi = -kSinTable1024[j];
103
104#ifdef WEBRTC_ARCH_ARM_V7
105                int32_t wri = 0;
106                __asm __volatile("pkhbt %0, %1, %2, lsl #16" : "=r"(wri) :
107                    "r"((int32_t)wr), "r"((int32_t)wi));
108#endif
109
110                for (i = m; i < n; i += istep)
111                {
112                    j = i + l;
113
114#ifdef WEBRTC_ARCH_ARM_V7
115                    register int32_t frfi_r;
116                    __asm __volatile(
117                        "pkhbt %[frfi_r], %[frfi_even], %[frfi_odd],"
118                        " lsl #16\n\t"
119                        "smlsd %[tr32], %[wri], %[frfi_r], %[cfftrnd]\n\t"
120                        "smladx %[ti32], %[wri], %[frfi_r], %[cfftrnd]\n\t"
121                        :[frfi_r]"=&r"(frfi_r),
122                         [tr32]"=&r"(tr32),
123                         [ti32]"=r"(ti32)
124                        :[frfi_even]"r"((int32_t)frfi[2*j]),
125                         [frfi_odd]"r"((int32_t)frfi[2*j +1]),
126                         [wri]"r"(wri),
127                         [cfftrnd]"r"(CFFTRND));
128#else
129                    tr32 = wr * frfi[2 * j] - wi * frfi[2 * j + 1] + CFFTRND;
130
131                    ti32 = wr * frfi[2 * j + 1] + wi * frfi[2 * j] + CFFTRND;
132#endif
133
134                    tr32 >>= 15 - CFFTSFT;
135                    ti32 >>= 15 - CFFTSFT;
136
137                    qr32 = ((int32_t)frfi[2 * i]) << CFFTSFT;
138                    qi32 = ((int32_t)frfi[2 * i + 1]) << CFFTSFT;
139
140                    frfi[2 * j] = (int16_t)(
141                        (qr32 - tr32 + CFFTRND2) >> (1 + CFFTSFT));
142                    frfi[2 * j + 1] = (int16_t)(
143                        (qi32 - ti32 + CFFTRND2) >> (1 + CFFTSFT));
144                    frfi[2 * i] = (int16_t)(
145                        (qr32 + tr32 + CFFTRND2) >> (1 + CFFTSFT));
146                    frfi[2 * i + 1] = (int16_t)(
147                        (qi32 + ti32 + CFFTRND2) >> (1 + CFFTSFT));
148                }
149            }
150
151            --k;
152            l = istep;
153        }
154    }
155    return 0;
156}
157
158int WebRtcSpl_ComplexIFFT(int16_t frfi[], int stages, int mode)
159{
160    size_t i, j, l, istep, n, m;
161    int k, scale, shift;
162    int16_t wr, wi;
163    int32_t tr32, ti32, qr32, qi32;
164    int32_t tmp32, round2;
165
166    /* The 1024-value is a constant given from the size of kSinTable1024[],
167     * and should not be changed depending on the input parameter 'stages'
168     */
169    n = 1 << stages;
170    if (n > 1024)
171        return -1;
172
173    scale = 0;
174
175    l = 1;
176    k = 10 - 1; /* Constant for given kSinTable1024[]. Do not change
177         depending on the input parameter 'stages' */
178
179    while (l < n)
180    {
181        // variable scaling, depending upon data
182        shift = 0;
183        round2 = 8192;
184
185        tmp32 = WebRtcSpl_MaxAbsValueW16(frfi, 2 * n);
186        if (tmp32 > 13573)
187        {
188            shift++;
189            scale++;
190            round2 <<= 1;
191        }
192        if (tmp32 > 27146)
193        {
194            shift++;
195            scale++;
196            round2 <<= 1;
197        }
198
199        istep = l << 1;
200
201        if (mode == 0)
202        {
203            // mode==0: Low-complexity and Low-accuracy mode
204            for (m = 0; m < l; ++m)
205            {
206                j = m << k;
207
208                /* The 256-value is a constant given as 1/4 of the size of
209                 * kSinTable1024[], and should not be changed depending on the input
210                 * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
211                 */
212                wr = kSinTable1024[j + 256];
213                wi = kSinTable1024[j];
214
215                for (i = m; i < n; i += istep)
216                {
217                    j = i + l;
218
219                    tr32 = (wr * frfi[2 * j] - wi * frfi[2 * j + 1]) >> 15;
220
221                    ti32 = (wr * frfi[2 * j + 1] + wi * frfi[2 * j]) >> 15;
222
223                    qr32 = (int32_t)frfi[2 * i];
224                    qi32 = (int32_t)frfi[2 * i + 1];
225                    frfi[2 * j] = (int16_t)((qr32 - tr32) >> shift);
226                    frfi[2 * j + 1] = (int16_t)((qi32 - ti32) >> shift);
227                    frfi[2 * i] = (int16_t)((qr32 + tr32) >> shift);
228                    frfi[2 * i + 1] = (int16_t)((qi32 + ti32) >> shift);
229                }
230            }
231        } else
232        {
233            // mode==1: High-complexity and High-accuracy mode
234
235            for (m = 0; m < l; ++m)
236            {
237                j = m << k;
238
239                /* The 256-value is a constant given as 1/4 of the size of
240                 * kSinTable1024[], and should not be changed depending on the input
241                 * parameter 'stages'. It will result in 0 <= j < N_SINE_WAVE/2
242                 */
243                wr = kSinTable1024[j + 256];
244                wi = kSinTable1024[j];
245
246#ifdef WEBRTC_ARCH_ARM_V7
247                int32_t wri = 0;
248                __asm __volatile("pkhbt %0, %1, %2, lsl #16" : "=r"(wri) :
249                    "r"((int32_t)wr), "r"((int32_t)wi));
250#endif
251
252                for (i = m; i < n; i += istep)
253                {
254                    j = i + l;
255
256#ifdef WEBRTC_ARCH_ARM_V7
257                    register int32_t frfi_r;
258                    __asm __volatile(
259                      "pkhbt %[frfi_r], %[frfi_even], %[frfi_odd], lsl #16\n\t"
260                      "smlsd %[tr32], %[wri], %[frfi_r], %[cifftrnd]\n\t"
261                      "smladx %[ti32], %[wri], %[frfi_r], %[cifftrnd]\n\t"
262                      :[frfi_r]"=&r"(frfi_r),
263                       [tr32]"=&r"(tr32),
264                       [ti32]"=r"(ti32)
265                      :[frfi_even]"r"((int32_t)frfi[2*j]),
266                       [frfi_odd]"r"((int32_t)frfi[2*j +1]),
267                       [wri]"r"(wri),
268                       [cifftrnd]"r"(CIFFTRND)
269                    );
270#else
271
272                    tr32 = wr * frfi[2 * j] - wi * frfi[2 * j + 1] + CIFFTRND;
273
274                    ti32 = wr * frfi[2 * j + 1] + wi * frfi[2 * j] + CIFFTRND;
275#endif
276                    tr32 >>= 15 - CIFFTSFT;
277                    ti32 >>= 15 - CIFFTSFT;
278
279                    qr32 = ((int32_t)frfi[2 * i]) << CIFFTSFT;
280                    qi32 = ((int32_t)frfi[2 * i + 1]) << CIFFTSFT;
281
282                    frfi[2 * j] = (int16_t)(
283                        (qr32 - tr32 + round2) >> (shift + CIFFTSFT));
284                    frfi[2 * j + 1] = (int16_t)(
285                        (qi32 - ti32 + round2) >> (shift + CIFFTSFT));
286                    frfi[2 * i] = (int16_t)(
287                        (qr32 + tr32 + round2) >> (shift + CIFFTSFT));
288                    frfi[2 * i + 1] = (int16_t)(
289                        (qi32 + ti32 + round2) >> (shift + CIFFTSFT));
290                }
291            }
292
293        }
294        --k;
295        l = istep;
296    }
297    return scale;
298}
299