libdspl-2.0
goertzel.c
1 /*
2 * Copyright (c) 2015-2018 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 #include <stdlib.h>
22 #include <string.h>
23 #include "dspl.h"
24 
25 
26 /******************************************************************************
27 Goertzel algorithm for real vector
28 *******************************************************************************/
29 int DSPL_API goertzel(double *x, int n, int *ind, int k, complex_t *y)
30 {
31 
32  int m, p;
33  double wR, wI;
34  double alpha;
35  double v[3];
36 
37  if(!x || !y || !ind)
38  return ERROR_PTR;
39 
40  if(n < 1 || k < 1)
41  return ERROR_SIZE;
42 
43  for(p = 0; p < k; p++)
44  {
45  wR = cos(M_2PI * (double)ind[p] / (double)n);
46  wI = sin(M_2PI * (double)ind[p] / (double)n);
47 
48  alpha = 2.0 * wR;
49  v[0] = v[1] = v[2] = 0.0;
50 
51  for(m = 0; m < n; m++)
52  {
53  v[2] = v[1];
54  v[1] = v[0];
55  v[0] = x[m]+alpha*v[1] - v[2];
56  }
57  RE(y[p]) = wR * v[0] - v[1];
58  IM(y[p]) = wI * v[0];
59  }
60  return RES_OK;
61 }
62 
63 
64 
65 
66 
67 /******************************************************************************
68 Goertzel algorithm for complex vector
69 *******************************************************************************/
70 int DSPL_API goertzel_cmplx(complex_t *x, int n, int *ind, int k, complex_t *y)
71 {
72 
73  int m, p;
74  complex_t w;
75  double alpha;
76  complex_t v[3];
77 
78  if(!x || !y || !ind)
79  return ERROR_PTR;
80 
81  if(n < 1 || k < 1)
82  return ERROR_SIZE;
83 
84  for(p = 0; p < k; p++)
85  {
86  RE(w) = cos(M_2PI * (double)ind[p] / (double)n);
87  IM(w) = sin(M_2PI * (double)ind[p] / (double)n);
88 
89  alpha = 2.0 * RE(w);
90  memset(v, 0, 3*sizeof(complex_t));
91 
92  for(m = 0; m < n; m++)
93  {
94  RE(v[2]) = RE(v[1]);
95  RE(v[1]) = RE(v[0]);
96  RE(v[0]) = RE(x[m]) + alpha * RE(v[1]) - RE(v[2]);
97 
98  IM(v[2]) = IM(v[1]);
99  IM(v[1]) = IM(v[0]);
100  IM(v[0]) = IM(x[m]) + alpha * IM(v[1]) - IM(v[2]);
101  }
102 
103  RE(y[p]) = CMRE(w, v[0]) - RE(v[1]);
104  IM(y[p]) = CMIM(w, v[0]) - IM(v[1]);
105  }
106 
107  return RES_OK;
108 }
109 
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:124
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:41
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:117
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:66
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:65
int goertzel_cmplx(complex_t *x, int n, int *ind, int k, complex_t *y)
Алгоритм Герцеля для расчета отдельных спектральных отсчетов дискретного преобразования Фурье комплек...
Definition: goertzel.c:70
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:81
int goertzel(double *x, int n, int *ind, int k, complex_t *y)
Алгоритм Герцеля для расчета отдельных спектральных отсчетов дискретного преобразования Фурье веществ...
Definition: goertzel.c:29