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 */
111
arm_fir_sparse_f32(arm_fir_sparse_instance_f32 * S,const float32_t * pSrc,float32_t * pDst,float32_t * pScratchIn,uint32_t blockSize)112 ARM_DSP_ATTRIBUTE void arm_fir_sparse_f32(
113 arm_fir_sparse_instance_f32 * S,
114 const float32_t * pSrc,
115 float32_t * pDst,
116 float32_t * pScratchIn,
117 uint32_t blockSize)
118 {
119 float32_t *pState = S->pState; /* State pointer */
120 const float32_t *pCoeffs = S->pCoeffs; /* Coefficient pointer */
121 float32_t *px; /* Scratch buffer pointer */
122 float32_t *py = pState; /* Temporary pointers for state buffer */
123 float32_t *pb = pScratchIn; /* Temporary pointers for scratch buffer */
124 float32_t *pOut; /* Destination pointer */
125 int32_t *pTapDelay = S->pTapDelay; /* Pointer to the array containing offset of the non-zero tap values. */
126 uint32_t delaySize = S->maxDelay + blockSize; /* state length */
127 uint16_t numTaps = S->numTaps; /* Number of filter coefficients in the filter */
128 int32_t readIndex; /* Read index of the state buffer */
129 uint32_t tapCnt, blkCnt; /* loop counters */
130 float32_t coeff = *pCoeffs++; /* Read the first coefficient value */
131
132
133 /* BlockSize of Input samples are copied into the state buffer */
134 /* StateIndex points to the starting position to write in the state buffer */
135 arm_circularWrite_f32((int32_t *) py, delaySize, &S->stateIndex, 1, (int32_t *) pSrc, 1, blockSize);
136
137 /* Read Index, from where the state buffer should be read, is calculated. */
138 readIndex = (int32_t) (S->stateIndex - blockSize) - *pTapDelay++;
139
140 /* Wraparound of readIndex */
141 if (readIndex < 0)
142 {
143 readIndex += (int32_t) delaySize;
144 }
145
146 /* Working pointer for state buffer is updated */
147 py = pState;
148
149 /* blockSize samples are read from the state buffer */
150 arm_circularRead_f32((int32_t *) py, delaySize, &readIndex, 1,
151 (int32_t *) pb, (int32_t *) pb, blockSize, 1, blockSize);
152
153 /* Working pointer for the scratch buffer of state values */
154 px = pb;
155
156 /* Working pointer for scratch buffer of output values */
157 pOut = pDst;
158
159
160 #if defined (ARM_MATH_LOOPUNROLL)
161
162 /* Loop unrolling: Compute 4 outputs at a time. */
163 blkCnt = blockSize >> 2U;
164
165 while (blkCnt > 0U)
166 {
167 /* Perform Multiplications and store in destination buffer */
168 *pOut++ = *px++ * coeff;
169
170 *pOut++ = *px++ * coeff;
171
172 *pOut++ = *px++ * coeff;
173
174 *pOut++ = *px++ * coeff;
175
176 /* Decrement loop counter */
177 blkCnt--;
178 }
179
180 /* Loop unrolling: Compute remaining outputs */
181 blkCnt = blockSize % 0x4U;
182
183 #else
184
185 /* Initialize blkCnt with number of samples */
186 blkCnt = blockSize;
187
188 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
189
190 while (blkCnt > 0U)
191 {
192 /* Perform Multiplication and store in destination buffer */
193 *pOut++ = *px++ * coeff;
194
195 /* Decrement loop counter */
196 blkCnt--;
197 }
198
199 /* Load the coefficient value and
200 * increment the coefficient buffer for the next set of state values */
201 coeff = *pCoeffs++;
202
203 /* Read Index, from where the state buffer should be read, is calculated. */
204 readIndex = (int32_t) (S->stateIndex - blockSize) - *pTapDelay++;
205
206 /* Wraparound of readIndex */
207 if (readIndex < 0)
208 {
209 readIndex += (int32_t) delaySize;
210 }
211
212 /* Loop over the number of taps. */
213 tapCnt = (uint32_t) numTaps - 2U;
214
215 while (tapCnt > 0U)
216 {
217 /* Working pointer for state buffer is updated */
218 py = pState;
219
220 /* blockSize samples are read from the state buffer */
221 arm_circularRead_f32((int32_t *) py, delaySize, &readIndex, 1,
222 (int32_t *) pb, (int32_t *) pb, blockSize, 1, blockSize);
223
224 /* Working pointer for the scratch buffer of state values */
225 px = pb;
226
227 /* Working pointer for scratch buffer of output values */
228 pOut = pDst;
229
230
231 #if defined (ARM_MATH_LOOPUNROLL)
232
233 /* Loop unrolling: Compute 4 outputs at a time. */
234 blkCnt = blockSize >> 2U;
235
236 while (blkCnt > 0U)
237 {
238 /* Perform Multiply-Accumulate */
239 *pOut++ += *px++ * coeff;
240
241 *pOut++ += *px++ * coeff;
242
243 *pOut++ += *px++ * coeff;
244
245 *pOut++ += *px++ * coeff;
246
247 /* Decrement loop counter */
248 blkCnt--;
249 }
250
251 /* Loop unrolling: Compute remaining outputs */
252 blkCnt = blockSize % 0x4U;
253
254 #else
255
256 /* Initialize blkCnt with number of samples */
257 blkCnt = blockSize;
258
259 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
260
261 while (blkCnt > 0U)
262 {
263 /* Perform Multiply-Accumulate */
264 *pOut++ += *px++ * coeff;
265
266 /* Decrement loop counter */
267 blkCnt--;
268 }
269
270 /* Load the coefficient value and
271 * increment the coefficient buffer for the next set of state values */
272 coeff = *pCoeffs++;
273
274 /* Read Index, from where the state buffer should be read, is calculated. */
275 readIndex = (int32_t) (S->stateIndex - blockSize) - *pTapDelay++;
276
277 /* Wraparound of readIndex */
278 if (readIndex < 0)
279 {
280 readIndex += (int32_t) delaySize;
281 }
282
283 /* Decrement tap loop counter */
284 tapCnt--;
285 }
286
287 /* Compute last tap without the final read of pTapDelay */
288
289 /* Working pointer for state buffer is updated */
290 py = pState;
291
292 /* blockSize samples are read from the state buffer */
293 arm_circularRead_f32((int32_t *) py, delaySize, &readIndex, 1,
294 (int32_t *) pb, (int32_t *) pb, blockSize, 1, blockSize);
295
296 /* Working pointer for the scratch buffer of state values */
297 px = pb;
298
299 /* Working pointer for scratch buffer of output values */
300 pOut = pDst;
301
302
303 #if defined (ARM_MATH_LOOPUNROLL)
304
305 /* Loop unrolling: Compute 4 outputs at a time. */
306 blkCnt = blockSize >> 2U;
307
308 while (blkCnt > 0U)
309 {
310 /* Perform Multiply-Accumulate */
311 *pOut++ += *px++ * coeff;
312 *pOut++ += *px++ * coeff;
313 *pOut++ += *px++ * coeff;
314 *pOut++ += *px++ * coeff;
315
316 /* Decrement loop counter */
317 blkCnt--;
318 }
319
320 /* Loop unrolling: Compute remaining outputs */
321 blkCnt = blockSize % 0x4U;
322
323 #else
324
325 /* Initialize blkCnt with number of samples */
326 blkCnt = blockSize;
327
328 #endif /* #if defined (ARM_MATH_LOOPUNROLL) */
329
330 while (blkCnt > 0U)
331 {
332 /* Perform Multiply-Accumulate */
333 *pOut++ += *px++ * coeff;
334
335 /* Decrement loop counter */
336 blkCnt--;
337 }
338
339 }
340
341 /**
342 @} end of FIR_Sparse group
343 */
344
