A bundled STM32F10x Std Periph and CMSIS library
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

306 lines
12 KiB

  1. /* ----------------------------------------------------------------------
  2. * Copyright (C) 2010-2014 ARM Limited. All rights reserved.
  3. *
  4. * $Date: 12. March 2014
  5. * $Revision: V1.4.4
  6. *
  7. * Project: CMSIS DSP Library
  8. * Title: arm_biquad_cascade_df1_fast_q31.c
  9. *
  10. * Description: Processing function for the
  11. * Q31 Fast Biquad cascade DirectFormI(DF1) filter.
  12. *
  13. * Target Processor: Cortex-M4/Cortex-M3
  14. *
  15. * Redistribution and use in source and binary forms, with or without
  16. * modification, are permitted provided that the following conditions
  17. * are met:
  18. * - Redistributions of source code must retain the above copyright
  19. * notice, this list of conditions and the following disclaimer.
  20. * - Redistributions in binary form must reproduce the above copyright
  21. * notice, this list of conditions and the following disclaimer in
  22. * the documentation and/or other materials provided with the
  23. * distribution.
  24. * - Neither the name of ARM LIMITED nor the names of its contributors
  25. * may be used to endorse or promote products derived from this
  26. * software without specific prior written permission.
  27. *
  28. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  29. * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  30. * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  31. * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  32. * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  33. * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  34. * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  35. * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  36. * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  37. * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  38. * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  39. * POSSIBILITY OF SUCH DAMAGE.
  40. * -------------------------------------------------------------------- */
  41. #include "arm_math.h"
  42. /**
  43. * @ingroup groupFilters
  44. */
  45. /**
  46. * @addtogroup BiquadCascadeDF1
  47. * @{
  48. */
  49. /**
  50. * @details
  51. *
  52. * @param[in] *S points to an instance of the Q31 Biquad cascade structure.
  53. * @param[in] *pSrc points to the block of input data.
  54. * @param[out] *pDst points to the block of output data.
  55. * @param[in] blockSize number of samples to process per call.
  56. * @return none.
  57. *
  58. * <b>Scaling and Overflow Behavior:</b>
  59. * \par
  60. * This function is optimized for speed at the expense of fixed-point precision and overflow protection.
  61. * The result of each 1.31 x 1.31 multiplication is truncated to 2.30 format.
  62. * These intermediate results are added to a 2.30 accumulator.
  63. * Finally, the accumulator is saturated and converted to a 1.31 result.
  64. * The fast version has the same overflow behavior as the standard version and provides less precision since it discards the low 32 bits of each multiplication result.
  65. * In order to avoid overflows completely the input signal must be scaled down by two bits and lie in the range [-0.25 +0.25). Use the intialization function
  66. * arm_biquad_cascade_df1_init_q31() to initialize filter structure.
  67. *
  68. * \par
  69. * Refer to the function <code>arm_biquad_cascade_df1_q31()</code> for a slower implementation of this function which uses 64-bit accumulation to provide higher precision. Both the slow and the fast versions use the same instance structure.
  70. * Use the function <code>arm_biquad_cascade_df1_init_q31()</code> to initialize the filter structure.
  71. */
  72. void arm_biquad_cascade_df1_fast_q31(
  73. const arm_biquad_casd_df1_inst_q31 * S,
  74. q31_t * pSrc,
  75. q31_t * pDst,
  76. uint32_t blockSize)
  77. {
  78. q31_t acc = 0; /* accumulator */
  79. q31_t Xn1, Xn2, Yn1, Yn2; /* Filter state variables */
  80. q31_t b0, b1, b2, a1, a2; /* Filter coefficients */
  81. q31_t *pIn = pSrc; /* input pointer initialization */
  82. q31_t *pOut = pDst; /* output pointer initialization */
  83. q31_t *pState = S->pState; /* pState pointer initialization */
  84. q31_t *pCoeffs = S->pCoeffs; /* coeff pointer initialization */
  85. q31_t Xn; /* temporary input */
  86. int32_t shift = (int32_t) S->postShift + 1; /* Shift to be applied to the output */
  87. uint32_t sample, stage = S->numStages; /* loop counters */
  88. do
  89. {
  90. /* Reading the coefficients */
  91. b0 = *pCoeffs++;
  92. b1 = *pCoeffs++;
  93. b2 = *pCoeffs++;
  94. a1 = *pCoeffs++;
  95. a2 = *pCoeffs++;
  96. /* Reading the state values */
  97. Xn1 = pState[0];
  98. Xn2 = pState[1];
  99. Yn1 = pState[2];
  100. Yn2 = pState[3];
  101. /* Apply loop unrolling and compute 4 output values simultaneously. */
  102. /* The variables acc ... acc3 hold output values that are being computed:
  103. *
  104. * acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]
  105. */
  106. sample = blockSize >> 2u;
  107. /* First part of the processing with loop unrolling. Compute 4 outputs at a time.
  108. ** a second loop below computes the remaining 1 to 3 samples. */
  109. while(sample > 0u)
  110. {
  111. /* Read the input */
  112. Xn = *pIn;
  113. /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
  114. /* acc = b0 * x[n] */
  115. //acc = (q31_t) (((q63_t) b1 * Xn1) >> 32);
  116. mult_32x32_keep32_R(acc, b1, Xn1);
  117. /* acc += b1 * x[n-1] */
  118. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b0 * (Xn))) >> 32);
  119. multAcc_32x32_keep32_R(acc, b0, Xn);
  120. /* acc += b[2] * x[n-2] */
  121. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn2))) >> 32);
  122. multAcc_32x32_keep32_R(acc, b2, Xn2);
  123. /* acc += a1 * y[n-1] */
  124. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn1))) >> 32);
  125. multAcc_32x32_keep32_R(acc, a1, Yn1);
  126. /* acc += a2 * y[n-2] */
  127. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn2))) >> 32);
  128. multAcc_32x32_keep32_R(acc, a2, Yn2);
  129. /* The result is converted to 1.31 , Yn2 variable is reused */
  130. Yn2 = acc << shift;
  131. /* Read the second input */
  132. Xn2 = *(pIn + 1u);
  133. /* Store the output in the destination buffer. */
  134. *pOut = Yn2;
  135. /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
  136. /* acc = b0 * x[n] */
  137. //acc = (q31_t) (((q63_t) b0 * (Xn2)) >> 32);
  138. mult_32x32_keep32_R(acc, b0, Xn2);
  139. /* acc += b1 * x[n-1] */
  140. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn))) >> 32);
  141. multAcc_32x32_keep32_R(acc, b1, Xn);
  142. /* acc += b[2] * x[n-2] */
  143. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn1))) >> 32);
  144. multAcc_32x32_keep32_R(acc, b2, Xn1);
  145. /* acc += a1 * y[n-1] */
  146. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn2))) >> 32);
  147. multAcc_32x32_keep32_R(acc, a1, Yn2);
  148. /* acc += a2 * y[n-2] */
  149. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn1))) >> 32);
  150. multAcc_32x32_keep32_R(acc, a2, Yn1);
  151. /* The result is converted to 1.31, Yn1 variable is reused */
  152. Yn1 = acc << shift;
  153. /* Read the third input */
  154. Xn1 = *(pIn + 2u);
  155. /* Store the output in the destination buffer. */
  156. *(pOut + 1u) = Yn1;
  157. /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
  158. /* acc = b0 * x[n] */
  159. //acc = (q31_t) (((q63_t) b0 * (Xn1)) >> 32);
  160. mult_32x32_keep32_R(acc, b0, Xn1);
  161. /* acc += b1 * x[n-1] */
  162. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn2))) >> 32);
  163. multAcc_32x32_keep32_R(acc, b1, Xn2);
  164. /* acc += b[2] * x[n-2] */
  165. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn))) >> 32);
  166. multAcc_32x32_keep32_R(acc, b2, Xn);
  167. /* acc += a1 * y[n-1] */
  168. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn1))) >> 32);
  169. multAcc_32x32_keep32_R(acc, a1, Yn1);
  170. /* acc += a2 * y[n-2] */
  171. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn2))) >> 32);
  172. multAcc_32x32_keep32_R(acc, a2, Yn2);
  173. /* The result is converted to 1.31, Yn2 variable is reused */
  174. Yn2 = acc << shift;
  175. /* Read the forth input */
  176. Xn = *(pIn + 3u);
  177. /* Store the output in the destination buffer. */
  178. *(pOut + 2u) = Yn2;
  179. pIn += 4u;
  180. /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
  181. /* acc = b0 * x[n] */
  182. //acc = (q31_t) (((q63_t) b0 * (Xn)) >> 32);
  183. mult_32x32_keep32_R(acc, b0, Xn);
  184. /* acc += b1 * x[n-1] */
  185. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn1))) >> 32);
  186. multAcc_32x32_keep32_R(acc, b1, Xn1);
  187. /* acc += b[2] * x[n-2] */
  188. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn2))) >> 32);
  189. multAcc_32x32_keep32_R(acc, b2, Xn2);
  190. /* acc += a1 * y[n-1] */
  191. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn2))) >> 32);
  192. multAcc_32x32_keep32_R(acc, a1, Yn2);
  193. /* acc += a2 * y[n-2] */
  194. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn1))) >> 32);
  195. multAcc_32x32_keep32_R(acc, a2, Yn1);
  196. /* Every time after the output is computed state should be updated. */
  197. /* The states should be updated as: */
  198. /* Xn2 = Xn1 */
  199. Xn2 = Xn1;
  200. /* The result is converted to 1.31, Yn1 variable is reused */
  201. Yn1 = acc << shift;
  202. /* Xn1 = Xn */
  203. Xn1 = Xn;
  204. /* Store the output in the destination buffer. */
  205. *(pOut + 3u) = Yn1;
  206. pOut += 4u;
  207. /* decrement the loop counter */
  208. sample--;
  209. }
  210. /* If the blockSize is not a multiple of 4, compute any remaining output samples here.
  211. ** No loop unrolling is used. */
  212. sample = (blockSize & 0x3u);
  213. while(sample > 0u)
  214. {
  215. /* Read the input */
  216. Xn = *pIn++;
  217. /* acc = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
  218. /* acc = b0 * x[n] */
  219. //acc = (q31_t) (((q63_t) b0 * (Xn)) >> 32);
  220. mult_32x32_keep32_R(acc, b0, Xn);
  221. /* acc += b1 * x[n-1] */
  222. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b1 * (Xn1))) >> 32);
  223. multAcc_32x32_keep32_R(acc, b1, Xn1);
  224. /* acc += b[2] * x[n-2] */
  225. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) b2 * (Xn2))) >> 32);
  226. multAcc_32x32_keep32_R(acc, b2, Xn2);
  227. /* acc += a1 * y[n-1] */
  228. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a1 * (Yn1))) >> 32);
  229. multAcc_32x32_keep32_R(acc, a1, Yn1);
  230. /* acc += a2 * y[n-2] */
  231. //acc = (q31_t) ((((q63_t) acc << 32) + ((q63_t) a2 * (Yn2))) >> 32);
  232. multAcc_32x32_keep32_R(acc, a2, Yn2);
  233. /* The result is converted to 1.31 */
  234. acc = acc << shift;
  235. /* Every time after the output is computed state should be updated. */
  236. /* The states should be updated as: */
  237. /* Xn2 = Xn1 */
  238. /* Xn1 = Xn */
  239. /* Yn2 = Yn1 */
  240. /* Yn1 = acc */
  241. Xn2 = Xn1;
  242. Xn1 = Xn;
  243. Yn2 = Yn1;
  244. Yn1 = acc;
  245. /* Store the output in the destination buffer. */
  246. *pOut++ = acc;
  247. /* decrement the loop counter */
  248. sample--;
  249. }
  250. /* The first stage goes from the input buffer to the output buffer. */
  251. /* Subsequent stages occur in-place in the output buffer */
  252. pIn = pDst;
  253. /* Reset to destination pointer */
  254. pOut = pDst;
  255. /* Store the updated state variables back into the pState array */
  256. *pState++ = Xn1;
  257. *pState++ = Xn2;
  258. *pState++ = Yn1;
  259. *pState++ = Yn2;
  260. } while(--stage);
  261. }
  262. /**
  263. * @} end of BiquadCascadeDF1 group
  264. */