1 /* -*- c++ -*- */
2 /*
3 * Copyright 2002,2007,2008,2012,2013,2018 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 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26
27 #include <gnuradio/fft/window.h>
28 #include <gnuradio/math.h>
29 #include <stdexcept>
30
31 namespace gr {
32 namespace fft {
33
34 #define IzeroEPSILON 1E-21 /* Max error acceptable in Izero */
35
Izero(double x)36 static double Izero(double x)
37 {
38 double sum, u, halfx, temp;
39 int n;
40
41 sum = u = n = 1;
42 halfx = x / 2.0;
43 do {
44 temp = halfx / (double)n;
45 n += 1;
46 temp *= temp;
47 u *= temp;
48 sum += u;
49 } while (u >= IzeroEPSILON * sum);
50 return (sum);
51 }
52
midn(int ntaps)53 double midn(int ntaps) { return ntaps / 2.0; }
54
midm1(int ntaps)55 double midm1(int ntaps) { return (ntaps - 1.0) / 2.0; }
56
midp1(int ntaps)57 double midp1(int ntaps) { return (ntaps + 1.0) / 2.0; }
58
freq(int ntaps)59 double freq(int ntaps) { return 2.0 * GR_M_PI / ntaps; }
60
rate(int ntaps)61 double rate(int ntaps) { return 1.0 / (ntaps >> 1); }
62
max_attenuation(win_type type,double beta)63 double window::max_attenuation(win_type type, double beta)
64 {
65 switch (type) {
66 case (WIN_HAMMING):
67 return 53;
68 break;
69 case (WIN_HANN):
70 return 44;
71 break;
72 case (WIN_BLACKMAN):
73 return 74;
74 break;
75 case (WIN_RECTANGULAR):
76 return 21;
77 break;
78 case (WIN_KAISER):
79 return (beta / 0.1102 + 8.7);
80 break;
81 case (WIN_BLACKMAN_hARRIS):
82 return 92;
83 break;
84 case (WIN_BARTLETT):
85 return 27;
86 break;
87 case (WIN_FLATTOP):
88 return 93;
89 break;
90 default:
91 throw std::out_of_range("window::max_attenuation: unknown window type provided.");
92 }
93 }
94
coswindow(int ntaps,float c0,float c1,float c2)95 std::vector<float> window::coswindow(int ntaps, float c0, float c1, float c2)
96 {
97 std::vector<float> taps(ntaps);
98 float M = static_cast<float>(ntaps - 1);
99
100 for (int n = 0; n < ntaps; n++)
101 taps[n] = c0 - c1 * cosf((2.0f * GR_M_PI * n) / M) +
102 c2 * cosf((4.0f * GR_M_PI * n) / M);
103 return taps;
104 }
105
coswindow(int ntaps,float c0,float c1,float c2,float c3)106 std::vector<float> window::coswindow(int ntaps, float c0, float c1, float c2, float c3)
107 {
108 std::vector<float> taps(ntaps);
109 float M = static_cast<float>(ntaps - 1);
110
111 for (int n = 0; n < ntaps; n++)
112 taps[n] = c0 - c1 * cosf((2.0f * GR_M_PI * n) / M) +
113 c2 * cosf((4.0f * GR_M_PI * n) / M) -
114 c3 * cosf((6.0f * GR_M_PI * n) / M);
115 return taps;
116 }
117
118 std::vector<float>
coswindow(int ntaps,float c0,float c1,float c2,float c3,float c4)119 window::coswindow(int ntaps, float c0, float c1, float c2, float c3, float c4)
120 {
121 std::vector<float> taps(ntaps);
122 float M = static_cast<float>(ntaps - 1);
123
124 for (int n = 0; n < ntaps; n++)
125 taps[n] = c0 - c1 * cosf((2.0f * GR_M_PI * n) / M) +
126 c2 * cosf((4.0f * GR_M_PI * n) / M) -
127 c3 * cosf((6.0f * GR_M_PI * n) / M) +
128 c4 * cosf((8.0f * GR_M_PI * n) / M);
129 return taps;
130 }
131
rectangular(int ntaps)132 std::vector<float> window::rectangular(int ntaps)
133 {
134 std::vector<float> taps(ntaps);
135 for (int n = 0; n < ntaps; n++)
136 taps[n] = 1;
137 return taps;
138 }
139
hamming(int ntaps)140 std::vector<float> window::hamming(int ntaps)
141 {
142 std::vector<float> taps(ntaps);
143 float M = static_cast<float>(ntaps - 1);
144
145 for (int n = 0; n < ntaps; n++)
146 taps[n] = 0.54 - 0.46 * cos((2 * GR_M_PI * n) / M);
147 return taps;
148 }
149
hann(int ntaps)150 std::vector<float> window::hann(int ntaps)
151 {
152 std::vector<float> taps(ntaps);
153 float M = static_cast<float>(ntaps - 1);
154
155 for (int n = 0; n < ntaps; n++)
156 taps[n] = 0.5 - 0.5 * cos((2 * GR_M_PI * n) / M);
157 return taps;
158 }
159
hanning(int ntaps)160 std::vector<float> window::hanning(int ntaps) { return hann(ntaps); }
161
blackman(int ntaps)162 std::vector<float> window::blackman(int ntaps)
163 {
164 return coswindow(ntaps, 0.42, 0.5, 0.08);
165 }
166
blackman2(int ntaps)167 std::vector<float> window::blackman2(int ntaps)
168 {
169 return coswindow(ntaps, 0.34401, 0.49755, 0.15844);
170 }
171
blackman3(int ntaps)172 std::vector<float> window::blackman3(int ntaps)
173 {
174 return coswindow(ntaps, 0.21747, 0.45325, 0.28256, 0.04672);
175 }
176
blackman4(int ntaps)177 std::vector<float> window::blackman4(int ntaps)
178 {
179 return coswindow(ntaps, 0.084037, 0.29145, 0.375696, 0.20762, 0.041194);
180 }
181
blackman_harris(int ntaps,int atten)182 std::vector<float> window::blackman_harris(int ntaps, int atten)
183 {
184 switch (atten) {
185 case (61):
186 return coswindow(ntaps, 0.42323, 0.49755, 0.07922);
187 case (67):
188 return coswindow(ntaps, 0.44959, 0.49364, 0.05677);
189 case (74):
190 return coswindow(ntaps, 0.40271, 0.49703, 0.09392, 0.00183);
191 case (92):
192 return coswindow(ntaps, 0.35875, 0.48829, 0.14128, 0.01168);
193 default:
194 throw std::out_of_range("window::blackman_harris: unknown attenuation value "
195 "(must be 61, 67, 74, or 92)");
196 }
197 }
198
blackmanharris(int ntaps,int atten)199 std::vector<float> window::blackmanharris(int ntaps, int atten)
200 {
201 return blackman_harris(ntaps, atten);
202 }
203
nuttal(int ntaps)204 std::vector<float> window::nuttal(int ntaps) { return nuttall(ntaps); }
205
nuttall(int ntaps)206 std::vector<float> window::nuttall(int ntaps)
207 {
208 return coswindow(ntaps, 0.3635819, 0.4891775, 0.1365995, 0.0106411);
209 }
210
blackman_nuttal(int ntaps)211 std::vector<float> window::blackman_nuttal(int ntaps) { return nuttall(ntaps); }
212
blackman_nuttall(int ntaps)213 std::vector<float> window::blackman_nuttall(int ntaps) { return nuttall(ntaps); }
214
nuttal_cfd(int ntaps)215 std::vector<float> window::nuttal_cfd(int ntaps) { return nuttall_cfd(ntaps); }
216
nuttall_cfd(int ntaps)217 std::vector<float> window::nuttall_cfd(int ntaps)
218 {
219 return coswindow(ntaps, 0.355768, 0.487396, 0.144232, 0.012604);
220 }
221
flattop(int ntaps)222 std::vector<float> window::flattop(int ntaps)
223 {
224 double scale = 4.63867;
225 return coswindow(
226 ntaps, 1.0 / scale, 1.93 / scale, 1.29 / scale, 0.388 / scale, 0.028 / scale);
227 }
228
kaiser(int ntaps,double beta)229 std::vector<float> window::kaiser(int ntaps, double beta)
230 {
231 if (beta < 0)
232 throw std::out_of_range("window::kaiser: beta must be >= 0");
233
234 std::vector<float> taps(ntaps);
235
236 double IBeta = 1.0 / Izero(beta);
237 double inm1 = 1.0 / ((double)(ntaps - 1));
238 double temp;
239
240 /* extracting first and last element out of the loop, since
241 sqrt(1.0-temp*temp) might trigger unexpected floating point behaviour
242 if |temp| = 1.0+epsilon, which can happen for i==0 and
243 1/i==1/(ntaps-1)==inm1 ; compare
244 https://github.com/gnuradio/gnuradio/issues/1348 .
245 In any case, the 0. Bessel function of first kind is 1 at point 0.
246 */
247 taps[0] = IBeta;
248 for (int i = 1; i < ntaps - 1; i++) {
249 temp = 2 * i * inm1 - 1;
250 taps[i] = Izero(beta * sqrt(1.0 - temp * temp)) * IBeta;
251 }
252 taps[ntaps - 1] = IBeta;
253 return taps;
254 }
255
bartlett(int ntaps)256 std::vector<float> window::bartlett(int ntaps)
257 {
258 std::vector<float> taps(ntaps);
259 float M = static_cast<float>(ntaps - 1);
260
261 for (int n = 0; n < ntaps / 2; n++)
262 taps[n] = 2 * n / M;
263 for (int n = ntaps / 2; n < ntaps; n++)
264 taps[n] = 2 - 2 * n / M;
265
266 return taps;
267 }
268
welch(int ntaps)269 std::vector<float> window::welch(int ntaps)
270 {
271 std::vector<float> taps(ntaps);
272 double m1 = midm1(ntaps);
273 double p1 = midp1(ntaps);
274 for (int i = 0; i < midn(ntaps) + 1; i++) {
275 taps[i] = 1.0 - pow((i - m1) / p1, 2);
276 taps[ntaps - i - 1] = taps[i];
277 }
278 return taps;
279 }
280
parzen(int ntaps)281 std::vector<float> window::parzen(int ntaps)
282 {
283 std::vector<float> taps(ntaps);
284 double m1 = midm1(ntaps);
285 double m = midn(ntaps);
286 int i;
287 for (i = ntaps / 4; i < 3 * ntaps / 4; i++) {
288 taps[i] = 1.0 - 6.0 * pow((i - m1) / m, 2.0) * (1.0 - fabs(i - m1) / m);
289 }
290 for (; i < ntaps; i++) {
291 taps[i] = 2.0 * pow(1.0 - fabs(i - m1) / m, 3.0);
292 taps[ntaps - i - 1] = taps[i];
293 }
294 return taps;
295 }
296
exponential(int ntaps,double d)297 std::vector<float> window::exponential(int ntaps, double d)
298 {
299 if (d < 0)
300 throw std::out_of_range("window::exponential: d must be >= 0");
301
302 double m1 = midm1(ntaps);
303 double p = 1.0 / (8.69 * ntaps / (2.0 * d));
304 std::vector<float> taps(ntaps);
305 for (int i = 0; i < midn(ntaps) + 1; i++) {
306 taps[i] = exp(-fabs(i - m1) * p);
307 taps[ntaps - i - 1] = taps[i];
308 }
309 return taps;
310 }
311
riemann(int ntaps)312 std::vector<float> window::riemann(int ntaps)
313 {
314 double cx;
315 double sr1 = freq(ntaps);
316 double mid = midn(ntaps);
317 std::vector<float> taps(ntaps);
318 for (int i = 0; i < mid; i++) {
319 if (i == midn(ntaps)) {
320 taps[i] = 1.0;
321 taps[ntaps - i - 1] = 1.0;
322 } else {
323 cx = sr1 * (i - mid);
324 taps[i] = sin(cx) / cx;
325 taps[ntaps - i - 1] = sin(cx) / cx;
326 }
327 }
328 return taps;
329 }
330
build(win_type type,int ntaps,double beta)331 std::vector<float> window::build(win_type type, int ntaps, double beta)
332 {
333 switch (type) {
334 case WIN_RECTANGULAR:
335 return rectangular(ntaps);
336 case WIN_HAMMING:
337 return hamming(ntaps);
338 case WIN_HANN:
339 return hann(ntaps);
340 case WIN_BLACKMAN:
341 return blackman(ntaps);
342 case WIN_BLACKMAN_hARRIS:
343 return blackman_harris(ntaps);
344 case WIN_KAISER:
345 return kaiser(ntaps, beta);
346 case WIN_BARTLETT:
347 return bartlett(ntaps);
348 case WIN_FLATTOP:
349 return flattop(ntaps);
350 default:
351 throw std::out_of_range("window::build: type out of range");
352 }
353 }
354
355 } /* namespace fft */
356 } /* namespace gr */
357