GNU Radio Manual and C++ API Reference  3.7.7
The Free & Open Software Radio Ecosystem
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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 
23 /*!
24  * \page volk_32f_s32f_calc_spectral_noise_floor_32f
25  *
26  * \b Overview
27  *
28  * Computes the spectral noise floor of an input power spectrum.
29  *
30  * Calculates the spectral noise floor of an input power spectrum by
31  * determining the mean of the input power spectrum, then
32  * recalculating the mean excluding any power spectrum values that
33  * exceed the mean by the spectralExclusionValue (in dB). Provides a
34  * rough estimation of the signal noise floor.
35  *
36  * <b>Dispatcher Prototype</b>
37  * \code
38  * void volk_32f_s32f_calc_spectral_noise_floor_32f(float* noiseFloorAmplitude, const float* realDataPoints, const float spectralExclusionValue, const unsigned int num_points)
39  * \endcode
40  *
41  * \b Inputs
42  * \li realDataPoints: The input power spectrum.
43  * \li spectralExclusionValue: The number of dB above the noise floor that a data point must be to be excluded from the noise floor calculation - default value is 20.
44  * \li num_points: The number of data points.
45  *
46  * \b Outputs
47  * \li noiseFloorAmplitude: The noise floor of the input spectrum, in dB.
48  *
49  * \b Example
50  * \code
51  * int N = 10000;
52  *
53  * volk_32f_s32f_calc_spectral_noise_floor_32f
54  *
55  * volk_free(x);
56  * \endcode
57  */
58 
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_SSE
67 #include <xmmintrin.h>
68 
69 static inline void
70 volk_32f_s32f_calc_spectral_noise_floor_32f_a_sse(float* noiseFloorAmplitude,
71  const float* realDataPoints,
72  const float spectralExclusionValue,
73  const unsigned int num_points)
74 {
75  unsigned int number = 0;
76  const unsigned int quarterPoints = num_points / 4;
77 
78  const float* dataPointsPtr = realDataPoints;
79  __VOLK_ATTR_ALIGNED(16) float avgPointsVector[4];
80 
81  __m128 dataPointsVal;
82  __m128 avgPointsVal = _mm_setzero_ps();
83  // Calculate the sum (for mean) for all points
84  for(; number < quarterPoints; number++){
85 
86  dataPointsVal = _mm_load_ps(dataPointsPtr);
87 
88  dataPointsPtr += 4;
89 
90  avgPointsVal = _mm_add_ps(avgPointsVal, dataPointsVal);
91  }
92 
93  _mm_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 
101  number = quarterPoints * 4;
102  for(;number < num_points; number++){
103  sumMean += realDataPoints[number];
104  }
105 
106  // calculate the spectral mean
107  // +20 because for the comparison below we only want to throw out bins
108  // that are significantly higher (and would, thus, affect the mean more
109  const float meanAmplitude = (sumMean / ((float)num_points)) + spectralExclusionValue;
110 
111  dataPointsPtr = realDataPoints; // Reset the dataPointsPtr
112  __m128 vMeanAmplitudeVector = _mm_set_ps1(meanAmplitude);
113  __m128 vOnesVector = _mm_set_ps1(1.0);
114  __m128 vValidBinCount = _mm_setzero_ps();
115  avgPointsVal = _mm_setzero_ps();
116  __m128 compareMask;
117  number = 0;
118  // Calculate the sum (for mean) for any points which do NOT exceed the mean amplitude
119  for(; number < quarterPoints; number++){
120 
121  dataPointsVal = _mm_load_ps(dataPointsPtr);
122 
123  dataPointsPtr += 4;
124 
125  // Identify which items do not exceed the mean amplitude
126  compareMask = _mm_cmple_ps(dataPointsVal, vMeanAmplitudeVector);
127 
128  // Mask off the items that exceed the mean amplitude and add the avg Points that do not exceed the mean amplitude
129  avgPointsVal = _mm_add_ps(avgPointsVal, _mm_and_ps(compareMask, dataPointsVal));
130 
131  // Count the number of bins which do not exceed the mean amplitude
132  vValidBinCount = _mm_add_ps(vValidBinCount, _mm_and_ps(compareMask, vOnesVector));
133  }
134 
135  // Calculate the mean from the remaining data points
136  _mm_store_ps(avgPointsVector, avgPointsVal);
137 
138  sumMean = 0.0;
139  sumMean += avgPointsVector[0];
140  sumMean += avgPointsVector[1];
141  sumMean += avgPointsVector[2];
142  sumMean += avgPointsVector[3];
143 
144  // Calculate the number of valid bins from the remaning count
145  __VOLK_ATTR_ALIGNED(16) float validBinCountVector[4];
146  _mm_store_ps(validBinCountVector, vValidBinCount);
147 
148  float validBinCount = 0;
149  validBinCount += validBinCountVector[0];
150  validBinCount += validBinCountVector[1];
151  validBinCount += validBinCountVector[2];
152  validBinCount += validBinCountVector[3];
153 
154  number = quarterPoints * 4;
155  for(;number < num_points; number++){
156  if(realDataPoints[number] <= meanAmplitude){
157  sumMean += realDataPoints[number];
158  validBinCount += 1.0;
159  }
160  }
161 
162  float localNoiseFloorAmplitude = 0;
163  if(validBinCount > 0.0){
164  localNoiseFloorAmplitude = sumMean / validBinCount;
165  }
166  else{
167  localNoiseFloorAmplitude = meanAmplitude; // For the odd case that all the amplitudes are equal...
168  }
169 
170  *noiseFloorAmplitude = localNoiseFloorAmplitude;
171 }
172 #endif /* LV_HAVE_SSE */
173 
174 
175 #ifdef LV_HAVE_GENERIC
176 
177 static inline void
178 volk_32f_s32f_calc_spectral_noise_floor_32f_generic(float* noiseFloorAmplitude,
179  const float* realDataPoints,
180  const float spectralExclusionValue,
181  const unsigned int num_points)
182 {
183  float sumMean = 0.0;
184  unsigned int number;
185  // find the sum (for mean), etc
186  for(number = 0; number < num_points; number++){
187  // sum (for mean)
188  sumMean += realDataPoints[number];
189  }
190 
191  // calculate the spectral mean
192  // +20 because for the comparison below we only want to throw out bins
193  // that are significantly higher (and would, thus, affect the mean more)
194  const float meanAmplitude = (sumMean / num_points) + spectralExclusionValue;
195 
196  // now throw out any bins higher than the mean
197  sumMean = 0.0;
198  unsigned int newNumDataPoints = num_points;
199  for(number = 0; number < num_points; number++){
200  if (realDataPoints[number] <= meanAmplitude)
201  sumMean += realDataPoints[number];
202  else
203  newNumDataPoints--;
204  }
205 
206  float localNoiseFloorAmplitude = 0.0;
207  if (newNumDataPoints == 0) // in the odd case that all
208  localNoiseFloorAmplitude = meanAmplitude; // amplitudes are equal!
209  else
210  localNoiseFloorAmplitude = sumMean / ((float)newNumDataPoints);
211 
212  *noiseFloorAmplitude = localNoiseFloorAmplitude;
213 }
214 #endif /* LV_HAVE_GENERIC */
215 
216 
217 #endif /* INCLUDED_volk_32f_s32f_calc_spectral_noise_floor_32f_a_H */
#define __VOLK_ATTR_ALIGNED(x)
Definition: volk_common.h:27