Instrument Neutral Distributed Interface INDI  1.9.5
filters.c
Go to the documentation of this file.
1 /*
2  * libDSPAU - a digital signal processing library for astronomy usage
3  * Copyright (C) 2017 Ilia Platone <info@iliaplatone.com>
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (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
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  */
18 
19 #include "dsp.h"
20 
22 {
23  dsp_t* in = stream->buf;
24  dsp_t *out = (dsp_t*)malloc(sizeof(dsp_t) * stream->len);
25  int len = stream->len;
26  dsp_t mean = dsp_stats_mean(stream->buf, stream->len);
27  int val = 0;
28  int i;
29  for(i = 0; i < len; i++) {
30  val = in [i] - mean;
31  out [i] = (abs (val) + mean);
32  }
33  memcpy(stream->buf, out, stream->len * sizeof(dsp_t));
34  free(out);
35 }
36 
37 void dsp_filter_calc_coefficients(double SamplingFrequency, double LowFrequency, double HighFrequency, double* CF, double* R, double *K)
38 {
39  double BW = (HighFrequency - LowFrequency) / SamplingFrequency;
40  *CF = 2.0 * cos((LowFrequency + HighFrequency) * M_PI / SamplingFrequency);
41  *R = 1.0 - 3.0 * BW;
42  *K = (1.0 - *R * *CF + *R * *R) / (2.0 - *CF);
43 }
44 
45 void dsp_filter_lowpass(dsp_stream_p stream, double SamplingFrequency, double Frequency)
46 {
47  double *out = (double*)malloc(sizeof(double) * stream->len);
48  double CF = cos(Frequency / 2.0 * M_PI / SamplingFrequency);
49  int dim = -1, i;
50  out[0] = stream->buf[0];
51  while (dim++ < stream->dims - 1) {
52  int size = (dim < 0 ? 1 : stream->sizes[dim]);
53  for(i = size; i < stream->len; i+=size) {
54  out[i] += stream->buf[i] + (out[i - size] - stream->buf[i]) * CF;
55  }
56  }
57  memcpy(stream->buf, out, stream->len * sizeof(double));
58  free(out);
59 }
60 
61 void dsp_filter_highpass(dsp_stream_p stream, double SamplingFrequency, double Frequency)
62 {
63  dsp_t *out = (dsp_t*)malloc(sizeof(dsp_t) * stream->len);
64  double CF = cos(Frequency / 2.0 * M_PI / SamplingFrequency);
65  int dim = -1, i;
66  out[0] = stream->buf[0];
67  while (dim++ < stream->dims - 1) {
68  int size = (dim < 0 ? 1 : stream->sizes[dim]);
69  for(i = size; i < stream->len; i+=size) {
70  out[i] += stream->buf[i] + (out[i - size] - stream->buf[i]) * CF;
71  }
72  }
73  dsp_buffer_sub(stream, out, stream->len);
74  free(out);
75 }
76 
77 void dsp_filter_bandreject(dsp_stream_p stream, double SamplingFrequency, double LowFrequency, double HighFrequency) {
78  dsp_t *out = (dsp_t*)malloc(sizeof(dsp_t) * stream->len);
79  double R, K, CF;
80  dsp_filter_calc_coefficients(SamplingFrequency, LowFrequency, HighFrequency, &CF, &R, &K);
81  double R2 = R*R;
82 
83  double a[3] = { K, -K*CF, K };
84  double b[2] = { R*CF, -R2 };
85 
86  int dim = -1, i, x;
87  while (dim++ < stream->dims - 1) {
88  for(i = 0; i < stream->len; i++) {
89  int size = (dim < 0 ? 1 : stream->sizes[dim]);
90  out[i] = 0;
91  if(i < stream->len - 2 * size) {
92  for(x = 0; x < 3; x++) {
93  out[i] += (double)stream->buf[i + x * size] * a[2 - x];
94  }
95  }
96  if(i > size) {
97  for(x = 0; x < 2; x++) {
98  out[i] -= out[i - 2 * size + x * size] * b[x];
99  }
100  }
101  }
102  }
103  memcpy(stream->buf, out, stream->len * sizeof(dsp_t));
104  free(out);
105 }
106 
107 void dsp_filter_bandpass(dsp_stream_p stream, double SamplingFrequency, double LowFrequency, double HighFrequency) {
108  dsp_t *out = (dsp_t*)malloc(sizeof(dsp_t) * stream->len);
109  double R, K, CF;
110  dsp_filter_calc_coefficients(SamplingFrequency, LowFrequency, HighFrequency, &CF, &R, &K);
111  double R2 = R*R;
112 
113  double a[3] = { 1 - K, (K-R)*CF, R2 - K };
114  double b[2] = { R*CF, -R2 };
115 
116  int dim = -1, i, x;
117  while (dim++ < stream->dims - 1) {
118  for(i = 0; i < stream->len; i++) {
119  int size = (dim < 0 ? 1 : stream->sizes[dim]);
120  out[i] = 0;
121  if(i < stream->len - 2 * size) {
122  for(x = 0; x < 3; x++) {
123  out[i] += (double)stream->buf[i + x * size] * a[2 - x];
124  }
125  }
126  if(i > size) {
127  for(x = 0; x < 2; x++) {
128  out[i] -= out[i - 2 * size + x * size] * b[x];
129  }
130  }
131  }
132  }
133  memcpy(stream->buf, out, stream->len * sizeof(dsp_t));
134  free(out);
135 }
dsp.h
dsp_filter_highpass
void dsp_filter_highpass(dsp_stream_p stream, double SamplingFrequency, double Frequency)
A high pass filter.
Definition: filters.c:61
dsp_filter_squarelaw
void dsp_filter_squarelaw(dsp_stream_p stream)
A square law filter.
Definition: filters.c:21
dsp_stream_t
Contains a set of informations and data relative to a buffer and how to use it.
Definition: dsp.h:175
dsp_filter_bandreject
void dsp_filter_bandreject(dsp_stream_p stream, double SamplingFrequency, double LowFrequency, double HighFrequency)
A band reject filter.
Definition: filters.c:77
dsp_filter_calc_coefficients
void dsp_filter_calc_coefficients(double SamplingFrequency, double LowFrequency, double HighFrequency, double *CF, double *R, double *K)
Definition: filters.c:37
dsp_stream_t::len
int len
The buffers length.
Definition: dsp.h:178
dsp_stats_mean
#define dsp_stats_mean(buf, len)
A mean calculator.
Definition: dsp.h:459
dsp_stream_t::sizes
int * sizes
Sizes of each dimension.
Definition: dsp.h:182
dsp_filter_lowpass
void dsp_filter_lowpass(dsp_stream_p stream, double SamplingFrequency, double Frequency)
A low pass filter.
Definition: filters.c:45
dsp_filter_bandpass
void dsp_filter_bandpass(dsp_stream_p stream, double SamplingFrequency, double LowFrequency, double HighFrequency)
A band pass filter.
Definition: filters.c:107
dsp_stream_t::dims
int dims
Number of dimensions of the buffers.
Definition: dsp.h:180
dsp_t
#define dsp_t
Definition: dsp.h:57
dsp_buffer_sub
void dsp_buffer_sub(dsp_stream_p stream, dsp_t *in, int inlen)
Subtract elements of one stream from another's.
Definition: buffer.c:54
dsp_stream_t::buf
dsp_t * buf
buffer
Definition: dsp.h:184