arm_cfft_q15.c 8.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345
  1. /* ----------------------------------------------------------------------
  2. * Project: CMSIS DSP Library
  3. * Title: arm_cfft_q15.c
  4. * Description: Combined Radix Decimation in Q15 Frequency CFFT processing function
  5. *
  6. * $Date: 27. January 2017
  7. * $Revision: V.1.5.1
  8. *
  9. * Target Processor: Cortex-M cores
  10. * -------------------------------------------------------------------- */
  11. /*
  12. * Copyright (C) 2010-2017 ARM Limited or its affiliates. All rights reserved.
  13. *
  14. * SPDX-License-Identifier: Apache-2.0
  15. *
  16. * Licensed under the Apache License, Version 2.0 (the License); you may
  17. * not use this file except in compliance with the License.
  18. * You may obtain a copy of the License at
  19. *
  20. * www.apache.org/licenses/LICENSE-2.0
  21. *
  22. * Unless required by applicable law or agreed to in writing, software
  23. * distributed under the License is distributed on an AS IS BASIS, WITHOUT
  24. * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  25. * See the License for the specific language governing permissions and
  26. * limitations under the License.
  27. */
  28. #include "arm_math.h"
  29. extern void arm_radix4_butterfly_q15(
  30. q15_t * pSrc,
  31. uint32_t fftLen,
  32. q15_t * pCoef,
  33. uint32_t twidCoefModifier);
  34. extern void arm_radix4_butterfly_inverse_q15(
  35. q15_t * pSrc,
  36. uint32_t fftLen,
  37. q15_t * pCoef,
  38. uint32_t twidCoefModifier);
  39. extern void arm_bitreversal_16(
  40. uint16_t * pSrc,
  41. const uint16_t bitRevLen,
  42. const uint16_t * pBitRevTable);
  43. void arm_cfft_radix4by2_q15(
  44. q15_t * pSrc,
  45. uint32_t fftLen,
  46. const q15_t * pCoef);
  47. void arm_cfft_radix4by2_inverse_q15(
  48. q15_t * pSrc,
  49. uint32_t fftLen,
  50. const q15_t * pCoef);
  51. /**
  52. * @ingroup groupTransforms
  53. */
  54. /**
  55. * @addtogroup ComplexFFT
  56. * @{
  57. */
  58. /**
  59. * @details
  60. * @brief Processing function for the Q15 complex FFT.
  61. * @param[in] *S points to an instance of the Q15 CFFT structure.
  62. * @param[in, out] *p1 points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.
  63. * @param[in] ifftFlag flag that selects forward (ifftFlag=0) or inverse (ifftFlag=1) transform.
  64. * @param[in] bitReverseFlag flag that enables (bitReverseFlag=1) or disables (bitReverseFlag=0) bit reversal of output.
  65. * @return none.
  66. */
  67. void arm_cfft_q15(
  68. const arm_cfft_instance_q15 * S,
  69. q15_t * p1,
  70. uint8_t ifftFlag,
  71. uint8_t bitReverseFlag)
  72. {
  73. uint32_t L = S->fftLen;
  74. if (ifftFlag == 1U)
  75. {
  76. switch (L)
  77. {
  78. case 16:
  79. case 64:
  80. case 256:
  81. case 1024:
  82. case 4096:
  83. arm_radix4_butterfly_inverse_q15 ( p1, L, (q15_t*)S->pTwiddle, 1 );
  84. break;
  85. case 32:
  86. case 128:
  87. case 512:
  88. case 2048:
  89. arm_cfft_radix4by2_inverse_q15 ( p1, L, S->pTwiddle );
  90. break;
  91. }
  92. }
  93. else
  94. {
  95. switch (L)
  96. {
  97. case 16:
  98. case 64:
  99. case 256:
  100. case 1024:
  101. case 4096:
  102. arm_radix4_butterfly_q15 ( p1, L, (q15_t*)S->pTwiddle, 1 );
  103. break;
  104. case 32:
  105. case 128:
  106. case 512:
  107. case 2048:
  108. arm_cfft_radix4by2_q15 ( p1, L, S->pTwiddle );
  109. break;
  110. }
  111. }
  112. if ( bitReverseFlag )
  113. arm_bitreversal_16((uint16_t*)p1,S->bitRevLength,S->pBitRevTable);
  114. }
  115. /**
  116. * @} end of ComplexFFT group
  117. */
  118. void arm_cfft_radix4by2_q15(
  119. q15_t * pSrc,
  120. uint32_t fftLen,
  121. const q15_t * pCoef)
  122. {
  123. uint32_t i;
  124. uint32_t n2;
  125. q15_t p0, p1, p2, p3;
  126. #if defined (ARM_MATH_DSP)
  127. q31_t T, S, R;
  128. q31_t coeff, out1, out2;
  129. const q15_t *pC = pCoef;
  130. q15_t *pSi = pSrc;
  131. q15_t *pSl = pSrc + fftLen;
  132. #else
  133. uint32_t ia, l;
  134. q15_t xt, yt, cosVal, sinVal;
  135. #endif
  136. n2 = fftLen >> 1;
  137. #if defined (ARM_MATH_DSP)
  138. for (i = n2; i > 0; i--)
  139. {
  140. coeff = _SIMD32_OFFSET(pC);
  141. pC += 2;
  142. T = _SIMD32_OFFSET(pSi);
  143. T = __SHADD16(T, 0); // this is just a SIMD arithmetic shift right by 1
  144. S = _SIMD32_OFFSET(pSl);
  145. S = __SHADD16(S, 0); // this is just a SIMD arithmetic shift right by 1
  146. R = __QSUB16(T, S);
  147. _SIMD32_OFFSET(pSi) = __SHADD16(T, S);
  148. pSi += 2;
  149. #ifndef ARM_MATH_BIG_ENDIAN
  150. out1 = __SMUAD(coeff, R) >> 16;
  151. out2 = __SMUSDX(coeff, R);
  152. #else
  153. out1 = __SMUSDX(R, coeff) >> 16U;
  154. out2 = __SMUAD(coeff, R);
  155. #endif // #ifndef ARM_MATH_BIG_ENDIAN
  156. _SIMD32_OFFSET(pSl) =
  157. (q31_t) ((out2) & 0xFFFF0000) | (out1 & 0x0000FFFF);
  158. pSl += 2;
  159. }
  160. #else // #if defined (ARM_MATH_DSP)
  161. ia = 0;
  162. for (i = 0; i < n2; i++)
  163. {
  164. cosVal = pCoef[ia * 2];
  165. sinVal = pCoef[(ia * 2) + 1];
  166. ia++;
  167. l = i + n2;
  168. xt = (pSrc[2 * i] >> 1U) - (pSrc[2 * l] >> 1U);
  169. pSrc[2 * i] = ((pSrc[2 * i] >> 1U) + (pSrc[2 * l] >> 1U)) >> 1U;
  170. yt = (pSrc[2 * i + 1] >> 1U) - (pSrc[2 * l + 1] >> 1U);
  171. pSrc[2 * i + 1] =
  172. ((pSrc[2 * l + 1] >> 1U) + (pSrc[2 * i + 1] >> 1U)) >> 1U;
  173. pSrc[2U * l] = (((int16_t) (((q31_t) xt * cosVal) >> 16)) +
  174. ((int16_t) (((q31_t) yt * sinVal) >> 16)));
  175. pSrc[2U * l + 1U] = (((int16_t) (((q31_t) yt * cosVal) >> 16)) -
  176. ((int16_t) (((q31_t) xt * sinVal) >> 16)));
  177. }
  178. #endif // #if defined (ARM_MATH_DSP)
  179. // first col
  180. arm_radix4_butterfly_q15( pSrc, n2, (q15_t*)pCoef, 2U);
  181. // second col
  182. arm_radix4_butterfly_q15( pSrc + fftLen, n2, (q15_t*)pCoef, 2U);
  183. for (i = 0; i < fftLen >> 1; i++)
  184. {
  185. p0 = pSrc[4*i+0];
  186. p1 = pSrc[4*i+1];
  187. p2 = pSrc[4*i+2];
  188. p3 = pSrc[4*i+3];
  189. p0 <<= 1;
  190. p1 <<= 1;
  191. p2 <<= 1;
  192. p3 <<= 1;
  193. pSrc[4*i+0] = p0;
  194. pSrc[4*i+1] = p1;
  195. pSrc[4*i+2] = p2;
  196. pSrc[4*i+3] = p3;
  197. }
  198. }
  199. void arm_cfft_radix4by2_inverse_q15(
  200. q15_t * pSrc,
  201. uint32_t fftLen,
  202. const q15_t * pCoef)
  203. {
  204. uint32_t i;
  205. uint32_t n2;
  206. q15_t p0, p1, p2, p3;
  207. #if defined (ARM_MATH_DSP)
  208. q31_t T, S, R;
  209. q31_t coeff, out1, out2;
  210. const q15_t *pC = pCoef;
  211. q15_t *pSi = pSrc;
  212. q15_t *pSl = pSrc + fftLen;
  213. #else
  214. uint32_t ia, l;
  215. q15_t xt, yt, cosVal, sinVal;
  216. #endif
  217. n2 = fftLen >> 1;
  218. #if defined (ARM_MATH_DSP)
  219. for (i = n2; i > 0; i--)
  220. {
  221. coeff = _SIMD32_OFFSET(pC);
  222. pC += 2;
  223. T = _SIMD32_OFFSET(pSi);
  224. T = __SHADD16(T, 0); // this is just a SIMD arithmetic shift right by 1
  225. S = _SIMD32_OFFSET(pSl);
  226. S = __SHADD16(S, 0); // this is just a SIMD arithmetic shift right by 1
  227. R = __QSUB16(T, S);
  228. _SIMD32_OFFSET(pSi) = __SHADD16(T, S);
  229. pSi += 2;
  230. #ifndef ARM_MATH_BIG_ENDIAN
  231. out1 = __SMUSD(coeff, R) >> 16;
  232. out2 = __SMUADX(coeff, R);
  233. #else
  234. out1 = __SMUADX(R, coeff) >> 16U;
  235. out2 = __SMUSD(__QSUB(0, coeff), R);
  236. #endif // #ifndef ARM_MATH_BIG_ENDIAN
  237. _SIMD32_OFFSET(pSl) =
  238. (q31_t) ((out2) & 0xFFFF0000) | (out1 & 0x0000FFFF);
  239. pSl += 2;
  240. }
  241. #else // #if defined (ARM_MATH_DSP)
  242. ia = 0;
  243. for (i = 0; i < n2; i++)
  244. {
  245. cosVal = pCoef[ia * 2];
  246. sinVal = pCoef[(ia * 2) + 1];
  247. ia++;
  248. l = i + n2;
  249. xt = (pSrc[2 * i] >> 1U) - (pSrc[2 * l] >> 1U);
  250. pSrc[2 * i] = ((pSrc[2 * i] >> 1U) + (pSrc[2 * l] >> 1U)) >> 1U;
  251. yt = (pSrc[2 * i + 1] >> 1U) - (pSrc[2 * l + 1] >> 1U);
  252. pSrc[2 * i + 1] =
  253. ((pSrc[2 * l + 1] >> 1U) + (pSrc[2 * i + 1] >> 1U)) >> 1U;
  254. pSrc[2U * l] = (((int16_t) (((q31_t) xt * cosVal) >> 16)) -
  255. ((int16_t) (((q31_t) yt * sinVal) >> 16)));
  256. pSrc[2U * l + 1U] = (((int16_t) (((q31_t) yt * cosVal) >> 16)) +
  257. ((int16_t) (((q31_t) xt * sinVal) >> 16)));
  258. }
  259. #endif // #if defined (ARM_MATH_DSP)
  260. // first col
  261. arm_radix4_butterfly_inverse_q15( pSrc, n2, (q15_t*)pCoef, 2U);
  262. // second col
  263. arm_radix4_butterfly_inverse_q15( pSrc + fftLen, n2, (q15_t*)pCoef, 2U);
  264. for (i = 0; i < fftLen >> 1; i++)
  265. {
  266. p0 = pSrc[4*i+0];
  267. p1 = pSrc[4*i+1];
  268. p2 = pSrc[4*i+2];
  269. p3 = pSrc[4*i+3];
  270. p0 <<= 1;
  271. p1 <<= 1;
  272. p2 <<= 1;
  273. p3 <<= 1;
  274. pSrc[4*i+0] = p0;
  275. pSrc[4*i+1] = p1;
  276. pSrc[4*i+2] = p2;
  277. pSrc[4*i+3] = p3;
  278. }
  279. }