SDRAngel  4.11.5
Developer docs for <a href="https://github.com/f4exb/sdrangel">SDRangel<\a>, an Open Source Qt5 / OpenGL 3.0+ SDR and signal analyzer frontend to various hardware.
wfir.cpp
Go to the documentation of this file.
1 /*
2  July 15, 2015
3  Iowa Hills Software LLC
4  http://www.iowahills.com
5  */
6 
7 #include <math.h>
8 #include <new>
9 #include <iostream>
10 
11 #include "wfir.h"
12 
13 #undef M_PI
14 #define M_PI 3.14159265358979323846
15 #define M_2PI 6.28318530717958647692
16 
17 // This first calculates the impulse response for a rectangular window.
18 // It then applies the windowing function of choice to the impulse response.
19 void WFIR::BasicFIR(double *FirCoeff, int NumTaps, TPassTypeName PassType,
20  double OmegaC, double BW, TWindowType WindowType, double WinBeta)
21 {
22  int j;
23  double Arg, OmegaLow, OmegaHigh;
24 
25  switch (PassType)
26  {
27  case LPF:
28  for (j = 0; j < NumTaps; j++)
29  {
30  Arg = (double) j - (double) (NumTaps - 1) / 2.0;
31  FirCoeff[j] = OmegaC * Sinc(OmegaC * Arg * M_PI);
32  }
33  break;
34 
35  case HPF:
36  if (NumTaps % 2 == 1) // Odd tap counts
37  {
38  for (j = 0; j < NumTaps; j++)
39  {
40  Arg = (double) j - (double) (NumTaps - 1) / 2.0;
41  FirCoeff[j] = Sinc(Arg * M_PI)
42  - OmegaC * Sinc(OmegaC * Arg * M_PI);
43  }
44  }
45 
46  else // Even tap counts
47  {
48  for (j = 0; j < NumTaps; j++)
49  {
50  Arg = (double) j - (double) (NumTaps - 1) / 2.0;
51  if (Arg == 0.0)
52  FirCoeff[j] = 0.0;
53  else
54  FirCoeff[j] = cos(OmegaC * Arg * M_PI) / M_PI / Arg
55  + cos(Arg * M_PI);
56  }
57  }
58  break;
59 
60  case BPF:
61  OmegaLow = OmegaC - BW / 2.0;
62  OmegaHigh = OmegaC + BW / 2.0;
63  for (j = 0; j < NumTaps; j++)
64  {
65  Arg = (double) j - (double) (NumTaps - 1) / 2.0;
66  if (Arg == 0.0)
67  FirCoeff[j] = 0.0;
68  else
69  FirCoeff[j] = (cos(OmegaLow * Arg * M_PI)
70  - cos(OmegaHigh * Arg * M_PI)) / M_PI / Arg;
71  }
72  break;
73 
74  case NOTCH: // If NumTaps is even for Notch filters, the response at Pi is attenuated.
75  OmegaLow = OmegaC - BW / 2.0;
76  OmegaHigh = OmegaC + BW / 2.0;
77  for (j = 0; j < NumTaps; j++)
78  {
79  Arg = (double) j - (double) (NumTaps - 1) / 2.0;
80  FirCoeff[j] = Sinc(Arg * M_PI)
81  - OmegaHigh * Sinc(OmegaHigh * Arg * M_PI)
82  - OmegaLow * Sinc(OmegaLow * Arg * M_PI);
83  }
84  break;
85  }
86 
87  // WindowData can be used to window data before an FFT. When used for FIR filters we set
88  // Alpha = 0.0 to prevent a flat top on the window and
89  // set UnityGain = false to prevent the window gain from getting set to unity.
90  WindowData(FirCoeff, NumTaps, WindowType, 0.0, WinBeta, false);
91 
92 }
93 
94 //---------------------------------------------------------------------------
95 
96 // This gets used with the Kaiser window.
97 double WFIR::Bessel(double x)
98 {
99  double Sum = 0.0, XtoIpower;
100  int i, j, Factorial;
101  for (i = 1; i < 10; i++)
102  {
103  XtoIpower = pow(x / 2.0, (double) i);
104  Factorial = 1;
105  for (j = 1; j <= i; j++)
106  Factorial *= j;
107  Sum += pow(XtoIpower / (double) Factorial, 2.0);
108  }
109  return (1.0 + Sum);
110 }
111 
112 //-----------------------------------------------------------------------------
113 
114 // This gets used with the Sinc window and various places in the BasicFIR function.
115 double WFIR::Sinc(double x)
116 {
117  if (x > -1.0E-5 && x < 1.0E-5)
118  return (1.0);
119  return (sin(x) / x);
120 }
121 
122 //---------------------------------------------------------------------------
123 
124 // These are the various windows definitions. These windows can be used for either
125 // FIR filter design or with an FFT for spectral analysis.
126 // Sourced verbatim from: ~MyDocs\Code\Common\FFTFunctions.cpp
127 // For definitions, see this article: http://en.wikipedia.org/wiki/Window_function
128 
129 // This function has 6 inputs
130 // Data is the array, of length N, containing the data to to be windowed.
131 // This data is either a FIR filter sinc pulse, or the data to be analyzed by an fft.
132 
133 // WindowType is an enum defined in the header file, which is at the bottom of this file.
134 // e.g. wtKAISER, wtSINC, wtHANNING, wtHAMMING, wtBLACKMAN, ...
135 
136 // Alpha sets the width of the flat top.
137 // Windows such as the Tukey and Trapezoid are defined to have a variably wide flat top.
138 // As can be seen by its definition, the Tukey is just a Hanning window with a flat top.
139 // Alpha can be used to give any of these windows a partial flat top, except the Flattop and Kaiser.
140 // Alpha = 0 gives the original window. (i.e. no flat top)
141 // To generate a Tukey window, use a Hanning with 0 < Alpha < 1
142 // To generate a Bartlett window (triangular), use a Trapezoid window with Alpha = 0.
143 // Alpha = 1 generates a rectangular window in all cases. (except the Flattop and Kaiser)
144 
145 // Beta is used with the Kaiser, Sinc, and Sine windows only.
146 // These three windows are primarily used for FIR filter design, not spectral analysis.
147 // In FIR filter design, Beta controls the filter's transition bandwidth and the sidelobe levels.
148 // The code ignores Beta except in the Kaiser, Sinc, and Sine window cases.
149 
150 // UnityGain controls whether the gain of these windows is set to unity.
151 // Only the Flattop window has unity gain by design. The Hanning window, for example, has a gain
152 // of 1/2. UnityGain = true will set the gain of all these windows to 1.
153 // Then, when the window is applied to a signal, the signal's energy content is preserved.
154 // Don't use this with FIR filter design however. Since most of the enegy in an FIR sinc pulse
155 // is in the middle of the window, the window needs a peak amplitude of one, not unity gain.
156 // Setting UnityGain = true will simply cause the resulting FIR filter to have excess gain.
157 
158 // If using these windows for FIR filters, start with the Kaiser, Sinc, or Sine windows and
159 // adjust Beta for the desired transition BW and sidelobe levels (set Alpha = 0).
160 // While the FlatTop is an excellent window for spectral analysis, don't use it for FIR filter design.
161 // It has a peak amplitude of ~ 4.7 which causes the resulting FIR filter to have about this much gain.
162 // It works poorly for FIR filters even if you adjust its peak amplitude.
163 // The Trapezoid also works poorly for FIR filter design.
164 
165 // If using these windows with an fft for spectral analysis, start with the Hanning, Gauss, or Flattop.
166 // When choosing a window for spectral analysis, you must trade off between resolution and amplitude accuracy.
167 // The Hanning has the best resolution while the Flatop has the best amplitude accuracy.
168 // The Gauss is midway between these two for both accuracy and resolution.
169 // These three were the only windows available in the HP 89410A Vector Signal Analyzer. Which is to say,
170 // unless you have specific windowing requirements, use one of these 3 for general purpose signal analysis.
171 // Set UnityGain = true when using any of these windows for spectral analysis to preserve the signal's enegy level.
172 
173 void WFIR::WindowData(double *Data, int N, TWindowType WindowType, double Alpha,
174  double Beta, bool UnityGain)
175 {
176  if (WindowType == wtNONE)
177  return;
178 
179  int j, M, TopWidth;
180  double dM, *WinCoeff;
181 
182  if (WindowType == wtKAISER || WindowType == wtFLATTOP)
183  Alpha = 0.0;
184 
185  if (Alpha < 0.0)
186  Alpha = 0.0;
187  if (Alpha > 1.0)
188  Alpha = 1.0;
189 
190  if (Beta < 0.0)
191  Beta = 0.0;
192  if (Beta > 10.0)
193  Beta = 10.0;
194 
195  WinCoeff = new (std::nothrow) double[N + 2];
196  if (WinCoeff == 0)
197  {
198  std::cerr
199  << "Failed to allocate memory in FFTFunctions::WindowFFTData() "
200  << std::endl;
201  return;
202  }
203 
204  TopWidth = (int) (Alpha * (double) N);
205  if (TopWidth % 2 != 0)
206  TopWidth++;
207  if (TopWidth > N)
208  TopWidth = N;
209  M = N - TopWidth;
210  dM = M + 1;
211 
212  // Calculate the window for N/2 points, then fold the window over (at the bottom).
213  // TopWidth points will be set to 1.
214  if (WindowType == wtKAISER)
215  {
216  double Arg;
217  for (j = 0; j < M; j++)
218  {
219  Arg = Beta * sqrt(1.0 - pow(((double) (2 * j + 2) - dM) / dM, 2.0));
220  WinCoeff[j] = Bessel(Arg) / Bessel(Beta);
221  }
222  }
223 
224  else if (WindowType == wtSINC) // Lanczos
225  {
226  for (j = 0; j < M; j++)
227  WinCoeff[j] = Sinc((double) (2 * j + 1 - M) / dM * M_PI);
228  for (j = 0; j < M; j++)
229  WinCoeff[j] = pow(WinCoeff[j], Beta);
230  }
231 
232  else if (WindowType == wtSINE) // Hanning if Beta = 2
233  {
234  for (j = 0; j < M / 2; j++)
235  WinCoeff[j] = sin((double) (j + 1) * M_PI / dM);
236  for (j = 0; j < M / 2; j++)
237  WinCoeff[j] = pow(WinCoeff[j], Beta);
238  }
239 
240  else if (WindowType == wtHANNING)
241  {
242  for (j = 0; j < M / 2; j++)
243  WinCoeff[j] = 0.5 - 0.5 * cos((double) (j + 1) * M_2PI / dM);
244  }
245 
246  else if (WindowType == wtHAMMING)
247  {
248  for (j = 0; j < M / 2; j++)
249  WinCoeff[j] = 0.54 - 0.46 * cos((double) (j + 1) * M_2PI / dM);
250  }
251 
252  else if (WindowType == wtBLACKMAN)
253  {
254  for (j = 0; j < M / 2; j++)
255  {
256  WinCoeff[j] = 0.42 - 0.50 * cos((double) (j + 1) * M_2PI / dM)
257  + 0.08 * cos((double) (j + 1) * M_2PI * 2.0 / dM);
258  }
259  }
260 
261  // See: http://www.bth.se/fou/forskinfo.nsf/0/130c0940c5e7ffcdc1256f7f0065ac60/$file/ICOTA_2004_ttr_icl_mdh.pdf
262  else if (WindowType == wtFLATTOP)
263  {
264  for (j = 0; j <= M / 2; j++)
265  {
266  WinCoeff[j] = 1.0
267  - 1.93293488969227 * cos((double) (j + 1) * M_2PI / dM)
268  + 1.28349769674027
269  * cos((double) (j + 1) * M_2PI * 2.0 / dM)
270  - 0.38130801681619
271  * cos((double) (j + 1) * M_2PI * 3.0 / dM)
272  + 0.02929730258511
273  * cos((double) (j + 1) * M_2PI * 4.0 / dM);
274  }
275  }
276 
277  else if (WindowType == wtBLACKMAN_HARRIS)
278  {
279  for (j = 0; j < M / 2; j++)
280  {
281  WinCoeff[j] = 0.35875 - 0.48829 * cos((double) (j + 1) * M_2PI / dM)
282  + 0.14128 * cos((double) (j + 1) * M_2PI * 2.0 / dM)
283  - 0.01168 * cos((double) (j + 1) * M_2PI * 3.0 / dM);
284  }
285  }
286 
287  else if (WindowType == wtBLACKMAN_NUTTALL)
288  {
289  for (j = 0; j < M / 2; j++)
290  {
291  WinCoeff[j] = 0.3535819
292  - 0.4891775 * cos((double) (j + 1) * M_2PI / dM)
293  + 0.1365995 * cos((double) (j + 1) * M_2PI * 2.0 / dM)
294  - 0.0106411 * cos((double) (j + 1) * M_2PI * 3.0 / dM);
295  }
296  }
297 
298  else if (WindowType == wtNUTTALL)
299  {
300  for (j = 0; j < M / 2; j++)
301  {
302  WinCoeff[j] = 0.355768
303  - 0.487396 * cos((double) (j + 1) * M_2PI / dM)
304  + 0.144232 * cos((double) (j + 1) * M_2PI * 2.0 / dM)
305  - 0.012604 * cos((double) (j + 1) * M_2PI * 3.0 / dM);
306  }
307  }
308 
309  else if (WindowType == wtKAISER_BESSEL)
310  {
311  for (j = 0; j <= M / 2; j++)
312  {
313  WinCoeff[j] = 0.402 - 0.498 * cos(M_2PI * (double) (j + 1) / dM)
314  + 0.098 * cos(2.0 * M_2PI * (double) (j + 1) / dM)
315  + 0.001 * cos(3.0 * M_2PI * (double) (j + 1) / dM);
316  }
317  }
318 
319  else if (WindowType == wtTRAPEZOID) // Rectangle for Alpha = 1 Triangle for Alpha = 0
320  {
321  int K = M / 2;
322  if (M % 2)
323  K++;
324  for (j = 0; j < K; j++)
325  WinCoeff[j] = (double) (j + 1) / (double) K;
326  }
327 
328  // This definition is from http://en.wikipedia.org/wiki/Window_function (Gauss Generalized normal window)
329  // We set their p = 2, and use Alpha in the numerator, instead of Sigma in the denominator, as most others do.
330  // Alpha = 2.718 puts the Gauss window response midway between the Hanning and the Flattop (basically what we want).
331  // It also gives the same BW as the Gauss window used in the HP 89410A Vector Signal Analyzer.
332  // Alpha = 1.8 puts it quite close to the Hanning.
333  else if (WindowType == wtGAUSS)
334  {
335  for (j = 0; j < M / 2; j++)
336  {
337  WinCoeff[j] = ((double) (j + 1) - dM / 2.0) / (dM / 2.0) * 2.7183;
338  WinCoeff[j] *= WinCoeff[j];
339  WinCoeff[j] = exp(-WinCoeff[j]);
340  }
341  }
342 
343  else // Error.
344  {
345  std::cerr << "Incorrect window type in WindowFFTData" << std::endl;
346  delete[] WinCoeff;
347  return;
348  }
349 
350  // Fold the coefficients over.
351  for (j = 0; j < M / 2; j++)
352  WinCoeff[N - j - 1] = WinCoeff[j];
353 
354  // This is the flat top if Alpha > 0. Cannot be applied to a Kaiser or Flat Top.
355  if (WindowType != wtKAISER && WindowType != wtFLATTOP)
356  {
357  for (j = M / 2; j < N - M / 2; j++)
358  WinCoeff[j] = 1.0;
359  }
360 
361  // This will set the gain of the window to 1. Only the Flattop window has unity gain by design.
362  if (UnityGain)
363  {
364  double Sum = 0.0;
365  for (j = 0; j < N; j++)
366  Sum += WinCoeff[j];
367  Sum /= (double) N;
368  if (Sum != 0.0)
369  for (j = 0; j < N; j++)
370  WinCoeff[j] /= Sum;
371  }
372 
373  // Apply the window to the data.
374  for (j = 0; j < N; j++)
375  Data[j] *= WinCoeff[j];
376 
377  delete[] WinCoeff;
378 
379 }
#define M_2PI
Definition: wfir.cpp:15
Fixed< IntType, IntBits > cos(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2271
#define M_PI
Definition: wfir.cpp:14
TWindowType
Definition: wfir.h:66
TPassTypeName
Definition: wfir.h:61
static void WindowData(double *Data, int N, TWindowType WindowType, double Alpha, double Beta, bool UnityGain)
Definition: wfir.cpp:173
Fixed< IntType, IntBits > exp(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2289
Definition: wfir.h:63
Definition: wfir.h:63
Definition: wfir.h:63
Fixed< IntType, IntBits > sin(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2265
int32_t i
Definition: decimators.h:244
static void BasicFIR(double *FirCoeff, int NumTaps, TPassTypeName PassType, double OmegaC, double BW, TWindowType WindowType, double WinBeta)
Definition: wfir.cpp:19
Fixed< IntType, IntBits > sqrt(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2283
static double Bessel(double x)
Definition: wfir.cpp:97
static double Sinc(double x)
Definition: wfir.cpp:115