Instrument Neutral Distributed Interface INDI  2.0.2
buffer.c
Go to the documentation of this file.
1 /*
2 * DSP API - a digital signal processing library for astronomy usage
3 * Copyright © 2017-2022 Ilia Platone
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 3 of the License, or (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with this program; if not, write to the Free Software Foundation,
17 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19 
20 #include "dsp.h"
21 #include <setjmp.h>
22 #include <signal.h>
23 
25 {
26  if(stream->dims == 0)
27  return;
28  dsp_t* tmp = (dsp_t*)malloc(sizeof(dsp_t) * stream->len);
29  int x, d;
30  for(x = 0; x < stream->len/2; x++) {
31  int* pos = dsp_stream_get_position(stream, x);
32  for(d = 0; d < stream->dims; d++) {
33  if(pos[d]<stream->sizes[d] / 2) {
34  pos[d] += stream->sizes[d] / 2;
35  } else {
36  pos[d] -= stream->sizes[d] / 2;
37  }
38  }
39  tmp[x] = stream->buf[dsp_stream_set_position(stream, pos)];
40  tmp[dsp_stream_set_position(stream, pos)] = stream->buf[x];
41  free(pos);
42  }
43  memcpy(stream->buf, tmp, stream->len * sizeof(dsp_t));
44  free(tmp);
45 }
46 
48 {
49  int k;
50 
51  dsp_t mean = dsp_stats_mean(stream->buf, stream->len);
52  for(k = 0; k < stream->len; k++)
53  stream->buf[k] = stream->buf[k] - mean;
54 
55 }
56 
57 void dsp_buffer_sub(dsp_stream_p stream, dsp_t* in, int inlen)
58 {
59  int len = Min(stream->len, inlen);
60 
61  int k;
62  for(k = 0; k < len; k++) {
63  stream->buf[k] = stream->buf[k] - in[k];
64  }
65 
66 }
67 
68 void dsp_buffer_sum(dsp_stream_p stream, dsp_t* in, int inlen)
69 {
70  int len = Min(stream->len, inlen);
71 
72  int k;
73  for(k = 0; k < len; k++) {
74  stream->buf[k] += in[k];
75  }
76 
77 }
78 
79 void dsp_buffer_max(dsp_stream_p stream, dsp_t* in, int inlen)
80 {
81  int len = Min(stream->len, inlen);
82 
83  int k;
84  for(k = 0; k < len; k++) {
85  stream->buf[k] = Max(stream->buf[k], in[k]);
86  }
87 
88 }
89 
90 void dsp_buffer_min(dsp_stream_p stream, dsp_t* in, int inlen)
91 {
92  int len = Min(stream->len, inlen);
93 
94  int k;
95  for(k = 0; k < len; k++) {
96  stream->buf[k] = Min(stream->buf[k], in[k]);
97  }
98 
99 }
100 
101 void dsp_buffer_div(dsp_stream_p stream, dsp_t* in, int inlen)
102 {
103  int len = Min(stream->len, inlen);
104 
105  int k;
106  for(k = 0; k < len; k++) {
107  stream->buf[k] = stream->buf[k] / in[k];
108  }
109 
110 }
111 
112 void dsp_buffer_mul(dsp_stream_p stream, dsp_t* in, int inlen)
113 {
114  int len = Min(stream->len, inlen);
115 
116  int k;
117  for(k = 0; k < len; k++) {
118  stream->buf[k] = stream->buf[k] * in[k];
119  }
120 
121 }
122 
123 void dsp_buffer_pow(dsp_stream_p stream, dsp_t* in, int inlen)
124 {
125  int len = Min(stream->len, inlen);
126 
127  int k;
128  for(k = 0; k < len; k++) {
129  stream->buf[k] = pow(stream->buf[k], in[k]);
130  }
131 
132 }
133 
134 void dsp_buffer_log(dsp_stream_p stream, dsp_t* in, int inlen)
135 {
136  int len = Min(stream->len, inlen);
137 
138  int k;
139  for(k = 0; k < len; k++) {
140  stream->buf[k] = Log(stream->buf[k], in[k]);
141  }
142 
143 }
144 
146 {
147  int k;
148 
149  for(k = 0; k < stream->len; k++) {
150  stream->buf[k] = val - stream->buf[k];
151  }
152 
153 }
154 
156 {
157  int k;
158 
159  for(k = 0; k < stream->len; k++) {
160  stream->buf[k] = stream->buf[k] - val;
161  }
162 
163 }
164 
166 {
167  int k;
168 
169  for(k = 0; k < stream->len; k++) {
170  stream->buf[k] += val;
171  }
172 
173 }
174 
175 void dsp_buffer_1div(dsp_stream_p stream, double val)
176 {
177  int k;
178 
179  for(k = 0; k < stream->len; k++) {
180  stream->buf[k] = val / stream->buf[k];
181  }
182 
183 }
184 
185 void dsp_buffer_div1(dsp_stream_p stream, double val)
186 {
187  int k;
188 
189  for(k = 0; k < stream->len; k++) {
190  stream->buf[k] /= val;
191  }
192 
193 }
194 
195 void dsp_buffer_mul1(dsp_stream_p stream, double val)
196 {
197  int k;
198 
199  for(k = 0; k < stream->len; k++) {
200  stream->buf[k] = stream->buf[k] * val;
201  }
202 
203 }
204 
205 void dsp_buffer_pow1(dsp_stream_p stream, double val)
206 {
207  int k;
208 
209  for(k = 0; k < stream->len; k++) {
210  stream->buf[k] = pow(stream->buf[k], val);
211  }
212 
213 }
214 
215 void dsp_buffer_log1(dsp_stream_p stream, double val)
216 {
217  int k;
218 
219  for(k = 0; k < stream->len; k++) {
220  stream->buf[k] = Log(stream->buf[k], val);
221  }
222 
223 }
224 
225 static int compare( const void* a, const void* b)
226 {
227  dsp_t int_a = * ( (dsp_t*) a );
228  dsp_t int_b = * ( (dsp_t*) b );
229 
230  if ( int_a == int_b ) return 0;
231  else if ( int_a < int_b ) return -1;
232  else return 1;
233 }
234 
235 static void* dsp_buffer_median_th(void* arg)
236 {
237  struct {
238  int cur_th;
239  int size;
240  int median;
241  dsp_stream_p stream;
242  dsp_stream_p box;
243  } *arguments = arg;
244  dsp_stream_p stream = arguments->stream;
245  dsp_stream_p box = arguments->box;
246  dsp_stream_p in = stream->parent;
247  int cur_th = arguments->cur_th;
248  int size = arguments->size;
249  int median = arguments->median;
250  int start = cur_th * stream->len / dsp_max_threads(0);
251  int end = start + stream->len / dsp_max_threads(0);
252  end = Min(stream->len, end);
253  int x, y, dim, idx;
254  dsp_t* sorted = (dsp_t*)malloc(pow(size, stream->dims) * sizeof(dsp_t));
255  int len = pow(size, in->dims);
256  for(x = start; x < end; x++) {
257  dsp_t* buf = sorted;
258  for(y = 0; y < box->len; y++) {
259  int *pos = dsp_stream_get_position(stream, x);
260  int *mat = dsp_stream_get_position(box, y);
261  for(dim = 0; dim < stream->dims; dim++) {
262  pos[dim] += mat[dim] - size / 2;
263  }
264  idx = dsp_stream_set_position(stream, pos);
265  if(idx >= 0 && idx < in->len) {
266  *buf++ = in->buf[idx];
267  }
268  free(pos);
269  free(mat);
270  }
271  qsort(sorted, len, sizeof(dsp_t), compare);
272  stream->buf[x] = sorted[median*box->len/size];
273  }
275  dsp_stream_free(box);
276  free(sorted);
277  return NULL;
278 }
279 
280 void dsp_buffer_median(dsp_stream_p in, int size, int median)
281 {
282  size_t y;
283  int d;
284  dsp_stream_p stream = dsp_stream_copy(in);
285  dsp_buffer_set(stream->buf, stream->len, 0);
286  stream->parent = in;
287  pthread_t *th = (pthread_t *)malloc(sizeof(pthread_t)*dsp_max_threads(0));
288  struct {
289  int cur_th;
290  int size;
291  int median;
292  dsp_stream_p stream;
293  dsp_stream_p box;
294  } thread_arguments[dsp_max_threads(0)];
295  for(y = 0; y < dsp_max_threads(0); y++)
296  {
297  thread_arguments[y].cur_th = y;
298  thread_arguments[y].size = size;
299  thread_arguments[y].median = median;
300  thread_arguments[y].stream = stream;
301  thread_arguments[y].box = dsp_stream_new();
302  for(d = 0; d < stream->dims; d++)
303  dsp_stream_add_dim(thread_arguments[y].box, size);
304  dsp_stream_alloc_buffer(thread_arguments[y].box, thread_arguments[y].box->len);
305  pthread_create(&th[y], NULL, dsp_buffer_median_th, &thread_arguments[y]);
306  }
307  for(y = 0; y < dsp_max_threads(0); y++)
308  pthread_join(th[y], NULL);
309  free(th);
310  stream->parent = NULL;
311  dsp_buffer_copy(stream->buf, in->buf, stream->len);
312  dsp_stream_free_buffer(stream);
313  dsp_stream_free(stream);
314 }
315 
316 static void* dsp_buffer_sigma_th(void* arg)
317 {
318  struct {
319  int cur_th;
320  int size;
321  dsp_stream_p stream;
322  dsp_stream_p box;
323  } *arguments = arg;
324  dsp_stream_p stream = arguments->stream;
325  dsp_stream_p in = stream->parent;
326  dsp_stream_p box = arguments->box;
327  int cur_th = arguments->cur_th;
328  int size = arguments->size;
329  int start = cur_th * stream->len / dsp_max_threads(0);
330  int end = start + stream->len / dsp_max_threads(0);
331  end = Min(stream->len, end);
332  int x, y, dim, idx;
333  dsp_t* sigma = (dsp_t*)malloc(pow(size, stream->dims) * sizeof(dsp_t));
334  int len = pow(size, in->dims);
335  for(x = start; x < end; x++) {
336  dsp_t* buf = sigma;
337  for(y = 0; y < box->len; y++) {
338  int *pos = dsp_stream_get_position(stream, x);
339  int *mat = dsp_stream_get_position(box, y);
340  for(dim = 0; dim < stream->dims; dim++) {
341  pos[dim] += mat[dim] - size / 2;
342  }
343  idx = dsp_stream_set_position(stream, pos);
344  if(idx >= 0 && idx < in->len) {
345  buf[y] = in->buf[idx];
346  }
347  free(pos);
348  free(mat);
349  }
350  stream->buf[x] = dsp_stats_stddev(buf, len);
351  }
353  dsp_stream_free(box);
354  free(sigma);
355  return NULL;
356 }
357 
358 void dsp_buffer_sigma(dsp_stream_p in, int size)
359 {
360  size_t y;
361  int d;
362  dsp_stream_p stream = dsp_stream_copy(in);
363  dsp_buffer_set(stream->buf, stream->len, 0);
364  stream->parent = in;
365  pthread_t *th = (pthread_t *)malloc(sizeof(pthread_t)*dsp_max_threads(0));
366  struct {
367  int cur_th;
368  int size;
369  dsp_stream_p stream;
370  dsp_stream_p box;
371  } thread_arguments[dsp_max_threads(0)];
372  for(y = 0; y < dsp_max_threads(0); y++)
373  {
374  thread_arguments[y].cur_th = y;
375  thread_arguments[y].size = size;
376  thread_arguments[y].stream = stream;
377  thread_arguments[y].box = dsp_stream_new();
378  for(d = 0; d < stream->dims; d++)
379  dsp_stream_add_dim(thread_arguments[y].box, size);
380  pthread_create(&th[y], NULL, dsp_buffer_sigma_th, &thread_arguments[y]);
381  }
382  for(y = 0; y < dsp_max_threads(0); y++)
383  pthread_join(th[y], NULL);
384  free(th);
385  stream->parent = NULL;
386  dsp_buffer_copy(stream->buf, in->buf, stream->len);
387  dsp_stream_free_buffer(stream);
388  dsp_stream_free(stream);
389 }
390 
391 void dsp_buffer_deviate(dsp_stream_p stream, dsp_t* deviation, dsp_t mindeviation, dsp_t maxdeviation)
392 {
393  dsp_stream_p tmp = dsp_stream_copy(stream);
394  int k;
395  for(k = 1; k < stream->len; k++) {
396  stream->buf[(int)Max(0, Min(stream->len, ((deviation[k] - mindeviation) * (maxdeviation - mindeviation) + mindeviation) + k))] = tmp->buf[k];
397  }
398  dsp_stream_free(tmp);
399 }
std::thread th
#define Max(a, b)
if max() is not present you can use this one
Definition: dsp.h:179
DLL_EXPORT unsigned long int dsp_max_threads(unsigned long value)
get/set the maximum number of threads allowed
Definition: stream.c:31
#define Log(a, b)
Logarithm of a with arbitrary base b.
Definition: dsp.h:186
double dsp_t
Definition: dsp.h:69
#define Min(a, b)
if min() is not present you can use this one
Definition: dsp.h:172
struct dsp_stream_t * parent
The parent stream.
Definition: dsp.h:381
int * sizes
Sizes of each dimension.
Definition: dsp.h:373
int dims
Number of dimensions of the buffers.
Definition: dsp.h:371
dsp_t * buf
buffer
Definition: dsp.h:375
int len
The buffers length.
Definition: dsp.h:369
#define dsp_buffer_copy(in, out, len)
Fill the output buffer with the values of the elements of the input stream by casting them to the out...
Definition: dsp.h:1038
void dsp_buffer_sum(dsp_stream_p stream, dsp_t *in, int inlen)
Sum elements of one stream to another's.
Definition: buffer.c:68
void dsp_buffer_log(dsp_stream_p stream, dsp_t *in, int inlen)
Logarithm elements of one stream using another's as base.
Definition: buffer.c:134
void dsp_buffer_median(dsp_stream_p in, int size, int median)
Median elements of the input stream.
Definition: buffer.c:280
void dsp_buffer_shift(dsp_stream_p stream)
Shift a stream on each dimension.
Definition: buffer.c:24
void dsp_buffer_deviate(dsp_stream_p stream, dsp_t *deviation, dsp_t mindeviation, dsp_t maxdeviation)
Deviate forward the first input stream using the second stream as indexing reference.
Definition: buffer.c:391
void dsp_buffer_div1(dsp_stream_p stream, double val)
Divide elements of the input stream to a value.
Definition: buffer.c:185
void dsp_buffer_removemean(dsp_stream_p stream)
Subtract mean from stream.
Definition: buffer.c:47
void dsp_buffer_1sub(dsp_stream_p stream, dsp_t val)
Subtract each element of the input stream a value.
Definition: buffer.c:145
void dsp_buffer_max(dsp_stream_p stream, dsp_t *in, int inlen)
Subtract elements of one stream from another's.
Definition: buffer.c:79
void dsp_buffer_sigma(dsp_stream_p in, int size)
Standard deviation of each element of the input stream within the given size.
Definition: buffer.c:358
#define dsp_buffer_set(buf, len, _val)
Fill the buffer with the passed value.
Definition: dsp.h:815
void dsp_buffer_div(dsp_stream_p stream, dsp_t *in, int inlen)
Divide elements of one stream to another's.
Definition: buffer.c:101
void dsp_buffer_sum1(dsp_stream_p stream, dsp_t val)
Sum elements of the input stream to a value.
Definition: buffer.c:165
void dsp_buffer_mul(dsp_stream_p stream, dsp_t *in, int inlen)
Multiply elements of one stream to another's.
Definition: buffer.c:112
void dsp_buffer_pow1(dsp_stream_p stream, double val)
Expose elements of the input stream to the given power.
Definition: buffer.c:205
void dsp_buffer_log1(dsp_stream_p stream, double val)
Logarithm elements of the input stream using the given base.
Definition: buffer.c:215
void dsp_buffer_pow(dsp_stream_p stream, dsp_t *in, int inlen)
Expose elements of one stream to another's.
Definition: buffer.c:123
void dsp_buffer_min(dsp_stream_p stream, dsp_t *in, int inlen)
Sum elements of one stream to another's.
Definition: buffer.c:90
void dsp_buffer_1div(dsp_stream_p stream, double val)
Divide a value to each element of the input stream.
Definition: buffer.c:175
void dsp_buffer_sub(dsp_stream_p stream, dsp_t *in, int inlen)
Subtract elements of one stream from another's.
Definition: buffer.c:57
void dsp_buffer_sub1(dsp_stream_p stream, dsp_t val)
Subtract a value from elements of the input stream.
Definition: buffer.c:155
void dsp_buffer_mul1(dsp_stream_p stream, double val)
Multiply elements of the input stream to a value.
Definition: buffer.c:195
DLL_EXPORT void dsp_stream_free(dsp_stream_p stream)
Free the DSP stream passed as argument.
Definition: stream.c:163
DLL_EXPORT int dsp_stream_set_position(dsp_stream_p stream, int *pos)
Obtain the position the DSP stream by parsing multidimensional indexes.
Definition: stream.c:440
DLL_EXPORT int * dsp_stream_get_position(dsp_stream_p stream, int index)
Return the multidimensional positional indexes of a DSP stream by specify a linear index.
Definition: stream.c:420
DLL_EXPORT dsp_stream_p dsp_stream_copy(dsp_stream_p stream)
Create a copy of the DSP stream passed as argument.
Definition: stream.c:192
DLL_EXPORT dsp_stream_p dsp_stream_new(void)
Allocate a new DSP stream type.
Definition: stream.c:124
DLL_EXPORT void dsp_stream_add_dim(dsp_stream_p stream, int len)
Add a dimension with length len to a DSP stream.
Definition: stream.c:227
DLL_EXPORT void dsp_stream_alloc_buffer(dsp_stream_p stream, int len)
Allocate a buffer with length len on the stream passed as argument.
Definition: stream.c:78
DLL_EXPORT void dsp_stream_free_buffer(dsp_stream_p stream)
Free the buffer of the DSP Stream passed as argument.
Definition: stream.c:112
#define dsp_stats_stddev(buf, len)
Standard deviation of the inut stream.
Definition: dsp.h:667
#define dsp_stats_mean(buf, len)
A mean calculator.
Definition: dsp.h:649
Contains a set of informations and data relative to a buffer and how to use it.
Definition: dsp.h:363