IT++ Logo
sigfun.cpp
Go to the documentation of this file.
1
29#include <itpp/signal/sigfun.h>
31#include <itpp/signal/window.h>
34#include <itpp/base/matfunc.h>
35#include <itpp/base/specmat.h>
36#include <itpp/base/itcompat.h>
37#include <itpp/stat/misc_stat.h>
38
39
40namespace itpp
41{
42
43vec xcorr_old(const vec &x, const int max_lag, const std::string scaleopt)
44{
45 vec out;
46 xcorr_old(x, x, out, max_lag, scaleopt);
47 return out;
48}
49
50vec xcorr(const vec &x, const int max_lag, const std::string scaleopt)
51{
52 cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
53 xcorr(to_cvec(x), to_cvec(x), out, max_lag, scaleopt, true);
54
55 return real(out);
56}
57
58cvec xcorr(const cvec &x, const int max_lag, const std::string scaleopt)
59{
60 cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
61 xcorr(x, x, out, max_lag, scaleopt, true);
62
63 return out;
64}
65
66vec xcorr(const vec &x, const vec &y, const int max_lag, const std::string scaleopt)
67{
68 cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
69 xcorr(to_cvec(x), to_cvec(y), out, max_lag, scaleopt, false);
70
71 return real(out);
72}
73
74cvec xcorr(const cvec &x, const cvec &y, const int max_lag, const std::string scaleopt)
75{
76 cvec out(2*x.length() - 1); //Initial size does ont matter, it will get adjusted
77 xcorr(x, y, out, max_lag, scaleopt, false);
78
79 return out;
80}
81
82void xcorr(const vec &x, const vec &y, vec &out, const int max_lag, const std::string scaleopt)
83{
84 cvec xx = to_cvec(x);
85 cvec yy = to_cvec(y);
86 cvec oo = to_cvec(out);
87 xcorr(xx, yy, oo, max_lag, scaleopt, false);
88
89 out = real(oo);
90}
91
92void xcorr_old(const vec &x, const vec &y, vec &out, const int max_lag, const std::string scaleopt)
93{
94 int m, n;
95 double s_plus, s_minus, M_double, coeff_scale = 0.0;
96 int M, N;
97
98 M = x.size();
99 M = std::max(x.size(), y.size());
100 M_double = double(M);
101
102 if (max_lag == -1) {
103 N = std::max(x.size(), y.size());
104 }
105 else {
106 N = max_lag + 1;
107 }
108
109 out.set_size(2*N - 1, false);
110
111 it_assert(N <= std::max(x.size(), y.size()), "max_lag cannot be as large as, or larger than, the maximum length of x and y.");
112
113 if (scaleopt == "coeff") {
114 coeff_scale = std::sqrt(energy(x)) * std::sqrt(energy(y));
115 }
116
117 for (m = 0; m < N; m++) {
118 s_plus = 0;
119 s_minus = 0;
120
121 for (n = 0;n < M - m;n++) {
122 s_minus += index_zero_pad(x, n) * index_zero_pad(y, n + m);
123 s_plus += index_zero_pad(x, n + m) * index_zero_pad(y, n);
124 }
125
126 if (scaleopt == "none") {
127 out(N + m - 1) = s_plus;
128 out(N - m - 1) = s_minus;
129 }
130 else if (scaleopt == "biased") {
131 out(N + m - 1) = s_plus / M_double;
132 out(N - m - 1) = s_minus / M_double;
133 }
134 else if (scaleopt == "unbiased") {
135 out(N + m - 1) = s_plus / double(M - m);
136 out(N - m - 1) = s_minus / double(M - m);
137 }
138 else if (scaleopt == "coeff") {
139 out(N + m - 1) = s_plus / coeff_scale;
140 out(N - m - 1) = s_minus / coeff_scale;
141 }
142 else
143 it_error("Incorrect scaleopt specified.");
144 }
145}
146
147
148vec xcorr_old(const vec &x, const vec &y, const int max_lag, const std::string scaleopt)
149{
150 vec out;
151 xcorr_old(x, y, out, max_lag, scaleopt);
152 return out;
153}
154
155//Correlation
156void xcorr(const cvec &x, const cvec &y, cvec &out, const int max_lag, const std::string scaleopt, bool autoflag)
157{
158 int N = std::max(x.length(), y.length());
159
160 //Compute the FFT size as the "next power of 2" of the input vector's length (max)
161 int b = ceil_i(::log2(2.0 * N - 1));
162 int fftsize = pow2i(b);
163
164 int end = fftsize - 1;
165
166 cvec temp2;
167 if (autoflag == true) {
168 //Take FFT of input vector
169 cvec X = fft(zero_pad(x, fftsize));
170
171 //Compute the abs(X).^2 and take the inverse FFT.
172 temp2 = ifft(elem_mult(X, conj(X)));
173 }
174 else {
175 //Take FFT of input vectors
176 cvec X = fft(zero_pad(x, fftsize));
177 cvec Y = fft(zero_pad(y, fftsize));
178
179 //Compute the crosscorrelation
180 temp2 = ifft(elem_mult(X, conj(Y)));
181 }
182
183 // Compute the total number of lags to keep. We truncate the maximum number of lags to N-1.
184 int maxlag;
185 if ((max_lag == -1) || (max_lag >= N))
186 maxlag = N - 1;
187 else
188 maxlag = max_lag;
189
190
191 //Move negative lags to the beginning of the vector. Drop extra values from the FFT/IFFt
192 if (maxlag == 0) {
193 out.set_size(1, false);
194 out = temp2(0);
195 }
196 else
197 out = concat(temp2(end - maxlag + 1, end), temp2(0, maxlag));
198
199
200 //Scale data
201 if (scaleopt == "biased")
202 //out = out / static_cast<double_complex>(N);
203 out = out / static_cast<std::complex<double> >(N);
204 else if (scaleopt == "unbiased") {
205 //Total lag vector
206 vec lags = linspace(-maxlag, maxlag, 2 * maxlag + 1);
207 cvec scale = to_cvec(static_cast<double>(N) - abs(lags));
208 out /= scale;
209 }
210 else if (scaleopt == "coeff") {
211 if (autoflag == true) // Normalize by Rxx(0)
212 out /= out(maxlag);
213 else { //Normalize by sqrt(Rxx(0)*Ryy(0))
214 double rxx0 = sum(abs(elem_mult(x, x)));
215 double ryy0 = sum(abs(elem_mult(y, y)));
216 out /= std::sqrt(rxx0 * ryy0);
217 }
218 }
219 else if (scaleopt == "none") {}
220 else
221 it_warning("Unknow scaling option in XCORR, defaulting to <none> ");
222
223}
224
225
226mat cov(const mat &X, bool is_zero_mean)
227{
228 int d = X.cols(), n = X.rows();
229 mat R(d, d), m2(n, d);
230 vec tmp;
231
232 R = 0.0;
233
234 if (!is_zero_mean) {
235 // Compute and remove mean
236 for (int i = 0; i < d; i++) {
237 tmp = X.get_col(i);
238 m2.set_col(i, tmp - mean(tmp));
239 }
240
241 // Calc corr matrix
242 for (int i = 0; i < d; i++) {
243 for (int j = 0; j <= i; j++) {
244 for (int k = 0; k < n; k++) {
245 R(i, j) += m2(k, i) * m2(k, j);
246 }
247 R(j, i) = R(i, j); // When i=j this is unnecassary work
248 }
249 }
250 }
251 else {
252 // Calc corr matrix
253 for (int i = 0; i < d; i++) {
254 for (int j = 0; j <= i; j++) {
255 for (int k = 0; k < n; k++) {
256 R(i, j) += X(k, i) * X(k, j);
257 }
258 R(j, i) = R(i, j); // When i=j this is unnecassary work
259 }
260 }
261 }
262 R /= n;
263
264 return R;
265}
266
267vec spectrum(const vec &v, int nfft, int noverlap)
268{
269 it_assert_debug(pow2i(levels2bits(nfft)) == nfft,
270 "nfft must be a power of two in spectrum()!");
271
272 vec P(nfft / 2 + 1), w(nfft), wd(nfft);
273
274 P = 0.0;
275 w = hanning(nfft);
276 double w_energy = nfft == 1 ? 1 : (nfft + 1) * .375; // Hanning energy
277
278 if (nfft > v.size()) {
279 P = sqr(abs(fft(to_cvec(elem_mult(zero_pad(v, nfft), w)))(0, nfft / 2)));
280 P /= w_energy;
281 }
282 else {
283 int k = (v.size() - noverlap) / (nfft - noverlap), idx = 0;
284 for (int i = 0; i < k; i++) {
285 wd = elem_mult(v(idx, idx + nfft - 1), w);
286 P += sqr(abs(fft(to_cvec(wd))(0, nfft / 2)));
287 idx += nfft - noverlap;
288 }
289 P /= k * w_energy;
290 }
291
292 P.set_size(nfft / 2 + 1, true);
293 return P;
294}
295
296vec spectrum(const vec &v, const vec &w, int noverlap)
297{
298 int nfft = w.size();
299 it_assert_debug(pow2i(levels2bits(nfft)) == nfft,
300 "The window size must be a power of two in spectrum()!");
301
302 vec P(nfft / 2 + 1), wd(nfft);
303
304 P = 0.0;
305 double w_energy = energy(w);
306
307 if (nfft > v.size()) {
308 P = sqr(abs(fft(to_cvec(elem_mult(zero_pad(v, nfft), w)))(0, nfft / 2)));
309 P /= w_energy;
310 }
311 else {
312 int k = (v.size() - noverlap) / (nfft - noverlap), idx = 0;
313 for (int i = 0; i < k; i++) {
314 wd = elem_mult(v(idx, idx + nfft - 1), w);
315 P += sqr(abs(fft(to_cvec(wd))(0, nfft / 2)));
316 idx += nfft - noverlap;
317 }
318 P /= k * w_energy;
319 }
320
321 P.set_size(nfft / 2 + 1, true);
322 return P;
323}
324
325vec filter_spectrum(const vec &a, int nfft)
326{
327 vec s = sqr(abs(fft(to_cvec(a), nfft)));
328 s.set_size(nfft / 2 + 1, true);
329 return s;
330}
331
332vec filter_spectrum(const vec &a, const vec &b, int nfft)
333{
334 vec s = sqr(abs(elem_div(fft(to_cvec(a), nfft), fft(to_cvec(b), nfft))));
335 s.set_size(nfft / 2 + 1, true);
336 return s;
337}
338
339} // namespace itpp
340
341
342
343
Definitions of converters between different vector and matrix types.
Elementary mathematical functions - header file.
#define it_error(s)
Abort unconditionally.
Definition: itassert.h:126
#define it_warning(s)
Display a warning message.
Definition: itassert.h:173
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition: itassert.h:107
#define it_assert(t, s)
Abort if t is not true.
Definition: itassert.h:94
ITPP_EXPORT void ifft(const cvec &in, cvec &out)
Inverse Fast Fourier Transform.
ITPP_EXPORT void fft(const cvec &in, cvec &out)
Fast Fourier Transform.
int pow2i(int x)
Calculate two to the power of x (2^x); x is integer.
Definition: log_exp.h:53
vec log2(const vec &x)
log-2 of the elements
Definition: log_exp.cpp:36
int levels2bits(int n)
Calculate the number of bits needed to represent n different values (levels).
Definition: log_exp.h:92
T index_zero_pad(const Vec< T > &v, const int index)
Definition: matfunc.h:297
Vec< T > zero_pad(const Vec< T > &v, int n)
Zero-pad a vector to size n.
Definition: matfunc.h:257
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
vec real(const cvec &data)
Real part of complex values.
Definition: elem_math.cpp:157
vec sqr(const cvec &data)
Absolute square of elements.
Definition: elem_math.cpp:36
cvec conj(const cvec &x)
Conjugate of complex value.
Definition: elem_math.h:226
vec filter_spectrum(const vec &a, int nfft)
Power spectrum calculation of a filter.
Definition: sigfun.cpp:325
vec spectrum(const vec &v, int nfft, int noverlap)
Power spectrum calculation.
Definition: sigfun.cpp:267
mat cov(const mat &X, bool is_zero_mean)
Covariance matrix calculation.
Definition: sigfun.cpp:226
vec xcorr_old(const vec &x, const int max_lag, const std::string scaleopt)
Auto-correlation calculation.
Definition: sigfun.cpp:43
vec xcorr(const vec &x, const int max_lag, const std::string scaleopt)
Auto-correlation calculation.
Definition: sigfun.cpp:50
vec linspace(double from, double to, int points)
linspace (works in the same way as the MATLAB version)
Definition: specmat.cpp:106
double energy(const Vec< T > &v)
Calculate the energy: squared 2-norm. energy(v)=sum(abs(v).^2)
Definition: misc_stat.h:253
double mean(const vec &v)
The mean value.
Definition: misc_stat.cpp:36
vec hanning(int n)
Hanning window.
Definition: window.cpp:56
IT++ compatibility types and functions.
Various functions on vectors and matrices - header file.
Miscellaneous statistics functions and classes - header file.
itpp namespace
Definition: itmex.h:37
cvec to_cvec(const Vec< T > &v)
Converts a Vec<T> to cvec.
Definition: converters.h:107
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.
Definition: array.h:486
int ceil_i(double x)
The nearest larger integer.
Definition: converters.h:339
Mat< Num_T > elem_div(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise division of two matrices.
Definition: mat.h:1688
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
Definition: mat.h:1582
bin abs(const bin &inbin)
absolute value of bin
Definition: binary.h:174
Definitions of signal processing functions.
Definitions of special vectors and matrices.
Fourier, Hadamard, Walsh-Hadamard, and 2D Hadamard transforms - header file.
Definitions of window functions.
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.6