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