libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
fourier_series.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 #include <stdlib.h>
23 #include <string.h>
24 #include <math.h>
25 #include "dspl.h"
26 
27 
28 /*******************************************************************************
29 Fourier Series Decomposition
30 *******************************************************************************/
31 int DSPL_API fourier_series_dec(double* t, double* s, int nt, double period,
32  int nw, double* w, complex_t* y)
33 {
34  int k, m;
35  double dw = M_2PI / period;
36  complex_t e[2];
37 
38  if(!t || !s || !w || !y)
39  return ERROR_PTR;
40  if(nt<1 || nw < 1)
41  return ERROR_SIZE;
42  if(period <= 0.0)
43  return ERROR_NEGATIVE;
44 
45  memset(y, 0 , nw*sizeof(complex_t));
46 
47  for(k = 0; k < nw; k++)
48  {
49  w[k] = (k - nw/2) * dw;
50  RE(e[1]) = s[0] * cos(w[k] * t[0]);
51  IM(e[1]) = -s[0] * sin(w[k] * t[0]);
52  for(m = 1; m < nt; m++)
53  {
54  RE(e[0]) = RE(e[1]);
55  IM(e[0]) = IM(e[1]);
56  RE(e[1]) = s[m] * cos(w[k] * t[m]);
57  IM(e[1]) = - s[m] * sin(w[k] * t[m]);
58  RE(y[k]) += 0.5 * (RE(e[0]) + RE(e[1]))*(t[m] - t[m-1]);
59  IM(y[k]) += 0.5 * (IM(e[0]) + IM(e[1]))*(t[m] - t[m-1]);
60  }
61  RE(y[k]) /= period;
62  IM(y[k]) /= period;
63  }
64 
65  if(!(nw%2))
66  RE(y[0]) = RE(y[1]) = 0.0;
67 
68  return RES_OK;
69 }
70 
71 
72 
73 /*******************************************************************************
74 Fourier Series Decomposition for complex input signal
75 *******************************************************************************/
76 int DSPL_API fourier_series_dec_cmplx(double* t, complex_t* s, int nt,
77  double period, int nw, double* w, complex_t* y)
78 {
79  int k, m;
80  double dw = M_2PI / period;
81  complex_t e[2];
82 
83  if(!t || !s || !w || !y)
84  return ERROR_PTR;
85  if(nt<1 || nw < 1)
86  return ERROR_SIZE;
87  if(period <= 0.0)
88  return ERROR_NEGATIVE;
89 
90  memset(y, 0 , nw*sizeof(complex_t));
91 
92  for(k = 0; k < nw; k++)
93  {
94  w[k] = (k - nw/2) * dw;
95  RE(e[1]) = RE(s[0]) * cos(w[k] * t[0]) +
96  IM(s[0]) * sin(w[k] * t[0]);
97  IM(e[1]) = -RE(s[0]) * sin(w[k] * t[0]) +
98  IM(s[0]) * cos(w[k] * t[0]);
99  for(m = 1; m < nt; m++)
100  {
101  RE(e[0]) = RE(e[1]);
102  IM(e[0]) = IM(e[1]);
103  RE(e[1]) = RE(s[m]) * cos(w[k] * t[m]) +
104  IM(s[m]) * sin(w[k] * t[m]);
105  IM(e[1]) = -RE(s[m]) * sin(w[k] * t[m]) +
106  IM(s[m]) * cos(w[k] * t[m]);
107  RE(y[k]) += 0.5 * (RE(e[0]) + RE(e[1]))*(t[m] - t[m-1]);
108  IM(y[k]) += 0.5 * (IM(e[0]) + IM(e[1]))*(t[m] - t[m-1]);
109  }
110  RE(y[k]) /= period;
111  IM(y[k]) /= period;
112  }
113 
114  if(!(nw%2))
115  RE(y[0]) = RE(y[1]) = 0.0;
116 
117  return RES_OK;
118 }
119 
120 
121 
122 
123 
124 
125 /*******************************************************************************
126 Fourier Transform
127 *******************************************************************************/
128 int DSPL_API fourier_integral_cmplx(double* t, complex_t* s, int nt,
129  int nw, double* w, complex_t* y)
130 {
131  int k, m;
132  complex_t e[2];
133 
134  if(!t || !s || !w || !y)
135  return ERROR_PTR;
136  if(nt<1 || nw < 1)
137  return ERROR_SIZE;
138 
139 
140  memset(y, 0 , nw*sizeof(complex_t));
141 
142  for(k = 0; k < nw; k++)
143  {
144  RE(e[1]) = RE(s[0]) * cos(w[k] * t[0]) +
145  IM(s[0]) * sin(w[k] * t[0]);
146  IM(e[1]) = -RE(s[0]) * sin(w[k] * t[0]) +
147  IM(s[0]) * cos(w[k] * t[0]);
148  for(m = 1; m < nt; m++)
149  {
150  RE(e[0]) = RE(e[1]);
151  IM(e[0]) = IM(e[1]);
152  RE(e[1]) = RE(s[m]) * cos(w[k] * t[m]) +
153  IM(s[m]) * sin(w[k] * t[m]);
154  IM(e[1]) = -RE(s[m]) * sin(w[k] * t[m]) +
155  IM(s[m]) * cos(w[k] * t[m]);
156  RE(y[k]) += 0.5 * (RE(e[0]) + RE(e[1]))*(t[m] - t[m-1]);
157  IM(y[k]) += 0.5 * (IM(e[0]) + IM(e[1]))*(t[m] - t[m-1]);
158  }
159  }
160 
161  return RES_OK;
162 }
163 
164 
165 
166 
167 
168 
169 /*******************************************************************************
170 Fourier Series Reconstruction
171 *******************************************************************************/
172 int DSPL_API fourier_series_rec(double* w, complex_t* s, int nw,
173  double *t, int nt, complex_t* y)
174 {
175  int k, m;
176  complex_t e;
177 
178  if(!t || !s || !w || !y)
179  return ERROR_PTR;
180  if(nt<1 || nw < 1)
181  return ERROR_SIZE;
182 
183  memset(y, 0, nt*sizeof(complex_t));
184 
185 
186  for(k = 0; k < nw; k++)
187  {
188  for(m = 0; m < nt; m++)
189  {
190  RE(e) = cos(w[k] * t[m]);
191  IM(e) = sin(w[k] * t[m]);
192 
193  RE(y[m]) += CMRE(s[k], e);
194  IM(y[m]) += CMIM(s[k], e);
195  }
196  }
197  return RES_OK;
198 }
199 
ERROR_SIZE
#define ERROR_SIZE
Ошибка при передаче размера массива.
Definition: dspl.h:125
complex_t
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:41
fourier_series_dec
int fourier_series_dec(double *t, double *s, int nt, double period, int nw, double *w, complex_t *y)
Расчет коэффициентов разложения в ряд Фурье
Definition: fourier_series.c:31
ERROR_PTR
#define ERROR_PTR
Ошибка указателя.
Definition: dspl.h:118
IM
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:65
RE
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:64
RES_OK
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:81
Документация по libdspl-2.0. Последние изменения: Вс 17 Фев 2019 16:20:10. Создано системой   doxygen 1.8.15
© Бахурин Сергей 2015 - 2019. Все права защищены. Любое копирование материалов сайте без разрешения автора запрещено.