libdspl-2.0
Digital Signal Processing Algorithm Library
conv.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 #include <stdlib.h>
22 #include <string.h>
23 #include "dspl.h"
24 
25 
26 
27 #ifdef DOXYGEN_ENGLISH
28 
90 #endif
91 #ifdef DOXYGEN_RUSSIAN
92 
156 #endif
157 int DSPL_API conv(double* a, int na, double* b, int nb, double* c)
158 {
159  int k;
160  int n;
161 
162  double *t;
163  size_t bufsize;
164 
165  if(!a || !b || !c)
166  return ERROR_PTR;
167  if(na < 1 || nb < 1)
168  return ERROR_SIZE;
169 
170 
171  bufsize = (na + nb - 1) * sizeof(double);
172 
173  if((a != c) && (b != c))
174  t = c;
175  else
176  t = (double*)malloc(bufsize);
177 
178  memset(t, 0, bufsize);
179 
180  for(k = 0; k < na; k++)
181  for(n = 0; n < nb; n++)
182  t[k+n] += a[k]*b[n];
183 
184  if(t!=c)
185  {
186  memcpy(c, t, bufsize);
187  free(t);
188  }
189  return RES_OK;
190 }
191 
192 
193 
194 
195 
196 
197 #ifdef DOXYGEN_ENGLISH
198 
260 #endif
261 #ifdef DOXYGEN_RUSSIAN
262 
326 #endif
327 int DSPL_API conv_cmplx(complex_t* a, int na, complex_t* b,
328  int nb, complex_t* c)
329 {
330  int k;
331  int n;
332 
333  complex_t *t;
334  size_t bufsize;
335 
336  if(!a || !b || !c)
337  return ERROR_PTR;
338  if(na < 1 || nb < 1)
339  return ERROR_SIZE;
340 
341  bufsize = (na + nb - 1) * sizeof(complex_t);
342 
343  if((a != c) && (b != c))
344  t = c;
345  else
346  t = (complex_t*)malloc(bufsize);
347 
348  memset(t, 0, bufsize);
349 
350  for(k = 0; k < na; k++)
351  {
352  for(n = 0; n < nb; n++)
353  {
354  RE(t[k+n]) += CMRE(a[k], b[n]);
355  IM(t[k+n]) += CMIM(a[k], b[n]);
356  }
357  }
358 
359  if(t!=c)
360  {
361  memcpy(c, t, bufsize);
362  free(t);
363  }
364 
365  return RES_OK;
366 }
367 
368 
369 
370 
371 
372 #ifdef DOXYGEN_ENGLISH
373 
451 #endif
452 #ifdef DOXYGEN_RUSSIAN
453 
542 #endif
543 int DSPL_API conv_fft(double* a, int na, double* b, int nb,
544  fft_t* pfft, int nfft, double* c)
545 {
546  complex_t *pa = NULL, *pb = NULL, *pc = NULL;
547  int err;
548 
549  if(!a || !b || !c || !pfft)
550  return ERROR_PTR;
551  if(na<1 || nb < 1)
552  return ERROR_SIZE;
553  if(nfft<2)
554  return ERROR_FFT_SIZE;
555 
556  pa = (complex_t*) malloc(na*sizeof(complex_t));
557  pb = (complex_t*) malloc(nb*sizeof(complex_t));
558  pc = (complex_t*) malloc((na+nb-1)*sizeof(complex_t));
559 
560  re2cmplx(a, na, pa);
561  re2cmplx(b, nb, pb);
562 
563  err = conv_fft_cmplx(pa, na, pb, nb, pfft, nfft, pc);
564  if(err != RES_OK)
565  goto exit_label;
566 
567  err = cmplx2re(pc, na+nb-1, c, NULL);
568 
569 exit_label:
570  if(pa) free(pa);
571  if(pb) free(pb);
572  if(pc) free(pc);
573 
574  return err;
575 }
576 
577 
578 
579 #ifdef DOXYGEN_ENGLISH
580 
656 #endif
657 #ifdef DOXYGEN_RUSSIAN
658 
744 #endif
745 int DSPL_API conv_fft_cmplx(complex_t* a, int na, complex_t* b, int nb,
746  fft_t* pfft, int nfft, complex_t* c)
747 {
748 
749  int La, Lb, Lc, Nz, n, p0, p1, ind, err;
750  complex_t *pa, *pb;
751  complex_t *pt, *pA, *pB, *pC;
752 
753  if(!a || !b || !c)
754  return ERROR_PTR;
755  if(na < 1 || nb < 1)
756  return ERROR_SIZE;
757 
758  if(na >= nb)
759  {
760  La = na;
761  Lb = nb;
762  pa = a;
763  pb = b;
764  }
765  else
766  {
767  La = nb;
768  pa = b;
769  Lb = na;
770  pb = a;
771  }
772 
773  Lc = La + Lb - 1;
774  Nz = nfft - Lb;
775 
776  if(Nz <= 0)
777  return ERROR_FFT_SIZE;
778 
779  pt = (complex_t*)malloc(nfft*sizeof(complex_t));
780  pB = (complex_t*)malloc(nfft*sizeof(complex_t));
781  pA = (complex_t*)malloc(nfft*sizeof(complex_t));
782  pC = (complex_t*)malloc(nfft*sizeof(complex_t));
783 
784  memset(pt, 0, nfft*sizeof(complex_t));
785  memcpy(pt+Nz, pb, Lb*sizeof(complex_t));
786 
787  err = fft_cmplx(pt, nfft, pfft, pB);
788  if(err != RES_OK)
789  goto exit_label;
790 
791  p0 = -Lb;
792  p1 = p0 + nfft;
793  ind = 0;
794  while(ind < Lc)
795  {
796  if(p0 >=0)
797  {
798  if(p1 < La)
799  err = fft_cmplx(pa + p0, nfft, pfft, pA);
800  else
801  {
802  memset(pt, 0, nfft*sizeof(complex_t));
803  memcpy(pt, pa+p0, (nfft+La-p1)*sizeof(complex_t));
804  err = fft_cmplx(pt, nfft, pfft, pA);
805  }
806  }
807  else
808  {
809  memset(pt, 0, nfft*sizeof(complex_t));
810  if(p1 < La)
811  memcpy(pt - p0, pa, (nfft+p0)*sizeof(complex_t));
812  else
813  memcpy(pt - p0, pa, La * sizeof(complex_t));
814  err = fft_cmplx(pt, nfft, pfft, pA);
815  }
816 
817  if(err != RES_OK)
818  goto exit_label;
819 
820  for(n = 0; n < nfft; n++)
821  {
822  RE(pC[n]) = CMRE(pA[n], pB[n]);
823  IM(pC[n]) = CMIM(pA[n], pB[n]);
824  }
825 
826 
827  if(ind+nfft < Lc)
828  err = ifft_cmplx(pC, nfft, pfft, c+ind);
829  else
830  {
831  err = ifft_cmplx(pC, nfft, pfft, pt);
832  memcpy(c+ind, pt, (Lc-ind)*sizeof(complex_t));
833  }
834  if(err != RES_OK)
835  goto exit_label;
836 
837  p0 += Nz;
838  p1 += Nz;
839  ind += Nz;
840  }
841 
842 exit_label:
843  if(pt) free(pt);
844  if(pB) free(pB);
845  if(pA) free(pA);
846  if(pC) free(pC);
847 
848  return err;
849 }
850 
851 
852 #ifdef DOXYGEN_ENGLISH
853 
924 #endif
925 #ifdef DOXYGEN_RUSSIAN
926 
992 #endif
993 int DSPL_API filter_iir(double* b, double* a, int ord,
994  double* x, int n, double* y)
995 {
996  double *buf = NULL;
997  double *an = NULL;
998  double *bn = NULL;
999  double u;
1000  int k;
1001  int m;
1002  int count;
1003 
1004  if(!b || !x || !y)
1005  return ERROR_PTR;
1006 
1007  if(ord < 1 || n < 1)
1008  return ERROR_SIZE;
1009 
1010  if(a && a[0]==0.0)
1011  return ERROR_FILTER_A0;
1012 
1013  count = ord + 1;
1014  buf = (double*) malloc(count*sizeof(double));
1015  an = (double*) malloc(count*sizeof(double));
1016 
1017  memset(buf, 0, count*sizeof(double));
1018 
1019  if(!a)
1020  {
1021  memset(an, 0, count*sizeof(double));
1022  bn = b;
1023  }
1024  else
1025  {
1026  bn = (double*) malloc(count*sizeof(double));
1027  for(k = 0; k < count; k++)
1028  {
1029  an[k] = a[k] / a[0];
1030  bn[k] = b[k] / a[0];
1031  }
1032  }
1033 
1034  for(k = 0; k < n; k++)
1035  {
1036  for(m = ord; m > 0; m--)
1037  buf[m] = buf[m-1];
1038  u = 0.0;
1039  for(m = ord; m > 0; m--)
1040  u += buf[m]*an[m];
1041 
1042  buf[0] = x[k] - u;
1043  y[k] = 0.0;
1044  for(m = 0; m < count; m++)
1045  y[k] += buf[m] * bn[m];
1046  }
1047 
1048  if(buf)
1049  free(buf);
1050  if(an)
1051  free(an);
1052  if(bn && (bn != b))
1053  free(bn);
1054  return RES_OK;
1055 }
1056 
int DSPL_API conv_fft(double *a, int na, double *b, int nb, fft_t *pfft, int nfft, double *c)
Real vectors fast linear convolution by using fast Fourier transform algorithms.
Definition: conv.c:543
int DSPL_API filter_iir(double *b, double *a, int ord, double *x, int n, double *y)
Real IIR filtration.
Definition: conv.c:993
int DSPL_API conv(double *a, int na, double *b, int nb, double *c)
Real vectors linear convolution.
Definition: conv.c:157
#define ERROR_FILTER_A0
Parameter cannot be zero digital for IIR filter transfer characteristic .
Definition: dspl.h:511
#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:549
#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:557
int DSPL_API fft_cmplx(complex_t *x, int n, fft_t *pfft, complex_t *y)
Fast Fourier transform for the complex vector.
Definition: fft.c:645
int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t *pfft, complex_t *y)
Inverse fast Fourier transform.
Definition: fft.c:179
int DSPL_API conv_fft_cmplx(complex_t *a, int na, complex_t *b, int nb, fft_t *pfft, int nfft, complex_t *c)
Complex vectors fast linear convolution by using fast Fourier transform algorithms.
Definition: conv.c:745
Fast Fourier Transform Object Data Structure.
Definition: dspl.h:227
int DSPL_API re2cmplx(double *x, int n, complex_t *y)
Convert real array to the complex array.
Definition: complex.c:246
int DSPL_API cmplx2re(complex_t *x, int n, double *re, double *im)
Separate complex vector to the real and image vectors.
Definition: complex.c:130
int DSPL_API conv_cmplx(complex_t *a, int na, complex_t *b, int nb, complex_t *c)
Complex vectors linear convolution.
Definition: conv.c:327
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 ERROR_FFT_SIZE
The FFT size is not set correctly. FFT size can be composite , where , and – an arbitrary prime fact...
Definition: dspl.h:510
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:417