libdspl-2.0
Библиотека алгоритмов цифровой обработки сигналов
resampling.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 
23 #include <stdlib.h>
24 #include <string.h>
25 #include <math.h>
26 #include "dspl.h"
27 #include "dspl_internal.h"
28 
29 
30 
31 
32 /******************************************************************************
33 Farrow resampler based on the cubic Lagrange polynomials
34 *******************************************************************************/
35 int DSPL_API farrow_lagrange(double *s, int n, double p, double q,
36  double frd, double **y, int *ny)
37 {
38  double a[4];
39  double t, x, dt;
40  int ind, k, res;
41  double g[4];
42  double *z;
43 
44  if(!s || !y)
45  return ERROR_PTR;
46 
47  if(n<1)
48  return ERROR_SIZE;
49 
50  if(p <= 0.0 || q <= 0.0)
51  return ERROR_RESAMPLE_RATIO;
52 
53  if(frd <= -1.0 || frd >= 1.0)
54  return ERROR_RESAMPLE_FRAC_DELAY;
55 
56  dt = q/p;
57 
58  if((*ny) != (int)((double)(n-1)/dt)+1 || !(*y))
59  {
60 
61  *ny = (int)((double)(n-1)/dt)+1;
62  (*y) = (double*)realloc((*y), (*ny)*sizeof(double));
63  }
64 
65  t = -frd;
66  k = 0;
67  while(k < (*ny))
68  {
69  ind = (int)floor(t)+1;
70  x = t - (double)ind;
71  ind-=2;
72  if(ind < 0)
73  {
74  memset(g, 0, 4*sizeof(double));
75  if(ind > (-3))
76  memcpy(g-ind, s, (4+ind)*sizeof(double));
77  z = g;
78  }
79  else
80  {
81  if(ind < n-3)
82  z = s+ind;
83  else
84  {
85  memset(g, 0, 4*sizeof(double));
86  if((n-ind)>0)
87  memcpy(g, s+ind, (n-ind)*sizeof(double));
88  z = g;
89  }
90  }
91  a[0] = z[2];
92  a[3] = DSPL_FARROW_LAGRANGE_COEFF*(z[3] -z[0]) + 0.5*(z[1] - z[2]);
93  a[1] = 0.5*(z[3] - z[1])-a[3];
94  a[2] = z[3] - z[2] -a[3]-a[1];
95 
96  res = polyval(a, 3, &x, 1, (*y)+k);
97 
98  if(res != RES_OK)
99  goto exit_label;
100  t+=dt;
101  k++;
102  }
103 
104 exit_label:
105  return res;
106 }
107 
108 
109 
110 
111 
112 /******************************************************************************
113 Farrow resampler based on the cubic splines
114 *******************************************************************************/
115 int DSPL_API farrow_spline(double *s, int n, double p, double q,
116  double frd, double **y, int *ny)
117 {
118  double a[4];
119  double t, x, dt;
120  int ind, k, res;
121  double g[4];
122  double *z;
123 
124  if(!s || !y)
125  return ERROR_PTR;
126 
127  if(n<1)
128  return ERROR_SIZE;
129 
130  if(p <= 0.0 || q <= 0.0)
131  return ERROR_RESAMPLE_RATIO;
132 
133  if(frd <= -1.0 || frd >= 1.0)
134  return ERROR_RESAMPLE_FRAC_DELAY;
135 
136  dt = q/p;
137 
138  if((*ny) != (int)((double)(n-1)/dt)+1 || !(*y))
139  {
140 
141  *ny = (int)((double)(n-1)/dt)+1;
142  (*y) = (double*)realloc((*y), (*ny)*sizeof(double));
143  }
144 
145  t = -frd;
146  k = 0;
147  while(k < (*ny))
148  {
149  ind = (int)floor(t)+1;
150  x = t - (double)ind;
151  ind-=2;
152  if(ind < 0)
153  {
154  memset(g, 0, 4*sizeof(double));
155  if(ind > (-3))
156  memcpy(g-ind, s, (4+ind)*sizeof(double));
157  z = g;
158  }
159  else
160  {
161  if(ind < n-3)
162  z = s+ind;
163  else
164  {
165  memset(g, 0, 4*sizeof(double));
166  if((n-ind)>0)
167  memcpy(g, s+ind, (n-ind)*sizeof(double));
168  z = g;
169  }
170  }
171  a[0] = z[2];
172  a[1] = 0.5*(z[3] - z[1]);
173  a[3] = 2.0*(z[1] - z[2]) + a[1] + 0.5*(z[2] - z[0]);
174  a[2] = z[1] - z[2] +a[3] + a[1];
175 
176  res = polyval(a, 3, &x, 1, (*y)+k);
177 
178  if(res != RES_OK)
179  goto exit_label;
180  t+=dt;
181  k++;
182  }
183 
184 exit_label:
185  return res;
186 }
187 
#define ERROR_SIZE
Ошибка при передаче размера массива. Данная ошибка возникает когда помимо указателя на массив входных...
Definition: dspl.h:148
#define ERROR_PTR
Ошибка указателя. Данная ошибка означает, что один из обязательных указателей (память под который дол...
Definition: dspl.h:140
int polyval(double *a, int ord, double *x, int n, double *y)
Расчет вещественного полинома
Definition: polyval.c:69
int farrow_lagrange(double *s, int n, double p, double q, double frd, double **y, int *ny)
Передискретизация вещественного сигнала на основе полиномиальной Лагранжевой интерполяции.
Definition: resampling.c:35
#define RES_OK
Функция завершилась корректно. Ошибки отсутствуют.
Definition: dspl.h:94
int farrow_spline(double *s, int n, double p, double q, double frd, double **y, int *ny)
Передискретизация вещественного сигнала на основе сплайн интерполяции.
Definition: resampling.c:115