Vector Optimized Library of Kernels  2.0
Architecture-tuned implementations of math kernels
volk_32f_s32f_calc_spectral_noise_floor_32f.h
Go to the documentation of this file.
1 /* -*- c++ -*- */
2 /*
3  * Copyright 2012, 2014 Free Software Foundation, Inc.
4  *
5  * This file is part of GNU Radio
6  *
7  * GNU Radio is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 3, or (at your option)
10  * any later version.
11  *
12  * GNU Radio is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with GNU Radio; see the file COPYING. If not, write to
19  * the Free Software Foundation, Inc., 51 Franklin Street,
20  * Boston, MA 02110-1301, USA.
21  */
22 
59 #ifndef INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_a_H
60 #define INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_a_H
61 
62 #include <volk/volk_common.h>
63 #include <inttypes.h>
64 #include <stdio.h>
65 
66 #ifdef LV_HAVE_AVX
67 #include <immintrin.h>
68 
69 static inline void
71  const float* realDataPoints,
72  const float spectralExclusionValue,
73  const unsigned int num_points)
74 {
75  unsigned int number = 0;
76  const unsigned int eighthPoints = num_points / 8;
77 
78  const float* dataPointsPtr = realDataPoints;
79  __VOLK_ATTR_ALIGNED(32) float avgPointsVector[8];
80 
81  __m256 dataPointsVal;
82  __m256 avgPointsVal = _mm256_setzero_ps();
83  // Calculate the sum (for mean) for all points
84  for(; number < eighthPoints; number++){
85 
86  dataPointsVal = _mm256_load_ps(dataPointsPtr);
87 
88  dataPointsPtr += 8;
89 
90  avgPointsVal = _mm256_add_ps(avgPointsVal, dataPointsVal);
91  }
92 
93  _mm256_store_ps(avgPointsVector, avgPointsVal);
94 
95  float sumMean = 0.0;
96  sumMean += avgPointsVector[0];
97  sumMean += avgPointsVector[1];
98  sumMean += avgPointsVector[2];
99  sumMean += avgPointsVector[3];
100  sumMean += avgPointsVector[4];
101  sumMean += avgPointsVector[5];
102  sumMean += avgPointsVector[6];
103  sumMean += avgPointsVector[7];
104 
105  number = eighthPoints * 8;
106  for(;number < num_points; number++){
107  sumMean += realDataPoints[number];
108  }
109 
110  // calculate the spectral mean
111  // +20 because for the comparison below we only want to throw out bins
112  // that are significantly higher (and would, thus, affect the mean more
113  const float meanAmplitude = (sumMean / ((float)num_points)) + spectralExclusionValue;
114 
115  dataPointsPtr = realDataPoints; // Reset the dataPointsPtr
116  __m256 vMeanAmplitudeVector = _mm256_set1_ps(meanAmplitude);
117  __m256 vOnesVector = _mm256_set1_ps(1.0);
118  __m256 vValidBinCount = _mm256_setzero_ps();
119  avgPointsVal = _mm256_setzero_ps();
120  __m256 compareMask;
121  number = 0;
122  // Calculate the sum (for mean) for any points which do NOT exceed the mean amplitude
123  for(; number < eighthPoints; number++){
124 
125  dataPointsVal = _mm256_load_ps(dataPointsPtr);
126 
127  dataPointsPtr += 8;
128 
129  // Identify which items do not exceed the mean amplitude
130  compareMask = _mm256_cmp_ps(dataPointsVal, vMeanAmplitudeVector, 18);
131 
132  // Mask off the items that exceed the mean amplitude and add the avg Points that do not exceed the mean amplitude
133  avgPointsVal = _mm256_add_ps(avgPointsVal, _mm256_and_ps(compareMask, dataPointsVal));
134 
135  // Count the number of bins which do not exceed the mean amplitude
136  vValidBinCount = _mm256_add_ps(vValidBinCount, _mm256_and_ps(compareMask, vOnesVector));
137  }
138 
139  // Calculate the mean from the remaining data points
140  _mm256_store_ps(avgPointsVector, avgPointsVal);
141 
142  sumMean = 0.0;
143  sumMean += avgPointsVector[0];
144  sumMean += avgPointsVector[1];
145  sumMean += avgPointsVector[2];
146  sumMean += avgPointsVector[3];
147  sumMean += avgPointsVector[4];
148  sumMean += avgPointsVector[5];
149  sumMean += avgPointsVector[6];
150  sumMean += avgPointsVector[7];
151 
152  // Calculate the number of valid bins from the remaning count
153  __VOLK_ATTR_ALIGNED(32) float validBinCountVector[8];
154  _mm256_store_ps(validBinCountVector, vValidBinCount);
155 
156  float validBinCount = 0;
157  validBinCount += validBinCountVector[0];
158  validBinCount += validBinCountVector[1];
159  validBinCount += validBinCountVector[2];
160  validBinCount += validBinCountVector[3];
161  validBinCount += validBinCountVector[4];
162  validBinCount += validBinCountVector[5];
163  validBinCount += validBinCountVector[6];
164  validBinCount += validBinCountVector[7];
165 
166  number = eighthPoints * 8;
167  for(;number < num_points; number++){
168  if(realDataPoints[number] <= meanAmplitude){
169  sumMean += realDataPoints[number];
170  validBinCount += 1.0;
171  }
172  }
173 
174  float localNoiseFloorAmplitude = 0;
175  if(validBinCount > 0.0){
176  localNoiseFloorAmplitude = sumMean / validBinCount;
177  }
178  else{
179  localNoiseFloorAmplitude = meanAmplitude; // For the odd case that all the amplitudes are equal...
180  }
181 
182  *noiseFloorAmplitude = localNoiseFloorAmplitude;
183 }
184 #endif /* LV_HAVE_AVX */
185 
186 #ifdef LV_HAVE_SSE
187 #include <xmmintrin.h>
188 
189 static inline void
191  const float* realDataPoints,
192  const float spectralExclusionValue,
193  const unsigned int num_points)
194 {
195  unsigned int number = 0;
196  const unsigned int quarterPoints = num_points / 4;
197 
198  const float* dataPointsPtr = realDataPoints;
199  __VOLK_ATTR_ALIGNED(16) float avgPointsVector[4];
200 
201  __m128 dataPointsVal;
202  __m128 avgPointsVal = _mm_setzero_ps();
203  // Calculate the sum (for mean) for all points
204  for(; number < quarterPoints; number++){
205 
206  dataPointsVal = _mm_load_ps(dataPointsPtr);
207 
208  dataPointsPtr += 4;
209 
210  avgPointsVal = _mm_add_ps(avgPointsVal, dataPointsVal);
211  }
212 
213  _mm_store_ps(avgPointsVector, avgPointsVal);
214 
215  float sumMean = 0.0;
216  sumMean += avgPointsVector[0];
217  sumMean += avgPointsVector[1];
218  sumMean += avgPointsVector[2];
219  sumMean += avgPointsVector[3];
220 
221  number = quarterPoints * 4;
222  for(;number < num_points; number++){
223  sumMean += realDataPoints[number];
224  }
225 
226  // calculate the spectral mean
227  // +20 because for the comparison below we only want to throw out bins
228  // that are significantly higher (and would, thus, affect the mean more
229  const float meanAmplitude = (sumMean / ((float)num_points)) + spectralExclusionValue;
230 
231  dataPointsPtr = realDataPoints; // Reset the dataPointsPtr
232  __m128 vMeanAmplitudeVector = _mm_set_ps1(meanAmplitude);
233  __m128 vOnesVector = _mm_set_ps1(1.0);
234  __m128 vValidBinCount = _mm_setzero_ps();
235  avgPointsVal = _mm_setzero_ps();
236  __m128 compareMask;
237  number = 0;
238  // Calculate the sum (for mean) for any points which do NOT exceed the mean amplitude
239  for(; number < quarterPoints; number++){
240 
241  dataPointsVal = _mm_load_ps(dataPointsPtr);
242 
243  dataPointsPtr += 4;
244 
245  // Identify which items do not exceed the mean amplitude
246  compareMask = _mm_cmple_ps(dataPointsVal, vMeanAmplitudeVector);
247 
248  // Mask off the items that exceed the mean amplitude and add the avg Points that do not exceed the mean amplitude
249  avgPointsVal = _mm_add_ps(avgPointsVal, _mm_and_ps(compareMask, dataPointsVal));
250 
251  // Count the number of bins which do not exceed the mean amplitude
252  vValidBinCount = _mm_add_ps(vValidBinCount, _mm_and_ps(compareMask, vOnesVector));
253  }
254 
255  // Calculate the mean from the remaining data points
256  _mm_store_ps(avgPointsVector, avgPointsVal);
257 
258  sumMean = 0.0;
259  sumMean += avgPointsVector[0];
260  sumMean += avgPointsVector[1];
261  sumMean += avgPointsVector[2];
262  sumMean += avgPointsVector[3];
263 
264  // Calculate the number of valid bins from the remaning count
265  __VOLK_ATTR_ALIGNED(16) float validBinCountVector[4];
266  _mm_store_ps(validBinCountVector, vValidBinCount);
267 
268  float validBinCount = 0;
269  validBinCount += validBinCountVector[0];
270  validBinCount += validBinCountVector[1];
271  validBinCount += validBinCountVector[2];
272  validBinCount += validBinCountVector[3];
273 
274  number = quarterPoints * 4;
275  for(;number < num_points; number++){
276  if(realDataPoints[number] <= meanAmplitude){
277  sumMean += realDataPoints[number];
278  validBinCount += 1.0;
279  }
280  }
281 
282  float localNoiseFloorAmplitude = 0;
283  if(validBinCount > 0.0){
284  localNoiseFloorAmplitude = sumMean / validBinCount;
285  }
286  else{
287  localNoiseFloorAmplitude = meanAmplitude; // For the odd case that all the amplitudes are equal...
288  }
289 
290  *noiseFloorAmplitude = localNoiseFloorAmplitude;
291 }
292 #endif /* LV_HAVE_SSE */
293 
294 
295 #ifdef LV_HAVE_GENERIC
296 
297 static inline void
299  const float* realDataPoints,
300  const float spectralExclusionValue,
301  const unsigned int num_points)
302 {
303  float sumMean = 0.0;
304  unsigned int number;
305  // find the sum (for mean), etc
306  for(number = 0; number < num_points; number++){
307  // sum (for mean)
308  sumMean += realDataPoints[number];
309  }
310 
311  // calculate the spectral mean
312  // +20 because for the comparison below we only want to throw out bins
313  // that are significantly higher (and would, thus, affect the mean more)
314  const float meanAmplitude = (sumMean / num_points) + spectralExclusionValue;
315 
316  // now throw out any bins higher than the mean
317  sumMean = 0.0;
318  unsigned int newNumDataPoints = num_points;
319  for(number = 0; number < num_points; number++){
320  if (realDataPoints[number] <= meanAmplitude)
321  sumMean += realDataPoints[number];
322  else
323  newNumDataPoints--;
324  }
325 
326  float localNoiseFloorAmplitude = 0.0;
327  if (newNumDataPoints == 0) // in the odd case that all
328  localNoiseFloorAmplitude = meanAmplitude; // amplitudes are equal!
329  else
330  localNoiseFloorAmplitude = sumMean / ((float)newNumDataPoints);
331 
332  *noiseFloorAmplitude = localNoiseFloorAmplitude;
333 }
334 #endif /* LV_HAVE_GENERIC */
335 
336 
337 #endif /* INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_a_H */
338 
339 #ifndef INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_u_H
340 #define INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_u_H
341 
342 #include <volk/volk_common.h>
343 #include <inttypes.h>
344 #include <stdio.h>
345 
346 #ifdef LV_HAVE_AVX
347 #include <immintrin.h>
348 
349 static inline void
351  const float* realDataPoints,
352  const float spectralExclusionValue,
353  const unsigned int num_points)
354 {
355  unsigned int number = 0;
356  const unsigned int eighthPoints = num_points / 8;
357 
358  const float* dataPointsPtr = realDataPoints;
359  __VOLK_ATTR_ALIGNED(16) float avgPointsVector[8];
360 
361  __m256 dataPointsVal;
362  __m256 avgPointsVal = _mm256_setzero_ps();
363  // Calculate the sum (for mean) for all points
364  for(; number < eighthPoints; number++){
365 
366  dataPointsVal = _mm256_loadu_ps(dataPointsPtr);
367 
368  dataPointsPtr += 8;
369 
370  avgPointsVal = _mm256_add_ps(avgPointsVal, dataPointsVal);
371  }
372 
373  _mm256_storeu_ps(avgPointsVector, avgPointsVal);
374 
375  float sumMean = 0.0;
376  sumMean += avgPointsVector[0];
377  sumMean += avgPointsVector[1];
378  sumMean += avgPointsVector[2];
379  sumMean += avgPointsVector[3];
380  sumMean += avgPointsVector[4];
381  sumMean += avgPointsVector[5];
382  sumMean += avgPointsVector[6];
383  sumMean += avgPointsVector[7];
384 
385  number = eighthPoints * 8;
386  for(;number < num_points; number++){
387  sumMean += realDataPoints[number];
388  }
389 
390  // calculate the spectral mean
391  // +20 because for the comparison below we only want to throw out bins
392  // that are significantly higher (and would, thus, affect the mean more
393  const float meanAmplitude = (sumMean / ((float)num_points)) + spectralExclusionValue;
394 
395  dataPointsPtr = realDataPoints; // Reset the dataPointsPtr
396  __m256 vMeanAmplitudeVector = _mm256_set1_ps(meanAmplitude);
397  __m256 vOnesVector = _mm256_set1_ps(1.0);
398  __m256 vValidBinCount = _mm256_setzero_ps();
399  avgPointsVal = _mm256_setzero_ps();
400  __m256 compareMask;
401  number = 0;
402  // Calculate the sum (for mean) for any points which do NOT exceed the mean amplitude
403  for(; number < eighthPoints; number++){
404 
405  dataPointsVal = _mm256_loadu_ps(dataPointsPtr);
406 
407  dataPointsPtr += 8;
408 
409  // Identify which items do not exceed the mean amplitude
410  compareMask = _mm256_cmp_ps(dataPointsVal, vMeanAmplitudeVector, 18);
411 
412  // Mask off the items that exceed the mean amplitude and add the avg Points that do not exceed the mean amplitude
413  avgPointsVal = _mm256_add_ps(avgPointsVal, _mm256_and_ps(compareMask, dataPointsVal));
414 
415  // Count the number of bins which do not exceed the mean amplitude
416  vValidBinCount = _mm256_add_ps(vValidBinCount, _mm256_and_ps(compareMask, vOnesVector));
417  }
418 
419  // Calculate the mean from the remaining data points
420  _mm256_storeu_ps(avgPointsVector, avgPointsVal);
421 
422  sumMean = 0.0;
423  sumMean += avgPointsVector[0];
424  sumMean += avgPointsVector[1];
425  sumMean += avgPointsVector[2];
426  sumMean += avgPointsVector[3];
427  sumMean += avgPointsVector[4];
428  sumMean += avgPointsVector[5];
429  sumMean += avgPointsVector[6];
430  sumMean += avgPointsVector[7];
431 
432  // Calculate the number of valid bins from the remaning count
433  __VOLK_ATTR_ALIGNED(16) float validBinCountVector[8];
434  _mm256_storeu_ps(validBinCountVector, vValidBinCount);
435 
436  float validBinCount = 0;
437  validBinCount += validBinCountVector[0];
438  validBinCount += validBinCountVector[1];
439  validBinCount += validBinCountVector[2];
440  validBinCount += validBinCountVector[3];
441  validBinCount += validBinCountVector[4];
442  validBinCount += validBinCountVector[5];
443  validBinCount += validBinCountVector[6];
444  validBinCount += validBinCountVector[7];
445 
446  number = eighthPoints * 8;
447  for(;number < num_points; number++){
448  if(realDataPoints[number] <= meanAmplitude){
449  sumMean += realDataPoints[number];
450  validBinCount += 1.0;
451  }
452  }
453 
454  float localNoiseFloorAmplitude = 0;
455  if(validBinCount > 0.0){
456  localNoiseFloorAmplitude = sumMean / validBinCount;
457  }
458  else{
459  localNoiseFloorAmplitude = meanAmplitude; // For the odd case that all the amplitudes are equal...
460  }
461 
462  *noiseFloorAmplitude = localNoiseFloorAmplitude;
463 }
464 #endif /* LV_HAVE_AVX */
465 #endif /* INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_u_H */
static void volk_32f_s32f_calc_spectral_noise_floor_32f_generic(float *noiseFloorAmplitude, const float *realDataPoints, const float spectralExclusionValue, const unsigned int num_points)
Definition: volk_32f_s32f_calc_spectral_noise_floor_32f.h:298
static void volk_32f_s32f_calc_spectral_noise_floor_32f_a_sse(float *noiseFloorAmplitude, const float *realDataPoints, const float spectralExclusionValue, const unsigned int num_points)
Definition: volk_32f_s32f_calc_spectral_noise_floor_32f.h:190
static void volk_32f_s32f_calc_spectral_noise_floor_32f_u_avx(float *noiseFloorAmplitude, const float *realDataPoints, const float spectralExclusionValue, const unsigned int num_points)
Definition: volk_32f_s32f_calc_spectral_noise_floor_32f.h:350
#define __VOLK_ATTR_ALIGNED(x)
Definition: volk_common.h:33
static void volk_32f_s32f_calc_spectral_noise_floor_32f_a_avx(float *noiseFloorAmplitude, const float *realDataPoints, const float spectralExclusionValue, const unsigned int num_points)
Definition: volk_32f_s32f_calc_spectral_noise_floor_32f.h:70