libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
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 /*******************************************************************************
31 Window function
32 *******************************************************************************/
33 int window(double* w, int n, int win_type, double param)
34 {
35 
36  switch(win_type & DSPL_WIN_MASK)
37  {
38  case DSPL_WIN_BARTLETT:
39  return win_bartlett(w, n, win_type);
40  case DSPL_WIN_BARTLETT_HANN:
41  return win_bartlett_hann(w, n, win_type);
42  case DSPL_WIN_BLACKMAN:
43  return win_blackman(w, n, win_type);
44  case DSPL_WIN_BLACKMAN_HARRIS:
45  return win_blackman_harris(w, n, win_type);
46  case DSPL_WIN_BLACKMAN_NUTTALL:
47  return win_blackman_nuttall(w, n, win_type);
48  case DSPL_WIN_CHEBY:
49  return win_cheby(w, n, param);
50  case DSPL_WIN_FLAT_TOP:
51  return win_flat_top(w, n, win_type);
52  case DSPL_WIN_GAUSSIAN:
53  return win_gaussian(w, n, win_type, param);
54  case DSPL_WIN_HAMMING:
55  return win_hamming(w, n, win_type);
56  case DSPL_WIN_HANN:
57  return win_hann(w, n, win_type);
58  case DSPL_WIN_KAISER:
59  return win_kaiser(w, n, win_type, param);
60  case DSPL_WIN_LANCZOS:
61  return win_lanczos(w, n, win_type);
62  case DSPL_WIN_NUTTALL:
63  return win_nuttall(w, n, win_type);
64  case DSPL_WIN_RECT:
65  return win_rect(w, n);
66  case DSPL_WIN_COS:
67  return win_cos(w, n, win_type);
68  default:
69  return ERROR_WIN_TYPE;
70  }
71  return RES_OK;
72 }
73 
74 
75 
76 /******************************************************************************
77 Barlett window function
78 *******************************************************************************/
79 int win_bartlett(double *w, int n, int win_type)
80 {
81  double x = 0.0;
82  int i;
83 
84  if(!w)
85  return ERROR_PTR;
86  if(n<2)
87  return ERROR_SIZE;
88 
89  switch(win_type & DSPL_WIN_SYM_MASK)
90  {
91  case DSPL_WIN_SYMMETRIC: x = (double)(n-1); break;
92  case DSPL_WIN_PERIODIC : x = (double)n; break;
93  default: return ERROR_WIN_SYM;
94  }
95 
96  for(i = 0; i < n; i++)
97  {
98  w[i] = 2.0 / x * (x * 0.5-fabs((double)i - x * 0.5));
99  }
100  return RES_OK;
101 }
102 
103 
104 
105 
106 
107 /******************************************************************************
108 Barlett - Hann window function
109 ******************************************************************************/
110 int win_bartlett_hann(double *w, int n, int win_type)
111 {
112  double y;
113  double x = 0.0;
114  int i;
115 
116  if(!w)
117  return ERROR_PTR;
118  if(n<2)
119  return ERROR_SIZE;
120 
121  switch(win_type & DSPL_WIN_SYM_MASK)
122  {
123  case DSPL_WIN_SYMMETRIC: x = 1.0/(double)(n-1); break;
124  case DSPL_WIN_PERIODIC : x = 1.0/(double)n; break;
125  default: return ERROR_WIN_SYM;
126  }
127 
128  y = 0.0;
129  for(i = 0; i<n; i++)
130  {
131  w[i] = 0.62 - 0.48 * fabs(y-0.5)-0.38*cos(M_2PI*y);
132  y += x;
133  }
134  return RES_OK;
135 }
136 
137 
138 
139 
140 
141 /******************************************************************************
142 Blackman window function
143 ******************************************************************************/
144 int win_blackman(double *w, int n, int win_type)
145 {
146  double y;
147  double x = 0.0;
148  int i;
149 
150  if(!w)
151  return ERROR_PTR;
152  if(n<2)
153  return ERROR_SIZE;
154 
155  switch(win_type & DSPL_WIN_SYM_MASK)
156  {
157  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
158  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
159  default: return ERROR_WIN_SYM;
160  }
161 
162  y = 0.0;
163  for(i = 0; i<n; i++)
164  {
165  w[i] = 0.42 - 0.5* cos(y)+0.08*cos(2.0*y);
166  y += x;
167  }
168  return RES_OK;
169 }
170 
171 
172 
173 
174 /******************************************************************************
175 Blackman - Harris window function
176 ******************************************************************************/
177 int win_blackman_harris(double *w, int n, int win_type)
178 {
179  double y;
180  double x = 0.0;
181  double a0 = 0.35875;
182  double a1 = 0.48829;
183  double a2 = 0.14128;
184  double a3 = 0.01168;
185  int i;
186 
187  if(!w)
188  return ERROR_PTR;
189  if(n<2)
190  return ERROR_SIZE;
191 
192  switch(win_type & DSPL_WIN_SYM_MASK)
193  {
194  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
195  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
196  default: return ERROR_WIN_SYM;
197  }
198 
199  y = 0.0;
200  for(i = 0; i<n; i++)
201  {
202  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
203  y += x;
204  }
205  return RES_OK;
206 }
207 
208 
209 
210 
211 /******************************************************************************
212 Blackman - Nuttull window function
213 ******************************************************************************/
214 int win_blackman_nuttall(double *w, int n, int win_type)
215 {
216  double y;
217  double x = 0.0;
218  double a0 = 0.3635819;
219  double a1 = 0.4891775;
220  double a2 = 0.1365995;
221  double a3 = 0.0106411;
222  int i;
223 
224 
225  if(!w)
226  return ERROR_PTR;
227  if(n<2)
228  return ERROR_SIZE;
229 
230  switch(win_type & DSPL_WIN_SYM_MASK)
231  {
232  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
233  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
234  default: return ERROR_WIN_SYM;
235  }
236 
237  y = 0.0;
238  for(i = 0; i<n; i++)
239  {
240  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
241  y += x;
242  }
243  return RES_OK;
244 }
245 
246 
247 
248 
249 
250 
251 /******************************************************************************
252 Chebyshev parametric window function
253 param sets spectrum sidelobes level in dB
254 ATTENTION! ONLY SYMMETRIC WINDOW
255 *******************************************************************************/
256 int win_cheby(double *w, int n, double param)
257 {
258  int k, i, m;
259  double z, dz, sum = 0, wmax=0, r1, x0, chx, chy, in;
260 
261  if(!w)
262  return ERROR_PTR;
263 
264  if(n<2)
265  return ERROR_SIZE;
266 
267  if(param <= 0.0)
268  return ERROR_WIN_PARAM;
269 
270  r1 = pow(10, param/20);
271  x0 = cosh((1.0/(double)(n-1)) * acosh(r1));
272 
273  /* check window length even or odd */
274  if(n%2==0)
275  {
276  dz = 0.5;
277  m = n/2-1;
278  }
279  else
280  {
281  m = (n-1)/2;
282  dz = 0.0;
283  }
284 
285  for(k = 0; k < m+2; k++)
286  {
287  z = (double)(k - m) - dz;
288  sum = 0;
289 
290  for(i = 1; i <= m; i++)
291  {
292  in = (double)i / (double)n;
293  chx = x0 * cos(M_PI * in);
294  cheby_poly1(&chx, 1, n-1, &chy);
295  sum += chy * cos(2.0 * z * M_PI * in);
296  }
297 
298  w[k] = r1 + 2.0 * sum;
299  w[n-1-k] = w[k];
300 
301  /* max value calculation */
302  if(w[k]>wmax)
303  wmax=w[k];
304  }
305 
306  /* normalization */
307  for(k=0; k < n; k++)
308  w[k] /= wmax;
309 
310  return RES_OK;
311 }
312 
313 
314 
315 /******************************************************************************
316 Cosine window function
317 ******************************************************************************/
318 int win_cos(double *w, int n, int win_type)
319 {
320  double y;
321  double x = 0.0;
322  int i;
323 
324  if(!w)
325  return ERROR_PTR;
326  if(n<2)
327  return ERROR_SIZE;
328 
329  switch(win_type & DSPL_WIN_SYM_MASK)
330  {
331  case DSPL_WIN_SYMMETRIC: x = M_PI/(double)(n-1); break;
332  case DSPL_WIN_PERIODIC : x = M_PI/(double)n; break;
333  default: return ERROR_WIN_SYM;
334  }
335 
336  y = 0.0;
337  for(i = 0; i<n; i++)
338  {
339  w[i] = sin(y);
340  y += x;
341  }
342  return RES_OK;
343 }
344 
345 
346 
347 
348 
349 
350 /******************************************************************************
351 Flat - Top window function
352 ******************************************************************************/
353 int win_flat_top(double *w, int n, int win_type)
354 {
355  double y;
356  double x = 0.0;
357  double a0 = 1.0;
358  double a1 = 1.93;
359  double a2 = 1.29;
360  double a3 = 0.388;
361  double a4 = 0.032;
362  int i;
363 
364  if(!w)
365  return ERROR_PTR;
366  if(n<2)
367  return ERROR_SIZE;
368 
369  switch(win_type & DSPL_WIN_SYM_MASK)
370  {
371  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
372  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
373  default: return ERROR_WIN_SYM;
374  }
375 
376  y = 0.0;
377  for(i = 0; i<n; i++)
378  {
379  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y)+a4*cos(4.0*y);
380  y += x;
381  }
382  return RES_OK;
383 }
384 
385 
386 
387 
388 
389 
390 /******************************************************************************
391 Gaussian window function
392 ******************************************************************************/
393 int win_gaussian(double *w, int n, int win_type, double alpha)
394 {
395  double x = 0.0;
396  double y;
397  double sigma;
398  int i;
399 
400  if(!w)
401  return ERROR_PTR;
402  if(n<2)
403  return ERROR_SIZE;
404 
405  switch(win_type & DSPL_WIN_SYM_MASK)
406  {
407  case DSPL_WIN_SYMMETRIC: x = (double)(n-1)*0.5; break;
408  case DSPL_WIN_PERIODIC : x = (double)(n)*0.5; break;
409  default: return ERROR_WIN_SYM;
410  }
411 
412 
413  sigma = alpha / x;
414  for(i = 0; i<n; i++)
415  {
416  y = ((double)i - x)*sigma;
417  w[i] = exp(-0.5*y*y);
418  }
419  return RES_OK;
420 }
421 
422 
423 
424 
425 
426 
427 /******************************************************************************
428 Hamming window function
429 ******************************************************************************/
430 int win_hamming(double *w, int n, int win_type)
431 {
432  double x = 0.0;
433  double y;
434  int i;
435 
436  if(!w)
437  return ERROR_PTR;
438  if(n<2)
439  return ERROR_SIZE;
440 
441  switch(win_type & DSPL_WIN_SYM_MASK)
442  {
443  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
444  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
445  default: return ERROR_WIN_SYM;
446  }
447 
448  y = 0.0;
449  for(i = 0; i < n; i++)
450  {
451  w[i] = 0.54-0.46*cos(y);
452  y += x;
453  }
454  return RES_OK;
455 }
456 
457 
458 
459 
460 /******************************************************************************
461 Hann window function
462 ******************************************************************************/
463 int win_hann(double *w, int n, int win_type)
464 {
465  double x = 0.0;
466  double y;
467  int i;
468 
469  if(!w)
470  return ERROR_PTR;
471  if(n<2)
472  return ERROR_SIZE;
473 
474  switch(win_type & DSPL_WIN_SYM_MASK)
475  {
476  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
477  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
478  default: return ERROR_WIN_SYM;
479  }
480 
481  y = 0.0;
482  for(i = 0; i < n; i++)
483  {
484  w[i] = 0.5*(1-cos(y));
485  y += x;
486  }
487  return RES_OK;
488 }
489 
490 
491 /******************************************************************************
492 Kaiser window function
493 ******************************************************************************/
494 int win_kaiser(double* w, int n, int win_type, double param)
495 {
496  double num, den, x, y;
497  int i, err;
498  if(!w)
499  return ERROR_PTR;
500  if(n<2)
501  return ERROR_SIZE;
502 
503  switch(win_type & DSPL_WIN_SYM_MASK)
504  {
505  case DSPL_WIN_SYMMETRIC: x = 1.0/(double)(n-1); break;
506  case DSPL_WIN_PERIODIC : x = 1.0/(double)n; break;
507  default: return ERROR_WIN_SYM;
508  }
509 
510  err = bessel_i0(&param, 1, &den);
511  if(err != RES_OK)
512  return err;
513  for(i = 0; i < n; i++)
514  {
515  y = (double)(2*i) / x - 1.0;
516  y = param * sqrt(1.0 - y*y);
517  err = bessel_i0(&y, 1, &num);
518  if(err != RES_OK)
519  return err;
520  w[i] = num / den;
521  }
522  return err;
523 }
524 
525 
526 
527 /******************************************************************************
528 Lanczos window function
529 ******************************************************************************/
530 int win_lanczos(double *w, int n, int win_type)
531 {
532  double y;
533  double x = 0.0;
534  int i;
535 
536  if(!w)
537  return ERROR_PTR;
538  if(n<2)
539  return ERROR_SIZE;
540 
541  switch(win_type & DSPL_WIN_SYM_MASK)
542  {
543  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
544  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
545  default: return ERROR_WIN_SYM;
546  }
547 
548  y = 0.0;
549  for(i = 0; i < n; i++)
550  {
551  if((y - M_PI)==0.0)
552  w[i] = 1.0;
553  else
554  w[i] = sin(y - M_PI)/(y - M_PI);
555  y += x;
556  }
557  return RES_OK;
558 
559 }
560 
561 
562 
563 /******************************************************************************
564 Nuttall window function
565 ******************************************************************************/
566 int win_nuttall(double *w, int n, int win_type)
567 {
568  double y;
569  double x = 0.0;
570  double a0 = 0.355768;
571  double a1 = 0.487396;
572  double a2 = 0.144232;
573  double a3 = 0.012604;
574  int i;
575 
576  if(!w)
577  return ERROR_PTR;
578  if(n<2)
579  return ERROR_SIZE;
580 
581  switch(win_type & DSPL_WIN_SYM_MASK)
582  {
583  case DSPL_WIN_SYMMETRIC: x = M_2PI/(double)(n-1); break;
584  case DSPL_WIN_PERIODIC : x = M_2PI/(double)n; break;
585  default: return ERROR_WIN_SYM;
586  }
587 
588  y = 0.0;
589  for(i = 0; i < n; i++)
590  {
591  w[i] = a0 - a1* cos(y)+a2*cos(2.0*y)-a3*cos(3.0*y);
592  y += x;
593  }
594  return RES_OK;
595 }
596 
597 
598 
599 
600 
601 /******************************************************************************
602 Rectangle window function
603 ******************************************************************************/
604 int win_rect(double *w, int n)
605 {
606  int i;
607 
608  if(!w)
609  return ERROR_PTR;
610  if(n<2)
611  return ERROR_SIZE;
612 
613  for(i = 0; i < n; i++)
614  w[i] = 1.0;
615  return RES_OK;
616 }
617 
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:148
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:140
int bessel_i0(double *x, int n, double *y)
Модифицированная функция Бесселя первого рода .
Definition: math.c:36
int window(double *w, int n, int win_type, double param)
Расчет функции оконного взвешивания
Definition: win.c:33
int cheby_poly1(double *x, int n, int ord, double *y)
Многочлен Чебышева первого рода порядка ord
Definition: cheby.c:78
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:94