libdspl-2.0
Digital Signal Processing Algorithm Library
win.c
1 /*
2 * Copyright (c) 2015-2019 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of libdspl-2.0.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * DSPL 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 Lesser General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 
23 #include <stdio.h>
24 #include <math.h>
25 #include "dspl.h"
26 #include "dspl_internal.h"
27 
28 
29 
30 #ifdef DOXYGEN_ENGLISH
31 
166 #endif
167 #ifdef DOXYGEN_RUSSIAN
168 
307 #endif
308 int DSPL_API window(double* w, int n, int win_type, double param)
309 {
310  switch(win_type & DSPL_WIN_MASK)
311  {
312  case DSPL_WIN_BARTLETT:
313  return win_bartlett(w, n, win_type);
314  case DSPL_WIN_BARTLETT_HANN:
315  return win_bartlett_hann(w, n, win_type);
316  case DSPL_WIN_BLACKMAN:
317  return win_blackman(w, n, win_type);
318  case DSPL_WIN_BLACKMAN_HARRIS:
319  return win_blackman_harris(w, n, win_type);
320  case DSPL_WIN_BLACKMAN_NUTTALL:
321  return win_blackman_nuttall(w, n, win_type);
322  case DSPL_WIN_CHEBY:
323  return win_cheby(w, n, param);
324  case DSPL_WIN_FLAT_TOP:
325  return win_flat_top(w, n, win_type);
326  case DSPL_WIN_GAUSSIAN:
327  return win_gaussian(w, n, win_type, param);
328  case DSPL_WIN_HAMMING:
329  return win_hamming(w, n, win_type);
330  case DSPL_WIN_HANN:
331  return win_hann(w, n, win_type);
332  case DSPL_WIN_KAISER:
333  return win_kaiser(w, n, win_type, param);
334  case DSPL_WIN_LANCZOS:
335  return win_lanczos(w, n, win_type);
336  case DSPL_WIN_NUTTALL:
337  return win_nuttall(w, n, win_type);
338  case DSPL_WIN_RECT:
339  return win_rect(w, n);
340  case DSPL_WIN_COS:
341  return win_cos(w, n, win_type);
342  default:
343  return ERROR_WIN_TYPE;
344  }
345  return RES_OK;
346 }
347 
348 
349 
350 /******************************************************************************
351 Barlett window function
352 *******************************************************************************/
353 int win_bartlett(double *w, int n, int win_type)
354 {
355  double x = 0.0;
356  int i;
357 
358  if(!w)
359  return ERROR_PTR;
360  if(n<2)
361  return ERROR_SIZE;
362 
363  switch(win_type & DSPL_WIN_SYM_MASK)
364  {
365  case DSPL_WIN_SYMMETRIC: x = (double)(n-1); break;
366  case DSPL_WIN_PERIODIC : x = (double)n; break;
367  default: return ERROR_WIN_SYM;
368  }
369 
370  for(i = 0; i < n; i++)
371  {
372  w[i] = 2.0 / x * (x * 0.5-fabs((double)i - x * 0.5));
373  }
374  return RES_OK;
375 }
376 
377 
378 
379 
380 
381 /******************************************************************************
382 Barlett - Hann window function
383 ******************************************************************************/
384 int win_bartlett_hann(double *w, int n, int win_type)
385 {
386  double y;
387  double x = 0.0;
388  int i;
389 
390  if(!w)
391  return ERROR_PTR;
392  if(n<2)
393  return ERROR_SIZE;
394 
395  switch(win_type & DSPL_WIN_SYM_MASK)
396  {
397  case DSPL_WIN_SYMMETRIC: x = 1.0/(double)(n-1); break;
398  case DSPL_WIN_PERIODIC : x = 1.0/(double)n; break;
399  default: return ERROR_WIN_SYM;
400  }
401 
402  y = 0.0;
403  for(i = 0; i<n; i++)
404  {
405  w[i] = 0.62 - 0.48 * fabs(y-0.5)-0.38*cos(M_2PI*y);
406  y += x;
407  }
408  return RES_OK;
409 }
410 
411 
412 
413 
414 
415 /******************************************************************************
416 Blackman window function
417 ******************************************************************************/
418 int win_blackman(double *w, int n, int win_type)
419 {
420  double y;
421  double x = 0.0;
422  int i;
423 
424  if(!w)
425  return ERROR_PTR;
426  if(n<2)
427  return ERROR_SIZE;
428 
429  switch(win_type & DSPL_WIN_SYM_MASK)
430  {
431  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
432  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
433  default: return ERROR_WIN_SYM;
434  }
435 
436  y = 0.0;
437  for(i = 0; i<n; i++)
438  {
439  w[i] = 0.42 - 0.5* cos(y)+0.08*cos(2.0*y);
440  y += x;
441  }
442  return RES_OK;
443 }
444 
445 
446 
447 
448 /******************************************************************************
449 Blackman - Harris window function
450 ******************************************************************************/
451 int win_blackman_harris(double *w, int n, int win_type)
452 {
453  double y;
454  double x = 0.0;
455  double a0 = 0.35875;
456  double a1 = 0.48829;
457  double a2 = 0.14128;
458  double a3 = 0.01168;
459  int i;
460 
461  if(!w)
462  return ERROR_PTR;
463  if(n<2)
464  return ERROR_SIZE;
465 
466  switch(win_type & DSPL_WIN_SYM_MASK)
467  {
468  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
469  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
470  default: return ERROR_WIN_SYM;
471  }
472 
473  y = 0.0;
474  for(i = 0; i<n; i++)
475  {
476  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
477  y += x;
478  }
479  return RES_OK;
480 }
481 
482 
483 
484 
485 /******************************************************************************
486 Blackman - Nuttull window function
487 ******************************************************************************/
488 int win_blackman_nuttall(double *w, int n, int win_type)
489 {
490  double y;
491  double x = 0.0;
492  double a0 = 0.3635819;
493  double a1 = 0.4891775;
494  double a2 = 0.1365995;
495  double a3 = 0.0106411;
496  int i;
497 
498 
499  if(!w)
500  return ERROR_PTR;
501  if(n<2)
502  return ERROR_SIZE;
503 
504  switch(win_type & DSPL_WIN_SYM_MASK)
505  {
506  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
507  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
508  default: return ERROR_WIN_SYM;
509  }
510 
511  y = 0.0;
512  for(i = 0; i<n; i++)
513  {
514  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
515  y += x;
516  }
517  return RES_OK;
518 }
519 
520 
521 
522 
523 
524 
525 /******************************************************************************
526 Chebyshev parametric window function
527 param sets spectrum sidelobes level in dB
528 ATTENTION! ONLY SYMMETRIC WINDOW
529 *******************************************************************************/
530 int win_cheby(double *w, int n, double param)
531 {
532  int k, i, m;
533  double z, dz, sum = 0, wmax=0, r1, x0, chx, chy, in;
534 
535  if(!w)
536  return ERROR_PTR;
537 
538  if(n<2)
539  return ERROR_SIZE;
540 
541  if(param <= 0.0)
542  return ERROR_WIN_PARAM;
543 
544  r1 = pow(10, param/20);
545  x0 = cosh((1.0/(double)(n-1)) * acosh(r1));
546 
547  /* check window length even or odd */
548  if(n%2==0)
549  {
550  dz = 0.5;
551  m = n/2-1;
552  }
553  else
554  {
555  m = (n-1)/2;
556  dz = 0.0;
557  }
558 
559  for(k = 0; k < m+2; k++)
560  {
561  z = (double)(k - m) - dz;
562  sum = 0;
563 
564  for(i = 1; i <= m; i++)
565  {
566  in = (double)i / (double)n;
567  chx = x0 * cos(M_PI * in);
568  cheby_poly1(&chx, 1, n-1, &chy);
569  sum += chy * cos(2.0 * z * M_PI * in);
570  }
571 
572  w[k] = r1 + 2.0 * sum;
573  w[n-1-k] = w[k];
574 
575  /* max value calculation */
576  if(w[k]>wmax)
577  wmax=w[k];
578  }
579 
580  /* normalization */
581  for(k=0; k < n; k++)
582  w[k] /= wmax;
583 
584  return RES_OK;
585 }
586 
587 
588 
589 /******************************************************************************
590 Cosine window function
591 ******************************************************************************/
592 int win_cos(double *w, int n, int win_type)
593 {
594  double y;
595  double x = 0.0;
596  int i;
597 
598  if(!w)
599  return ERROR_PTR;
600  if(n<2)
601  return ERROR_SIZE;
602 
603  switch(win_type & DSPL_WIN_SYM_MASK)
604  {
605  case DSPL_WIN_SYMMETRIC: x = M_PI/(double)(n-1); break;
606  case DSPL_WIN_PERIODIC : x = M_PI/(double)n; break;
607  default: return ERROR_WIN_SYM;
608  }
609 
610  y = 0.0;
611  for(i = 0; i<n; i++)
612  {
613  w[i] = sin(y);
614  y += x;
615  }
616  return RES_OK;
617 }
618 
619 
620 
621 
622 
623 
624 /******************************************************************************
625 Flat - Top window function
626 ******************************************************************************/
627 int win_flat_top(double *w, int n, int win_type)
628 {
629  double y;
630  double x = 0.0;
631  double a0 = 1.0;
632  double a1 = 1.93;
633  double a2 = 1.29;
634  double a3 = 0.388;
635  double a4 = 0.032;
636  int i;
637 
638  if(!w)
639  return ERROR_PTR;
640  if(n<2)
641  return ERROR_SIZE;
642 
643  switch(win_type & DSPL_WIN_SYM_MASK)
644  {
645  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
646  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
647  default: return ERROR_WIN_SYM;
648  }
649 
650  y = 0.0;
651  for(i = 0; i<n; i++)
652  {
653  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y)+a4*cos(4.0*y);
654  y += x;
655  }
656  return RES_OK;
657 }
658 
659 
660 
661 
662 
663 
664 /******************************************************************************
665 Gaussian window function
666 ******************************************************************************/
667 int win_gaussian(double *w, int n, int win_type, double alpha)
668 {
669  double a = 0.0;
670  double y;
671  double sigma;
672  int i;
673 
674  if(!w)
675  return ERROR_PTR;
676  if(n<2)
677  return ERROR_SIZE;
678 
679  switch(win_type & DSPL_WIN_SYM_MASK)
680  {
681  case DSPL_WIN_SYMMETRIC: a = (double)(n-1)*0.5; break;
682  case DSPL_WIN_PERIODIC : a = (double)(n)*0.5; break;
683  default: return ERROR_WIN_SYM;
684  }
685 
686 
687  sigma = 1.0 / (alpha * a);
688  for(i = 0; i<n; i++)
689  {
690  y = ((double)i - a)*sigma;
691  w[i] = exp(-0.5*y*y);
692  }
693  return RES_OK;
694 }
695 
696 
697 
698 
699 
700 
701 /******************************************************************************
702 Hamming window function
703 ******************************************************************************/
704 int win_hamming(double *w, int n, int win_type)
705 {
706  double x = 0.0;
707  double y;
708  int i;
709 
710  if(!w)
711  return ERROR_PTR;
712  if(n<2)
713  return ERROR_SIZE;
714 
715  switch(win_type & DSPL_WIN_SYM_MASK)
716  {
717  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
718  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
719  default: return ERROR_WIN_SYM;
720  }
721 
722  y = 0.0;
723  for(i = 0; i < n; i++)
724  {
725  w[i] = 0.54-0.46*cos(y);
726  y += x;
727  }
728  return RES_OK;
729 }
730 
731 
732 
733 
734 /******************************************************************************
735 Hann window function
736 ******************************************************************************/
737 int win_hann(double *w, int n, int win_type)
738 {
739  double x = 0.0;
740  double y;
741  int i;
742 
743  if(!w)
744  return ERROR_PTR;
745  if(n<2)
746  return ERROR_SIZE;
747 
748  switch(win_type & DSPL_WIN_SYM_MASK)
749  {
750  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
751  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
752  default: return ERROR_WIN_SYM;
753  }
754 
755  y = 0.0;
756  for(i = 0; i < n; i++)
757  {
758  w[i] = 0.5*(1-cos(y));
759  y += x;
760  }
761  return RES_OK;
762 }
763 
764 
765 /******************************************************************************
766 Kaiser window function
767 ******************************************************************************/
768 int win_kaiser(double* w, int n, int win_type, double param)
769 {
770  double num, den, x, y, L;
771  int i, err;
772  if(!w)
773  return ERROR_PTR;
774  if(n<2)
775  return ERROR_SIZE;
776 
777  switch(win_type & DSPL_WIN_SYM_MASK)
778  {
779  case DSPL_WIN_SYMMETRIC: L = (double)(n-1) / 2.0; break;
780  case DSPL_WIN_PERIODIC : L = (double)n / 2.0; break;
781  default: return ERROR_WIN_SYM;
782  }
783 
784  err = bessel_i0(&param, 1, &den);
785  if(err != RES_OK)
786  return err;
787  for(i = 0; i < n; i++)
788  {
789  x = 2.0*((double)i - L) / (double)n;
790  y = param * sqrt(1.0 - x*x);
791  err = bessel_i0(&y, 1, &num);
792  if(err != RES_OK)
793  return err;
794  w[i] = num / den;
795  }
796  return err;
797 }
798 
799 
800 
801 /******************************************************************************
802 Lanczos window function
803 ******************************************************************************/
804 int win_lanczos(double *w, int n, int win_type)
805 {
806  double y;
807  double x = 0.0;
808  int i;
809 
810  if(!w)
811  return ERROR_PTR;
812  if(n<2)
813  return ERROR_SIZE;
814 
815  switch(win_type & DSPL_WIN_SYM_MASK)
816  {
817  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
818  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
819  default: return ERROR_WIN_SYM;
820  }
821 
822  y = 0.0;
823  for(i = 0; i < n; i++)
824  {
825  if((y - M_PI)==0.0)
826  w[i] = 1.0;
827  else
828  w[i] = sin(y - M_PI)/(y - M_PI);
829  y += x;
830  }
831  return RES_OK;
832 
833 }
834 
835 
836 
837 /******************************************************************************
838 Nuttall window function
839 ******************************************************************************/
840 int win_nuttall(double *w, int n, int win_type)
841 {
842  double y;
843  double x = 0.0;
844  double a0 = 0.355768;
845  double a1 = 0.487396;
846  double a2 = 0.144232;
847  double a3 = 0.012604;
848  int i;
849 
850  if(!w)
851  return ERROR_PTR;
852  if(n<2)
853  return ERROR_SIZE;
854 
855  switch(win_type & DSPL_WIN_SYM_MASK)
856  {
857  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
858  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
859  default: return ERROR_WIN_SYM;
860  }
861 
862  y = 0.0;
863  for(i = 0; i < n; i++)
864  {
865  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
866  y += x;
867  }
868  return RES_OK;
869 }
870 
871 
872 
873 
874 
875 /******************************************************************************
876 Rectangle window function
877 ******************************************************************************/
878 int win_rect(double *w, int n)
879 {
880  int i;
881 
882  if(!w)
883  return ERROR_PTR;
884  if(n<2)
885  return ERROR_SIZE;
886 
887  for(i = 0; i < n; i++)
888  w[i] = 1.0;
889  return RES_OK;
890 }
891 
int DSPL_API window(double *w, int n, int win_type, double param)
Window function calculation.
Definition: win.c:308
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:549
#define ERROR_SIZE
Error array size. This error occurs when in addition to the pointer the wrong input is passed to the ...
Definition: dspl.h:557
#define ERROR_WIN_TYPE
Unknown window function type.
Definition: dspl.h:566
int DSPL_API cheby_poly1(double *x, int n, int ord, double *y)
Chebyshev polynomial of the first kind order ord
Definition: cheby.c:134
#define ERROR_WIN_PARAM
window function parameter error. Valid parameter values exist for each parametric window function.
Definition: dspl.h:564
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:497
int DSPL_API bessel_i0(double *x, int n, double *y)
Modified Bessel Function of the First Kind – [1].
Definition: math.c:386
#define ERROR_WIN_SYM
Symmetry or periodicity of a given window is not supported.
Definition: dspl.h:565