Line data Source code
1 : /** \file signalWindows.hpp
2 : * \brief Procedures to calculate window functions for signal processing
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup signal_processing_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2015 - 2020 Jared R. Males (jaredmales@gmail.com)
12 : //
13 : // This file is part of mxlib.
14 : //
15 : // mxlib is free software: you can redistribute it and/or modify
16 : // it under the terms of the GNU General Public License as published by
17 : // the Free Software Foundation, either version 3 of the License, or
18 : // (at your option) any later version.
19 : //
20 : // mxlib is distributed in the hope that it will be useful,
21 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 : // GNU General Public License for more details.
24 : //
25 : // You should have received a copy of the GNU General Public License
26 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
27 : //***********************************************************************//
28 :
29 : #ifndef signalWindows_hpp
30 : #define signalWindows_hpp
31 :
32 : #include <cmath>
33 : #include "../math/constants.hpp"
34 :
35 : namespace mx
36 : {
37 : namespace sigproc
38 : {
39 : namespace window
40 : {
41 :
42 : /// The Tukey Window
43 : /**
44 : * The width of the window is controlled by alpha. alpha = 0 is a square wave, alpha=1.0 is the Hann window.
45 : *
46 : * See https://en.wikipedia.org/wiki/Window_function
47 : *
48 : * \tparam realT a floating point type
49 : *
50 : * \ingroup signal_windows1D
51 : */
52 : template <typename realT>
53 2 : void tukey( realT *filt, ///< [out] the pre-allocated array to hold the filter
54 : int N, ///< [in] the size of the filter
55 : realT alpha ///< [in] the width parameter
56 : )
57 : {
58 2 : constexpr realT pi = math::pi<realT>();
59 :
60 2 : realT lim1 = alpha * ( N - 1.0 ) / 2.0;
61 2 : realT lim2 = ( N - 1.0 ) * ( 1.0 - 0.5 * alpha );
62 :
63 514 : for( int ii = 0; ii < N; ++ii )
64 : {
65 :
66 512 : if( ii < lim1 && alpha > 0. )
67 : {
68 192 : filt[ii] = 0.5 * ( 1.0 + cos( pi * ( 2. * ( ii ) / ( alpha * ( N - 1 ) ) - 1.0 ) ) );
69 : }
70 320 : else if( ii > lim2 && alpha > 0. )
71 : {
72 192 : filt[ii] = 0.5 * ( 1.0 + cos( pi * ( 2. * ( ii ) / ( alpha * ( N - 1 ) ) - 2. / alpha + 1.0 ) ) );
73 : }
74 : else
75 : {
76 128 : filt[ii] = 1.0;
77 : }
78 : }
79 2 : }
80 :
81 : /// The Tukey Window
82 : /**
83 : * The width of the window is controlled by alpha. alpha = 0 is a square wave, alpha=1.0 is the Hann window.
84 : *
85 : * See https://en.wikipedia.org/wiki/Window_function
86 : *
87 : * \overload
88 : *
89 : * \tparam realT a floating point type
90 : *
91 : * \ingroup signal_windows1D
92 : */
93 : template <typename realT>
94 2 : void tukey( std::vector<realT> &filt, ///< [out] the pre-allocated vector to hold the filter
95 : realT alpha ///< [in] the width parameter
96 : )
97 : {
98 2 : tukey( filt.data(), filt.size(), alpha );
99 2 : }
100 :
101 : /// The generalized 2 parameter cosine Window
102 : /** See https://en.wikipedia.org/wiki/Window_function
103 : *
104 : * \tparam realT a real floating point type
105 : *
106 : * \ingroup signal_windows1D
107 : */
108 : template <typename realT>
109 : void genCosine( realT *filt, ///< [out] The pre-allocated vector which will store the filter
110 : size_t N, ///< [in] the size of the filter vector
111 : realT a0, ///< [in] parameter of the generalized cosine window
112 : realT a1 ///< [in] parameter of the generalized cosine window
113 : )
114 : {
115 : constexpr realT pi = math::pi<realT>();
116 :
117 : for( size_t n = 0; n < N; ++n )
118 : {
119 : filt[n] = a0 - a1 * cos( 2 * pi * n / N );
120 : }
121 : }
122 :
123 : /// The generalized 3 parameter cosine Window
124 : /** See https://en.wikipedia.org/wiki/Window_function
125 : *
126 : * \tparam realT a real floating point type
127 : *
128 : * \ingroup signal_windows1D
129 : */
130 : template <typename realT>
131 : void genCosine( realT *filt, ///< [out] The pre-allocated vector which will store the filter
132 : size_t N, ///< [in] the size of the filter vector
133 : realT a0, ///< [in] parameter of the generalized cosine window
134 : realT a1, ///< [in] parameter of the generalized cosine window
135 : realT a2 ///< [in] parameter of the generalized cosine window
136 : )
137 : {
138 : constexpr realT pi = math::pi<realT>();
139 :
140 : for( size_t n = 0; n < N; ++n )
141 : {
142 : filt[n] = a0 - a1 * cos( 2 * pi * n / N ) + a2 * cos( 4 * pi * n / N );
143 : }
144 : }
145 :
146 : /// The generalized 4 parameter cosine Window
147 : /** See https://en.wikipedia.org/wiki/Window_function
148 : *
149 : * \tparam realT a real floating point type
150 : *
151 : * \ingroup signal_windows1D
152 : */
153 : template <typename realT>
154 : void genCosine( realT *filt, ///< [out] The pre-allocated vector which will store the filter
155 : size_t N, ///< [in] the size of the filter vector
156 : realT a0, ///< [in] parameter of the generalized cosine window
157 : realT a1, ///< [in] parameter of the generalized cosine window
158 : realT a2, ///< [in] parameter of the generalized cosine window
159 : realT a3 ///< [in] parameter of the generalized cosine window
160 : )
161 : {
162 : constexpr realT pi = math::pi<realT>();
163 :
164 : for( size_t n = 0; n < N; ++n )
165 : {
166 : filt[n] = a0 - a1 * cos( 2 * pi * n / N ) + a2 * cos( 4 * pi * n / N ) - a3 * cos( 6 * pi * n / N );
167 : }
168 : }
169 :
170 : /// The Hann Window
171 : /**
172 : * See https://en.wikipedia.org/wiki/Window_function
173 : *
174 : * \tparam realT a floating point type
175 : *
176 : * \ingroup signal_windows1D
177 : */
178 : template <typename realT>
179 : void hann( realT *filt, ///< [out] the pre-allocated array to hold the filter
180 : int N ///< [in] the size of the filter
181 : )
182 : {
183 : genCosine<realT>( filt, N, 0.5, 0.5 );
184 : }
185 :
186 : /// The Hann Window
187 : /**
188 : * See https://en.wikipedia.org/wiki/Window_function
189 : *
190 : * \overload
191 : *
192 : * \tparam realT a floating point type
193 : *
194 : * \ingroup signal_windows1D
195 : */
196 : template <typename realT>
197 : void hann( std::vector<realT> &filt /**< [out] the pre-allocated vector to hold the filter */ )
198 : {
199 : hann( filt.data(), filt.size() );
200 : }
201 :
202 : /// The Hamming Window
203 : /**
204 : * See https://en.wikipedia.org/wiki/Window_function
205 : *
206 : * \tparam realT a floating point type
207 : *
208 : * \ingroup signal_windows1D
209 : */
210 : template <typename realT>
211 : void hamming( realT *filt, ///< [out] the pre-allocated array to hold the filter
212 : int N ///< [in] the size of the filter
213 : )
214 : {
215 : genCosine<realT>( filt, N, 25.0 / 46.0, 1.0 - 25.0 / 46.0 );
216 : }
217 :
218 : /// The Hamming Window
219 : /**
220 : * See https://en.wikipedia.org/wiki/Window_function
221 : *
222 : * \overload
223 : *
224 : * \tparam realT a floating point type
225 : *
226 : * \ingroup signal_windows1D
227 : */
228 : template <typename realT>
229 : void hamming( std::vector<realT> &filt /**< [out] the pre-allocated vector to hold the filter */ )
230 : {
231 : hamming( filt.data(), filt.size() );
232 : }
233 :
234 : /// The Blackman Window
235 : /** See https://en.wikipedia.org/wiki/Window_function
236 : *
237 : * \tparam realT a real floating point type
238 : *
239 : * \ingroup signal_windows1D
240 : */
241 : template <typename realT>
242 : void blackman( realT *filt, ///< [out] The pre-allocated vector which will store the filter
243 : size_t N ///< [in] the size of the filter vector
244 : )
245 : {
246 : genCosine<realT>( filt, N, 0.42, 0.5, 0.08 );
247 : }
248 :
249 : /// The Blackman Window
250 : /** See https://en.wikipedia.org/wiki/Window_function
251 : *
252 : * \overload
253 : *
254 : * \tparam realT a real floating point type
255 : *
256 : * \ingroup signal_windows1D
257 : */
258 : template <typename realT>
259 : void blackman( std::vector<realT> &filt /**< [out] The pre-allocated vector which will store the filter */ )
260 : {
261 : blackman( filt.data(), filt.size() );
262 : }
263 :
264 : /// The Exact Blackman Window
265 : /** See https://en.wikipedia.org/wiki/Window_function
266 : *
267 : * \tparam realT a real floating point type
268 : *
269 : * \ingroup signal_windows1D
270 : */
271 : template <typename realT>
272 : void exactBlackman( realT *filt, ///< [out] The pre-allocated vector which will store the filter
273 : size_t N ///< [in] the size of the filter vector
274 : )
275 : {
276 : genCosine<realT>( filt, N, 0.42659, 0.49656, 0.076849 );
277 : }
278 :
279 : /// The Exact Blackman Windwo
280 : /** See https://en.wikipedia.org/wiki/Window_function
281 : *
282 : * \overload
283 : *
284 : * \tparam realT a real floating point type
285 : */
286 : template <typename realT>
287 : void exactBlackman( std::vector<realT> &filt /**< [out] The pre-allocated vector which will store the filter */ )
288 : {
289 : exactBlackman( filt.data(), filt.size() );
290 : }
291 :
292 : /// The Nuttal Window
293 : /** See https://en.wikipedia.org/wiki/Window_function
294 : *
295 : * \tparam realT a real floating point type
296 : *
297 : * \ingroup signal_windows1D
298 : */
299 : template <typename realT>
300 : void nuttal( realT *filt, ///< [out] The pre-allocated vector which will store the filter
301 : size_t N ///< [in] the size of the filter vector
302 : )
303 : {
304 : genCosine<realT>( filt, N, 0.355768, 0.487396, 0.144232, 0.012604 );
305 : }
306 :
307 : /// The Nuttal Window
308 : /** See https://en.wikipedia.org/wiki/Window_function
309 : *
310 : * \overload
311 : *
312 : * \tparam realT a real floating point type
313 : *
314 : * \ingroup signal_windows1D
315 : */
316 : template <typename realT>
317 : void nuttal( std::vector<realT> &filt /**< [out] The pre-allocated vector which will store the filter */ )
318 : {
319 : nuttal( filt.data(), filt.size() );
320 : }
321 :
322 : /// The Blackman-Nuttal Windwo
323 : /** See https://en.wikipedia.org/wiki/Window_function
324 : *
325 : * \tparam realT a real floating point type
326 : *
327 : * \ingroup signal_windows1D
328 : */
329 : template <typename realT>
330 : void blackmanNuttal( realT *filt, ///< [out] The pre-allocated vector which will store the filter
331 : size_t N ///< [in] the size of the filter vector
332 : )
333 : {
334 : genCosine<realT>( filt, N, 0.3635819, 0.4891775, 0.1365995, 0.0106411 );
335 : }
336 :
337 : /// The Blackman-Nuttal Windwo
338 : /** See https://en.wikipedia.org/wiki/Window_function
339 : *
340 : * \overload
341 : *
342 : * \tparam realT a real floating point type
343 : *
344 : * \ingroup signal_windows1D
345 : */
346 : template <typename realT>
347 : void blackmanNuttal( std::vector<realT> &filt /**< [out] The pre-allocated vector which will store the filter */ )
348 : {
349 : blackmanNuttal( filt.data(), filt.size() );
350 : }
351 :
352 : /** \brief Create a 2-D Tukey window
353 : *
354 : * Function to create a 2-D Tukey window.
355 : *
356 : * The shape of the window is controlled by alpha. alpha = 0 is a cylinder, alpha=1.0 is the Hann window.
357 : *
358 : * See https://en.wikipedia.org/wiki/Window_function
359 : *
360 : * \ingroup signal_windows2D
361 : */
362 : template <typename realT>
363 : void tukey2d(
364 : realT *filt, ///< [out] a pre-allocated array of size \p rows x \p cols (column major)
365 : int rows, ///< [in] the number of rows in filt
366 : int cols, ///< [in] the number of cols in filt
367 : realT D, ///< [in] the diameter of the window
368 : realT alpha, ///< [in] controls the window shape. 1.0 gives a Hann window, 0.0 gives a cylinder (a.k.a. no window)
369 : realT xc, ///< [in] the desired x center of the window.
370 : realT yc ///< [in] the desired y center of the window.
371 : )
372 : {
373 :
374 : int ii, jj;
375 :
376 : realT rad = 0.5 * ( D - 1.0 );
377 :
378 : realT pi = math::pi<realT>();
379 :
380 : for( int cc = 0; cc < cols; ++cc )
381 : {
382 : realT y = ( (realT)cc ) - yc;
383 : for( int rr = 0; rr < rows; ++rr )
384 : {
385 : realT x = ( (realT)rr ) - xc;
386 :
387 : realT r = sqrt( x * x + y * y );
388 :
389 : // Following mxlib convention of including half pixels
390 : if( r > rad + 0.5 )
391 : {
392 : filt[cc * rows + rr] = 0.0;
393 : }
394 : else if( rad + r > ( D - 1 ) * ( 1 - 0.5 * alpha ) && alpha > 0. )
395 : {
396 : // Have to prevent going greater than N-1 due to half pixel inclusion.
397 : realT dr = rad + r;
398 : if( dr > D - 1 )
399 : dr = D - 1;
400 :
401 : filt[cc * rows + rr] =
402 : 0.5 * ( 1.0 + cos( pi * ( 2. * ( dr ) / ( alpha * ( D - 1 ) ) - 2. / alpha + 1.0 ) ) );
403 : }
404 : else
405 : {
406 : filt[cc * rows + rr] = 1.0;
407 : }
408 : }
409 : }
410 : }
411 :
412 : /** \brief Create a 2-D Tukey window
413 : *
414 : * Function to create a 2-D Tukey window.
415 : *
416 : * The shape of the window is controlled by alpha. alpha = 0 is a cylinder, alpha=1.0 is the Hann window.
417 : *
418 : * See https://en.wikipedia.org/wiki/Window_function
419 : *
420 : * \tparam arrT an Eigen-like array type
421 : *
422 : * \overload
423 : *
424 : * \ingroup signal_windows2D
425 : */
426 : template <typename arrT>
427 : void tukey2d( arrT &filt, ///< [in.out] a pre-allocated array
428 : typename arrT::Scalar D, ///< [in] the diameter of the window
429 : typename arrT::Scalar alpha, ///< [in] controls the window shape. 1.0 gives a Hann window, 0.0 gives a
430 : ///< cylinder (a.k.a. no window)
431 : typename arrT::Scalar xc, ///< [in] the desired x center of the window.
432 : typename arrT::Scalar yc ///< [in] the desired y center of the window.
433 : )
434 : {
435 : tukey2d( filt.data(), filt.rows(), filt.cols(), D, alpha, xc, yc );
436 : }
437 :
438 : /** \brief Create a 2-D Tukey window on an annulus
439 : *
440 : * Function to create a 2-D Tukey window on an annulus.
441 : *
442 : * The shape of the window is controlled by alpha. alpha = 0 is a cylinder, alpha=1.0 is the Hann window.
443 : *
444 : * See https://en.wikipedia.org/wiki/Window_function
445 : *
446 : * \tparam realT is a floating point type
447 : *
448 : * \ingroup signal_windows2D
449 : */
450 : template <typename realT>
451 : void tukey2dAnnulus(
452 : realT *filt, ///< [out] a pre-allocated array of size \p rows x \p cols (column major)
453 : int rows, ///< [in] the number of rows in filt
454 : int cols, ///< [in] the number of cols in filt
455 : realT D, ///< [in] the outer diameter of the window
456 : realT eps, ///< [in] the ratio of inner diameter to the outer diameter
457 : realT alpha, ///< [in] controls the window shape. 1.0 gives a Hann window, 0.0 gives a cylinder (a.k.a. no window)
458 : realT xc, ///< [in] the desired x center of the window.
459 : realT yc ///< [in] the desired y center of the window.
460 : )
461 : {
462 : realT rad = 0.5 * ( D - 1.0 );
463 :
464 : int Z = ( 1 - eps ) * rad + 1.0; // floor((1.0-eps)*(rad)) + 1.0;
465 :
466 : realT pi = math::pi<realT>();
467 :
468 : for( int cc = 0; cc < cols; ++cc )
469 : {
470 : realT y = ( (realT)cc ) - yc;
471 : for( int rr = 0; rr < rows; ++rr )
472 : {
473 : realT x = ( (realT)rr ) - xc;
474 :
475 : realT r = sqrt( x * x + y * y );
476 :
477 : realT z = ( r - eps * ( rad ) );
478 :
479 : // Following mxlib convention of including half pixels
480 : if( r > rad + 0.5 || r < eps * rad )
481 : {
482 : filt[cc * rows + rr] = 0.0;
483 : }
484 : else if( z <= 0.5 * alpha * ( Z - 1 ) && alpha > 0 )
485 : {
486 : filt[cc * rows + rr] = 0.5 * ( 1.0 + cos( pi * ( 2. * z / ( alpha * ( Z - 1 ) ) - 1.0 ) ) );
487 : }
488 : else if( z > ( Z - 1 ) * ( 1. - 0.5 * alpha ) && alpha > 0 )
489 : {
490 : z = z * ( ( Z - 0.5 ) / Z ); // Stretch a little to help with the half pixel
491 : if( z > Z )
492 : z = Z - 1;
493 : filt[cc * rows + rr] =
494 : 0.5 * ( 1.0 + cos( pi * ( 2. * ( z ) / ( alpha * ( Z - 1 ) ) - 2. / alpha + 1.0 ) ) );
495 : }
496 : else
497 : {
498 : filt[cc * rows + rr] = 1.0;
499 : }
500 : }
501 : }
502 : }
503 :
504 : /** \brief Create a 2-D Tukey window on an annulus
505 : *
506 : * Function to create a 2-D Tukey window on an annulus.
507 : *
508 : * The shape of the window is controlled by alpha. alpha = 0 is a cylinder, alpha=1.0 is the Hann window.
509 : *
510 : * See https://en.wikipedia.org/wiki/Window_function
511 : *
512 : * \overload
513 : *
514 : * \tparam realT is a floating point type
515 : *
516 : * \ingroup signal_windows2D
517 : */
518 : template <typename arrT>
519 : void tukey2dAnnulus( arrT &filt, ///< [in,out] a pre-allocated array
520 : typename arrT::Scalar D, ///< [in] the outer diameter of the window
521 : typename arrT::Scalar eps, ///< [in] the ratio of inner diameter to the outer diameter
522 : typename arrT::Scalar alpha, ///< [in] controls the window shape. 1.0 gives a Hann window, 0.0
523 : ///< gives a cylinder (a.k.a. no window)
524 : typename arrT::Scalar xc, ///< [in] the desired x center of the window.
525 : typename arrT::Scalar yc ///< [in] the desired y center of the window.
526 : )
527 : {
528 : return tukey2dAnnulus( filt.data(), filt.rows(), filt.cols(), D, eps, alpha, xc, yc );
529 : }
530 :
531 : /** \brief Create a 2-D Tukey window on a rectangle
532 : *
533 : * Function to create a 2-D Tukey window on a rectangle.
534 : *
535 : * The shape of the window is controlled by alpha. alpha = 0 is a rectangle window, alpha=1.0 is the Hann window.
536 : *
537 : * See https://en.wikipedia.org/wiki/Window_function
538 : *
539 : * \ingroup signal_windows2D
540 : */
541 : template <typename realT>
542 3 : void tukey2dSquare(
543 : realT *filt, ///< [out] a pre-allocated array of size \p rows x \p cols (column major)
544 : int rows, ///< [in] the number of rows in filt
545 : int cols, ///< [in] the number of cols in filt
546 : realT W, ///< [in] the width of the window, corresponds to rows
547 : realT H, ///< [in] the height of the window, corresponds to cols
548 : realT alpha, ///< [in] controls the window shape. 1.0 gives a Hann window, 0.0 gives a cylinder (a.k.a. no window)
549 : realT xc, ///< [in] the desired x center of the window.
550 : realT yc ///< [in] the desired y center of the window.
551 : )
552 : {
553 :
554 : int ii, jj;
555 :
556 3 : realT W2 = 0.5 * ( W - 1.0 );
557 3 : realT H2 = 0.5 * ( H - 1.0 );
558 :
559 3 : realT pi = math::pi<realT>();
560 :
561 771 : for( int cc = 0; cc < cols; ++cc )
562 : {
563 768 : realT y = fabs( ( (realT)cc ) - yc );
564 :
565 768 : if( y > H2 + 0.5 )
566 : {
567 0 : for( int rr = 0; rr < rows; ++rr )
568 : {
569 0 : filt[cc * rows + rr] = 0;
570 : }
571 0 : continue;
572 0 : }
573 :
574 197376 : for( int rr = 0; rr < rows; ++rr )
575 : {
576 196608 : realT x = fabs( ( (realT)rr ) - xc );
577 :
578 196608 : if( x > W2 + 0.5 )
579 : {
580 0 : filt[cc * rows + rr] = 0;
581 0 : continue;
582 : }
583 196608 : else if( ( W2 + x > ( W - 1 ) * ( 1 - 0.5 * alpha ) ) && ( H2 + x > ( H - 1 ) * ( 1 - 0.5 * alpha ) ) &&
584 : alpha > 0. )
585 : {
586 : // Have to prevent going greater than N-1 due to half pixel inclusion.
587 98304 : realT dx = W2 + x;
588 98304 : if( dx > W - 1 )
589 0 : dx = W - 1;
590 :
591 98304 : realT dy = H2 + y;
592 98304 : if( dy > H - 1 )
593 0 : dy = H - 1;
594 :
595 98304 : filt[cc * rows + rr] =
596 98304 : 0.5 * ( 1.0 + cos( pi * ( 2. * ( dx ) / ( alpha * ( W - 1 ) ) - 2. / alpha + 1.0 ) ) );
597 98304 : filt[cc * rows + rr] *=
598 98304 : 0.5 * ( 1.0 + cos( pi * ( 2. * ( dy ) / ( alpha * ( H - 1 ) ) - 2. / alpha + 1.0 ) ) );
599 98304 : }
600 : else
601 : {
602 98304 : filt[cc * rows + rr] = 1.0;
603 : }
604 : }
605 : }
606 3 : }
607 :
608 : } // namespace window
609 : } // namespace sigproc
610 : } // namespace mx
611 :
612 : #endif // signalWindows_hpp
|
