1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
|
#ifndef CONVOLUTIONS_H
#define CONVOLUTIONS_H
#include "../../structures/types.h"
#include "../../util/rng.h"
#include <cmath>
class Convolutions {
public:
/**
* This function will convolve data with the specified kernel. It will
* place the element at position (kernelSize/2) of the kernel in the
* centre, i.e. data[0] := sum over i in kS : data[i - kS + _kS/2_] *
* kernel[i]. Therefore, it makes most sense to specify an odd kernelSize if
* the kernel consists of a peak / symmetric function.
*
* This function assumes the function represented by @c data is zero outside
* its domain [0..dataSize}.
*
* This function does not reverse one of the input functions, therefore this
* function is not adhering to the strict definition of a convolution.
*
* This function does not use an FFT to calculate the convolution, and is
* therefore O(dataSize x kernelSize).
*
* @param data The data to be convolved (will also be the output)
* @param dataSize Number of samples in data
* @param kernel The kernel to be convolved, central sample = kernelSize/2
* @param kernelSize Number of samples in the kernel, probably desired to be
* odd if the kernel is symmetric.
*/
static void OneDimensionalConvolutionBorderZero(num_t *data,
unsigned dataSize,
const num_t *kernel,
unsigned kernelSize) {
num_t *tmp = new num_t[dataSize];
for (unsigned i = 0; i < dataSize; ++i) {
unsigned kStart, kEnd;
const int offset = i - kernelSize / 2;
// Make sure that kStart >= 0
if (offset < 0)
kStart = -offset;
else
kStart = 0;
// Make sure that kEnd + offset <= dataSize
if (offset + kernelSize > dataSize)
kEnd = dataSize - offset;
else
kEnd = kernelSize;
num_t sum = 0.0;
for (unsigned k = kStart; k < kEnd; ++k)
sum += data[k + offset] * kernel[k];
tmp[i] = sum;
}
for (unsigned i = 0; i < dataSize; ++i) data[i] = tmp[i];
delete[] tmp;
}
/**
* This function will convolve data with the specified kernel. It will
* place the element at position (kernelSize/2) of the kernel in the
* centre, i.e. data[0] := sum over i in kS : data[i - kS + _kS/2_] *
* kernel[i]. Therefore, it makes most sense to specify an odd kernelSize if
* the kernel consists of a peak / symmetric function.
*
* This function assumes the function represented by @c data is not zero
* outside its domain [0..dataSize}, but is "unknown". It assumes that the
* integral over @c kernel is one. If this is not the case, you have to
* multiply all output values with the kernels integral after convolution.
*
* This function does not reverse one of the input functions, therefore this
* function is not adhering to the strict definition of a convolution.
*
* This function does not use an FFT to calculate the convolution, and is
* therefore O(dataSize x kernelSize).
*
* @param data The data to be convolved (will also be the output)
* @param dataSize Number of samples in data
* @param kernel The kernel to be convolved, central sample = kernelSize/2
* @param kernelSize Number of samples in the kernel, probably desired to be
* odd if the kernel is symmetric.
*/
static void OneDimensionalConvolutionBorderInterp(num_t *data,
unsigned dataSize,
const num_t *kernel,
unsigned kernelSize) {
num_t *tmp = new num_t[dataSize];
for (unsigned i = 0; i < dataSize; ++i) {
unsigned kStart, kEnd;
const int offset = i - kernelSize / 2;
// Make sure that kStart >= 0
if (offset < 0)
kStart = -offset;
else
kStart = 0;
// Make sure that kEnd + offset <= dataSize
if (offset + kernelSize > dataSize)
kEnd = dataSize - offset;
else
kEnd = kernelSize;
num_t sum = 0.0;
num_t weight = 0.0;
for (unsigned k = kStart; k < kEnd; ++k) {
sum += data[k + offset] * kernel[k];
weight += kernel[k];
}
// Weighting is performed to correct for the "missing data" outside the
// domain of data. The actual correct value should be sum * (sumover
// kernel / sumover kernel-used ). We however assume "sumover kernel" to
// be 1. sumover kernel-used is the weight.
tmp[i] = sum / weight;
}
for (unsigned i = 0; i < dataSize; ++i) data[i] = tmp[i];
delete[] tmp;
}
static void OneDimensionalGausConvolution(num_t *data, unsigned dataSize,
num_t sigma) {
unsigned kernelSize = (unsigned)roundn(sigma * 3.0L);
if (kernelSize % 2 == 0) ++kernelSize;
if (kernelSize > dataSize * 2) {
if (dataSize == 0) return;
kernelSize = dataSize * 2 - 1;
}
unsigned centreElement = kernelSize / 2;
num_t *kernel = new num_t[kernelSize];
for (unsigned i = 0; i < kernelSize; ++i) {
num_t x = ((num_t)i - (num_t)centreElement);
kernel[i] = RNG::EvaluateGaussian(x, sigma);
}
OneDimensionalConvolutionBorderInterp(data, dataSize, kernel, kernelSize);
delete[] kernel;
}
/**
* Perform a sinc convolution. This is equivalent with a low-pass filter,
* where @c frequency specifies the frequency in index units.
*
* With frequency=1, the data will be convolved with the
* delta function and have no effect. With frequency = 0.25,
* all fringes faster than 0.25 fringes/index units will be
* filtered.
*
* The function is O(dataSize^2).
*/
static void OneDimensionalSincConvolution(num_t *data, unsigned dataSize,
num_t frequency) {
if (dataSize == 0) return;
const unsigned kernelSize = dataSize * 2 - 1;
const unsigned centreElement = kernelSize / 2;
num_t *kernel = new num_t[kernelSize];
const numl_t factor = 2.0 * frequency * M_PInl;
numl_t sum = 0.0;
for (unsigned i = 0; i < kernelSize; ++i) {
const numl_t x = (((numl_t)i - (numl_t)centreElement) * factor);
if (x != 0.0)
kernel[i] = (num_t)(sinnl(x) / x);
else
kernel[i] = 1.0;
sum += kernel[i];
}
for (unsigned i = 0; i < kernelSize; ++i) {
kernel[i] /= sum;
}
OneDimensionalConvolutionBorderZero(data, dataSize, kernel, kernelSize);
delete[] kernel;
}
/**
* Perform a sinc convolution, with a convolution kernel that has been Hamming
* windowed. Because of the hamming window, the filter has less of a steep
* cutting edge, but less ripples.
*
* See OneDimensionalSincConvolution() for more info.
*
* The function is O(dataSize^2).
*/
static void OneDimensionalSincConvolutionHammingWindow(num_t *data,
unsigned dataSize,
num_t frequency) {
if (dataSize == 0) return;
const unsigned kernelSize = dataSize * 2 - 1;
const unsigned centreElement = kernelSize / 2;
num_t *kernel = new num_t[kernelSize];
const numl_t sincFactor = 2.0 * frequency * M_PInl;
const numl_t hammingFactor = 2.0 * M_PInl * (numl_t)(kernelSize - 1);
numl_t sum = 0.0;
for (unsigned i = 0; i < kernelSize; ++i) {
const numl_t hamming = 0.54 - 0.46 * cosnl(hammingFactor * (numl_t)i);
const numl_t x = (((numl_t)i - (numl_t)centreElement) * sincFactor);
if (x != 0.0)
kernel[i] = (num_t)(sinnl(x) / x) * hamming;
else
kernel[i] = 1.0;
sum += kernel[i];
}
for (unsigned i = 0; i < kernelSize; ++i) {
kernel[i] /= sum;
}
OneDimensionalConvolutionBorderZero(data, dataSize, kernel, kernelSize);
delete[] kernel;
}
};
#endif
|
