libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
complex.c
1 /*
2 * Copyright (c) 2015-2019 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of DSPL.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU 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 General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include "dspl.h"
25 
26 
27 
28 /******************************************************************************
29 \ingroup SPEC_MATH_TRIG_GROUP
30 \fn int acos_cmplx(complex_t* x, int n, complex_t *y)
31 \brief The inverse of the cosine function the complex vector argument `x`
32 
33 Function calculates the inverse of the cosine function as: \n
34 
35 \f[
36 \textrm{Arccos}(x) = \frac{\pi}{2} - \textrm{Arcsin}(x) =
37 \frac{\pi}{2} -j \textrm{Ln}\left( j x + \sqrt{1 - x^2} \right)
38 \f]
39 
40 
41 \param[in] x Pointer to the argument vector `x`. \n
42  Vector size is `[n x 1]`. \n \n
43 
44 \param[in] n Input vector `x` and the inverse cosine vector `y` size. \n \n
45 
46 
47 \param[out] y Pointer to the output complex vector `y`,
48  corresponds to the input vector `x`. \n
49  Vector size is `[n x 1]`. \n
50  Memory must be allocated. \n \n
51 
52 \return
53 `RES_OK` if function calculated successfully. \n
54 Else \ref ERROR_CODE_GROUP "code error". \n
55 
56 Example: \n
57 \code{.cpp}
58  complex_t x[3] = {{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
59  complex_t y[3];
60  int k;
61 
62  acos_cmplx(x, 3, y);
63 
64  for(k = 0; k < 3; k++)
65  printf("acos_cmplx(%.1f%+.1fj) = %.3f%+.3fj\n",
66  RE(x[k]), IM(x[k]), RE(y[k]), IM(y[k]));
67 \endcode
68  \n
69 
70 Output is: \n
71 \verbatim
72 acos_cmplx(1.0+2.0j) = 1.144-1.529j
73 acos_cmplx(3.0+4.0j) = 0.937-2.306j
74 acos_cmplx(5.0+6.0j) = 0.880-2.749j
75 \endverbatim
76 
77 \author
78 Sergey Bakhurin www.dsplib.org
79 *******************************************************************************/
80 int DSPL_API acos_cmplx(complex_t* x, int n, complex_t *y)
81 {
82  int k, res;
83  double pi2 = 0.5 * M_PI;
84 
85  res = asin_cmplx(x, n, y);
86  if(res != RES_OK)
87  return res;
88 
89  for(k = 0; k < n; k++)
90  {
91  RE(y[k]) = pi2 - RE(y[k]);
92  IM(y[k]) = - IM(y[k]);
93  }
94  return RES_OK;
95 }
96 
97 
98 
99 
100 /******************************************************************************
101 \ingroup SPEC_MATH_TRIG_GROUP
102 \fn int asin_cmplx(complex_t* x, int n, complex_t *y)
103 \brief The inverse of the sine function the complex vector argument `x`
104 
105 Function calculates the inverse of the sine function as: \n
106 
107 \f[
108  \textrm{Arcsin}(x) = j \textrm{Ln}\left( j x + \sqrt{1 - x^2} \right)
109 \f]
110 
111 
112 \param[in] x Pointer to the argument vector `x`. \n
113  Vector size is `[n x 1]`. \n \n
114 
115 \param[in] n Input vector `x` and the inverse sine vector `y` size. \n \n
116 
117 
118 \param[out] y Pointer to the output complex vector `y`,
119  corresponds to the input vector `x`. \n
120  Vector size is `[n x 1]`. \n
121  Memory must be allocated. \n \n
122 
123 \return
124 `RES_OK` if function calculated successfully. \n
125 Else \ref ERROR_CODE_GROUP "code error". \n
126 
127 Example: \n
128 \code{.cpp}
129  complex_t x[3] = {{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
130  complex_t y[3];
131  int k;
132 
133  asin_cmplx(x, 3, y);
134  for(k = 0; k < 3; k++)
135  printf("asin_cmplx(%.1f%+.1fj) = %.3f%+.3fj\n",
136  RE(x[k]), IM(x[k]), RE(y[k]), IM(y[k]));
137 
138 \endcode
139  \n
140 
141 Output is: \n
142 \verbatim
143 asin_cmplx(1.0+2.0j) = 0.427+1.529j
144 asin_cmplx(3.0+4.0j) = 0.634+2.306j
145 asin_cmplx(5.0+6.0j) = 0.691+2.749j
146 \endverbatim
147 
148 \author
149 Sergey Bakhurin www.dsplib.org
150 *******************************************************************************/
151 int DSPL_API asin_cmplx(complex_t* x, int n, complex_t *y)
152 {
153  int k;
154  complex_t tmp;
155  if(!x || !y)
156  return ERROR_PTR;
157  if(n < 1)
158  return ERROR_SIZE;
159 
160  for(k = 0; k < n; k++)
161  {
162  RE(tmp) = 1.0 - CMRE(x[k], x[k]); /* 1-x[k]^2 */
163  IM(tmp) = - CMIM(x[k], x[k]); /* 1-x[k]^2 */
164  sqrt_cmplx(&tmp, 1, y+k); /* sqrt(1 - x[k]^2) */
165  RE(y[k]) -= IM(x[k]); /* j * x[k] + sqrt(1 - x[k]^2) */
166  IM(y[k]) += RE(x[k]); /* j * x[k] + sqrt(1 - x[k]^2) */
167  log_cmplx(y+k, 1, &tmp); /* log( j * x[k] + sqrt(1 - x[k]^2) ) */
168  RE(y[k]) = IM(tmp); /* -j * log( j * x[k] + sqrt(1 - x[k]^2) ) */
169  IM(y[k]) = -RE(tmp); /* -j * log( j * x[k] + sqrt(1 - x[k]^2) ) */
170  }
171  return RES_OK;
172 }
173 
174 
175 
176 
177 /******************************************************************************
178 \ingroup TYPES_GROUP
179 \fn int cmplx2re(complex_t* x, int n, double* re, double* im)
180 \brief Separate complex vector to the real and image vectors
181 
182 Function fills `re` and `im` vectors corresponds to real and image
183 parts of the input complex array `x`. \n
184 
185 
186 \param[in] x Pointer to the real complex vector. \n
187  Vector size is `[n x 1]`. \n \n
188 
189 \param[in] n Size of the input complex vector `x` and real and image
190  vectors `re` and `im`. \n \n
191 
192 \param[out] re Pointer to the real part vector. \n
193  Vector size is `[n x 1]`. \n
194  Memory must be allocated. \n \n
195 
196 \param[out] im Pointer to the image part vector. \n
197  Vector size is `[n x 1]`. \n
198  Memory must be allocated. \n \n
199 
200 \return
201 `RES_OK` if function converts complex vector successfully. \n
202 Else \ref ERROR_CODE_GROUP "code error". \n
203 
204 Example: \n
205 \code{.cpp}
206  complex_t x[3] = {{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
207  double re[3], im[3];
208 
209  cmplx2re(x, 3, re, im);
210 \endcode
211 
212 Vectors `re` and `im` will contains:
213 
214 \verbatim
215 re[0] = 1.0; im[0] = 2.0;
216 re[1] = 3.0; im[1] = 4.0;
217 re[2] = 5.0; im[2] = 6.0;
218 \endverbatim
219 
220 \author Sergey Bakhurin. www.dsplib.org
221 *******************************************************************************/
222 int DSPL_API cmplx2re(complex_t* x, int n, double* re, double* im)
223 {
224  int k;
225  if(!x)
226  return ERROR_PTR;
227  if(n < 1)
228  return ERROR_SIZE;
229 
230  if(re)
231  {
232  for(k = 0; k < n; k++)
233  re[k] = RE(x[k]);
234  }
235  if(im)
236  {
237  for(k = 0; k < n; k++)
238  im[k] = IM(x[k]);
239  }
240  return RES_OK;
241 }
242 
243 
244 
245 
246 
247 /******************************************************************************
248 \ingroup SPEC_MATH_TRIG_GROUP
249 \fn int cos_cmplx(complex_t* x, int n, complex_t *y)
250 \brief The cosine function the complex vector argument `x`
251 
252 Function calculates the cosine function as: \n
253 
254 \f[
255 \textrm{cos}(x) = \frac{\exp(jx) + \exp(-jx)}{2}
256 \f]
257 
258 
259 \param[in] x Pointer to the argument vector `x`. \n
260  Vector size is `[n x 1]`. \n \n
261 
262 \param[in] n Input vector `x` and the cosine vector `y` size. \n \n
263 
264 
265 \param[out] y Pointer to the output complex vector `y`,
266  corresponds to the input vector `x`. \n
267  Vector size is `[n x 1]`. \n
268  Memory must be allocated. \n \n
269 
270 \return
271 `RES_OK` if function calculated successfully. \n
272 Else \ref ERROR_CODE_GROUP "code error". \n
273 
274 Example: \n
275 \code{.cpp}
276  complex_t x[3] = {{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
277  complex_t y[3];
278  int k;
279 
280  cos_cmplx(x, 3, y);
281 
282  for(k = 0; k < 3; k++)
283  printf("cos_cmplx(%.1f%+.1fj) = %9.3f%+9.3fj\n",
284  RE(x[k]), IM(x[k]), RE(y[k]), IM(y[k]));
285 
286 \endcode
287  \n
288 
289 Output is: \n
290 \verbatim
291 cos_cmplx(1.0+2.0j) = 2.033 -3.052j
292 cos_cmplx(3.0+4.0j) = -27.035 -3.851j
293 cos_cmplx(5.0+6.0j) = 57.219 +193.428j
294 \endverbatim
295 
296 \author
297 Sergey Bakhurin www.dsplib.org
298 *******************************************************************************/
299 int DSPL_API cos_cmplx(complex_t* x, int n, complex_t *y)
300 {
301  int k;
302  double ep, em, sx, cx;
303  if(!x || !y)
304  return ERROR_PTR;
305  if(n < 1)
306  return ERROR_SIZE;
307 
308  for(k = 0; k < n; k++)
309  {
310  ep = exp( IM(x[k]));
311  em = exp(-IM(x[k]));
312  sx = 0.5 * sin(RE(x[k]));
313  cx = 0.5 * cos(RE(x[k]));
314  RE(y[k]) = cx * (em + ep);
315  IM(y[k]) = sx * (em - ep);
316  }
317  return RES_OK;
318 }
319 
320 
321 
322 
323 /******************************************************************************
324 \ingroup SPEC_MATH_COMMON_GROUP
325 \fn int log_cmplx(complex_t* x, int n, complex_t *y)
326 \brief The logarithm function the complex vector argument `x`
327 
328 Function calculates the logarithm function as: \n
329 
330 \f[
331 \textrm{Ln}(x) = j \varphi + \ln(|x|),
332 \f]
333 here \f$\varphi\f$ - the complex number phase.
334 
335 \param[in] x Pointer to the argument vector `x`. \n
336  Vector size is `[n x 1]`. \n \n
337 
338 \param[in] n Input vector `x` and the logarithm vector `y` size. \n \n
339 
340 
341 \param[out] y Pointer to the output complex vector `y`,
342  corresponds to the input vector `x`. \n
343  Vector size is `[n x 1]`. \n
344  Memory must be allocated. \n \n
345 
346 \return
347 `RES_OK` if function calculated successfully. \n
348 Else \ref ERROR_CODE_GROUP "code error". \n
349 
350 Example: \n
351 \code{.cpp}
352  complex_t x[3] = {{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
353  complex_t y[3];
354  int k;
355 
356  log_cmplx(x, 3, y);
357 
358  for(k = 0; k < 3; k++)
359  printf("log_cmplx(%.1f%+.1fj) = %.3f%+.3fj\n",
360  RE(x[k]), IM(x[k]), RE(y[k]), IM(y[k]));
361 \endcode
362  \n
363 
364 Output is: \n
365 \verbatim
366 log_cmplx(1.0+2.0j) = 0.805+1.107j
367 log_cmplx(3.0+4.0j) = 1.609+0.927j
368 log_cmplx(5.0+6.0j) = 2.055+0.876j
369 \endverbatim
370 
371 \author
372 Sergey Bakhurin www.dsplib.org
373 *******************************************************************************/
374 int DSPL_API log_cmplx(complex_t* x, int n, complex_t *y)
375 {
376  int k;
377  if(!x || !y)
378  return ERROR_PTR;
379  if(n < 1)
380  return ERROR_SIZE;
381 
382  for(k = 0; k < n; k++)
383  {
384  RE(y[k]) = 0.5 * log(ABSSQR(x[k]));
385  IM(y[k]) = atan2(IM(x[k]), RE(x[k]));
386  }
387  return RES_OK;
388 }
389 
390 
391 
392 
393 
394 
395 
396 /******************************************************************************
397 \ingroup TYPES_GROUP
398 \fn int re2cmplx(double* x, int n, complex_t *y)
399 \brief Convert real array to the complex array.
400 
401 Function copies the vector `x` to the real part of vector `y`.
402 Image part of the vector `y` sets as zero. \n
403 So complex vector contains data: \n
404 `y[i] = x[i] + j0, here i = 0,1,2 ... n-1`
405 
406 
407 \param[in] x Pointer to the real vector `x`. \n
408  Vector size is `[n x 1]`. \n \n
409 
410 \param[in] n Size of the real vector `x` and complex vector `y`. \n \n
411 
412 \param[out] y Pointer to the complex vector `y`. \n
413  Vector size is `[n x 1]`. \n
414  Memory must be allocated. \n \n
415 
416 
417 \return
418 `RES_OK` if function returns successfully. \n
419 Else \ref ERROR_CODE_GROUP "code error": \n
420 
421 
422 
423 Например при выполнении следующего кода
424 \code{.cpp}
425  double x[3] = {1.0, 2.0, 3.0};
426  complex_t y[3];
427 
428  re2cmplx(x, 3, y);
429 \endcode
430 
431 Vector `y` will keep:
432 
433 \verbatim
434  y[0] = 1+0j;
435  y[1] = 2+0j;
436  y[2] = 3+0j.
437 \endverbatim
438 
439 \author Sergey Bakhurin. www.dsplib.org
440 *******************************************************************************/
441 int DSPL_API re2cmplx(double* x, int n, complex_t* y)
442 {
443  int k;
444  if(!x || !y)
445  return ERROR_PTR;
446  if(n < 1)
447  return ERROR_SIZE;
448 
449  for(k = 0; k < n; k++)
450  {
451  RE(y[k]) = x[k];
452  IM(y[k]) = 0.0;
453  }
454  return RES_OK;
455 }
456 
457 
458 
459 
460 
461 /******************************************************************************
462 \ingroup SPEC_MATH_TRIG_GROUP
463 \fn int sin_cmplx(complex_t* x, int n, complex_t *y)
464 \brief The sine function the complex vector argument `x`
465 
466 Function calculates the sine function as: \n
467 
468 \f[
469 \textrm{cos}(x) = \frac{\exp(jx) - \exp(-jx)}{2j}
470 \f]
471 
472 
473 \param[in] x Pointer to the argument vector `x`. \n
474  Vector size is `[n x 1]`. \n \n
475 
476 \param[in] n Input vector `x` and the sine vector `y` size. \n \n
477 
478 
479 \param[out] y Pointer to the output complex vector `y`,
480  corresponds to the input vector `x`. \n
481  Vector size is `[n x 1]`. \n
482  Memory must be allocated. \n \n
483 
484 \return
485 `RES_OK` if function calculated successfully. \n
486 Else \ref ERROR_CODE_GROUP "code error". \n
487 
488 Example: \n
489 \code{.cpp}
490  complex_t x[3] = {{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
491  complex_t y[3];
492  int k;
493 
494  sin_cmplx(x, 3, y);
495 
496  for(k = 0; k < 3; k++)
497  printf("sin_cmplx(%.1f%+.1fj) = %9.3f%+9.3fj\n",
498  RE(x[k]), IM(x[k]), RE(y[k]), IM(y[k]));
499 
500 \endcode
501  \n
502 
503 Output is: \n
504 \verbatim
505 sin_cmplx(1.0+2.0j) = 3.166 +1.960j
506 sin_cmplx(3.0+4.0j) = 3.854 -27.017j
507 sin_cmplx(5.0+6.0j) = -193.430 +57.218j
508 \endverbatim
509 
510 \author
511 Sergey Bakhurin www.dsplib.org
512 *******************************************************************************/
513 int DSPL_API sin_cmplx(complex_t* x, int n, complex_t *y)
514 {
515  int k;
516  double ep, em, sx, cx;
517  if(!x || !y)
518  return ERROR_PTR;
519  if(n < 1)
520  return ERROR_SIZE;
521 
522  for(k = 0; k < n; k++)
523  {
524  ep = exp( IM(x[k]));
525  em = exp(-IM(x[k]));
526  sx = 0.5 * sin(RE(x[k]));
527  cx = 0.5 * cos(RE(x[k]));
528  RE(y[k]) = sx * (em + ep);
529  IM(y[k]) = cx * (ep - em);
530  }
531  return RES_OK;
532 }
533 
534 
535 
536 
537 
538 
539 /******************************************************************************
540 \ingroup SPEC_MATH_COMMON_GROUP
541 \fn int sqrt_cmplx(complex_t* x, int n, complex_t *y)
542 \brief Square root of the complex vector argguument `x`.
543 
544 Function calculates square root value of vector `x` length `n`: \n
545 \f[
546 y(k) = \sqrt{x(k)}, \qquad k = 0 \ldots n-1.
547 \f]
548 
549 
550 \param[in] x Pointer to the input complex vector `x`. \n
551  Vector size is `[n x 1]`. \n \n
552 
553 \param[in] n Size of input and output vectors `x` and `y`. \n \n
554 
555 
556 \param[out] y Pointer to the square root vector `y`. \n
557  Vector size is `[n x 1]`. \n
558  Memory must be allocated. \n \n
559 
560 \return `RES_OK` if function is calculated successfully. \n
561 Else \ref ERROR_CODE_GROUP "code error". \n
562 
563 Example
564 \code{.cpp}
565  complex_t x[3] = {{1.0, 2.0}, {3.0, 4.0}, {5.0, 6.0}};
566  complex_t y[3]
567  int k;
568 
569  sqrt_cmplx(x, 3, y);
570 
571  for(k = 0; k < 3; k++)
572  printf("sqrt_cmplx(%.1f%+.1fj) = %.3f%+.3fj\n",
573  RE(x[k]), IM(x[k]), RE(y[k]), IM(y[k]));
574 
575  \endcode
576  \n
577 
578 Результатом работы будет
579 
580 \verbatim
581 sqrt_cmplx(1.0+2.0j) = 1.272+0.786j
582 sqrt_cmplx(3.0+4.0j) = 2.000+1.000j
583 sqrt_cmplx(5.0+6.0j) = 2.531+1.185j
584 \endverbatim
585 
586 \author Sergey Bakhurin www.dsplib.org
587 *******************************************************************************/
588 int DSPL_API sqrt_cmplx(complex_t* x, int n, complex_t *y)
589 {
590  int k;
591  double r, zr, at;
592  complex_t t;
593  if(!x || !y)
594  return ERROR_PTR;
595  if(n < 1)
596  return ERROR_SIZE;
597 
598  for(k = 0; k < n; k++)
599  {
600  r = ABS(x[k]);
601  if(r == 0.0)
602  {
603  RE(y[k]) = 0.0;
604  IM(y[k]) = 0.0;
605  }
606  else
607  {
608  RE(t) = RE(x[k]) + r;
609  IM(t) = IM(x[k]);
610  at = ABS(t);
611  if(at == 0.0)
612  {
613  RE(y[k]) = 0.0;
614  IM(y[k]) = sqrt(r);
615  }
616  else
617  {
618  zr = 1.0 / ABS(t);
619  r = sqrt(r);
620  RE(y[k]) = RE(t) * zr * r;
621  IM(y[k]) = IM(t) * zr * r;
622  }
623  }
624  }
625  return RES_OK;
626 }
627 
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:148
double complex_t[2]
Описание комплексного типа данных.
Definition: dspl.h:41
int log_cmplx(complex_t *x, int n, complex_t *y)
Натуральный логарифм комплексного аргумента x
Definition: complex.c:374
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:140
int sin_cmplx(complex_t *x, int n, complex_t *y)
Синус комплексного аргумента x
Definition: complex.c:513
#define IM(x)
Макрос определяющий мнимую часть комплексного числа.
Definition: dspl.h:78
#define ABSSQR(x)
Макрос возвращает квадрат модуля комплексного числа x.
Definition: dspl.h:82
#define RE(x)
Макрос определяющий реальную часть комплексного числа.
Definition: dspl.h:77
int re2cmplx(double *x, int n, complex_t *y)
Преобразование массива вещественных данных в массив комплексных данных.
Definition: complex.c:441
int cmplx2re(complex_t *x, int n, double *re, double *im)
Преобразование массива комплексных данных в два массива вещественных данных, содержащих реальную и мн...
Definition: complex.c:222
int sqrt_cmplx(complex_t *x, int n, complex_t *y)
Квадратный корень из комплексного вектора x (поэлементный).
Definition: complex.c:588
int cos_cmplx(complex_t *x, int n, complex_t *y)
Косинус комплексного аргумента x
Definition: complex.c:299
int acos_cmplx(complex_t *x, int n, complex_t *y)
Арккосинус комплексного аргумента x
Definition: complex.c:80
int asin_cmplx(complex_t *x, int n, complex_t *y)
Арксинус комплексного аргумента x
Definition: complex.c:151
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:94