libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
win.c
1 /*
2 * Copyright (c) 2015-2024 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 
27 
28 
29 
30 /* Window functions */
31 int win_bartlett (double *w, int n, int win_type);
32 int win_bartlett_hann (double *w, int n, int win_type);
33 int win_blackman (double *w, int n, int win_type);
34 int win_blackman_harris (double *w, int n, int win_type);
35 int win_blackman_nuttall(double *w, int n, int win_type);
36 int win_cheby (double *w, int n, double param);
37 int win_cos (double *w, int n, int win_type);
38 int win_flat_top (double *w, int n, int win_type);
39 int win_gaussian (double *w, int n, int win_type, double sigma);
40 int win_hamming (double *w, int n, int win_type);
41 int win_hann (double *w, int n, int win_type);
42 int win_kaiser (double* w, int n, int win_type, double param);
43 int win_lanczos (double *w, int n, int win_type);
44 int win_nuttall (double *w, int n, int win_type);
45 int win_rect (double *w, int n);
46 
47 
48 
49 
50 
51 #ifdef DOXYGEN_ENGLISH
52 
187 #endif
188 #ifdef DOXYGEN_RUSSIAN
189 
328 #endif
329 int DSPL_API window(double* w, int n, int win_type, double param)
330 {
331  switch(win_type & DSPL_WIN_MASK)
332  {
333  case DSPL_WIN_BARTLETT:
334  return win_bartlett(w, n, win_type);
335  case DSPL_WIN_BARTLETT_HANN:
336  return win_bartlett_hann(w, n, win_type);
337  case DSPL_WIN_BLACKMAN:
338  return win_blackman(w, n, win_type);
339  case DSPL_WIN_BLACKMAN_HARRIS:
340  return win_blackman_harris(w, n, win_type);
341  case DSPL_WIN_BLACKMAN_NUTTALL:
342  return win_blackman_nuttall(w, n, win_type);
343  case DSPL_WIN_CHEBY:
344  return win_cheby(w, n, param);
345  case DSPL_WIN_FLAT_TOP:
346  return win_flat_top(w, n, win_type);
347  case DSPL_WIN_GAUSSIAN:
348  return win_gaussian(w, n, win_type, param);
349  case DSPL_WIN_HAMMING:
350  return win_hamming(w, n, win_type);
351  case DSPL_WIN_HANN:
352  return win_hann(w, n, win_type);
353  case DSPL_WIN_KAISER:
354  return win_kaiser(w, n, win_type, param);
355  case DSPL_WIN_LANCZOS:
356  return win_lanczos(w, n, win_type);
357  case DSPL_WIN_NUTTALL:
358  return win_nuttall(w, n, win_type);
359  case DSPL_WIN_RECT:
360  return win_rect(w, n);
361  case DSPL_WIN_COS:
362  return win_cos(w, n, win_type);
363  default:
364  return ERROR_WIN_TYPE;
365  }
366  return RES_OK;
367 }
368 
369 
370 
371 /******************************************************************************
372 Barlett window function
373 *******************************************************************************/
374 int win_bartlett(double *w, int n, int win_type)
375 {
376  double x = 0.0;
377  int i;
378 
379  if(!w)
380  return ERROR_PTR;
381  if(n<2)
382  return ERROR_SIZE;
383 
384  switch(win_type & DSPL_WIN_SYM_MASK)
385  {
386  case DSPL_WIN_SYMMETRIC: x = (double)(n-1); break;
387  case DSPL_WIN_PERIODIC : x = (double)n; break;
388  default: return ERROR_WIN_SYM;
389  }
390 
391  for(i = 0; i < n; i++)
392  {
393  w[i] = 2.0 / x * (x * 0.5-fabs((double)i - x * 0.5));
394  }
395  return RES_OK;
396 }
397 
398 
399 
400 
401 
402 /******************************************************************************
403 Barlett - Hann window function
404 ******************************************************************************/
405 int win_bartlett_hann(double *w, int n, int win_type)
406 {
407  double y;
408  double x = 0.0;
409  int i;
410 
411  if(!w)
412  return ERROR_PTR;
413  if(n<2)
414  return ERROR_SIZE;
415 
416  switch(win_type & DSPL_WIN_SYM_MASK)
417  {
418  case DSPL_WIN_SYMMETRIC: x = 1.0/(double)(n-1); break;
419  case DSPL_WIN_PERIODIC : x = 1.0/(double)n; break;
420  default: return ERROR_WIN_SYM;
421  }
422 
423  y = 0.0;
424  for(i = 0; i<n; i++)
425  {
426  w[i] = 0.62 - 0.48 * fabs(y-0.5)-0.38*cos(M_2PI*y);
427  y += x;
428  }
429  return RES_OK;
430 }
431 
432 
433 
434 
435 
436 /******************************************************************************
437 Blackman window function
438 ******************************************************************************/
439 int win_blackman(double *w, int n, int win_type)
440 {
441  double y;
442  double x = 0.0;
443  int i;
444 
445  if(!w)
446  return ERROR_PTR;
447  if(n<2)
448  return ERROR_SIZE;
449 
450  switch(win_type & DSPL_WIN_SYM_MASK)
451  {
452  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
453  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
454  default: return ERROR_WIN_SYM;
455  }
456 
457  y = 0.0;
458  for(i = 0; i<n; i++)
459  {
460  w[i] = 0.42 - 0.5* cos(y)+0.08*cos(2.0*y);
461  y += x;
462  }
463  return RES_OK;
464 }
465 
466 
467 
468 
469 /******************************************************************************
470 Blackman - Harris window function
471 ******************************************************************************/
472 int win_blackman_harris(double *w, int n, int win_type)
473 {
474  double y;
475  double x = 0.0;
476  double a0 = 0.35875;
477  double a1 = 0.48829;
478  double a2 = 0.14128;
479  double a3 = 0.01168;
480  int i;
481 
482  if(!w)
483  return ERROR_PTR;
484  if(n<2)
485  return ERROR_SIZE;
486 
487  switch(win_type & DSPL_WIN_SYM_MASK)
488  {
489  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
490  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
491  default: return ERROR_WIN_SYM;
492  }
493 
494  y = 0.0;
495  for(i = 0; i<n; i++)
496  {
497  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
498  y += x;
499  }
500  return RES_OK;
501 }
502 
503 
504 
505 
506 /******************************************************************************
507 Blackman - Nuttull window function
508 ******************************************************************************/
509 int win_blackman_nuttall(double *w, int n, int win_type)
510 {
511  double y;
512  double x = 0.0;
513  double a0 = 0.3635819;
514  double a1 = 0.4891775;
515  double a2 = 0.1365995;
516  double a3 = 0.0106411;
517  int i;
518 
519 
520  if(!w)
521  return ERROR_PTR;
522  if(n<2)
523  return ERROR_SIZE;
524 
525  switch(win_type & DSPL_WIN_SYM_MASK)
526  {
527  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
528  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
529  default: return ERROR_WIN_SYM;
530  }
531 
532  y = 0.0;
533  for(i = 0; i<n; i++)
534  {
535  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
536  y += x;
537  }
538  return RES_OK;
539 }
540 
541 
542 
543 
544 
545 
546 /******************************************************************************
547 Chebyshev parametric window function
548 param sets spectrum sidelobes level in dB
549 ATTENTION! ONLY SYMMETRIC WINDOW
550 *******************************************************************************/
551 int win_cheby(double *w, int n, double param)
552 {
553  int k, i, m;
554  double z, dz, sum = 0, wmax=0, r1, x0, chx, chy, in;
555 
556  if(!w)
557  return ERROR_PTR;
558 
559  if(n<2)
560  return ERROR_SIZE;
561 
562  if(param <= 0.0)
563  return ERROR_WIN_PARAM;
564 
565  r1 = pow(10, param/20);
566  x0 = cosh((1.0/(double)(n-1)) * acosh(r1));
567 
568  /* check window length even or odd */
569  if(n%2==0)
570  {
571  dz = 0.5;
572  m = n/2-1;
573  }
574  else
575  {
576  m = (n-1)/2;
577  dz = 0.0;
578  }
579 
580  for(k = 0; k < m+2; k++)
581  {
582  z = (double)(k - m) - dz;
583  sum = 0;
584 
585  for(i = 1; i <= m; i++)
586  {
587  in = (double)i / (double)n;
588  chx = x0 * cos(M_PI * in);
589  cheby_poly1(&chx, 1, n-1, &chy);
590  sum += chy * cos(2.0 * z * M_PI * in);
591  }
592 
593  w[k] = r1 + 2.0 * sum;
594  w[n-1-k] = w[k];
595 
596  /* max value calculation */
597  if(w[k]>wmax)
598  wmax=w[k];
599  }
600 
601  /* normalization */
602  for(k=0; k < n; k++)
603  w[k] /= wmax;
604 
605  return RES_OK;
606 }
607 
608 
609 
610 /******************************************************************************
611 Cosine window function
612 ******************************************************************************/
613 int win_cos(double *w, int n, int win_type)
614 {
615  double y;
616  double x = 0.0;
617  int i;
618 
619  if(!w)
620  return ERROR_PTR;
621  if(n<2)
622  return ERROR_SIZE;
623 
624  switch(win_type & DSPL_WIN_SYM_MASK)
625  {
626  case DSPL_WIN_SYMMETRIC: x = M_PI/(double)(n-1); break;
627  case DSPL_WIN_PERIODIC : x = M_PI/(double)n; break;
628  default: return ERROR_WIN_SYM;
629  }
630 
631  y = 0.0;
632  for(i = 0; i<n; i++)
633  {
634  w[i] = sin(y);
635  y += x;
636  }
637  return RES_OK;
638 }
639 
640 
641 
642 
643 
644 
645 /******************************************************************************
646 Flat - Top window function
647 ******************************************************************************/
648 int win_flat_top(double *w, int n, int win_type)
649 {
650  double y;
651  double x = 0.0;
652  double a0 = 1.0;
653  double a1 = 1.93;
654  double a2 = 1.29;
655  double a3 = 0.388;
656  double a4 = 0.032;
657  int i;
658 
659  if(!w)
660  return ERROR_PTR;
661  if(n<2)
662  return ERROR_SIZE;
663 
664  switch(win_type & DSPL_WIN_SYM_MASK)
665  {
666  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
667  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
668  default: return ERROR_WIN_SYM;
669  }
670 
671  y = 0.0;
672  for(i = 0; i<n; i++)
673  {
674  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y)+a4*cos(4.0*y);
675  y += x;
676  }
677  return RES_OK;
678 }
679 
680 
681 
682 
683 
684 
685 /******************************************************************************
686 Gaussian window function
687 ******************************************************************************/
688 int win_gaussian(double *w, int n, int win_type, double alpha)
689 {
690  double a = 0.0;
691  double y;
692  double sigma;
693  int i;
694 
695  if(!w)
696  return ERROR_PTR;
697  if(n<2)
698  return ERROR_SIZE;
699 
700  switch(win_type & DSPL_WIN_SYM_MASK)
701  {
702  case DSPL_WIN_SYMMETRIC: a = (double)(n-1)*0.5; break;
703  case DSPL_WIN_PERIODIC : a = (double)(n)*0.5; break;
704  default: return ERROR_WIN_SYM;
705  }
706 
707 
708  sigma = 1.0 / (alpha * a);
709  for(i = 0; i<n; i++)
710  {
711  y = ((double)i - a)*sigma;
712  w[i] = exp(-0.5*y*y);
713  }
714  return RES_OK;
715 }
716 
717 
718 
719 
720 
721 
722 /******************************************************************************
723 Hamming window function
724 ******************************************************************************/
725 int win_hamming(double *w, int n, int win_type)
726 {
727  double x = 0.0;
728  double y;
729  int i;
730 
731  if(!w)
732  return ERROR_PTR;
733  if(n<2)
734  return ERROR_SIZE;
735 
736  switch(win_type & DSPL_WIN_SYM_MASK)
737  {
738  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
739  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
740  default: return ERROR_WIN_SYM;
741  }
742 
743  y = 0.0;
744  for(i = 0; i < n; i++)
745  {
746  w[i] = 0.54-0.46*cos(y);
747  y += x;
748  }
749  return RES_OK;
750 }
751 
752 
753 
754 
755 /******************************************************************************
756 Hann window function
757 ******************************************************************************/
758 int win_hann(double *w, int n, int win_type)
759 {
760  double x = 0.0;
761  double y;
762  int i;
763 
764  if(!w)
765  return ERROR_PTR;
766  if(n<2)
767  return ERROR_SIZE;
768 
769  switch(win_type & DSPL_WIN_SYM_MASK)
770  {
771  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
772  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
773  default: return ERROR_WIN_SYM;
774  }
775 
776  y = 0.0;
777  for(i = 0; i < n; i++)
778  {
779  w[i] = 0.5*(1-cos(y));
780  y += x;
781  }
782  return RES_OK;
783 }
784 
785 
786 /******************************************************************************
787 Kaiser window function
788 ******************************************************************************/
789 int win_kaiser(double* w, int n, int win_type, double param)
790 {
791  double num, den, x, y, L;
792  int i, err;
793  if(!w)
794  return ERROR_PTR;
795  if(n<2)
796  return ERROR_SIZE;
797 
798  switch(win_type & DSPL_WIN_SYM_MASK)
799  {
800  case DSPL_WIN_SYMMETRIC: L = (double)(n-1) / 2.0; break;
801  case DSPL_WIN_PERIODIC : L = (double)n / 2.0; break;
802  default: return ERROR_WIN_SYM;
803  }
804 
805  err = bessel_i0(&param, 1, &den);
806  if(err != RES_OK)
807  return err;
808  for(i = 0; i < n; i++)
809  {
810  x = 2.0*((double)i - L) / (double)n;
811  y = param * sqrt(1.0 - x*x);
812  err = bessel_i0(&y, 1, &num);
813  if(err != RES_OK)
814  return err;
815  w[i] = num / den;
816  }
817  return err;
818 }
819 
820 
821 
822 /******************************************************************************
823 Lanczos window function
824 ******************************************************************************/
825 int win_lanczos(double *w, int n, int win_type)
826 {
827  double y;
828  double x = 0.0;
829  int i;
830 
831  if(!w)
832  return ERROR_PTR;
833  if(n<2)
834  return ERROR_SIZE;
835 
836  switch(win_type & DSPL_WIN_SYM_MASK)
837  {
838  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
839  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
840  default: return ERROR_WIN_SYM;
841  }
842 
843  y = 0.0;
844  for(i = 0; i < n; i++)
845  {
846  if((y - M_PI)==0.0)
847  w[i] = 1.0;
848  else
849  w[i] = sin(y - M_PI)/(y - M_PI);
850  y += x;
851  }
852  return RES_OK;
853 
854 }
855 
856 
857 
858 /******************************************************************************
859 Nuttall window function
860 ******************************************************************************/
861 int win_nuttall(double *w, int n, int win_type)
862 {
863  double y;
864  double x = 0.0;
865  double a0 = 0.355768;
866  double a1 = 0.487396;
867  double a2 = 0.144232;
868  double a3 = 0.012604;
869  int i;
870 
871  if(!w)
872  return ERROR_PTR;
873  if(n<2)
874  return ERROR_SIZE;
875 
876  switch(win_type & DSPL_WIN_SYM_MASK)
877  {
878  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
879  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
880  default: return ERROR_WIN_SYM;
881  }
882 
883  y = 0.0;
884  for(i = 0; i < n; i++)
885  {
886  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
887  y += x;
888  }
889  return RES_OK;
890 }
891 
892 
893 
894 
895 
896 /******************************************************************************
897 Rectangle window function
898 ******************************************************************************/
899 int win_rect(double *w, int n)
900 {
901  int i;
902 
903  if(!w)
904  return ERROR_PTR;
905  if(n<2)
906  return ERROR_SIZE;
907 
908  for(i = 0; i < n; i++)
909  w[i] = 1.0;
910  return RES_OK;
911 }
912 
window
int window(double *w, int n, int win_type, double param)
Расчет функции оконного взвешивания
Definition: win.c:329
ERROR_PTR
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:610
ERROR_SIZE
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:618
ERROR_WIN_TYPE
#define ERROR_WIN_TYPE
Неизвестный тип оконной функции.
Definition: dspl.h:627
cheby_poly1
int cheby_poly1(double *x, int n, int ord, double *y)
Многочлен Чебышева первого рода порядка ord
Definition: cheby_poly1.c:136
ERROR_WIN_PARAM
#define ERROR_WIN_PARAM
Ошибка значения параметра оконной функции. Для каждой параметрической оконной функции существуют допу...
Definition: dspl.h:625
RES_OK
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:558
sum
int sum(double *x, int n, double *s)
Расчет суммы элеметов массива
Definition: sum.c:76
bessel_i0
int bessel_i0(double *x, int n, double *y)
Модифицированная функция Бесселя первого рода .
Definition: bessel_i0.c:116
ERROR_WIN_SYM
#define ERROR_WIN_SYM
Симметричность или периодичность заданного окна не поддерживается.
Definition: dspl.h:626
Создано системой   doxygen 1.8.17
Документация libdspl-2.0 доступна по лицензии Creative Commons «Attribution-ShareAlike» («Атрибуция-СохранениеУсловий») 4.0 Всемирная.