GCC Code Coverage Report


Directory: ./
File: kernels/volk/volk_32f_s32f_calc_spectral_noise_floor_32f.h
Date: 2023-10-23 23:10:04
Exec Total Coverage
Lines: 189 194 97.4%
Functions: 4 4 100.0%
Branches: 36 44 81.8%

Line Branch Exec Source
1 /* -*- c++ -*- */
2 /*
3 * Copyright 2012, 2014 Free Software Foundation, Inc.
4 *
5 * This file is part of VOLK
6 *
7 * SPDX-License-Identifier: LGPL-3.0-or-later
8 */
9
10 /*!
11 * \page volk_32f_s32f_calc_spectral_noise_floor_32f
12 *
13 * \b Overview
14 *
15 * Computes the spectral noise floor of an input power spectrum.
16 *
17 * Calculates the spectral noise floor of an input power spectrum by
18 * determining the mean of the input power spectrum, then
19 * recalculating the mean excluding any power spectrum values that
20 * exceed the mean by the spectralExclusionValue (in dB). Provides a
21 * rough estimation of the signal noise floor.
22 *
23 * <b>Dispatcher Prototype</b>
24 * \code
25 * void volk_32f_s32f_calc_spectral_noise_floor_32f(float* noiseFloorAmplitude, const
26 * float* realDataPoints, const float spectralExclusionValue, const unsigned int
27 * num_points) \endcode
28 *
29 * \b Inputs
30 * \li realDataPoints: The input power spectrum.
31 * \li spectralExclusionValue: The number of dB above the noise floor that a data point
32 * must be to be excluded from the noise floor calculation - default value is 20. \li
33 * num_points: The number of data points.
34 *
35 * \b Outputs
36 * \li noiseFloorAmplitude: The noise floor of the input spectrum, in dB.
37 *
38 * \b Example
39 * \code
40 * int N = 10000;
41 *
42 * volk_32f_s32f_calc_spectral_noise_floor_32f
43 *
44 * volk_free(x);
45 * \endcode
46 */
47
48 #ifndef INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_a_H
49 #define INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_a_H
50
51 #include <inttypes.h>
52 #include <stdio.h>
53 #include <volk/volk_common.h>
54
55 #ifdef LV_HAVE_AVX
56 #include <immintrin.h>
57
58 static inline void
59 2 volk_32f_s32f_calc_spectral_noise_floor_32f_a_avx(float* noiseFloorAmplitude,
60 const float* realDataPoints,
61 const float spectralExclusionValue,
62 const unsigned int num_points)
63 {
64 2 unsigned int number = 0;
65 2 const unsigned int eighthPoints = num_points / 8;
66
67 2 const float* dataPointsPtr = realDataPoints;
68 __VOLK_ATTR_ALIGNED(32) float avgPointsVector[8];
69
70 __m256 dataPointsVal;
71 2 __m256 avgPointsVal = _mm256_setzero_ps();
72 // Calculate the sum (for mean) for all points
73
2/2
✓ Branch 0 taken 32766 times.
✓ Branch 1 taken 2 times.
32768 for (; number < eighthPoints; number++) {
74
75 32766 dataPointsVal = _mm256_load_ps(dataPointsPtr);
76
77 32766 dataPointsPtr += 8;
78
79 32766 avgPointsVal = _mm256_add_ps(avgPointsVal, dataPointsVal);
80 }
81
82 _mm256_store_ps(avgPointsVector, avgPointsVal);
83
84 2 float sumMean = 0.0;
85 2 sumMean += avgPointsVector[0];
86 2 sumMean += avgPointsVector[1];
87 2 sumMean += avgPointsVector[2];
88 2 sumMean += avgPointsVector[3];
89 2 sumMean += avgPointsVector[4];
90 2 sumMean += avgPointsVector[5];
91 2 sumMean += avgPointsVector[6];
92 2 sumMean += avgPointsVector[7];
93
94 2 number = eighthPoints * 8;
95
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 2 times.
16 for (; number < num_points; number++) {
96 14 sumMean += realDataPoints[number];
97 }
98
99 // calculate the spectral mean
100 // +20 because for the comparison below we only want to throw out bins
101 // that are significantly higher (and would, thus, affect the mean more
102 2 const float meanAmplitude = (sumMean / ((float)num_points)) + spectralExclusionValue;
103
104 2 dataPointsPtr = realDataPoints; // Reset the dataPointsPtr
105 2 __m256 vMeanAmplitudeVector = _mm256_set1_ps(meanAmplitude);
106 2 __m256 vOnesVector = _mm256_set1_ps(1.0);
107 2 __m256 vValidBinCount = _mm256_setzero_ps();
108 2 avgPointsVal = _mm256_setzero_ps();
109 __m256 compareMask;
110 2 number = 0;
111 // Calculate the sum (for mean) for any points which do NOT exceed the mean amplitude
112
2/2
✓ Branch 0 taken 32766 times.
✓ Branch 1 taken 2 times.
32768 for (; number < eighthPoints; number++) {
113
114 32766 dataPointsVal = _mm256_load_ps(dataPointsPtr);
115
116 32766 dataPointsPtr += 8;
117
118 // Identify which items do not exceed the mean amplitude
119 32766 compareMask = _mm256_cmp_ps(dataPointsVal, vMeanAmplitudeVector, _CMP_LE_OQ);
120
121 // Mask off the items that exceed the mean amplitude and add the avg Points that
122 // do not exceed the mean amplitude
123 avgPointsVal =
124 65532 _mm256_add_ps(avgPointsVal, _mm256_and_ps(compareMask, dataPointsVal));
125
126 // Count the number of bins which do not exceed the mean amplitude
127 vValidBinCount =
128 65532 _mm256_add_ps(vValidBinCount, _mm256_and_ps(compareMask, vOnesVector));
129 }
130
131 // Calculate the mean from the remaining data points
132 _mm256_store_ps(avgPointsVector, avgPointsVal);
133
134 2 sumMean = 0.0;
135 2 sumMean += avgPointsVector[0];
136 2 sumMean += avgPointsVector[1];
137 2 sumMean += avgPointsVector[2];
138 2 sumMean += avgPointsVector[3];
139 2 sumMean += avgPointsVector[4];
140 2 sumMean += avgPointsVector[5];
141 2 sumMean += avgPointsVector[6];
142 2 sumMean += avgPointsVector[7];
143
144 // Calculate the number of valid bins from the remaining count
145 __VOLK_ATTR_ALIGNED(32) float validBinCountVector[8];
146 _mm256_store_ps(validBinCountVector, vValidBinCount);
147
148 2 float validBinCount = 0;
149 2 validBinCount += validBinCountVector[0];
150 2 validBinCount += validBinCountVector[1];
151 2 validBinCount += validBinCountVector[2];
152 2 validBinCount += validBinCountVector[3];
153 2 validBinCount += validBinCountVector[4];
154 2 validBinCount += validBinCountVector[5];
155 2 validBinCount += validBinCountVector[6];
156 2 validBinCount += validBinCountVector[7];
157
158 2 number = eighthPoints * 8;
159
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 2 times.
16 for (; number < num_points; number++) {
160
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
14 if (realDataPoints[number] <= meanAmplitude) {
161 14 sumMean += realDataPoints[number];
162 14 validBinCount += 1.0;
163 }
164 }
165
166 2 float localNoiseFloorAmplitude = 0;
167
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 if (validBinCount > 0.0) {
168 2 localNoiseFloorAmplitude = sumMean / validBinCount;
169 } else {
170 localNoiseFloorAmplitude =
171 meanAmplitude; // For the odd case that all the amplitudes are equal...
172 }
173
174 2 *noiseFloorAmplitude = localNoiseFloorAmplitude;
175 2 }
176 #endif /* LV_HAVE_AVX */
177
178 #ifdef LV_HAVE_SSE
179 #include <xmmintrin.h>
180
181 static inline void
182 2 volk_32f_s32f_calc_spectral_noise_floor_32f_a_sse(float* noiseFloorAmplitude,
183 const float* realDataPoints,
184 const float spectralExclusionValue,
185 const unsigned int num_points)
186 {
187 2 unsigned int number = 0;
188 2 const unsigned int quarterPoints = num_points / 4;
189
190 2 const float* dataPointsPtr = realDataPoints;
191 __VOLK_ATTR_ALIGNED(16) float avgPointsVector[4];
192
193 __m128 dataPointsVal;
194 2 __m128 avgPointsVal = _mm_setzero_ps();
195 // Calculate the sum (for mean) for all points
196
2/2
✓ Branch 0 taken 65534 times.
✓ Branch 1 taken 2 times.
65536 for (; number < quarterPoints; number++) {
197
198 65534 dataPointsVal = _mm_load_ps(dataPointsPtr);
199
200 65534 dataPointsPtr += 4;
201
202 65534 avgPointsVal = _mm_add_ps(avgPointsVal, dataPointsVal);
203 }
204
205 _mm_store_ps(avgPointsVector, avgPointsVal);
206
207 2 float sumMean = 0.0;
208 2 sumMean += avgPointsVector[0];
209 2 sumMean += avgPointsVector[1];
210 2 sumMean += avgPointsVector[2];
211 2 sumMean += avgPointsVector[3];
212
213 2 number = quarterPoints * 4;
214
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
8 for (; number < num_points; number++) {
215 6 sumMean += realDataPoints[number];
216 }
217
218 // calculate the spectral mean
219 // +20 because for the comparison below we only want to throw out bins
220 // that are significantly higher (and would, thus, affect the mean more
221 2 const float meanAmplitude = (sumMean / ((float)num_points)) + spectralExclusionValue;
222
223 2 dataPointsPtr = realDataPoints; // Reset the dataPointsPtr
224 2 __m128 vMeanAmplitudeVector = _mm_set_ps1(meanAmplitude);
225 2 __m128 vOnesVector = _mm_set_ps1(1.0);
226 2 __m128 vValidBinCount = _mm_setzero_ps();
227 2 avgPointsVal = _mm_setzero_ps();
228 __m128 compareMask;
229 2 number = 0;
230 // Calculate the sum (for mean) for any points which do NOT exceed the mean amplitude
231
2/2
✓ Branch 0 taken 65534 times.
✓ Branch 1 taken 2 times.
65536 for (; number < quarterPoints; number++) {
232
233 65534 dataPointsVal = _mm_load_ps(dataPointsPtr);
234
235 65534 dataPointsPtr += 4;
236
237 // Identify which items do not exceed the mean amplitude
238 65534 compareMask = _mm_cmple_ps(dataPointsVal, vMeanAmplitudeVector);
239
240 // Mask off the items that exceed the mean amplitude and add the avg Points that
241 // do not exceed the mean amplitude
242 131068 avgPointsVal = _mm_add_ps(avgPointsVal, _mm_and_ps(compareMask, dataPointsVal));
243
244 // Count the number of bins which do not exceed the mean amplitude
245 131068 vValidBinCount = _mm_add_ps(vValidBinCount, _mm_and_ps(compareMask, vOnesVector));
246 }
247
248 // Calculate the mean from the remaining data points
249 _mm_store_ps(avgPointsVector, avgPointsVal);
250
251 2 sumMean = 0.0;
252 2 sumMean += avgPointsVector[0];
253 2 sumMean += avgPointsVector[1];
254 2 sumMean += avgPointsVector[2];
255 2 sumMean += avgPointsVector[3];
256
257 // Calculate the number of valid bins from the remaining count
258 __VOLK_ATTR_ALIGNED(16) float validBinCountVector[4];
259 _mm_store_ps(validBinCountVector, vValidBinCount);
260
261 2 float validBinCount = 0;
262 2 validBinCount += validBinCountVector[0];
263 2 validBinCount += validBinCountVector[1];
264 2 validBinCount += validBinCountVector[2];
265 2 validBinCount += validBinCountVector[3];
266
267 2 number = quarterPoints * 4;
268
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
8 for (; number < num_points; number++) {
269
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
6 if (realDataPoints[number] <= meanAmplitude) {
270 6 sumMean += realDataPoints[number];
271 6 validBinCount += 1.0;
272 }
273 }
274
275 2 float localNoiseFloorAmplitude = 0;
276
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 if (validBinCount > 0.0) {
277 2 localNoiseFloorAmplitude = sumMean / validBinCount;
278 } else {
279 localNoiseFloorAmplitude =
280 meanAmplitude; // For the odd case that all the amplitudes are equal...
281 }
282
283 2 *noiseFloorAmplitude = localNoiseFloorAmplitude;
284 2 }
285 #endif /* LV_HAVE_SSE */
286
287
288 #ifdef LV_HAVE_GENERIC
289
290 static inline void
291 2 volk_32f_s32f_calc_spectral_noise_floor_32f_generic(float* noiseFloorAmplitude,
292 const float* realDataPoints,
293 const float spectralExclusionValue,
294 const unsigned int num_points)
295 {
296 2 float sumMean = 0.0;
297 unsigned int number;
298 // find the sum (for mean), etc
299
2/2
✓ Branch 0 taken 262142 times.
✓ Branch 1 taken 2 times.
262144 for (number = 0; number < num_points; number++) {
300 // sum (for mean)
301 262142 sumMean += realDataPoints[number];
302 }
303
304 // calculate the spectral mean
305 // +20 because for the comparison below we only want to throw out bins
306 // that are significantly higher (and would, thus, affect the mean more)
307 2 const float meanAmplitude = (sumMean / num_points) + spectralExclusionValue;
308
309 // now throw out any bins higher than the mean
310 2 sumMean = 0.0;
311 2 unsigned int newNumDataPoints = num_points;
312
2/2
✓ Branch 0 taken 262142 times.
✓ Branch 1 taken 2 times.
262144 for (number = 0; number < num_points; number++) {
313
1/2
✓ Branch 0 taken 262142 times.
✗ Branch 1 not taken.
262142 if (realDataPoints[number] <= meanAmplitude)
314 262142 sumMean += realDataPoints[number];
315 else
316 newNumDataPoints--;
317 }
318
319 2 float localNoiseFloorAmplitude = 0.0;
320
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if (newNumDataPoints == 0) // in the odd case that all
321 localNoiseFloorAmplitude = meanAmplitude; // amplitudes are equal!
322 else
323 2 localNoiseFloorAmplitude = sumMean / ((float)newNumDataPoints);
324
325 2 *noiseFloorAmplitude = localNoiseFloorAmplitude;
326 2 }
327 #endif /* LV_HAVE_GENERIC */
328
329
330 #endif /* INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_a_H */
331
332 #ifndef INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_u_H
333 #define INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_u_H
334
335 #include <inttypes.h>
336 #include <stdio.h>
337 #include <volk/volk_common.h>
338
339 #ifdef LV_HAVE_AVX
340 #include <immintrin.h>
341
342 static inline void
343 2 volk_32f_s32f_calc_spectral_noise_floor_32f_u_avx(float* noiseFloorAmplitude,
344 const float* realDataPoints,
345 const float spectralExclusionValue,
346 const unsigned int num_points)
347 {
348 2 unsigned int number = 0;
349 2 const unsigned int eighthPoints = num_points / 8;
350
351 2 const float* dataPointsPtr = realDataPoints;
352 __VOLK_ATTR_ALIGNED(16) float avgPointsVector[8];
353
354 __m256 dataPointsVal;
355 2 __m256 avgPointsVal = _mm256_setzero_ps();
356 // Calculate the sum (for mean) for all points
357
2/2
✓ Branch 0 taken 32766 times.
✓ Branch 1 taken 2 times.
32768 for (; number < eighthPoints; number++) {
358
359 32766 dataPointsVal = _mm256_loadu_ps(dataPointsPtr);
360
361 32766 dataPointsPtr += 8;
362
363 32766 avgPointsVal = _mm256_add_ps(avgPointsVal, dataPointsVal);
364 }
365
366 _mm256_storeu_ps(avgPointsVector, avgPointsVal);
367
368 2 float sumMean = 0.0;
369 2 sumMean += avgPointsVector[0];
370 2 sumMean += avgPointsVector[1];
371 2 sumMean += avgPointsVector[2];
372 2 sumMean += avgPointsVector[3];
373 2 sumMean += avgPointsVector[4];
374 2 sumMean += avgPointsVector[5];
375 2 sumMean += avgPointsVector[6];
376 2 sumMean += avgPointsVector[7];
377
378 2 number = eighthPoints * 8;
379
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 2 times.
16 for (; number < num_points; number++) {
380 14 sumMean += realDataPoints[number];
381 }
382
383 // calculate the spectral mean
384 // +20 because for the comparison below we only want to throw out bins
385 // that are significantly higher (and would, thus, affect the mean more
386 2 const float meanAmplitude = (sumMean / ((float)num_points)) + spectralExclusionValue;
387
388 2 dataPointsPtr = realDataPoints; // Reset the dataPointsPtr
389 2 __m256 vMeanAmplitudeVector = _mm256_set1_ps(meanAmplitude);
390 2 __m256 vOnesVector = _mm256_set1_ps(1.0);
391 2 __m256 vValidBinCount = _mm256_setzero_ps();
392 2 avgPointsVal = _mm256_setzero_ps();
393 __m256 compareMask;
394 2 number = 0;
395 // Calculate the sum (for mean) for any points which do NOT exceed the mean amplitude
396
2/2
✓ Branch 0 taken 32766 times.
✓ Branch 1 taken 2 times.
32768 for (; number < eighthPoints; number++) {
397
398 32766 dataPointsVal = _mm256_loadu_ps(dataPointsPtr);
399
400 32766 dataPointsPtr += 8;
401
402 // Identify which items do not exceed the mean amplitude
403 32766 compareMask = _mm256_cmp_ps(dataPointsVal, vMeanAmplitudeVector, _CMP_LE_OQ);
404
405 // Mask off the items that exceed the mean amplitude and add the avg Points that
406 // do not exceed the mean amplitude
407 avgPointsVal =
408 65532 _mm256_add_ps(avgPointsVal, _mm256_and_ps(compareMask, dataPointsVal));
409
410 // Count the number of bins which do not exceed the mean amplitude
411 vValidBinCount =
412 65532 _mm256_add_ps(vValidBinCount, _mm256_and_ps(compareMask, vOnesVector));
413 }
414
415 // Calculate the mean from the remaining data points
416 _mm256_storeu_ps(avgPointsVector, avgPointsVal);
417
418 2 sumMean = 0.0;
419 2 sumMean += avgPointsVector[0];
420 2 sumMean += avgPointsVector[1];
421 2 sumMean += avgPointsVector[2];
422 2 sumMean += avgPointsVector[3];
423 2 sumMean += avgPointsVector[4];
424 2 sumMean += avgPointsVector[5];
425 2 sumMean += avgPointsVector[6];
426 2 sumMean += avgPointsVector[7];
427
428 // Calculate the number of valid bins from the remaining count
429 __VOLK_ATTR_ALIGNED(16) float validBinCountVector[8];
430 _mm256_storeu_ps(validBinCountVector, vValidBinCount);
431
432 2 float validBinCount = 0;
433 2 validBinCount += validBinCountVector[0];
434 2 validBinCount += validBinCountVector[1];
435 2 validBinCount += validBinCountVector[2];
436 2 validBinCount += validBinCountVector[3];
437 2 validBinCount += validBinCountVector[4];
438 2 validBinCount += validBinCountVector[5];
439 2 validBinCount += validBinCountVector[6];
440 2 validBinCount += validBinCountVector[7];
441
442 2 number = eighthPoints * 8;
443
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 2 times.
16 for (; number < num_points; number++) {
444
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
14 if (realDataPoints[number] <= meanAmplitude) {
445 14 sumMean += realDataPoints[number];
446 14 validBinCount += 1.0;
447 }
448 }
449
450 2 float localNoiseFloorAmplitude = 0;
451
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
2 if (validBinCount > 0.0) {
452 2 localNoiseFloorAmplitude = sumMean / validBinCount;
453 } else {
454 localNoiseFloorAmplitude =
455 meanAmplitude; // For the odd case that all the amplitudes are equal...
456 }
457
458 2 *noiseFloorAmplitude = localNoiseFloorAmplitude;
459 2 }
460 #endif /* LV_HAVE_AVX */
461 #endif /* INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_u_H */
462