libdspl-2.0
Digital Signal Processing Algorithm Library
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 #ifdef DOXYGEN_ENGLISH
29 
85 #endif
86 #ifdef DOXYGEN_RUSSIAN
87 
150 #endif
151 int DSPL_API fourier_series_dec(double* t, double* s, int nt, double period,
152  int nw, double* w, complex_t* y)
153 {
154  int k, m;
155  double dw = M_2PI / period;
156  complex_t e[2];
157 
158  if(!t || !s || !w || !y)
159  return ERROR_PTR;
160  if(nt<1 || nw < 1)
161  return ERROR_SIZE;
162  if(period <= 0.0)
163  return ERROR_NEGATIVE;
164 
165  memset(y, 0 , nw*sizeof(complex_t));
166 
167  for(k = 0; k < nw; k++)
168  {
169  w[k] = (k - nw/2) * dw;
170  RE(e[1]) = s[0] * cos(w[k] * t[0]);
171  IM(e[1]) = -s[0] * sin(w[k] * t[0]);
172  for(m = 1; m < nt; m++)
173  {
174  RE(e[0]) = RE(e[1]);
175  IM(e[0]) = IM(e[1]);
176  RE(e[1]) = s[m] * cos(w[k] * t[m]);
177  IM(e[1]) = - s[m] * sin(w[k] * t[m]);
178  RE(y[k]) += 0.5 * (RE(e[0]) + RE(e[1]))*(t[m] - t[m-1]);
179  IM(y[k]) += 0.5 * (IM(e[0]) + IM(e[1]))*(t[m] - t[m-1]);
180  }
181  RE(y[k]) /= period;
182  IM(y[k]) /= period;
183  }
184 
185  if(!(nw%2))
186  RE(y[0]) = RE(y[1]) = 0.0;
187 
188  return RES_OK;
189 }
190 
191 
192 
193 #ifdef DOXYGEN_ENGLISH
194 
195 #endif
196 #ifdef DOXYGEN_RUSSIAN
197 
198 #endif
199 int DSPL_API fourier_series_dec_cmplx(double* t, complex_t* s, int nt,
200  double period, int nw, double* w, complex_t* y)
201 {
202  int k, m;
203  double dw = M_2PI / period;
204  complex_t e[2];
205 
206  if(!t || !s || !w || !y)
207  return ERROR_PTR;
208  if(nt<1 || nw < 1)
209  return ERROR_SIZE;
210  if(period <= 0.0)
211  return ERROR_NEGATIVE;
212 
213  memset(y, 0 , nw*sizeof(complex_t));
214 
215  for(k = 0; k < nw; k++)
216  {
217  w[k] = (k - nw/2) * dw;
218  RE(e[1]) = RE(s[0]) * cos(w[k] * t[0]) +
219  IM(s[0]) * sin(w[k] * t[0]);
220  IM(e[1]) = -RE(s[0]) * sin(w[k] * t[0]) +
221  IM(s[0]) * cos(w[k] * t[0]);
222  for(m = 1; m < nt; m++)
223  {
224  RE(e[0]) = RE(e[1]);
225  IM(e[0]) = IM(e[1]);
226  RE(e[1]) = RE(s[m]) * cos(w[k] * t[m]) +
227  IM(s[m]) * sin(w[k] * t[m]);
228  IM(e[1]) = -RE(s[m]) * sin(w[k] * t[m]) +
229  IM(s[m]) * cos(w[k] * t[m]);
230  RE(y[k]) += 0.5 * (RE(e[0]) + RE(e[1]))*(t[m] - t[m-1]);
231  IM(y[k]) += 0.5 * (IM(e[0]) + IM(e[1]))*(t[m] - t[m-1]);
232  }
233  RE(y[k]) /= period;
234  IM(y[k]) /= period;
235  }
236 
237  if(!(nw%2))
238  RE(y[0]) = RE(y[1]) = 0.0;
239 
240  return RES_OK;
241 }
242 
243 
244 
245 
246 
247 #ifdef DOXYGEN_ENGLISH
248 
249 #endif
250 #ifdef DOXYGEN_RUSSIAN
251 
252 #endif
253 int DSPL_API fourier_integral_cmplx(double* t, complex_t* s, int nt,
254  int nw, double* w, complex_t* y)
255 {
256  int k, m;
257  complex_t e[2];
258 
259  if(!t || !s || !w || !y)
260  return ERROR_PTR;
261  if(nt<1 || nw < 1)
262  return ERROR_SIZE;
263 
264 
265  memset(y, 0 , nw*sizeof(complex_t));
266 
267  for(k = 0; k < nw; k++)
268  {
269  RE(e[1]) = RE(s[0]) * cos(w[k] * t[0]) +
270  IM(s[0]) * sin(w[k] * t[0]);
271  IM(e[1]) = -RE(s[0]) * sin(w[k] * t[0]) +
272  IM(s[0]) * cos(w[k] * t[0]);
273  for(m = 1; m < nt; m++)
274  {
275  RE(e[0]) = RE(e[1]);
276  IM(e[0]) = IM(e[1]);
277  RE(e[1]) = RE(s[m]) * cos(w[k] * t[m]) +
278  IM(s[m]) * sin(w[k] * t[m]);
279  IM(e[1]) = -RE(s[m]) * sin(w[k] * t[m]) +
280  IM(s[m]) * cos(w[k] * t[m]);
281  RE(y[k]) += 0.5 * (RE(e[0]) + RE(e[1]))*(t[m] - t[m-1]);
282  IM(y[k]) += 0.5 * (IM(e[0]) + IM(e[1]))*(t[m] - t[m-1]);
283  }
284  }
285 
286  return RES_OK;
287 }
288 
289 
290 
291 
292 
293 
294 #ifdef DOXYGEN_ENGLISH
295 
352 #endif
353 #ifdef DOXYGEN_RUSSIAN
354 
414 #endif
415 int DSPL_API fourier_series_rec(double* w, complex_t* s, int nw,
416  double* t, int nt, complex_t* y)
417 {
418  int k, m;
419  complex_t e;
420 
421  if(!t || !s || !w || !y)
422  return ERROR_PTR;
423  if(nt<1 || nw < 1)
424  return ERROR_SIZE;
425 
426  memset(y, 0, nt*sizeof(complex_t));
427 
428 
429  for(k = 0; k < nw; k++)
430  {
431  for(m = 0; m < nt; m++)
432  {
433  RE(e) = cos(w[k] * t[m]);
434  IM(e) = sin(w[k] * t[m]);
435 
436  RE(y[m]) += CMRE(s[k], e);
437  IM(y[m]) += CMIM(s[k], e);
438  }
439  }
440  return RES_OK;
441 }
442 
int DSPL_API fourier_series_rec(double *w, complex_t *s, int nw, double *t, int nt, complex_t *y)
Time signal reconstruction from Fourier series coefficients.
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:359
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:545
#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:553
int DSPL_API fourier_series_dec(double *t, double *s, int nt, double period, int nw, double *w, complex_t *y)
Fourier series coefficient calculation for periodic signal.
double complex_t[2]
Complex data type.
Definition: dspl.h:86
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:497
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:417
#define ERROR_NEGATIVE
Negative parameter. The function returns the given error code when it takes a negative parameter....
Definition: dspl.h:540