libdspl-2.0
fourier_series.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 
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 /*******************************************************************************
127 Fourier Series Reconstruction
128 *******************************************************************************/
129 int DSPL_API fourier_series_rec(double* w, complex_t* s, int nw,
130  double *t, int nt, complex_t* y)
131 {
132  int k, m;
133  complex_t e;
134 
135  if(!t || !s || !w || !y)
136  return ERROR_PTR;
137  if(nt<1 || nw < 1)
138  return ERROR_SIZE;
139 
140  memset(y, 0, nt*sizeof(complex_t));
141 
142 
143  for(k = 0; k < nw; k++)
144  {
145  for(m = 0; m < nt; m++)
146  {
147  RE(e) = cos(w[k] * t[m]);
148  IM(e) = sin(w[k] * t[m]);
149 
150  RE(y[m]) += CMRE(s[k], e);
151  IM(y[m]) += CMIM(s[k], e);
152  }
153  }
154  return RES_OK;
155 }
156 
ERROR_SIZE
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:123
complex_t
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:40
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:116
IM
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:65
RE
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:64
fourier_series_rec
int fourier_series_rec(double *w, complex_t *s, int nw, double *t, int nt, complex_t *y)
Восстановление сигнала при усечении ряда Фурье
Definition: fourier_series.c:129
RES_OK
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:80
Документация по libdspl-2.0. Последние изменения: Вт 16 Окт 2018 19:20:41. Создано системой   doxygen 1.8.13