1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_fir_sparse_f32.c
4 * Description: Floating-point sparse FIR filter processing function
5 *
6 * $Date: 23 April 2021
7 * $Revision: V1.9.0
8 *
9 * Target Processor: Cortex-M and Cortex-A cores
10 * -------------------------------------------------------------------- */
11 /*
12 * Copyright (C) 2010-2021 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
29 #include "dsp/filtering_functions.h"
30
31 /**
32 @ingroup groupFilters
33 */
34
35 /**
36 @defgroup FIR_Sparse Finite Impulse Response (FIR) Sparse Filters
37
38 @deprecated Those functions are no more tested nor maintained and will be removed in
39 a future version.
40
41 This group of functions implements sparse FIR filters.
42 Sparse FIR filters are equivalent to standard FIR filters except that most of the coefficients are equal to zero.
43 Sparse filters are used for simulating reflections in communications and audio applications.
44
45 There are separate functions for Q7, Q15, Q31, and floating-point data types.
46 The functions operate on blocks of input and output data and each call to the function processes
47 <code>blockSize</code> samples through the filter. <code>pSrc</code> and
48 <code>pDst</code> points to input and output arrays respectively containing <code>blockSize</code> values.
49
50 @par Algorithm
51 The sparse filter instant structure contains an array of tap indices <code>pTapDelay</code> which specifies the locations of the non-zero coefficients.
52 This is in addition to the coefficient array <code>b</code>.
53 The implementation essentially skips the multiplications by zero and leads to an efficient realization.
54 <pre>
55 y[n] = b[0] * x[n-pTapDelay[0]] + b[1] * x[n-pTapDelay[1]] + b[2] * x[n-pTapDelay[2]] + ...+ b[numTaps-1] * x[n-pTapDelay[numTaps-1]]
56 </pre>
57 @par
58 \image html FIRSparse.gif "Sparse FIR filter. b[n] represents the filter coefficients"
59 @par
60 <code>pCoeffs</code> points to a coefficient array of size <code>numTaps</code>;
61 <code>pTapDelay</code> points to an array of nonzero indices and is also of size <code>numTaps</code>;
62 <code>pState</code> points to a state array of size <code>maxDelay + blockSize</code>, where
63 <code>maxDelay</code> is the largest offset value that is ever used in the <code>pTapDelay</code> array.
64 Some of the processing functions also require temporary working buffers.
65
66 @par Instance Structure
67 The coefficients and state variables for a filter are stored together in an instance data structure.
68 A separate instance structure must be defined for each filter.
69 Coefficient and offset arrays may be shared among several instances while state variable arrays cannot be shared.
70 There are separate instance structure declarations for each of the 4 supported data types.
71
72 @par Initialization Functions
73 There is also an associated initialization function for each data type.
74 The initialization function performs the following operations:
75 - Sets the values of the internal structure fields.
76 - Zeros out the values in the state buffer.
77 To do this manually without calling the init function, assign the follow subfields of the instance structure:
78 numTaps, pCoeffs, pTapDelay, maxDelay, stateIndex, pState. Also set all of the values in pState to zero.
79 @par
80 Use of the initialization function is optional.
81 However, if the initialization function is used, then the instance structure cannot be placed into a const data section.
82 To place an instance structure into a const data section, the instance structure must be manually initialized.
83 Set the values in the state buffer to zeros before static initialization.
84 The code below statically initializes each of the 4 different data type filter instance structures
85 <pre>
86 arm_fir_sparse_instance_f32 S = {numTaps, 0, pState, pCoeffs, maxDelay, pTapDelay};
87 arm_fir_sparse_instance_q31 S = {numTaps, 0, pState, pCoeffs, maxDelay, pTapDelay};
88 arm_fir_sparse_instance_q15 S = {numTaps, 0, pState, pCoeffs, maxDelay, pTapDelay};
89 arm_fir_sparse_instance_q7 S = {numTaps, 0, pState, pCoeffs, maxDelay, pTapDelay};
90 </pre>
91
92 @par Fixed-Point Behavior
93 Care must be taken when using the fixed-point versions of the sparse FIR filter functions.
94 In particular, the overflow and saturation behavior of the accumulator used in each function must be considered.
95 Refer to the function specific documentation below for usage guidelines.
96 */
97
98 /**
99 @addtogroup FIR_Sparse
100 @{
101 */
102
103 /**
104 @brief Processing function for the floating-point sparse FIR filter.
105 @param[in] S points to an instance of the floating-point sparse FIR structure
106 @param[in] pSrc points to the block of input data
107 @param[out] pDst points to the block of output data
108 @param[in] pScratchIn points to a temporary buffer of size blockSize
109 @param[in] blockSize number of input samples to process
110 @return none
111 */
112
arm_fir_sparse_f32(arm_fir_sparse_instance_f32 * S,const float32_t * pSrc,float32_t * pDst,float32_t * pScratchIn,uint32_t blockSize)113 void arm_fir_sparse_f32(
114 arm_fir_sparse_instance_f32 * S,
115 const float32_t * pSrc,
116 float32_t * pDst,
117 float32_t * pScratchIn,
118 uint32_t blockSize)
119 {
120 float32_t *pState = S->pState; /* State pointer */
121 const float32_t *pCoeffs = S->pCoeffs; /* Coefficient pointer */
122 float32_t *px; /* Scratch buffer pointer */
123 float32_t *py = pState; /* Temporary pointers for state buffer */
124 float32_t *pb = pScratchIn; /* Temporary pointers for scratch buffer */
125 float32_t *pOut; /* Destination pointer */
126 int32_t *pTapDelay = S->pTapDelay; /* Pointer to the array containing offset of the non-zero tap values. */
127 uint32_t delaySize = S->maxDelay + blockSize; /* state length */
128 uint16_t numTaps = S->numTaps; /* Number of filter coefficients in the filter */
129 int32_t readIndex; /* Read index of the state buffer */
130 uint32_t tapCnt, blkCnt; /* loop counters */
131 float32_t coeff = *pCoeffs++; /* Read the first coefficient value */
132
133
134 /* BlockSize of Input samples are copied into the state buffer */
135 /* StateIndex points to the starting position to write in the state buffer */
136 arm_circularWrite_f32((int32_t *) py, delaySize, &S->stateIndex, 1, (int32_t *) pSrc, 1, blockSize);
137
138 /* Read Index, from where the state buffer should be read, is calculated. */
139 readIndex = (int32_t) (S->stateIndex - blockSize) - *pTapDelay++;
140
141 /* Wraparound of readIndex */
142 if (readIndex < 0)
143 {
144 readIndex += (int32_t) delaySize;
145 }
146
147 /* Working pointer for state buffer is updated */
148 py = pState;
149
150 /* blockSize samples are read from the state buffer */
151 arm_circularRead_f32((int32_t *) py, delaySize, &readIndex, 1,
152 (int32_t *) pb, (int32_t *) pb, blockSize, 1, blockSize);
153
154 /* Working pointer for the scratch buffer of state values */
155 px = pb;
156
157 /* Working pointer for scratch buffer of output values */
158 pOut = pDst;
159
160
161 #if defined (ARM_MATH_LOOPUNROLL)
162
163 /* Loop unrolling: Compute 4 outputs at a time. */
164 blkCnt = blockSize >> 2U;
165
166 while (blkCnt > 0U)
167 {
168 /* Perform Multiplications and store in destination buffer */
169 *pOut++ = *px++ * coeff;
170
171 *pOut++ = *px++ * coeff;
172
173 *pOut++ = *px++ * coeff;
174
175 *pOut++ = *px++ * coeff;
176
177 /* Decrement loop counter */
178 blkCnt--;
179 }
180
181 /* Loop unrolling: Compute remaining outputs */
182 blkCnt = blockSize % 0x4U;
183
184 #else
185
186 /* Initialize blkCnt with number of samples */
187 blkCnt = blockSize;
188
189 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
190
191 while (blkCnt > 0U)
192 {
193 /* Perform Multiplication and store in destination buffer */
194 *pOut++ = *px++ * coeff;
195
196 /* Decrement loop counter */
197 blkCnt--;
198 }
199
200 /* Load the coefficient value and
201 * increment the coefficient buffer for the next set of state values */
202 coeff = *pCoeffs++;
203
204 /* Read Index, from where the state buffer should be read, is calculated. */
205 readIndex = (int32_t) (S->stateIndex - blockSize) - *pTapDelay++;
206
207 /* Wraparound of readIndex */
208 if (readIndex < 0)
209 {
210 readIndex += (int32_t) delaySize;
211 }
212
213 /* Loop over the number of taps. */
214 tapCnt = (uint32_t) numTaps - 2U;
215
216 while (tapCnt > 0U)
217 {
218 /* Working pointer for state buffer is updated */
219 py = pState;
220
221 /* blockSize samples are read from the state buffer */
222 arm_circularRead_f32((int32_t *) py, delaySize, &readIndex, 1,
223 (int32_t *) pb, (int32_t *) pb, blockSize, 1, blockSize);
224
225 /* Working pointer for the scratch buffer of state values */
226 px = pb;
227
228 /* Working pointer for scratch buffer of output values */
229 pOut = pDst;
230
231
232 #if defined (ARM_MATH_LOOPUNROLL)
233
234 /* Loop unrolling: Compute 4 outputs at a time. */
235 blkCnt = blockSize >> 2U;
236
237 while (blkCnt > 0U)
238 {
239 /* Perform Multiply-Accumulate */
240 *pOut++ += *px++ * coeff;
241
242 *pOut++ += *px++ * coeff;
243
244 *pOut++ += *px++ * coeff;
245
246 *pOut++ += *px++ * coeff;
247
248 /* Decrement loop counter */
249 blkCnt--;
250 }
251
252 /* Loop unrolling: Compute remaining outputs */
253 blkCnt = blockSize % 0x4U;
254
255 #else
256
257 /* Initialize blkCnt with number of samples */
258 blkCnt = blockSize;
259
260 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
261
262 while (blkCnt > 0U)
263 {
264 /* Perform Multiply-Accumulate */
265 *pOut++ += *px++ * coeff;
266
267 /* Decrement loop counter */
268 blkCnt--;
269 }
270
271 /* Load the coefficient value and
272 * increment the coefficient buffer for the next set of state values */
273 coeff = *pCoeffs++;
274
275 /* Read Index, from where the state buffer should be read, is calculated. */
276 readIndex = (int32_t) (S->stateIndex - blockSize) - *pTapDelay++;
277
278 /* Wraparound of readIndex */
279 if (readIndex < 0)
280 {
281 readIndex += (int32_t) delaySize;
282 }
283
284 /* Decrement tap loop counter */
285 tapCnt--;
286 }
287
288 /* Compute last tap without the final read of pTapDelay */
289
290 /* Working pointer for state buffer is updated */
291 py = pState;
292
293 /* blockSize samples are read from the state buffer */
294 arm_circularRead_f32((int32_t *) py, delaySize, &readIndex, 1,
295 (int32_t *) pb, (int32_t *) pb, blockSize, 1, blockSize);
296
297 /* Working pointer for the scratch buffer of state values */
298 px = pb;
299
300 /* Working pointer for scratch buffer of output values */
301 pOut = pDst;
302
303
304 #if defined (ARM_MATH_LOOPUNROLL)
305
306 /* Loop unrolling: Compute 4 outputs at a time. */
307 blkCnt = blockSize >> 2U;
308
309 while (blkCnt > 0U)
310 {
311 /* Perform Multiply-Accumulate */
312 *pOut++ += *px++ * coeff;
313 *pOut++ += *px++ * coeff;
314 *pOut++ += *px++ * coeff;
315 *pOut++ += *px++ * coeff;
316
317 /* Decrement loop counter */
318 blkCnt--;
319 }
320
321 /* Loop unrolling: Compute remaining outputs */
322 blkCnt = blockSize % 0x4U;
323
324 #else
325
326 /* Initialize blkCnt with number of samples */
327 blkCnt = blockSize;
328
329 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
330
331 while (blkCnt > 0U)
332 {
333 /* Perform Multiply-Accumulate */
334 *pOut++ += *px++ * coeff;
335
336 /* Decrement loop counter */
337 blkCnt--;
338 }
339
340 }
341
342 /**
343 @} end of FIR_Sparse group
344 */
345