1 /*************************************************************************
2 ALGLIB 3.18.0 (source code generated 2021-10-25)
3 Copyright (c) Sergey Bochkanov (ALGLIB project).
4
5 >>> SOURCE LICENSE >>>
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation (www.fsf.org); either version 2 of the
9 License, or (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 A copy of the GNU General Public License is available at
17 http://www.fsf.org/licensing/licenses
18 >>> END OF LICENSE >>>
19 *************************************************************************/
20 #ifdef _MSC_VER
21 #define _CRT_SECURE_NO_WARNINGS
22 #endif
23 #include "stdafx.h"
24 #include "fasttransforms.h"
25
26 // disable some irrelevant warnings
27 #if (AE_COMPILER==AE_MSVC) && !defined(AE_ALL_WARNINGS)
28 #pragma warning(disable:4100)
29 #pragma warning(disable:4127)
30 #pragma warning(disable:4611)
31 #pragma warning(disable:4702)
32 #pragma warning(disable:4996)
33 #endif
34
35 /////////////////////////////////////////////////////////////////////////
36 //
37 // THIS SECTION CONTAINS IMPLEMENTATION OF C++ INTERFACE
38 //
39 /////////////////////////////////////////////////////////////////////////
40 namespace alglib
41 {
42
43 #if defined(AE_COMPILE_FFT) || !defined(AE_PARTIAL_BUILD)
44
45 #endif
46
47 #if defined(AE_COMPILE_FHT) || !defined(AE_PARTIAL_BUILD)
48
49 #endif
50
51 #if defined(AE_COMPILE_CONV) || !defined(AE_PARTIAL_BUILD)
52
53 #endif
54
55 #if defined(AE_COMPILE_CORR) || !defined(AE_PARTIAL_BUILD)
56
57 #endif
58
59 #if defined(AE_COMPILE_FFT) || !defined(AE_PARTIAL_BUILD)
60 /*************************************************************************
61 1-dimensional complex FFT.
62
63 Array size N may be arbitrary number (composite or prime). Composite N's
64 are handled with cache-oblivious variation of a Cooley-Tukey algorithm.
65 Small prime-factors are transformed using hard coded codelets (similar to
66 FFTW codelets, but without low-level optimization), large prime-factors
67 are handled with Bluestein's algorithm.
68
69 Fastests transforms are for smooth N's (prime factors are 2, 3, 5 only),
70 most fast for powers of 2. When N have prime factors larger than these,
71 but orders of magnitude smaller than N, computations will be about 4 times
72 slower than for nearby highly composite N's. When N itself is prime, speed
73 will be 6 times lower.
74
75 Algorithm has O(N*logN) complexity for any N (composite or prime).
76
77 INPUT PARAMETERS
78 A - array[0..N-1] - complex function to be transformed
79 N - problem size
80
81 OUTPUT PARAMETERS
82 A - DFT of a input array, array[0..N-1]
83 A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
84
85
86 -- ALGLIB --
87 Copyright 29.05.2009 by Bochkanov Sergey
88 *************************************************************************/
fftc1d(complex_1d_array & a,const ae_int_t n,const xparams _xparams)89 void fftc1d(complex_1d_array &a, const ae_int_t n, const xparams _xparams)
90 {
91 jmp_buf _break_jump;
92 alglib_impl::ae_state _alglib_env_state;
93 alglib_impl::ae_state_init(&_alglib_env_state);
94 if( setjmp(_break_jump) )
95 {
96 #if !defined(AE_NO_EXCEPTIONS)
97 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
98 #else
99 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
100 return;
101 #endif
102 }
103 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
104 if( _xparams.flags!=0x0 )
105 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
106 alglib_impl::fftc1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
107 alglib_impl::ae_state_clear(&_alglib_env_state);
108 return;
109 }
110
111 /*************************************************************************
112 1-dimensional complex FFT.
113
114 Array size N may be arbitrary number (composite or prime). Composite N's
115 are handled with cache-oblivious variation of a Cooley-Tukey algorithm.
116 Small prime-factors are transformed using hard coded codelets (similar to
117 FFTW codelets, but without low-level optimization), large prime-factors
118 are handled with Bluestein's algorithm.
119
120 Fastests transforms are for smooth N's (prime factors are 2, 3, 5 only),
121 most fast for powers of 2. When N have prime factors larger than these,
122 but orders of magnitude smaller than N, computations will be about 4 times
123 slower than for nearby highly composite N's. When N itself is prime, speed
124 will be 6 times lower.
125
126 Algorithm has O(N*logN) complexity for any N (composite or prime).
127
128 INPUT PARAMETERS
129 A - array[0..N-1] - complex function to be transformed
130 N - problem size
131
132 OUTPUT PARAMETERS
133 A - DFT of a input array, array[0..N-1]
134 A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
135
136
137 -- ALGLIB --
138 Copyright 29.05.2009 by Bochkanov Sergey
139 *************************************************************************/
140 #if !defined(AE_NO_EXCEPTIONS)
fftc1d(complex_1d_array & a,const xparams _xparams)141 void fftc1d(complex_1d_array &a, const xparams _xparams)
142 {
143 jmp_buf _break_jump;
144 alglib_impl::ae_state _alglib_env_state;
145 ae_int_t n;
146
147 n = a.length();
148 alglib_impl::ae_state_init(&_alglib_env_state);
149 if( setjmp(_break_jump) )
150 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
151 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
152 if( _xparams.flags!=0x0 )
153 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
154 alglib_impl::fftc1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
155
156 alglib_impl::ae_state_clear(&_alglib_env_state);
157 return;
158 }
159 #endif
160
161 /*************************************************************************
162 1-dimensional complex inverse FFT.
163
164 Array size N may be arbitrary number (composite or prime). Algorithm has
165 O(N*logN) complexity for any N (composite or prime).
166
167 See FFTC1D() description for more information about algorithm performance.
168
169 INPUT PARAMETERS
170 A - array[0..N-1] - complex array to be transformed
171 N - problem size
172
173 OUTPUT PARAMETERS
174 A - inverse DFT of a input array, array[0..N-1]
175 A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
176
177
178 -- ALGLIB --
179 Copyright 29.05.2009 by Bochkanov Sergey
180 *************************************************************************/
fftc1dinv(complex_1d_array & a,const ae_int_t n,const xparams _xparams)181 void fftc1dinv(complex_1d_array &a, const ae_int_t n, const xparams _xparams)
182 {
183 jmp_buf _break_jump;
184 alglib_impl::ae_state _alglib_env_state;
185 alglib_impl::ae_state_init(&_alglib_env_state);
186 if( setjmp(_break_jump) )
187 {
188 #if !defined(AE_NO_EXCEPTIONS)
189 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
190 #else
191 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
192 return;
193 #endif
194 }
195 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
196 if( _xparams.flags!=0x0 )
197 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
198 alglib_impl::fftc1dinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
199 alglib_impl::ae_state_clear(&_alglib_env_state);
200 return;
201 }
202
203 /*************************************************************************
204 1-dimensional complex inverse FFT.
205
206 Array size N may be arbitrary number (composite or prime). Algorithm has
207 O(N*logN) complexity for any N (composite or prime).
208
209 See FFTC1D() description for more information about algorithm performance.
210
211 INPUT PARAMETERS
212 A - array[0..N-1] - complex array to be transformed
213 N - problem size
214
215 OUTPUT PARAMETERS
216 A - inverse DFT of a input array, array[0..N-1]
217 A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
218
219
220 -- ALGLIB --
221 Copyright 29.05.2009 by Bochkanov Sergey
222 *************************************************************************/
223 #if !defined(AE_NO_EXCEPTIONS)
fftc1dinv(complex_1d_array & a,const xparams _xparams)224 void fftc1dinv(complex_1d_array &a, const xparams _xparams)
225 {
226 jmp_buf _break_jump;
227 alglib_impl::ae_state _alglib_env_state;
228 ae_int_t n;
229
230 n = a.length();
231 alglib_impl::ae_state_init(&_alglib_env_state);
232 if( setjmp(_break_jump) )
233 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
234 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
235 if( _xparams.flags!=0x0 )
236 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
237 alglib_impl::fftc1dinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
238
239 alglib_impl::ae_state_clear(&_alglib_env_state);
240 return;
241 }
242 #endif
243
244 /*************************************************************************
245 1-dimensional real FFT.
246
247 Algorithm has O(N*logN) complexity for any N (composite or prime).
248
249 INPUT PARAMETERS
250 A - array[0..N-1] - real function to be transformed
251 N - problem size
252
253 OUTPUT PARAMETERS
254 F - DFT of a input array, array[0..N-1]
255 F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
256
257 NOTE:
258 F[] satisfies symmetry property F[k] = conj(F[N-k]), so just one half
259 of array is usually needed. But for convinience subroutine returns full
260 complex array (with frequencies above N/2), so its result may be used by
261 other FFT-related subroutines.
262
263
264 -- ALGLIB --
265 Copyright 01.06.2009 by Bochkanov Sergey
266 *************************************************************************/
fftr1d(const real_1d_array & a,const ae_int_t n,complex_1d_array & f,const xparams _xparams)267 void fftr1d(const real_1d_array &a, const ae_int_t n, complex_1d_array &f, const xparams _xparams)
268 {
269 jmp_buf _break_jump;
270 alglib_impl::ae_state _alglib_env_state;
271 alglib_impl::ae_state_init(&_alglib_env_state);
272 if( setjmp(_break_jump) )
273 {
274 #if !defined(AE_NO_EXCEPTIONS)
275 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
276 #else
277 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
278 return;
279 #endif
280 }
281 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
282 if( _xparams.flags!=0x0 )
283 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
284 alglib_impl::fftr1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(f.c_ptr()), &_alglib_env_state);
285 alglib_impl::ae_state_clear(&_alglib_env_state);
286 return;
287 }
288
289 /*************************************************************************
290 1-dimensional real FFT.
291
292 Algorithm has O(N*logN) complexity for any N (composite or prime).
293
294 INPUT PARAMETERS
295 A - array[0..N-1] - real function to be transformed
296 N - problem size
297
298 OUTPUT PARAMETERS
299 F - DFT of a input array, array[0..N-1]
300 F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
301
302 NOTE:
303 F[] satisfies symmetry property F[k] = conj(F[N-k]), so just one half
304 of array is usually needed. But for convinience subroutine returns full
305 complex array (with frequencies above N/2), so its result may be used by
306 other FFT-related subroutines.
307
308
309 -- ALGLIB --
310 Copyright 01.06.2009 by Bochkanov Sergey
311 *************************************************************************/
312 #if !defined(AE_NO_EXCEPTIONS)
fftr1d(const real_1d_array & a,complex_1d_array & f,const xparams _xparams)313 void fftr1d(const real_1d_array &a, complex_1d_array &f, const xparams _xparams)
314 {
315 jmp_buf _break_jump;
316 alglib_impl::ae_state _alglib_env_state;
317 ae_int_t n;
318
319 n = a.length();
320 alglib_impl::ae_state_init(&_alglib_env_state);
321 if( setjmp(_break_jump) )
322 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
323 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
324 if( _xparams.flags!=0x0 )
325 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
326 alglib_impl::fftr1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(f.c_ptr()), &_alglib_env_state);
327
328 alglib_impl::ae_state_clear(&_alglib_env_state);
329 return;
330 }
331 #endif
332
333 /*************************************************************************
334 1-dimensional real inverse FFT.
335
336 Algorithm has O(N*logN) complexity for any N (composite or prime).
337
338 INPUT PARAMETERS
339 F - array[0..floor(N/2)] - frequencies from forward real FFT
340 N - problem size
341
342 OUTPUT PARAMETERS
343 A - inverse DFT of a input array, array[0..N-1]
344
345 NOTE:
346 F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just one
347 half of frequencies array is needed - elements from 0 to floor(N/2). F[0]
348 is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd, then
349 F[floor(N/2)] has no special properties.
350
351 Relying on properties noted above, FFTR1DInv subroutine uses only elements
352 from 0th to floor(N/2)-th. It ignores imaginary part of F[0], and in case
353 N is even it ignores imaginary part of F[floor(N/2)] too.
354
355 When you call this function using full arguments list - "FFTR1DInv(F,N,A)"
356 - you can pass either either frequencies array with N elements or reduced
357 array with roughly N/2 elements - subroutine will successfully transform
358 both.
359
360 If you call this function using reduced arguments list - "FFTR1DInv(F,A)"
361 - you must pass FULL array with N elements (although higher N/2 are still
362 not used) because array size is used to automatically determine FFT length
363
364
365 -- ALGLIB --
366 Copyright 01.06.2009 by Bochkanov Sergey
367 *************************************************************************/
fftr1dinv(const complex_1d_array & f,const ae_int_t n,real_1d_array & a,const xparams _xparams)368 void fftr1dinv(const complex_1d_array &f, const ae_int_t n, real_1d_array &a, const xparams _xparams)
369 {
370 jmp_buf _break_jump;
371 alglib_impl::ae_state _alglib_env_state;
372 alglib_impl::ae_state_init(&_alglib_env_state);
373 if( setjmp(_break_jump) )
374 {
375 #if !defined(AE_NO_EXCEPTIONS)
376 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
377 #else
378 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
379 return;
380 #endif
381 }
382 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
383 if( _xparams.flags!=0x0 )
384 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
385 alglib_impl::fftr1dinv(const_cast<alglib_impl::ae_vector*>(f.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(a.c_ptr()), &_alglib_env_state);
386 alglib_impl::ae_state_clear(&_alglib_env_state);
387 return;
388 }
389
390 /*************************************************************************
391 1-dimensional real inverse FFT.
392
393 Algorithm has O(N*logN) complexity for any N (composite or prime).
394
395 INPUT PARAMETERS
396 F - array[0..floor(N/2)] - frequencies from forward real FFT
397 N - problem size
398
399 OUTPUT PARAMETERS
400 A - inverse DFT of a input array, array[0..N-1]
401
402 NOTE:
403 F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just one
404 half of frequencies array is needed - elements from 0 to floor(N/2). F[0]
405 is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd, then
406 F[floor(N/2)] has no special properties.
407
408 Relying on properties noted above, FFTR1DInv subroutine uses only elements
409 from 0th to floor(N/2)-th. It ignores imaginary part of F[0], and in case
410 N is even it ignores imaginary part of F[floor(N/2)] too.
411
412 When you call this function using full arguments list - "FFTR1DInv(F,N,A)"
413 - you can pass either either frequencies array with N elements or reduced
414 array with roughly N/2 elements - subroutine will successfully transform
415 both.
416
417 If you call this function using reduced arguments list - "FFTR1DInv(F,A)"
418 - you must pass FULL array with N elements (although higher N/2 are still
419 not used) because array size is used to automatically determine FFT length
420
421
422 -- ALGLIB --
423 Copyright 01.06.2009 by Bochkanov Sergey
424 *************************************************************************/
425 #if !defined(AE_NO_EXCEPTIONS)
fftr1dinv(const complex_1d_array & f,real_1d_array & a,const xparams _xparams)426 void fftr1dinv(const complex_1d_array &f, real_1d_array &a, const xparams _xparams)
427 {
428 jmp_buf _break_jump;
429 alglib_impl::ae_state _alglib_env_state;
430 ae_int_t n;
431
432 n = f.length();
433 alglib_impl::ae_state_init(&_alglib_env_state);
434 if( setjmp(_break_jump) )
435 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
436 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
437 if( _xparams.flags!=0x0 )
438 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
439 alglib_impl::fftr1dinv(const_cast<alglib_impl::ae_vector*>(f.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(a.c_ptr()), &_alglib_env_state);
440
441 alglib_impl::ae_state_clear(&_alglib_env_state);
442 return;
443 }
444 #endif
445 #endif
446
447 #if defined(AE_COMPILE_FHT) || !defined(AE_PARTIAL_BUILD)
448 /*************************************************************************
449 1-dimensional Fast Hartley Transform.
450
451 Algorithm has O(N*logN) complexity for any N (composite or prime).
452
453 INPUT PARAMETERS
454 A - array[0..N-1] - real function to be transformed
455 N - problem size
456
457 OUTPUT PARAMETERS
458 A - FHT of a input array, array[0..N-1],
459 A_out[k] = sum(A_in[j]*(cos(2*pi*j*k/N)+sin(2*pi*j*k/N)), j=0..N-1)
460
461
462 -- ALGLIB --
463 Copyright 04.06.2009 by Bochkanov Sergey
464 *************************************************************************/
fhtr1d(real_1d_array & a,const ae_int_t n,const xparams _xparams)465 void fhtr1d(real_1d_array &a, const ae_int_t n, const xparams _xparams)
466 {
467 jmp_buf _break_jump;
468 alglib_impl::ae_state _alglib_env_state;
469 alglib_impl::ae_state_init(&_alglib_env_state);
470 if( setjmp(_break_jump) )
471 {
472 #if !defined(AE_NO_EXCEPTIONS)
473 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
474 #else
475 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
476 return;
477 #endif
478 }
479 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
480 if( _xparams.flags!=0x0 )
481 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
482 alglib_impl::fhtr1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
483 alglib_impl::ae_state_clear(&_alglib_env_state);
484 return;
485 }
486
487 /*************************************************************************
488 1-dimensional inverse FHT.
489
490 Algorithm has O(N*logN) complexity for any N (composite or prime).
491
492 INPUT PARAMETERS
493 A - array[0..N-1] - complex array to be transformed
494 N - problem size
495
496 OUTPUT PARAMETERS
497 A - inverse FHT of a input array, array[0..N-1]
498
499
500 -- ALGLIB --
501 Copyright 29.05.2009 by Bochkanov Sergey
502 *************************************************************************/
fhtr1dinv(real_1d_array & a,const ae_int_t n,const xparams _xparams)503 void fhtr1dinv(real_1d_array &a, const ae_int_t n, const xparams _xparams)
504 {
505 jmp_buf _break_jump;
506 alglib_impl::ae_state _alglib_env_state;
507 alglib_impl::ae_state_init(&_alglib_env_state);
508 if( setjmp(_break_jump) )
509 {
510 #if !defined(AE_NO_EXCEPTIONS)
511 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
512 #else
513 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
514 return;
515 #endif
516 }
517 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
518 if( _xparams.flags!=0x0 )
519 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
520 alglib_impl::fhtr1dinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
521 alglib_impl::ae_state_clear(&_alglib_env_state);
522 return;
523 }
524 #endif
525
526 #if defined(AE_COMPILE_CONV) || !defined(AE_PARTIAL_BUILD)
527 /*************************************************************************
528 1-dimensional complex convolution.
529
530 For given A/B returns conv(A,B) (non-circular). Subroutine can automatically
531 choose between three implementations: straightforward O(M*N) formula for
532 very small N (or M), overlap-add algorithm for cases where max(M,N) is
533 significantly larger than min(M,N), but O(M*N) algorithm is too slow, and
534 general FFT-based formula for cases where two previois algorithms are too
535 slow.
536
537 Algorithm has max(M,N)*log(max(M,N)) complexity for any M/N.
538
539 INPUT PARAMETERS
540 A - array[0..M-1] - complex function to be transformed
541 M - problem size
542 B - array[0..N-1] - complex function to be transformed
543 N - problem size
544
545 OUTPUT PARAMETERS
546 R - convolution: A*B. array[0..N+M-2].
547
548 NOTE:
549 It is assumed that A is zero at T<0, B is zero too. If one or both
550 functions have non-zero values at negative T's, you can still use this
551 subroutine - just shift its result correspondingly.
552
553 -- ALGLIB --
554 Copyright 21.07.2009 by Bochkanov Sergey
555 *************************************************************************/
convc1d(const complex_1d_array & a,const ae_int_t m,const complex_1d_array & b,const ae_int_t n,complex_1d_array & r,const xparams _xparams)556 void convc1d(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r, const xparams _xparams)
557 {
558 jmp_buf _break_jump;
559 alglib_impl::ae_state _alglib_env_state;
560 alglib_impl::ae_state_init(&_alglib_env_state);
561 if( setjmp(_break_jump) )
562 {
563 #if !defined(AE_NO_EXCEPTIONS)
564 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
565 #else
566 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
567 return;
568 #endif
569 }
570 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
571 if( _xparams.flags!=0x0 )
572 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
573 alglib_impl::convc1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), &_alglib_env_state);
574 alglib_impl::ae_state_clear(&_alglib_env_state);
575 return;
576 }
577
578 /*************************************************************************
579 1-dimensional complex non-circular deconvolution (inverse of ConvC1D()).
580
581 Algorithm has M*log(M)) complexity for any M (composite or prime).
582
583 INPUT PARAMETERS
584 A - array[0..M-1] - convolved signal, A = conv(R, B)
585 M - convolved signal length
586 B - array[0..N-1] - response
587 N - response length, N<=M
588
589 OUTPUT PARAMETERS
590 R - deconvolved signal. array[0..M-N].
591
592 NOTE:
593 deconvolution is unstable process and may result in division by zero
594 (if your response function is degenerate, i.e. has zero Fourier coefficient).
595
596 NOTE:
597 It is assumed that A is zero at T<0, B is zero too. If one or both
598 functions have non-zero values at negative T's, you can still use this
599 subroutine - just shift its result correspondingly.
600
601 -- ALGLIB --
602 Copyright 21.07.2009 by Bochkanov Sergey
603 *************************************************************************/
convc1dinv(const complex_1d_array & a,const ae_int_t m,const complex_1d_array & b,const ae_int_t n,complex_1d_array & r,const xparams _xparams)604 void convc1dinv(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r, const xparams _xparams)
605 {
606 jmp_buf _break_jump;
607 alglib_impl::ae_state _alglib_env_state;
608 alglib_impl::ae_state_init(&_alglib_env_state);
609 if( setjmp(_break_jump) )
610 {
611 #if !defined(AE_NO_EXCEPTIONS)
612 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
613 #else
614 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
615 return;
616 #endif
617 }
618 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
619 if( _xparams.flags!=0x0 )
620 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
621 alglib_impl::convc1dinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), &_alglib_env_state);
622 alglib_impl::ae_state_clear(&_alglib_env_state);
623 return;
624 }
625
626 /*************************************************************************
627 1-dimensional circular complex convolution.
628
629 For given S/R returns conv(S,R) (circular). Algorithm has linearithmic
630 complexity for any M/N.
631
632 IMPORTANT: normal convolution is commutative, i.e. it is symmetric -
633 conv(A,B)=conv(B,A). Cyclic convolution IS NOT. One function - S - is a
634 signal, periodic function, and another - R - is a response, non-periodic
635 function with limited length.
636
637 INPUT PARAMETERS
638 S - array[0..M-1] - complex periodic signal
639 M - problem size
640 B - array[0..N-1] - complex non-periodic response
641 N - problem size
642
643 OUTPUT PARAMETERS
644 R - convolution: A*B. array[0..M-1].
645
646 NOTE:
647 It is assumed that B is zero at T<0. If it has non-zero values at
648 negative T's, you can still use this subroutine - just shift its result
649 correspondingly.
650
651 -- ALGLIB --
652 Copyright 21.07.2009 by Bochkanov Sergey
653 *************************************************************************/
convc1dcircular(const complex_1d_array & s,const ae_int_t m,const complex_1d_array & r,const ae_int_t n,complex_1d_array & c,const xparams _xparams)654 void convc1dcircular(const complex_1d_array &s, const ae_int_t m, const complex_1d_array &r, const ae_int_t n, complex_1d_array &c, const xparams _xparams)
655 {
656 jmp_buf _break_jump;
657 alglib_impl::ae_state _alglib_env_state;
658 alglib_impl::ae_state_init(&_alglib_env_state);
659 if( setjmp(_break_jump) )
660 {
661 #if !defined(AE_NO_EXCEPTIONS)
662 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
663 #else
664 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
665 return;
666 #endif
667 }
668 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
669 if( _xparams.flags!=0x0 )
670 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
671 alglib_impl::convc1dcircular(const_cast<alglib_impl::ae_vector*>(s.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(c.c_ptr()), &_alglib_env_state);
672 alglib_impl::ae_state_clear(&_alglib_env_state);
673 return;
674 }
675
676 /*************************************************************************
677 1-dimensional circular complex deconvolution (inverse of ConvC1DCircular()).
678
679 Algorithm has M*log(M)) complexity for any M (composite or prime).
680
681 INPUT PARAMETERS
682 A - array[0..M-1] - convolved periodic signal, A = conv(R, B)
683 M - convolved signal length
684 B - array[0..N-1] - non-periodic response
685 N - response length
686
687 OUTPUT PARAMETERS
688 R - deconvolved signal. array[0..M-1].
689
690 NOTE:
691 deconvolution is unstable process and may result in division by zero
692 (if your response function is degenerate, i.e. has zero Fourier coefficient).
693
694 NOTE:
695 It is assumed that B is zero at T<0. If it has non-zero values at
696 negative T's, you can still use this subroutine - just shift its result
697 correspondingly.
698
699 -- ALGLIB --
700 Copyright 21.07.2009 by Bochkanov Sergey
701 *************************************************************************/
convc1dcircularinv(const complex_1d_array & a,const ae_int_t m,const complex_1d_array & b,const ae_int_t n,complex_1d_array & r,const xparams _xparams)702 void convc1dcircularinv(const complex_1d_array &a, const ae_int_t m, const complex_1d_array &b, const ae_int_t n, complex_1d_array &r, const xparams _xparams)
703 {
704 jmp_buf _break_jump;
705 alglib_impl::ae_state _alglib_env_state;
706 alglib_impl::ae_state_init(&_alglib_env_state);
707 if( setjmp(_break_jump) )
708 {
709 #if !defined(AE_NO_EXCEPTIONS)
710 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
711 #else
712 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
713 return;
714 #endif
715 }
716 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
717 if( _xparams.flags!=0x0 )
718 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
719 alglib_impl::convc1dcircularinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), &_alglib_env_state);
720 alglib_impl::ae_state_clear(&_alglib_env_state);
721 return;
722 }
723
724 /*************************************************************************
725 1-dimensional real convolution.
726
727 Analogous to ConvC1D(), see ConvC1D() comments for more details.
728
729 INPUT PARAMETERS
730 A - array[0..M-1] - real function to be transformed
731 M - problem size
732 B - array[0..N-1] - real function to be transformed
733 N - problem size
734
735 OUTPUT PARAMETERS
736 R - convolution: A*B. array[0..N+M-2].
737
738 NOTE:
739 It is assumed that A is zero at T<0, B is zero too. If one or both
740 functions have non-zero values at negative T's, you can still use this
741 subroutine - just shift its result correspondingly.
742
743 -- ALGLIB --
744 Copyright 21.07.2009 by Bochkanov Sergey
745 *************************************************************************/
convr1d(const real_1d_array & a,const ae_int_t m,const real_1d_array & b,const ae_int_t n,real_1d_array & r,const xparams _xparams)746 void convr1d(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r, const xparams _xparams)
747 {
748 jmp_buf _break_jump;
749 alglib_impl::ae_state _alglib_env_state;
750 alglib_impl::ae_state_init(&_alglib_env_state);
751 if( setjmp(_break_jump) )
752 {
753 #if !defined(AE_NO_EXCEPTIONS)
754 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
755 #else
756 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
757 return;
758 #endif
759 }
760 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
761 if( _xparams.flags!=0x0 )
762 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
763 alglib_impl::convr1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), &_alglib_env_state);
764 alglib_impl::ae_state_clear(&_alglib_env_state);
765 return;
766 }
767
768 /*************************************************************************
769 1-dimensional real deconvolution (inverse of ConvC1D()).
770
771 Algorithm has M*log(M)) complexity for any M (composite or prime).
772
773 INPUT PARAMETERS
774 A - array[0..M-1] - convolved signal, A = conv(R, B)
775 M - convolved signal length
776 B - array[0..N-1] - response
777 N - response length, N<=M
778
779 OUTPUT PARAMETERS
780 R - deconvolved signal. array[0..M-N].
781
782 NOTE:
783 deconvolution is unstable process and may result in division by zero
784 (if your response function is degenerate, i.e. has zero Fourier coefficient).
785
786 NOTE:
787 It is assumed that A is zero at T<0, B is zero too. If one or both
788 functions have non-zero values at negative T's, you can still use this
789 subroutine - just shift its result correspondingly.
790
791 -- ALGLIB --
792 Copyright 21.07.2009 by Bochkanov Sergey
793 *************************************************************************/
convr1dinv(const real_1d_array & a,const ae_int_t m,const real_1d_array & b,const ae_int_t n,real_1d_array & r,const xparams _xparams)794 void convr1dinv(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r, const xparams _xparams)
795 {
796 jmp_buf _break_jump;
797 alglib_impl::ae_state _alglib_env_state;
798 alglib_impl::ae_state_init(&_alglib_env_state);
799 if( setjmp(_break_jump) )
800 {
801 #if !defined(AE_NO_EXCEPTIONS)
802 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
803 #else
804 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
805 return;
806 #endif
807 }
808 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
809 if( _xparams.flags!=0x0 )
810 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
811 alglib_impl::convr1dinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), &_alglib_env_state);
812 alglib_impl::ae_state_clear(&_alglib_env_state);
813 return;
814 }
815
816 /*************************************************************************
817 1-dimensional circular real convolution.
818
819 Analogous to ConvC1DCircular(), see ConvC1DCircular() comments for more details.
820
821 INPUT PARAMETERS
822 S - array[0..M-1] - real signal
823 M - problem size
824 B - array[0..N-1] - real response
825 N - problem size
826
827 OUTPUT PARAMETERS
828 R - convolution: A*B. array[0..M-1].
829
830 NOTE:
831 It is assumed that B is zero at T<0. If it has non-zero values at
832 negative T's, you can still use this subroutine - just shift its result
833 correspondingly.
834
835 -- ALGLIB --
836 Copyright 21.07.2009 by Bochkanov Sergey
837 *************************************************************************/
convr1dcircular(const real_1d_array & s,const ae_int_t m,const real_1d_array & r,const ae_int_t n,real_1d_array & c,const xparams _xparams)838 void convr1dcircular(const real_1d_array &s, const ae_int_t m, const real_1d_array &r, const ae_int_t n, real_1d_array &c, const xparams _xparams)
839 {
840 jmp_buf _break_jump;
841 alglib_impl::ae_state _alglib_env_state;
842 alglib_impl::ae_state_init(&_alglib_env_state);
843 if( setjmp(_break_jump) )
844 {
845 #if !defined(AE_NO_EXCEPTIONS)
846 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
847 #else
848 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
849 return;
850 #endif
851 }
852 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
853 if( _xparams.flags!=0x0 )
854 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
855 alglib_impl::convr1dcircular(const_cast<alglib_impl::ae_vector*>(s.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(c.c_ptr()), &_alglib_env_state);
856 alglib_impl::ae_state_clear(&_alglib_env_state);
857 return;
858 }
859
860 /*************************************************************************
861 1-dimensional complex deconvolution (inverse of ConvC1D()).
862
863 Algorithm has M*log(M)) complexity for any M (composite or prime).
864
865 INPUT PARAMETERS
866 A - array[0..M-1] - convolved signal, A = conv(R, B)
867 M - convolved signal length
868 B - array[0..N-1] - response
869 N - response length
870
871 OUTPUT PARAMETERS
872 R - deconvolved signal. array[0..M-N].
873
874 NOTE:
875 deconvolution is unstable process and may result in division by zero
876 (if your response function is degenerate, i.e. has zero Fourier coefficient).
877
878 NOTE:
879 It is assumed that B is zero at T<0. If it has non-zero values at
880 negative T's, you can still use this subroutine - just shift its result
881 correspondingly.
882
883 -- ALGLIB --
884 Copyright 21.07.2009 by Bochkanov Sergey
885 *************************************************************************/
convr1dcircularinv(const real_1d_array & a,const ae_int_t m,const real_1d_array & b,const ae_int_t n,real_1d_array & r,const xparams _xparams)886 void convr1dcircularinv(const real_1d_array &a, const ae_int_t m, const real_1d_array &b, const ae_int_t n, real_1d_array &r, const xparams _xparams)
887 {
888 jmp_buf _break_jump;
889 alglib_impl::ae_state _alglib_env_state;
890 alglib_impl::ae_state_init(&_alglib_env_state);
891 if( setjmp(_break_jump) )
892 {
893 #if !defined(AE_NO_EXCEPTIONS)
894 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
895 #else
896 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
897 return;
898 #endif
899 }
900 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
901 if( _xparams.flags!=0x0 )
902 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
903 alglib_impl::convr1dcircularinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), &_alglib_env_state);
904 alglib_impl::ae_state_clear(&_alglib_env_state);
905 return;
906 }
907 #endif
908
909 #if defined(AE_COMPILE_CORR) || !defined(AE_PARTIAL_BUILD)
910 /*************************************************************************
911 1-dimensional complex cross-correlation.
912
913 For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).
914
915 Correlation is calculated using reduction to convolution. Algorithm with
916 max(N,N)*log(max(N,N)) complexity is used (see ConvC1D() for more info
917 about performance).
918
919 IMPORTANT:
920 for historical reasons subroutine accepts its parameters in reversed
921 order: CorrC1D(Signal, Pattern) = Pattern x Signal (using traditional
922 definition of cross-correlation, denoting cross-correlation as "x").
923
924 INPUT PARAMETERS
925 Signal - array[0..N-1] - complex function to be transformed,
926 signal containing pattern
927 N - problem size
928 Pattern - array[0..M-1] - complex function to be transformed,
929 pattern to search withing signal
930 M - problem size
931
932 OUTPUT PARAMETERS
933 R - cross-correlation, array[0..N+M-2]:
934 * positive lags are stored in R[0..N-1],
935 R[i] = sum(conj(pattern[j])*signal[i+j]
936 * negative lags are stored in R[N..N+M-2],
937 R[N+M-1-i] = sum(conj(pattern[j])*signal[-i+j]
938
939 NOTE:
940 It is assumed that pattern domain is [0..M-1]. If Pattern is non-zero
941 on [-K..M-1], you can still use this subroutine, just shift result by K.
942
943 -- ALGLIB --
944 Copyright 21.07.2009 by Bochkanov Sergey
945 *************************************************************************/
corrc1d(const complex_1d_array & signal,const ae_int_t n,const complex_1d_array & pattern,const ae_int_t m,complex_1d_array & r,const xparams _xparams)946 void corrc1d(const complex_1d_array &signal, const ae_int_t n, const complex_1d_array &pattern, const ae_int_t m, complex_1d_array &r, const xparams _xparams)
947 {
948 jmp_buf _break_jump;
949 alglib_impl::ae_state _alglib_env_state;
950 alglib_impl::ae_state_init(&_alglib_env_state);
951 if( setjmp(_break_jump) )
952 {
953 #if !defined(AE_NO_EXCEPTIONS)
954 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
955 #else
956 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
957 return;
958 #endif
959 }
960 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
961 if( _xparams.flags!=0x0 )
962 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
963 alglib_impl::corrc1d(const_cast<alglib_impl::ae_vector*>(signal.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(pattern.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), &_alglib_env_state);
964 alglib_impl::ae_state_clear(&_alglib_env_state);
965 return;
966 }
967
968 /*************************************************************************
969 1-dimensional circular complex cross-correlation.
970
971 For given Pattern/Signal returns corr(Pattern,Signal) (circular).
972 Algorithm has linearithmic complexity for any M/N.
973
974 IMPORTANT:
975 for historical reasons subroutine accepts its parameters in reversed
976 order: CorrC1DCircular(Signal, Pattern) = Pattern x Signal (using
977 traditional definition of cross-correlation, denoting cross-correlation
978 as "x").
979
980 INPUT PARAMETERS
981 Signal - array[0..N-1] - complex function to be transformed,
982 periodic signal containing pattern
983 N - problem size
984 Pattern - array[0..M-1] - complex function to be transformed,
985 non-periodic pattern to search withing signal
986 M - problem size
987
988 OUTPUT PARAMETERS
989 R - convolution: A*B. array[0..M-1].
990
991
992 -- ALGLIB --
993 Copyright 21.07.2009 by Bochkanov Sergey
994 *************************************************************************/
corrc1dcircular(const complex_1d_array & signal,const ae_int_t m,const complex_1d_array & pattern,const ae_int_t n,complex_1d_array & c,const xparams _xparams)995 void corrc1dcircular(const complex_1d_array &signal, const ae_int_t m, const complex_1d_array &pattern, const ae_int_t n, complex_1d_array &c, const xparams _xparams)
996 {
997 jmp_buf _break_jump;
998 alglib_impl::ae_state _alglib_env_state;
999 alglib_impl::ae_state_init(&_alglib_env_state);
1000 if( setjmp(_break_jump) )
1001 {
1002 #if !defined(AE_NO_EXCEPTIONS)
1003 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1004 #else
1005 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1006 return;
1007 #endif
1008 }
1009 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1010 if( _xparams.flags!=0x0 )
1011 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1012 alglib_impl::corrc1dcircular(const_cast<alglib_impl::ae_vector*>(signal.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(pattern.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(c.c_ptr()), &_alglib_env_state);
1013 alglib_impl::ae_state_clear(&_alglib_env_state);
1014 return;
1015 }
1016
1017 /*************************************************************************
1018 1-dimensional real cross-correlation.
1019
1020 For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).
1021
1022 Correlation is calculated using reduction to convolution. Algorithm with
1023 max(N,N)*log(max(N,N)) complexity is used (see ConvC1D() for more info
1024 about performance).
1025
1026 IMPORTANT:
1027 for historical reasons subroutine accepts its parameters in reversed
1028 order: CorrR1D(Signal, Pattern) = Pattern x Signal (using traditional
1029 definition of cross-correlation, denoting cross-correlation as "x").
1030
1031 INPUT PARAMETERS
1032 Signal - array[0..N-1] - real function to be transformed,
1033 signal containing pattern
1034 N - problem size
1035 Pattern - array[0..M-1] - real function to be transformed,
1036 pattern to search withing signal
1037 M - problem size
1038
1039 OUTPUT PARAMETERS
1040 R - cross-correlation, array[0..N+M-2]:
1041 * positive lags are stored in R[0..N-1],
1042 R[i] = sum(pattern[j]*signal[i+j]
1043 * negative lags are stored in R[N..N+M-2],
1044 R[N+M-1-i] = sum(pattern[j]*signal[-i+j]
1045
1046 NOTE:
1047 It is assumed that pattern domain is [0..M-1]. If Pattern is non-zero
1048 on [-K..M-1], you can still use this subroutine, just shift result by K.
1049
1050 -- ALGLIB --
1051 Copyright 21.07.2009 by Bochkanov Sergey
1052 *************************************************************************/
corrr1d(const real_1d_array & signal,const ae_int_t n,const real_1d_array & pattern,const ae_int_t m,real_1d_array & r,const xparams _xparams)1053 void corrr1d(const real_1d_array &signal, const ae_int_t n, const real_1d_array &pattern, const ae_int_t m, real_1d_array &r, const xparams _xparams)
1054 {
1055 jmp_buf _break_jump;
1056 alglib_impl::ae_state _alglib_env_state;
1057 alglib_impl::ae_state_init(&_alglib_env_state);
1058 if( setjmp(_break_jump) )
1059 {
1060 #if !defined(AE_NO_EXCEPTIONS)
1061 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1062 #else
1063 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1064 return;
1065 #endif
1066 }
1067 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1068 if( _xparams.flags!=0x0 )
1069 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1070 alglib_impl::corrr1d(const_cast<alglib_impl::ae_vector*>(signal.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(pattern.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(r.c_ptr()), &_alglib_env_state);
1071 alglib_impl::ae_state_clear(&_alglib_env_state);
1072 return;
1073 }
1074
1075 /*************************************************************************
1076 1-dimensional circular real cross-correlation.
1077
1078 For given Pattern/Signal returns corr(Pattern,Signal) (circular).
1079 Algorithm has linearithmic complexity for any M/N.
1080
1081 IMPORTANT:
1082 for historical reasons subroutine accepts its parameters in reversed
1083 order: CorrR1DCircular(Signal, Pattern) = Pattern x Signal (using
1084 traditional definition of cross-correlation, denoting cross-correlation
1085 as "x").
1086
1087 INPUT PARAMETERS
1088 Signal - array[0..N-1] - real function to be transformed,
1089 periodic signal containing pattern
1090 N - problem size
1091 Pattern - array[0..M-1] - real function to be transformed,
1092 non-periodic pattern to search withing signal
1093 M - problem size
1094
1095 OUTPUT PARAMETERS
1096 R - convolution: A*B. array[0..M-1].
1097
1098
1099 -- ALGLIB --
1100 Copyright 21.07.2009 by Bochkanov Sergey
1101 *************************************************************************/
corrr1dcircular(const real_1d_array & signal,const ae_int_t m,const real_1d_array & pattern,const ae_int_t n,real_1d_array & c,const xparams _xparams)1102 void corrr1dcircular(const real_1d_array &signal, const ae_int_t m, const real_1d_array &pattern, const ae_int_t n, real_1d_array &c, const xparams _xparams)
1103 {
1104 jmp_buf _break_jump;
1105 alglib_impl::ae_state _alglib_env_state;
1106 alglib_impl::ae_state_init(&_alglib_env_state);
1107 if( setjmp(_break_jump) )
1108 {
1109 #if !defined(AE_NO_EXCEPTIONS)
1110 _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1111 #else
1112 _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1113 return;
1114 #endif
1115 }
1116 ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1117 if( _xparams.flags!=0x0 )
1118 ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1119 alglib_impl::corrr1dcircular(const_cast<alglib_impl::ae_vector*>(signal.c_ptr()), m, const_cast<alglib_impl::ae_vector*>(pattern.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(c.c_ptr()), &_alglib_env_state);
1120 alglib_impl::ae_state_clear(&_alglib_env_state);
1121 return;
1122 }
1123 #endif
1124 }
1125
1126 /////////////////////////////////////////////////////////////////////////
1127 //
1128 // THIS SECTION CONTAINS IMPLEMENTATION OF COMPUTATIONAL CORE
1129 //
1130 /////////////////////////////////////////////////////////////////////////
1131 namespace alglib_impl
1132 {
1133 #if defined(AE_COMPILE_FFT) || !defined(AE_PARTIAL_BUILD)
1134
1135
1136 #endif
1137 #if defined(AE_COMPILE_FHT) || !defined(AE_PARTIAL_BUILD)
1138
1139
1140 #endif
1141 #if defined(AE_COMPILE_CONV) || !defined(AE_PARTIAL_BUILD)
1142
1143
1144 #endif
1145 #if defined(AE_COMPILE_CORR) || !defined(AE_PARTIAL_BUILD)
1146
1147
1148 #endif
1149
1150 #if defined(AE_COMPILE_FFT) || !defined(AE_PARTIAL_BUILD)
1151
1152
1153 /*************************************************************************
1154 1-dimensional complex FFT.
1155
1156 Array size N may be arbitrary number (composite or prime). Composite N's
1157 are handled with cache-oblivious variation of a Cooley-Tukey algorithm.
1158 Small prime-factors are transformed using hard coded codelets (similar to
1159 FFTW codelets, but without low-level optimization), large prime-factors
1160 are handled with Bluestein's algorithm.
1161
1162 Fastests transforms are for smooth N's (prime factors are 2, 3, 5 only),
1163 most fast for powers of 2. When N have prime factors larger than these,
1164 but orders of magnitude smaller than N, computations will be about 4 times
1165 slower than for nearby highly composite N's. When N itself is prime, speed
1166 will be 6 times lower.
1167
1168 Algorithm has O(N*logN) complexity for any N (composite or prime).
1169
1170 INPUT PARAMETERS
1171 A - array[0..N-1] - complex function to be transformed
1172 N - problem size
1173
1174 OUTPUT PARAMETERS
1175 A - DFT of a input array, array[0..N-1]
1176 A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
1177
1178
1179 -- ALGLIB --
1180 Copyright 29.05.2009 by Bochkanov Sergey
1181 *************************************************************************/
fftc1d(ae_vector * a,ae_int_t n,ae_state * _state)1182 void fftc1d(/* Complex */ ae_vector* a, ae_int_t n, ae_state *_state)
1183 {
1184 ae_frame _frame_block;
1185 fasttransformplan plan;
1186 ae_int_t i;
1187 ae_vector buf;
1188
1189 ae_frame_make(_state, &_frame_block);
1190 memset(&plan, 0, sizeof(plan));
1191 memset(&buf, 0, sizeof(buf));
1192 _fasttransformplan_init(&plan, _state, ae_true);
1193 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
1194
1195 ae_assert(n>0, "FFTC1D: incorrect N!", _state);
1196 ae_assert(a->cnt>=n, "FFTC1D: Length(A)<N!", _state);
1197 ae_assert(isfinitecvector(a, n, _state), "FFTC1D: A contains infinite or NAN values!", _state);
1198
1199 /*
1200 * Special case: N=1, FFT is just identity transform.
1201 * After this block we assume that N is strictly greater than 1.
1202 */
1203 if( n==1 )
1204 {
1205 ae_frame_leave(_state);
1206 return;
1207 }
1208
1209 /*
1210 * convert input array to the more convinient format
1211 */
1212 ae_vector_set_length(&buf, 2*n, _state);
1213 for(i=0; i<=n-1; i++)
1214 {
1215 buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
1216 buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
1217 }
1218
1219 /*
1220 * Generate plan and execute it.
1221 *
1222 * Plan is a combination of a successive factorizations of N and
1223 * precomputed data. It is much like a FFTW plan, but is not stored
1224 * between subroutine calls and is much simpler.
1225 */
1226 ftcomplexfftplan(n, 1, &plan, _state);
1227 ftapplyplan(&plan, &buf, 0, 1, _state);
1228
1229 /*
1230 * result
1231 */
1232 for(i=0; i<=n-1; i++)
1233 {
1234 a->ptr.p_complex[i].x = buf.ptr.p_double[2*i+0];
1235 a->ptr.p_complex[i].y = buf.ptr.p_double[2*i+1];
1236 }
1237 ae_frame_leave(_state);
1238 }
1239
1240
1241 /*************************************************************************
1242 1-dimensional complex inverse FFT.
1243
1244 Array size N may be arbitrary number (composite or prime). Algorithm has
1245 O(N*logN) complexity for any N (composite or prime).
1246
1247 See FFTC1D() description for more information about algorithm performance.
1248
1249 INPUT PARAMETERS
1250 A - array[0..N-1] - complex array to be transformed
1251 N - problem size
1252
1253 OUTPUT PARAMETERS
1254 A - inverse DFT of a input array, array[0..N-1]
1255 A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
1256
1257
1258 -- ALGLIB --
1259 Copyright 29.05.2009 by Bochkanov Sergey
1260 *************************************************************************/
fftc1dinv(ae_vector * a,ae_int_t n,ae_state * _state)1261 void fftc1dinv(/* Complex */ ae_vector* a, ae_int_t n, ae_state *_state)
1262 {
1263 ae_int_t i;
1264
1265
1266 ae_assert(n>0, "FFTC1DInv: incorrect N!", _state);
1267 ae_assert(a->cnt>=n, "FFTC1DInv: Length(A)<N!", _state);
1268 ae_assert(isfinitecvector(a, n, _state), "FFTC1DInv: A contains infinite or NAN values!", _state);
1269
1270 /*
1271 * Inverse DFT can be expressed in terms of the DFT as
1272 *
1273 * invfft(x) = fft(x')'/N
1274 *
1275 * here x' means conj(x).
1276 */
1277 for(i=0; i<=n-1; i++)
1278 {
1279 a->ptr.p_complex[i].y = -a->ptr.p_complex[i].y;
1280 }
1281 fftc1d(a, n, _state);
1282 for(i=0; i<=n-1; i++)
1283 {
1284 a->ptr.p_complex[i].x = a->ptr.p_complex[i].x/n;
1285 a->ptr.p_complex[i].y = -a->ptr.p_complex[i].y/n;
1286 }
1287 }
1288
1289
1290 /*************************************************************************
1291 1-dimensional real FFT.
1292
1293 Algorithm has O(N*logN) complexity for any N (composite or prime).
1294
1295 INPUT PARAMETERS
1296 A - array[0..N-1] - real function to be transformed
1297 N - problem size
1298
1299 OUTPUT PARAMETERS
1300 F - DFT of a input array, array[0..N-1]
1301 F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
1302
1303 NOTE:
1304 F[] satisfies symmetry property F[k] = conj(F[N-k]), so just one half
1305 of array is usually needed. But for convinience subroutine returns full
1306 complex array (with frequencies above N/2), so its result may be used by
1307 other FFT-related subroutines.
1308
1309
1310 -- ALGLIB --
1311 Copyright 01.06.2009 by Bochkanov Sergey
1312 *************************************************************************/
fftr1d(ae_vector * a,ae_int_t n,ae_vector * f,ae_state * _state)1313 void fftr1d(/* Real */ ae_vector* a,
1314 ae_int_t n,
1315 /* Complex */ ae_vector* f,
1316 ae_state *_state)
1317 {
1318 ae_frame _frame_block;
1319 ae_int_t i;
1320 ae_int_t n2;
1321 ae_int_t idx;
1322 ae_complex hn;
1323 ae_complex hmnc;
1324 ae_complex v;
1325 ae_vector buf;
1326 fasttransformplan plan;
1327
1328 ae_frame_make(_state, &_frame_block);
1329 memset(&buf, 0, sizeof(buf));
1330 memset(&plan, 0, sizeof(plan));
1331 ae_vector_clear(f);
1332 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
1333 _fasttransformplan_init(&plan, _state, ae_true);
1334
1335 ae_assert(n>0, "FFTR1D: incorrect N!", _state);
1336 ae_assert(a->cnt>=n, "FFTR1D: Length(A)<N!", _state);
1337 ae_assert(isfinitevector(a, n, _state), "FFTR1D: A contains infinite or NAN values!", _state);
1338
1339 /*
1340 * Special cases:
1341 * * N=1, FFT is just identity transform.
1342 * * N=2, FFT is simple too
1343 *
1344 * After this block we assume that N is strictly greater than 2
1345 */
1346 if( n==1 )
1347 {
1348 ae_vector_set_length(f, 1, _state);
1349 f->ptr.p_complex[0] = ae_complex_from_d(a->ptr.p_double[0]);
1350 ae_frame_leave(_state);
1351 return;
1352 }
1353 if( n==2 )
1354 {
1355 ae_vector_set_length(f, 2, _state);
1356 f->ptr.p_complex[0].x = a->ptr.p_double[0]+a->ptr.p_double[1];
1357 f->ptr.p_complex[0].y = (double)(0);
1358 f->ptr.p_complex[1].x = a->ptr.p_double[0]-a->ptr.p_double[1];
1359 f->ptr.p_complex[1].y = (double)(0);
1360 ae_frame_leave(_state);
1361 return;
1362 }
1363
1364 /*
1365 * Choose between odd-size and even-size FFTs
1366 */
1367 if( n%2==0 )
1368 {
1369
1370 /*
1371 * even-size real FFT, use reduction to the complex task
1372 */
1373 n2 = n/2;
1374 ae_vector_set_length(&buf, n, _state);
1375 ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,n-1));
1376 ftcomplexfftplan(n2, 1, &plan, _state);
1377 ftapplyplan(&plan, &buf, 0, 1, _state);
1378 ae_vector_set_length(f, n, _state);
1379 for(i=0; i<=n2; i++)
1380 {
1381 idx = 2*(i%n2);
1382 hn.x = buf.ptr.p_double[idx+0];
1383 hn.y = buf.ptr.p_double[idx+1];
1384 idx = 2*((n2-i)%n2);
1385 hmnc.x = buf.ptr.p_double[idx+0];
1386 hmnc.y = -buf.ptr.p_double[idx+1];
1387 v.x = -ae_sin(-2*ae_pi*i/n, _state);
1388 v.y = ae_cos(-2*ae_pi*i/n, _state);
1389 f->ptr.p_complex[i] = ae_c_sub(ae_c_add(hn,hmnc),ae_c_mul(v,ae_c_sub(hn,hmnc)));
1390 f->ptr.p_complex[i].x = 0.5*f->ptr.p_complex[i].x;
1391 f->ptr.p_complex[i].y = 0.5*f->ptr.p_complex[i].y;
1392 }
1393 for(i=n2+1; i<=n-1; i++)
1394 {
1395 f->ptr.p_complex[i] = ae_c_conj(f->ptr.p_complex[n-i], _state);
1396 }
1397 }
1398 else
1399 {
1400
1401 /*
1402 * use complex FFT
1403 */
1404 ae_vector_set_length(f, n, _state);
1405 for(i=0; i<=n-1; i++)
1406 {
1407 f->ptr.p_complex[i] = ae_complex_from_d(a->ptr.p_double[i]);
1408 }
1409 fftc1d(f, n, _state);
1410 }
1411 ae_frame_leave(_state);
1412 }
1413
1414
1415 /*************************************************************************
1416 1-dimensional real inverse FFT.
1417
1418 Algorithm has O(N*logN) complexity for any N (composite or prime).
1419
1420 INPUT PARAMETERS
1421 F - array[0..floor(N/2)] - frequencies from forward real FFT
1422 N - problem size
1423
1424 OUTPUT PARAMETERS
1425 A - inverse DFT of a input array, array[0..N-1]
1426
1427 NOTE:
1428 F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just one
1429 half of frequencies array is needed - elements from 0 to floor(N/2). F[0]
1430 is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd, then
1431 F[floor(N/2)] has no special properties.
1432
1433 Relying on properties noted above, FFTR1DInv subroutine uses only elements
1434 from 0th to floor(N/2)-th. It ignores imaginary part of F[0], and in case
1435 N is even it ignores imaginary part of F[floor(N/2)] too.
1436
1437 When you call this function using full arguments list - "FFTR1DInv(F,N,A)"
1438 - you can pass either either frequencies array with N elements or reduced
1439 array with roughly N/2 elements - subroutine will successfully transform
1440 both.
1441
1442 If you call this function using reduced arguments list - "FFTR1DInv(F,A)"
1443 - you must pass FULL array with N elements (although higher N/2 are still
1444 not used) because array size is used to automatically determine FFT length
1445
1446
1447 -- ALGLIB --
1448 Copyright 01.06.2009 by Bochkanov Sergey
1449 *************************************************************************/
fftr1dinv(ae_vector * f,ae_int_t n,ae_vector * a,ae_state * _state)1450 void fftr1dinv(/* Complex */ ae_vector* f,
1451 ae_int_t n,
1452 /* Real */ ae_vector* a,
1453 ae_state *_state)
1454 {
1455 ae_frame _frame_block;
1456 ae_int_t i;
1457 ae_vector h;
1458 ae_vector fh;
1459
1460 ae_frame_make(_state, &_frame_block);
1461 memset(&h, 0, sizeof(h));
1462 memset(&fh, 0, sizeof(fh));
1463 ae_vector_clear(a);
1464 ae_vector_init(&h, 0, DT_REAL, _state, ae_true);
1465 ae_vector_init(&fh, 0, DT_COMPLEX, _state, ae_true);
1466
1467 ae_assert(n>0, "FFTR1DInv: incorrect N!", _state);
1468 ae_assert(f->cnt>=ae_ifloor((double)n/(double)2, _state)+1, "FFTR1DInv: Length(F)<Floor(N/2)+1!", _state);
1469 ae_assert(ae_isfinite(f->ptr.p_complex[0].x, _state), "FFTR1DInv: F contains infinite or NAN values!", _state);
1470 for(i=1; i<=ae_ifloor((double)n/(double)2, _state)-1; i++)
1471 {
1472 ae_assert(ae_isfinite(f->ptr.p_complex[i].x, _state)&&ae_isfinite(f->ptr.p_complex[i].y, _state), "FFTR1DInv: F contains infinite or NAN values!", _state);
1473 }
1474 ae_assert(ae_isfinite(f->ptr.p_complex[ae_ifloor((double)n/(double)2, _state)].x, _state), "FFTR1DInv: F contains infinite or NAN values!", _state);
1475 if( n%2!=0 )
1476 {
1477 ae_assert(ae_isfinite(f->ptr.p_complex[ae_ifloor((double)n/(double)2, _state)].y, _state), "FFTR1DInv: F contains infinite or NAN values!", _state);
1478 }
1479
1480 /*
1481 * Special case: N=1, FFT is just identity transform.
1482 * After this block we assume that N is strictly greater than 1.
1483 */
1484 if( n==1 )
1485 {
1486 ae_vector_set_length(a, 1, _state);
1487 a->ptr.p_double[0] = f->ptr.p_complex[0].x;
1488 ae_frame_leave(_state);
1489 return;
1490 }
1491
1492 /*
1493 * inverse real FFT is reduced to the inverse real FHT,
1494 * which is reduced to the forward real FHT,
1495 * which is reduced to the forward real FFT.
1496 *
1497 * Don't worry, it is really compact and efficient reduction :)
1498 */
1499 ae_vector_set_length(&h, n, _state);
1500 ae_vector_set_length(a, n, _state);
1501 h.ptr.p_double[0] = f->ptr.p_complex[0].x;
1502 for(i=1; i<=ae_ifloor((double)n/(double)2, _state)-1; i++)
1503 {
1504 h.ptr.p_double[i] = f->ptr.p_complex[i].x-f->ptr.p_complex[i].y;
1505 h.ptr.p_double[n-i] = f->ptr.p_complex[i].x+f->ptr.p_complex[i].y;
1506 }
1507 if( n%2==0 )
1508 {
1509 h.ptr.p_double[ae_ifloor((double)n/(double)2, _state)] = f->ptr.p_complex[ae_ifloor((double)n/(double)2, _state)].x;
1510 }
1511 else
1512 {
1513 h.ptr.p_double[ae_ifloor((double)n/(double)2, _state)] = f->ptr.p_complex[ae_ifloor((double)n/(double)2, _state)].x-f->ptr.p_complex[ae_ifloor((double)n/(double)2, _state)].y;
1514 h.ptr.p_double[ae_ifloor((double)n/(double)2, _state)+1] = f->ptr.p_complex[ae_ifloor((double)n/(double)2, _state)].x+f->ptr.p_complex[ae_ifloor((double)n/(double)2, _state)].y;
1515 }
1516 fftr1d(&h, n, &fh, _state);
1517 for(i=0; i<=n-1; i++)
1518 {
1519 a->ptr.p_double[i] = (fh.ptr.p_complex[i].x-fh.ptr.p_complex[i].y)/n;
1520 }
1521 ae_frame_leave(_state);
1522 }
1523
1524
1525 /*************************************************************************
1526 Internal subroutine. Never call it directly!
1527
1528
1529 -- ALGLIB --
1530 Copyright 01.06.2009 by Bochkanov Sergey
1531 *************************************************************************/
fftr1dinternaleven(ae_vector * a,ae_int_t n,ae_vector * buf,fasttransformplan * plan,ae_state * _state)1532 void fftr1dinternaleven(/* Real */ ae_vector* a,
1533 ae_int_t n,
1534 /* Real */ ae_vector* buf,
1535 fasttransformplan* plan,
1536 ae_state *_state)
1537 {
1538 double x;
1539 double y;
1540 ae_int_t i;
1541 ae_int_t n2;
1542 ae_int_t idx;
1543 ae_complex hn;
1544 ae_complex hmnc;
1545 ae_complex v;
1546
1547
1548 ae_assert(n>0&&n%2==0, "FFTR1DEvenInplace: incorrect N!", _state);
1549
1550 /*
1551 * Special cases:
1552 * * N=2
1553 *
1554 * After this block we assume that N is strictly greater than 2
1555 */
1556 if( n==2 )
1557 {
1558 x = a->ptr.p_double[0]+a->ptr.p_double[1];
1559 y = a->ptr.p_double[0]-a->ptr.p_double[1];
1560 a->ptr.p_double[0] = x;
1561 a->ptr.p_double[1] = y;
1562 return;
1563 }
1564
1565 /*
1566 * even-size real FFT, use reduction to the complex task
1567 */
1568 n2 = n/2;
1569 ae_v_move(&buf->ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,n-1));
1570 ftapplyplan(plan, buf, 0, 1, _state);
1571 a->ptr.p_double[0] = buf->ptr.p_double[0]+buf->ptr.p_double[1];
1572 for(i=1; i<=n2-1; i++)
1573 {
1574 idx = 2*(i%n2);
1575 hn.x = buf->ptr.p_double[idx+0];
1576 hn.y = buf->ptr.p_double[idx+1];
1577 idx = 2*(n2-i);
1578 hmnc.x = buf->ptr.p_double[idx+0];
1579 hmnc.y = -buf->ptr.p_double[idx+1];
1580 v.x = -ae_sin(-2*ae_pi*i/n, _state);
1581 v.y = ae_cos(-2*ae_pi*i/n, _state);
1582 v = ae_c_sub(ae_c_add(hn,hmnc),ae_c_mul(v,ae_c_sub(hn,hmnc)));
1583 a->ptr.p_double[2*i+0] = 0.5*v.x;
1584 a->ptr.p_double[2*i+1] = 0.5*v.y;
1585 }
1586 a->ptr.p_double[1] = buf->ptr.p_double[0]-buf->ptr.p_double[1];
1587 }
1588
1589
1590 /*************************************************************************
1591 Internal subroutine. Never call it directly!
1592
1593
1594 -- ALGLIB --
1595 Copyright 01.06.2009 by Bochkanov Sergey
1596 *************************************************************************/
fftr1dinvinternaleven(ae_vector * a,ae_int_t n,ae_vector * buf,fasttransformplan * plan,ae_state * _state)1597 void fftr1dinvinternaleven(/* Real */ ae_vector* a,
1598 ae_int_t n,
1599 /* Real */ ae_vector* buf,
1600 fasttransformplan* plan,
1601 ae_state *_state)
1602 {
1603 double x;
1604 double y;
1605 double t;
1606 ae_int_t i;
1607 ae_int_t n2;
1608
1609
1610 ae_assert(n>0&&n%2==0, "FFTR1DInvInternalEven: incorrect N!", _state);
1611
1612 /*
1613 * Special cases:
1614 * * N=2
1615 *
1616 * After this block we assume that N is strictly greater than 2
1617 */
1618 if( n==2 )
1619 {
1620 x = 0.5*(a->ptr.p_double[0]+a->ptr.p_double[1]);
1621 y = 0.5*(a->ptr.p_double[0]-a->ptr.p_double[1]);
1622 a->ptr.p_double[0] = x;
1623 a->ptr.p_double[1] = y;
1624 return;
1625 }
1626
1627 /*
1628 * inverse real FFT is reduced to the inverse real FHT,
1629 * which is reduced to the forward real FHT,
1630 * which is reduced to the forward real FFT.
1631 *
1632 * Don't worry, it is really compact and efficient reduction :)
1633 */
1634 n2 = n/2;
1635 buf->ptr.p_double[0] = a->ptr.p_double[0];
1636 for(i=1; i<=n2-1; i++)
1637 {
1638 x = a->ptr.p_double[2*i+0];
1639 y = a->ptr.p_double[2*i+1];
1640 buf->ptr.p_double[i] = x-y;
1641 buf->ptr.p_double[n-i] = x+y;
1642 }
1643 buf->ptr.p_double[n2] = a->ptr.p_double[1];
1644 fftr1dinternaleven(buf, n, a, plan, _state);
1645 a->ptr.p_double[0] = buf->ptr.p_double[0]/n;
1646 t = (double)1/(double)n;
1647 for(i=1; i<=n2-1; i++)
1648 {
1649 x = buf->ptr.p_double[2*i+0];
1650 y = buf->ptr.p_double[2*i+1];
1651 a->ptr.p_double[i] = t*(x-y);
1652 a->ptr.p_double[n-i] = t*(x+y);
1653 }
1654 a->ptr.p_double[n2] = buf->ptr.p_double[1]/n;
1655 }
1656
1657
1658 #endif
1659 #if defined(AE_COMPILE_FHT) || !defined(AE_PARTIAL_BUILD)
1660
1661
1662 /*************************************************************************
1663 1-dimensional Fast Hartley Transform.
1664
1665 Algorithm has O(N*logN) complexity for any N (composite or prime).
1666
1667 INPUT PARAMETERS
1668 A - array[0..N-1] - real function to be transformed
1669 N - problem size
1670
1671 OUTPUT PARAMETERS
1672 A - FHT of a input array, array[0..N-1],
1673 A_out[k] = sum(A_in[j]*(cos(2*pi*j*k/N)+sin(2*pi*j*k/N)), j=0..N-1)
1674
1675
1676 -- ALGLIB --
1677 Copyright 04.06.2009 by Bochkanov Sergey
1678 *************************************************************************/
fhtr1d(ae_vector * a,ae_int_t n,ae_state * _state)1679 void fhtr1d(/* Real */ ae_vector* a, ae_int_t n, ae_state *_state)
1680 {
1681 ae_frame _frame_block;
1682 ae_int_t i;
1683 ae_vector fa;
1684
1685 ae_frame_make(_state, &_frame_block);
1686 memset(&fa, 0, sizeof(fa));
1687 ae_vector_init(&fa, 0, DT_COMPLEX, _state, ae_true);
1688
1689 ae_assert(n>0, "FHTR1D: incorrect N!", _state);
1690
1691 /*
1692 * Special case: N=1, FHT is just identity transform.
1693 * After this block we assume that N is strictly greater than 1.
1694 */
1695 if( n==1 )
1696 {
1697 ae_frame_leave(_state);
1698 return;
1699 }
1700
1701 /*
1702 * Reduce FHt to real FFT
1703 */
1704 fftr1d(a, n, &fa, _state);
1705 for(i=0; i<=n-1; i++)
1706 {
1707 a->ptr.p_double[i] = fa.ptr.p_complex[i].x-fa.ptr.p_complex[i].y;
1708 }
1709 ae_frame_leave(_state);
1710 }
1711
1712
1713 /*************************************************************************
1714 1-dimensional inverse FHT.
1715
1716 Algorithm has O(N*logN) complexity for any N (composite or prime).
1717
1718 INPUT PARAMETERS
1719 A - array[0..N-1] - complex array to be transformed
1720 N - problem size
1721
1722 OUTPUT PARAMETERS
1723 A - inverse FHT of a input array, array[0..N-1]
1724
1725
1726 -- ALGLIB --
1727 Copyright 29.05.2009 by Bochkanov Sergey
1728 *************************************************************************/
fhtr1dinv(ae_vector * a,ae_int_t n,ae_state * _state)1729 void fhtr1dinv(/* Real */ ae_vector* a, ae_int_t n, ae_state *_state)
1730 {
1731 ae_int_t i;
1732
1733
1734 ae_assert(n>0, "FHTR1DInv: incorrect N!", _state);
1735
1736 /*
1737 * Special case: N=1, iFHT is just identity transform.
1738 * After this block we assume that N is strictly greater than 1.
1739 */
1740 if( n==1 )
1741 {
1742 return;
1743 }
1744
1745 /*
1746 * Inverse FHT can be expressed in terms of the FHT as
1747 *
1748 * invfht(x) = fht(x)/N
1749 */
1750 fhtr1d(a, n, _state);
1751 for(i=0; i<=n-1; i++)
1752 {
1753 a->ptr.p_double[i] = a->ptr.p_double[i]/n;
1754 }
1755 }
1756
1757
1758 #endif
1759 #if defined(AE_COMPILE_CONV) || !defined(AE_PARTIAL_BUILD)
1760
1761
1762 /*************************************************************************
1763 1-dimensional complex convolution.
1764
1765 For given A/B returns conv(A,B) (non-circular). Subroutine can automatically
1766 choose between three implementations: straightforward O(M*N) formula for
1767 very small N (or M), overlap-add algorithm for cases where max(M,N) is
1768 significantly larger than min(M,N), but O(M*N) algorithm is too slow, and
1769 general FFT-based formula for cases where two previois algorithms are too
1770 slow.
1771
1772 Algorithm has max(M,N)*log(max(M,N)) complexity for any M/N.
1773
1774 INPUT PARAMETERS
1775 A - array[0..M-1] - complex function to be transformed
1776 M - problem size
1777 B - array[0..N-1] - complex function to be transformed
1778 N - problem size
1779
1780 OUTPUT PARAMETERS
1781 R - convolution: A*B. array[0..N+M-2].
1782
1783 NOTE:
1784 It is assumed that A is zero at T<0, B is zero too. If one or both
1785 functions have non-zero values at negative T's, you can still use this
1786 subroutine - just shift its result correspondingly.
1787
1788 -- ALGLIB --
1789 Copyright 21.07.2009 by Bochkanov Sergey
1790 *************************************************************************/
convc1d(ae_vector * a,ae_int_t m,ae_vector * b,ae_int_t n,ae_vector * r,ae_state * _state)1791 void convc1d(/* Complex */ ae_vector* a,
1792 ae_int_t m,
1793 /* Complex */ ae_vector* b,
1794 ae_int_t n,
1795 /* Complex */ ae_vector* r,
1796 ae_state *_state)
1797 {
1798
1799 ae_vector_clear(r);
1800
1801 ae_assert(n>0&&m>0, "ConvC1D: incorrect N or M!", _state);
1802
1803 /*
1804 * normalize task: make M>=N,
1805 * so A will be longer that B.
1806 */
1807 if( m<n )
1808 {
1809 convc1d(b, n, a, m, r, _state);
1810 return;
1811 }
1812 convc1dx(a, m, b, n, ae_false, -1, 0, r, _state);
1813 }
1814
1815
1816 /*************************************************************************
1817 1-dimensional complex non-circular deconvolution (inverse of ConvC1D()).
1818
1819 Algorithm has M*log(M)) complexity for any M (composite or prime).
1820
1821 INPUT PARAMETERS
1822 A - array[0..M-1] - convolved signal, A = conv(R, B)
1823 M - convolved signal length
1824 B - array[0..N-1] - response
1825 N - response length, N<=M
1826
1827 OUTPUT PARAMETERS
1828 R - deconvolved signal. array[0..M-N].
1829
1830 NOTE:
1831 deconvolution is unstable process and may result in division by zero
1832 (if your response function is degenerate, i.e. has zero Fourier coefficient).
1833
1834 NOTE:
1835 It is assumed that A is zero at T<0, B is zero too. If one or both
1836 functions have non-zero values at negative T's, you can still use this
1837 subroutine - just shift its result correspondingly.
1838
1839 -- ALGLIB --
1840 Copyright 21.07.2009 by Bochkanov Sergey
1841 *************************************************************************/
convc1dinv(ae_vector * a,ae_int_t m,ae_vector * b,ae_int_t n,ae_vector * r,ae_state * _state)1842 void convc1dinv(/* Complex */ ae_vector* a,
1843 ae_int_t m,
1844 /* Complex */ ae_vector* b,
1845 ae_int_t n,
1846 /* Complex */ ae_vector* r,
1847 ae_state *_state)
1848 {
1849 ae_frame _frame_block;
1850 ae_int_t i;
1851 ae_int_t p;
1852 ae_vector buf;
1853 ae_vector buf2;
1854 fasttransformplan plan;
1855 ae_complex c1;
1856 ae_complex c2;
1857 ae_complex c3;
1858 double t;
1859
1860 ae_frame_make(_state, &_frame_block);
1861 memset(&buf, 0, sizeof(buf));
1862 memset(&buf2, 0, sizeof(buf2));
1863 memset(&plan, 0, sizeof(plan));
1864 ae_vector_clear(r);
1865 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
1866 ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
1867 _fasttransformplan_init(&plan, _state, ae_true);
1868
1869 ae_assert((n>0&&m>0)&&n<=m, "ConvC1DInv: incorrect N or M!", _state);
1870 p = ftbasefindsmooth(m, _state);
1871 ftcomplexfftplan(p, 1, &plan, _state);
1872 ae_vector_set_length(&buf, 2*p, _state);
1873 for(i=0; i<=m-1; i++)
1874 {
1875 buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
1876 buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
1877 }
1878 for(i=m; i<=p-1; i++)
1879 {
1880 buf.ptr.p_double[2*i+0] = (double)(0);
1881 buf.ptr.p_double[2*i+1] = (double)(0);
1882 }
1883 ae_vector_set_length(&buf2, 2*p, _state);
1884 for(i=0; i<=n-1; i++)
1885 {
1886 buf2.ptr.p_double[2*i+0] = b->ptr.p_complex[i].x;
1887 buf2.ptr.p_double[2*i+1] = b->ptr.p_complex[i].y;
1888 }
1889 for(i=n; i<=p-1; i++)
1890 {
1891 buf2.ptr.p_double[2*i+0] = (double)(0);
1892 buf2.ptr.p_double[2*i+1] = (double)(0);
1893 }
1894 ftapplyplan(&plan, &buf, 0, 1, _state);
1895 ftapplyplan(&plan, &buf2, 0, 1, _state);
1896 for(i=0; i<=p-1; i++)
1897 {
1898 c1.x = buf.ptr.p_double[2*i+0];
1899 c1.y = buf.ptr.p_double[2*i+1];
1900 c2.x = buf2.ptr.p_double[2*i+0];
1901 c2.y = buf2.ptr.p_double[2*i+1];
1902 c3 = ae_c_div(c1,c2);
1903 buf.ptr.p_double[2*i+0] = c3.x;
1904 buf.ptr.p_double[2*i+1] = -c3.y;
1905 }
1906 ftapplyplan(&plan, &buf, 0, 1, _state);
1907 t = (double)1/(double)p;
1908 ae_vector_set_length(r, m-n+1, _state);
1909 for(i=0; i<=m-n; i++)
1910 {
1911 r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
1912 r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
1913 }
1914 ae_frame_leave(_state);
1915 }
1916
1917
1918 /*************************************************************************
1919 1-dimensional circular complex convolution.
1920
1921 For given S/R returns conv(S,R) (circular). Algorithm has linearithmic
1922 complexity for any M/N.
1923
1924 IMPORTANT: normal convolution is commutative, i.e. it is symmetric -
1925 conv(A,B)=conv(B,A). Cyclic convolution IS NOT. One function - S - is a
1926 signal, periodic function, and another - R - is a response, non-periodic
1927 function with limited length.
1928
1929 INPUT PARAMETERS
1930 S - array[0..M-1] - complex periodic signal
1931 M - problem size
1932 B - array[0..N-1] - complex non-periodic response
1933 N - problem size
1934
1935 OUTPUT PARAMETERS
1936 R - convolution: A*B. array[0..M-1].
1937
1938 NOTE:
1939 It is assumed that B is zero at T<0. If it has non-zero values at
1940 negative T's, you can still use this subroutine - just shift its result
1941 correspondingly.
1942
1943 -- ALGLIB --
1944 Copyright 21.07.2009 by Bochkanov Sergey
1945 *************************************************************************/
convc1dcircular(ae_vector * s,ae_int_t m,ae_vector * r,ae_int_t n,ae_vector * c,ae_state * _state)1946 void convc1dcircular(/* Complex */ ae_vector* s,
1947 ae_int_t m,
1948 /* Complex */ ae_vector* r,
1949 ae_int_t n,
1950 /* Complex */ ae_vector* c,
1951 ae_state *_state)
1952 {
1953 ae_frame _frame_block;
1954 ae_vector buf;
1955 ae_int_t i1;
1956 ae_int_t i2;
1957 ae_int_t j2;
1958
1959 ae_frame_make(_state, &_frame_block);
1960 memset(&buf, 0, sizeof(buf));
1961 ae_vector_clear(c);
1962 ae_vector_init(&buf, 0, DT_COMPLEX, _state, ae_true);
1963
1964 ae_assert(n>0&&m>0, "ConvC1DCircular: incorrect N or M!", _state);
1965
1966 /*
1967 * normalize task: make M>=N,
1968 * so A will be longer (at least - not shorter) that B.
1969 */
1970 if( m<n )
1971 {
1972 ae_vector_set_length(&buf, m, _state);
1973 for(i1=0; i1<=m-1; i1++)
1974 {
1975 buf.ptr.p_complex[i1] = ae_complex_from_i(0);
1976 }
1977 i1 = 0;
1978 while(i1<n)
1979 {
1980 i2 = ae_minint(i1+m-1, n-1, _state);
1981 j2 = i2-i1;
1982 ae_v_cadd(&buf.ptr.p_complex[0], 1, &r->ptr.p_complex[i1], 1, "N", ae_v_len(0,j2));
1983 i1 = i1+m;
1984 }
1985 convc1dcircular(s, m, &buf, m, c, _state);
1986 ae_frame_leave(_state);
1987 return;
1988 }
1989 convc1dx(s, m, r, n, ae_true, -1, 0, c, _state);
1990 ae_frame_leave(_state);
1991 }
1992
1993
1994 /*************************************************************************
1995 1-dimensional circular complex deconvolution (inverse of ConvC1DCircular()).
1996
1997 Algorithm has M*log(M)) complexity for any M (composite or prime).
1998
1999 INPUT PARAMETERS
2000 A - array[0..M-1] - convolved periodic signal, A = conv(R, B)
2001 M - convolved signal length
2002 B - array[0..N-1] - non-periodic response
2003 N - response length
2004
2005 OUTPUT PARAMETERS
2006 R - deconvolved signal. array[0..M-1].
2007
2008 NOTE:
2009 deconvolution is unstable process and may result in division by zero
2010 (if your response function is degenerate, i.e. has zero Fourier coefficient).
2011
2012 NOTE:
2013 It is assumed that B is zero at T<0. If it has non-zero values at
2014 negative T's, you can still use this subroutine - just shift its result
2015 correspondingly.
2016
2017 -- ALGLIB --
2018 Copyright 21.07.2009 by Bochkanov Sergey
2019 *************************************************************************/
convc1dcircularinv(ae_vector * a,ae_int_t m,ae_vector * b,ae_int_t n,ae_vector * r,ae_state * _state)2020 void convc1dcircularinv(/* Complex */ ae_vector* a,
2021 ae_int_t m,
2022 /* Complex */ ae_vector* b,
2023 ae_int_t n,
2024 /* Complex */ ae_vector* r,
2025 ae_state *_state)
2026 {
2027 ae_frame _frame_block;
2028 ae_int_t i;
2029 ae_int_t i1;
2030 ae_int_t i2;
2031 ae_int_t j2;
2032 ae_vector buf;
2033 ae_vector buf2;
2034 ae_vector cbuf;
2035 fasttransformplan plan;
2036 ae_complex c1;
2037 ae_complex c2;
2038 ae_complex c3;
2039 double t;
2040
2041 ae_frame_make(_state, &_frame_block);
2042 memset(&buf, 0, sizeof(buf));
2043 memset(&buf2, 0, sizeof(buf2));
2044 memset(&cbuf, 0, sizeof(cbuf));
2045 memset(&plan, 0, sizeof(plan));
2046 ae_vector_clear(r);
2047 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2048 ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
2049 ae_vector_init(&cbuf, 0, DT_COMPLEX, _state, ae_true);
2050 _fasttransformplan_init(&plan, _state, ae_true);
2051
2052 ae_assert(n>0&&m>0, "ConvC1DCircularInv: incorrect N or M!", _state);
2053
2054 /*
2055 * normalize task: make M>=N,
2056 * so A will be longer (at least - not shorter) that B.
2057 */
2058 if( m<n )
2059 {
2060 ae_vector_set_length(&cbuf, m, _state);
2061 for(i=0; i<=m-1; i++)
2062 {
2063 cbuf.ptr.p_complex[i] = ae_complex_from_i(0);
2064 }
2065 i1 = 0;
2066 while(i1<n)
2067 {
2068 i2 = ae_minint(i1+m-1, n-1, _state);
2069 j2 = i2-i1;
2070 ae_v_cadd(&cbuf.ptr.p_complex[0], 1, &b->ptr.p_complex[i1], 1, "N", ae_v_len(0,j2));
2071 i1 = i1+m;
2072 }
2073 convc1dcircularinv(a, m, &cbuf, m, r, _state);
2074 ae_frame_leave(_state);
2075 return;
2076 }
2077
2078 /*
2079 * Task is normalized
2080 */
2081 ftcomplexfftplan(m, 1, &plan, _state);
2082 ae_vector_set_length(&buf, 2*m, _state);
2083 for(i=0; i<=m-1; i++)
2084 {
2085 buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
2086 buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
2087 }
2088 ae_vector_set_length(&buf2, 2*m, _state);
2089 for(i=0; i<=n-1; i++)
2090 {
2091 buf2.ptr.p_double[2*i+0] = b->ptr.p_complex[i].x;
2092 buf2.ptr.p_double[2*i+1] = b->ptr.p_complex[i].y;
2093 }
2094 for(i=n; i<=m-1; i++)
2095 {
2096 buf2.ptr.p_double[2*i+0] = (double)(0);
2097 buf2.ptr.p_double[2*i+1] = (double)(0);
2098 }
2099 ftapplyplan(&plan, &buf, 0, 1, _state);
2100 ftapplyplan(&plan, &buf2, 0, 1, _state);
2101 for(i=0; i<=m-1; i++)
2102 {
2103 c1.x = buf.ptr.p_double[2*i+0];
2104 c1.y = buf.ptr.p_double[2*i+1];
2105 c2.x = buf2.ptr.p_double[2*i+0];
2106 c2.y = buf2.ptr.p_double[2*i+1];
2107 c3 = ae_c_div(c1,c2);
2108 buf.ptr.p_double[2*i+0] = c3.x;
2109 buf.ptr.p_double[2*i+1] = -c3.y;
2110 }
2111 ftapplyplan(&plan, &buf, 0, 1, _state);
2112 t = (double)1/(double)m;
2113 ae_vector_set_length(r, m, _state);
2114 for(i=0; i<=m-1; i++)
2115 {
2116 r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
2117 r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
2118 }
2119 ae_frame_leave(_state);
2120 }
2121
2122
2123 /*************************************************************************
2124 1-dimensional real convolution.
2125
2126 Analogous to ConvC1D(), see ConvC1D() comments for more details.
2127
2128 INPUT PARAMETERS
2129 A - array[0..M-1] - real function to be transformed
2130 M - problem size
2131 B - array[0..N-1] - real function to be transformed
2132 N - problem size
2133
2134 OUTPUT PARAMETERS
2135 R - convolution: A*B. array[0..N+M-2].
2136
2137 NOTE:
2138 It is assumed that A is zero at T<0, B is zero too. If one or both
2139 functions have non-zero values at negative T's, you can still use this
2140 subroutine - just shift its result correspondingly.
2141
2142 -- ALGLIB --
2143 Copyright 21.07.2009 by Bochkanov Sergey
2144 *************************************************************************/
convr1d(ae_vector * a,ae_int_t m,ae_vector * b,ae_int_t n,ae_vector * r,ae_state * _state)2145 void convr1d(/* Real */ ae_vector* a,
2146 ae_int_t m,
2147 /* Real */ ae_vector* b,
2148 ae_int_t n,
2149 /* Real */ ae_vector* r,
2150 ae_state *_state)
2151 {
2152
2153 ae_vector_clear(r);
2154
2155 ae_assert(n>0&&m>0, "ConvR1D: incorrect N or M!", _state);
2156
2157 /*
2158 * normalize task: make M>=N,
2159 * so A will be longer that B.
2160 */
2161 if( m<n )
2162 {
2163 convr1d(b, n, a, m, r, _state);
2164 return;
2165 }
2166 convr1dx(a, m, b, n, ae_false, -1, 0, r, _state);
2167 }
2168
2169
2170 /*************************************************************************
2171 1-dimensional real deconvolution (inverse of ConvC1D()).
2172
2173 Algorithm has M*log(M)) complexity for any M (composite or prime).
2174
2175 INPUT PARAMETERS
2176 A - array[0..M-1] - convolved signal, A = conv(R, B)
2177 M - convolved signal length
2178 B - array[0..N-1] - response
2179 N - response length, N<=M
2180
2181 OUTPUT PARAMETERS
2182 R - deconvolved signal. array[0..M-N].
2183
2184 NOTE:
2185 deconvolution is unstable process and may result in division by zero
2186 (if your response function is degenerate, i.e. has zero Fourier coefficient).
2187
2188 NOTE:
2189 It is assumed that A is zero at T<0, B is zero too. If one or both
2190 functions have non-zero values at negative T's, you can still use this
2191 subroutine - just shift its result correspondingly.
2192
2193 -- ALGLIB --
2194 Copyright 21.07.2009 by Bochkanov Sergey
2195 *************************************************************************/
convr1dinv(ae_vector * a,ae_int_t m,ae_vector * b,ae_int_t n,ae_vector * r,ae_state * _state)2196 void convr1dinv(/* Real */ ae_vector* a,
2197 ae_int_t m,
2198 /* Real */ ae_vector* b,
2199 ae_int_t n,
2200 /* Real */ ae_vector* r,
2201 ae_state *_state)
2202 {
2203 ae_frame _frame_block;
2204 ae_int_t i;
2205 ae_int_t p;
2206 ae_vector buf;
2207 ae_vector buf2;
2208 ae_vector buf3;
2209 fasttransformplan plan;
2210 ae_complex c1;
2211 ae_complex c2;
2212 ae_complex c3;
2213
2214 ae_frame_make(_state, &_frame_block);
2215 memset(&buf, 0, sizeof(buf));
2216 memset(&buf2, 0, sizeof(buf2));
2217 memset(&buf3, 0, sizeof(buf3));
2218 memset(&plan, 0, sizeof(plan));
2219 ae_vector_clear(r);
2220 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2221 ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
2222 ae_vector_init(&buf3, 0, DT_REAL, _state, ae_true);
2223 _fasttransformplan_init(&plan, _state, ae_true);
2224
2225 ae_assert((n>0&&m>0)&&n<=m, "ConvR1DInv: incorrect N or M!", _state);
2226 p = ftbasefindsmootheven(m, _state);
2227 ae_vector_set_length(&buf, p, _state);
2228 ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1));
2229 for(i=m; i<=p-1; i++)
2230 {
2231 buf.ptr.p_double[i] = (double)(0);
2232 }
2233 ae_vector_set_length(&buf2, p, _state);
2234 ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2235 for(i=n; i<=p-1; i++)
2236 {
2237 buf2.ptr.p_double[i] = (double)(0);
2238 }
2239 ae_vector_set_length(&buf3, p, _state);
2240 ftcomplexfftplan(p/2, 1, &plan, _state);
2241 fftr1dinternaleven(&buf, p, &buf3, &plan, _state);
2242 fftr1dinternaleven(&buf2, p, &buf3, &plan, _state);
2243 buf.ptr.p_double[0] = buf.ptr.p_double[0]/buf2.ptr.p_double[0];
2244 buf.ptr.p_double[1] = buf.ptr.p_double[1]/buf2.ptr.p_double[1];
2245 for(i=1; i<=p/2-1; i++)
2246 {
2247 c1.x = buf.ptr.p_double[2*i+0];
2248 c1.y = buf.ptr.p_double[2*i+1];
2249 c2.x = buf2.ptr.p_double[2*i+0];
2250 c2.y = buf2.ptr.p_double[2*i+1];
2251 c3 = ae_c_div(c1,c2);
2252 buf.ptr.p_double[2*i+0] = c3.x;
2253 buf.ptr.p_double[2*i+1] = c3.y;
2254 }
2255 fftr1dinvinternaleven(&buf, p, &buf3, &plan, _state);
2256 ae_vector_set_length(r, m-n+1, _state);
2257 ae_v_move(&r->ptr.p_double[0], 1, &buf.ptr.p_double[0], 1, ae_v_len(0,m-n));
2258 ae_frame_leave(_state);
2259 }
2260
2261
2262 /*************************************************************************
2263 1-dimensional circular real convolution.
2264
2265 Analogous to ConvC1DCircular(), see ConvC1DCircular() comments for more details.
2266
2267 INPUT PARAMETERS
2268 S - array[0..M-1] - real signal
2269 M - problem size
2270 B - array[0..N-1] - real response
2271 N - problem size
2272
2273 OUTPUT PARAMETERS
2274 R - convolution: A*B. array[0..M-1].
2275
2276 NOTE:
2277 It is assumed that B is zero at T<0. If it has non-zero values at
2278 negative T's, you can still use this subroutine - just shift its result
2279 correspondingly.
2280
2281 -- ALGLIB --
2282 Copyright 21.07.2009 by Bochkanov Sergey
2283 *************************************************************************/
convr1dcircular(ae_vector * s,ae_int_t m,ae_vector * r,ae_int_t n,ae_vector * c,ae_state * _state)2284 void convr1dcircular(/* Real */ ae_vector* s,
2285 ae_int_t m,
2286 /* Real */ ae_vector* r,
2287 ae_int_t n,
2288 /* Real */ ae_vector* c,
2289 ae_state *_state)
2290 {
2291 ae_frame _frame_block;
2292 ae_vector buf;
2293 ae_int_t i1;
2294 ae_int_t i2;
2295 ae_int_t j2;
2296
2297 ae_frame_make(_state, &_frame_block);
2298 memset(&buf, 0, sizeof(buf));
2299 ae_vector_clear(c);
2300 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2301
2302 ae_assert(n>0&&m>0, "ConvC1DCircular: incorrect N or M!", _state);
2303
2304 /*
2305 * normalize task: make M>=N,
2306 * so A will be longer (at least - not shorter) that B.
2307 */
2308 if( m<n )
2309 {
2310 ae_vector_set_length(&buf, m, _state);
2311 for(i1=0; i1<=m-1; i1++)
2312 {
2313 buf.ptr.p_double[i1] = (double)(0);
2314 }
2315 i1 = 0;
2316 while(i1<n)
2317 {
2318 i2 = ae_minint(i1+m-1, n-1, _state);
2319 j2 = i2-i1;
2320 ae_v_add(&buf.ptr.p_double[0], 1, &r->ptr.p_double[i1], 1, ae_v_len(0,j2));
2321 i1 = i1+m;
2322 }
2323 convr1dcircular(s, m, &buf, m, c, _state);
2324 ae_frame_leave(_state);
2325 return;
2326 }
2327
2328 /*
2329 * reduce to usual convolution
2330 */
2331 convr1dx(s, m, r, n, ae_true, -1, 0, c, _state);
2332 ae_frame_leave(_state);
2333 }
2334
2335
2336 /*************************************************************************
2337 1-dimensional complex deconvolution (inverse of ConvC1D()).
2338
2339 Algorithm has M*log(M)) complexity for any M (composite or prime).
2340
2341 INPUT PARAMETERS
2342 A - array[0..M-1] - convolved signal, A = conv(R, B)
2343 M - convolved signal length
2344 B - array[0..N-1] - response
2345 N - response length
2346
2347 OUTPUT PARAMETERS
2348 R - deconvolved signal. array[0..M-N].
2349
2350 NOTE:
2351 deconvolution is unstable process and may result in division by zero
2352 (if your response function is degenerate, i.e. has zero Fourier coefficient).
2353
2354 NOTE:
2355 It is assumed that B is zero at T<0. If it has non-zero values at
2356 negative T's, you can still use this subroutine - just shift its result
2357 correspondingly.
2358
2359 -- ALGLIB --
2360 Copyright 21.07.2009 by Bochkanov Sergey
2361 *************************************************************************/
convr1dcircularinv(ae_vector * a,ae_int_t m,ae_vector * b,ae_int_t n,ae_vector * r,ae_state * _state)2362 void convr1dcircularinv(/* Real */ ae_vector* a,
2363 ae_int_t m,
2364 /* Real */ ae_vector* b,
2365 ae_int_t n,
2366 /* Real */ ae_vector* r,
2367 ae_state *_state)
2368 {
2369 ae_frame _frame_block;
2370 ae_int_t i;
2371 ae_int_t i1;
2372 ae_int_t i2;
2373 ae_int_t j2;
2374 ae_vector buf;
2375 ae_vector buf2;
2376 ae_vector buf3;
2377 ae_vector cbuf;
2378 ae_vector cbuf2;
2379 fasttransformplan plan;
2380 ae_complex c1;
2381 ae_complex c2;
2382 ae_complex c3;
2383
2384 ae_frame_make(_state, &_frame_block);
2385 memset(&buf, 0, sizeof(buf));
2386 memset(&buf2, 0, sizeof(buf2));
2387 memset(&buf3, 0, sizeof(buf3));
2388 memset(&cbuf, 0, sizeof(cbuf));
2389 memset(&cbuf2, 0, sizeof(cbuf2));
2390 memset(&plan, 0, sizeof(plan));
2391 ae_vector_clear(r);
2392 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2393 ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
2394 ae_vector_init(&buf3, 0, DT_REAL, _state, ae_true);
2395 ae_vector_init(&cbuf, 0, DT_COMPLEX, _state, ae_true);
2396 ae_vector_init(&cbuf2, 0, DT_COMPLEX, _state, ae_true);
2397 _fasttransformplan_init(&plan, _state, ae_true);
2398
2399 ae_assert(n>0&&m>0, "ConvR1DCircularInv: incorrect N or M!", _state);
2400
2401 /*
2402 * normalize task: make M>=N,
2403 * so A will be longer (at least - not shorter) that B.
2404 */
2405 if( m<n )
2406 {
2407 ae_vector_set_length(&buf, m, _state);
2408 for(i=0; i<=m-1; i++)
2409 {
2410 buf.ptr.p_double[i] = (double)(0);
2411 }
2412 i1 = 0;
2413 while(i1<n)
2414 {
2415 i2 = ae_minint(i1+m-1, n-1, _state);
2416 j2 = i2-i1;
2417 ae_v_add(&buf.ptr.p_double[0], 1, &b->ptr.p_double[i1], 1, ae_v_len(0,j2));
2418 i1 = i1+m;
2419 }
2420 convr1dcircularinv(a, m, &buf, m, r, _state);
2421 ae_frame_leave(_state);
2422 return;
2423 }
2424
2425 /*
2426 * Task is normalized
2427 */
2428 if( m%2==0 )
2429 {
2430
2431 /*
2432 * size is even, use fast even-size FFT
2433 */
2434 ae_vector_set_length(&buf, m, _state);
2435 ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1));
2436 ae_vector_set_length(&buf2, m, _state);
2437 ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2438 for(i=n; i<=m-1; i++)
2439 {
2440 buf2.ptr.p_double[i] = (double)(0);
2441 }
2442 ae_vector_set_length(&buf3, m, _state);
2443 ftcomplexfftplan(m/2, 1, &plan, _state);
2444 fftr1dinternaleven(&buf, m, &buf3, &plan, _state);
2445 fftr1dinternaleven(&buf2, m, &buf3, &plan, _state);
2446 buf.ptr.p_double[0] = buf.ptr.p_double[0]/buf2.ptr.p_double[0];
2447 buf.ptr.p_double[1] = buf.ptr.p_double[1]/buf2.ptr.p_double[1];
2448 for(i=1; i<=m/2-1; i++)
2449 {
2450 c1.x = buf.ptr.p_double[2*i+0];
2451 c1.y = buf.ptr.p_double[2*i+1];
2452 c2.x = buf2.ptr.p_double[2*i+0];
2453 c2.y = buf2.ptr.p_double[2*i+1];
2454 c3 = ae_c_div(c1,c2);
2455 buf.ptr.p_double[2*i+0] = c3.x;
2456 buf.ptr.p_double[2*i+1] = c3.y;
2457 }
2458 fftr1dinvinternaleven(&buf, m, &buf3, &plan, _state);
2459 ae_vector_set_length(r, m, _state);
2460 ae_v_move(&r->ptr.p_double[0], 1, &buf.ptr.p_double[0], 1, ae_v_len(0,m-1));
2461 }
2462 else
2463 {
2464
2465 /*
2466 * odd-size, use general real FFT
2467 */
2468 fftr1d(a, m, &cbuf, _state);
2469 ae_vector_set_length(&buf2, m, _state);
2470 ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2471 for(i=n; i<=m-1; i++)
2472 {
2473 buf2.ptr.p_double[i] = (double)(0);
2474 }
2475 fftr1d(&buf2, m, &cbuf2, _state);
2476 for(i=0; i<=ae_ifloor((double)m/(double)2, _state); i++)
2477 {
2478 cbuf.ptr.p_complex[i] = ae_c_div(cbuf.ptr.p_complex[i],cbuf2.ptr.p_complex[i]);
2479 }
2480 fftr1dinv(&cbuf, m, r, _state);
2481 }
2482 ae_frame_leave(_state);
2483 }
2484
2485
2486 /*************************************************************************
2487 1-dimensional complex convolution.
2488
2489 Extended subroutine which allows to choose convolution algorithm.
2490 Intended for internal use, ALGLIB users should call ConvC1D()/ConvC1DCircular().
2491
2492 INPUT PARAMETERS
2493 A - array[0..M-1] - complex function to be transformed
2494 M - problem size
2495 B - array[0..N-1] - complex function to be transformed
2496 N - problem size, N<=M
2497 Alg - algorithm type:
2498 *-2 auto-select Q for overlap-add
2499 *-1 auto-select algorithm and parameters
2500 * 0 straightforward formula for small N's
2501 * 1 general FFT-based code
2502 * 2 overlap-add with length Q
2503 Q - length for overlap-add
2504
2505 OUTPUT PARAMETERS
2506 R - convolution: A*B. array[0..N+M-1].
2507
2508 -- ALGLIB --
2509 Copyright 21.07.2009 by Bochkanov Sergey
2510 *************************************************************************/
convc1dx(ae_vector * a,ae_int_t m,ae_vector * b,ae_int_t n,ae_bool circular,ae_int_t alg,ae_int_t q,ae_vector * r,ae_state * _state)2511 void convc1dx(/* Complex */ ae_vector* a,
2512 ae_int_t m,
2513 /* Complex */ ae_vector* b,
2514 ae_int_t n,
2515 ae_bool circular,
2516 ae_int_t alg,
2517 ae_int_t q,
2518 /* Complex */ ae_vector* r,
2519 ae_state *_state)
2520 {
2521 ae_frame _frame_block;
2522 ae_int_t i;
2523 ae_int_t j;
2524 ae_int_t p;
2525 ae_int_t ptotal;
2526 ae_int_t i1;
2527 ae_int_t i2;
2528 ae_int_t j1;
2529 ae_int_t j2;
2530 ae_vector bbuf;
2531 ae_complex v;
2532 double ax;
2533 double ay;
2534 double bx;
2535 double by;
2536 double t;
2537 double tx;
2538 double ty;
2539 double flopcand;
2540 double flopbest;
2541 ae_int_t algbest;
2542 fasttransformplan plan;
2543 ae_vector buf;
2544 ae_vector buf2;
2545
2546 ae_frame_make(_state, &_frame_block);
2547 memset(&bbuf, 0, sizeof(bbuf));
2548 memset(&plan, 0, sizeof(plan));
2549 memset(&buf, 0, sizeof(buf));
2550 memset(&buf2, 0, sizeof(buf2));
2551 ae_vector_clear(r);
2552 ae_vector_init(&bbuf, 0, DT_COMPLEX, _state, ae_true);
2553 _fasttransformplan_init(&plan, _state, ae_true);
2554 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2555 ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
2556
2557 ae_assert(n>0&&m>0, "ConvC1DX: incorrect N or M!", _state);
2558 ae_assert(n<=m, "ConvC1DX: N<M assumption is false!", _state);
2559
2560 /*
2561 * Auto-select
2562 */
2563 if( alg==-1||alg==-2 )
2564 {
2565
2566 /*
2567 * Initial candidate: straightforward implementation.
2568 *
2569 * If we want to use auto-fitted overlap-add,
2570 * flop count is initialized by large real number - to force
2571 * another algorithm selection
2572 */
2573 algbest = 0;
2574 if( alg==-1 )
2575 {
2576 flopbest = (double)(2*m*n);
2577 }
2578 else
2579 {
2580 flopbest = ae_maxrealnumber;
2581 }
2582
2583 /*
2584 * Another candidate - generic FFT code
2585 */
2586 if( alg==-1 )
2587 {
2588 if( circular&&ftbaseissmooth(m, _state) )
2589 {
2590
2591 /*
2592 * special code for circular convolution of a sequence with a smooth length
2593 */
2594 flopcand = 3*ftbasegetflopestimate(m, _state)+6*m;
2595 if( ae_fp_less(flopcand,flopbest) )
2596 {
2597 algbest = 1;
2598 flopbest = flopcand;
2599 }
2600 }
2601 else
2602 {
2603
2604 /*
2605 * general cyclic/non-cyclic convolution
2606 */
2607 p = ftbasefindsmooth(m+n-1, _state);
2608 flopcand = 3*ftbasegetflopestimate(p, _state)+6*p;
2609 if( ae_fp_less(flopcand,flopbest) )
2610 {
2611 algbest = 1;
2612 flopbest = flopcand;
2613 }
2614 }
2615 }
2616
2617 /*
2618 * Another candidate - overlap-add
2619 */
2620 q = 1;
2621 ptotal = 1;
2622 while(ptotal<n)
2623 {
2624 ptotal = ptotal*2;
2625 }
2626 while(ptotal<=m+n-1)
2627 {
2628 p = ptotal-n+1;
2629 flopcand = ae_iceil((double)m/(double)p, _state)*(2*ftbasegetflopestimate(ptotal, _state)+8*ptotal);
2630 if( ae_fp_less(flopcand,flopbest) )
2631 {
2632 flopbest = flopcand;
2633 algbest = 2;
2634 q = p;
2635 }
2636 ptotal = ptotal*2;
2637 }
2638 alg = algbest;
2639 convc1dx(a, m, b, n, circular, alg, q, r, _state);
2640 ae_frame_leave(_state);
2641 return;
2642 }
2643
2644 /*
2645 * straightforward formula for
2646 * circular and non-circular convolutions.
2647 *
2648 * Very simple code, no further comments needed.
2649 */
2650 if( alg==0 )
2651 {
2652
2653 /*
2654 * Special case: N=1
2655 */
2656 if( n==1 )
2657 {
2658 ae_vector_set_length(r, m, _state);
2659 v = b->ptr.p_complex[0];
2660 ae_v_cmovec(&r->ptr.p_complex[0], 1, &a->ptr.p_complex[0], 1, "N", ae_v_len(0,m-1), v);
2661 ae_frame_leave(_state);
2662 return;
2663 }
2664
2665 /*
2666 * use straightforward formula
2667 */
2668 if( circular )
2669 {
2670
2671 /*
2672 * circular convolution
2673 */
2674 ae_vector_set_length(r, m, _state);
2675 v = b->ptr.p_complex[0];
2676 ae_v_cmovec(&r->ptr.p_complex[0], 1, &a->ptr.p_complex[0], 1, "N", ae_v_len(0,m-1), v);
2677 for(i=1; i<=n-1; i++)
2678 {
2679 v = b->ptr.p_complex[i];
2680 i1 = 0;
2681 i2 = i-1;
2682 j1 = m-i;
2683 j2 = m-1;
2684 ae_v_caddc(&r->ptr.p_complex[i1], 1, &a->ptr.p_complex[j1], 1, "N", ae_v_len(i1,i2), v);
2685 i1 = i;
2686 i2 = m-1;
2687 j1 = 0;
2688 j2 = m-i-1;
2689 ae_v_caddc(&r->ptr.p_complex[i1], 1, &a->ptr.p_complex[j1], 1, "N", ae_v_len(i1,i2), v);
2690 }
2691 }
2692 else
2693 {
2694
2695 /*
2696 * non-circular convolution
2697 */
2698 ae_vector_set_length(r, m+n-1, _state);
2699 for(i=0; i<=m+n-2; i++)
2700 {
2701 r->ptr.p_complex[i] = ae_complex_from_i(0);
2702 }
2703 for(i=0; i<=n-1; i++)
2704 {
2705 v = b->ptr.p_complex[i];
2706 ae_v_caddc(&r->ptr.p_complex[i], 1, &a->ptr.p_complex[0], 1, "N", ae_v_len(i,i+m-1), v);
2707 }
2708 }
2709 ae_frame_leave(_state);
2710 return;
2711 }
2712
2713 /*
2714 * general FFT-based code for
2715 * circular and non-circular convolutions.
2716 *
2717 * First, if convolution is circular, we test whether M is smooth or not.
2718 * If it is smooth, we just use M-length FFT to calculate convolution.
2719 * If it is not, we calculate non-circular convolution and wrap it arount.
2720 *
2721 * IF convolution is non-circular, we use zero-padding + FFT.
2722 */
2723 if( alg==1 )
2724 {
2725 if( circular&&ftbaseissmooth(m, _state) )
2726 {
2727
2728 /*
2729 * special code for circular convolution with smooth M
2730 */
2731 ftcomplexfftplan(m, 1, &plan, _state);
2732 ae_vector_set_length(&buf, 2*m, _state);
2733 for(i=0; i<=m-1; i++)
2734 {
2735 buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
2736 buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
2737 }
2738 ae_vector_set_length(&buf2, 2*m, _state);
2739 for(i=0; i<=n-1; i++)
2740 {
2741 buf2.ptr.p_double[2*i+0] = b->ptr.p_complex[i].x;
2742 buf2.ptr.p_double[2*i+1] = b->ptr.p_complex[i].y;
2743 }
2744 for(i=n; i<=m-1; i++)
2745 {
2746 buf2.ptr.p_double[2*i+0] = (double)(0);
2747 buf2.ptr.p_double[2*i+1] = (double)(0);
2748 }
2749 ftapplyplan(&plan, &buf, 0, 1, _state);
2750 ftapplyplan(&plan, &buf2, 0, 1, _state);
2751 for(i=0; i<=m-1; i++)
2752 {
2753 ax = buf.ptr.p_double[2*i+0];
2754 ay = buf.ptr.p_double[2*i+1];
2755 bx = buf2.ptr.p_double[2*i+0];
2756 by = buf2.ptr.p_double[2*i+1];
2757 tx = ax*bx-ay*by;
2758 ty = ax*by+ay*bx;
2759 buf.ptr.p_double[2*i+0] = tx;
2760 buf.ptr.p_double[2*i+1] = -ty;
2761 }
2762 ftapplyplan(&plan, &buf, 0, 1, _state);
2763 t = (double)1/(double)m;
2764 ae_vector_set_length(r, m, _state);
2765 for(i=0; i<=m-1; i++)
2766 {
2767 r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
2768 r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
2769 }
2770 }
2771 else
2772 {
2773
2774 /*
2775 * M is non-smooth, general code (circular/non-circular):
2776 * * first part is the same for circular and non-circular
2777 * convolutions. zero padding, FFTs, inverse FFTs
2778 * * second part differs:
2779 * * for non-circular convolution we just copy array
2780 * * for circular convolution we add array tail to its head
2781 */
2782 p = ftbasefindsmooth(m+n-1, _state);
2783 ftcomplexfftplan(p, 1, &plan, _state);
2784 ae_vector_set_length(&buf, 2*p, _state);
2785 for(i=0; i<=m-1; i++)
2786 {
2787 buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
2788 buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
2789 }
2790 for(i=m; i<=p-1; i++)
2791 {
2792 buf.ptr.p_double[2*i+0] = (double)(0);
2793 buf.ptr.p_double[2*i+1] = (double)(0);
2794 }
2795 ae_vector_set_length(&buf2, 2*p, _state);
2796 for(i=0; i<=n-1; i++)
2797 {
2798 buf2.ptr.p_double[2*i+0] = b->ptr.p_complex[i].x;
2799 buf2.ptr.p_double[2*i+1] = b->ptr.p_complex[i].y;
2800 }
2801 for(i=n; i<=p-1; i++)
2802 {
2803 buf2.ptr.p_double[2*i+0] = (double)(0);
2804 buf2.ptr.p_double[2*i+1] = (double)(0);
2805 }
2806 ftapplyplan(&plan, &buf, 0, 1, _state);
2807 ftapplyplan(&plan, &buf2, 0, 1, _state);
2808 for(i=0; i<=p-1; i++)
2809 {
2810 ax = buf.ptr.p_double[2*i+0];
2811 ay = buf.ptr.p_double[2*i+1];
2812 bx = buf2.ptr.p_double[2*i+0];
2813 by = buf2.ptr.p_double[2*i+1];
2814 tx = ax*bx-ay*by;
2815 ty = ax*by+ay*bx;
2816 buf.ptr.p_double[2*i+0] = tx;
2817 buf.ptr.p_double[2*i+1] = -ty;
2818 }
2819 ftapplyplan(&plan, &buf, 0, 1, _state);
2820 t = (double)1/(double)p;
2821 if( circular )
2822 {
2823
2824 /*
2825 * circular, add tail to head
2826 */
2827 ae_vector_set_length(r, m, _state);
2828 for(i=0; i<=m-1; i++)
2829 {
2830 r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
2831 r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
2832 }
2833 for(i=m; i<=m+n-2; i++)
2834 {
2835 r->ptr.p_complex[i-m].x = r->ptr.p_complex[i-m].x+t*buf.ptr.p_double[2*i+0];
2836 r->ptr.p_complex[i-m].y = r->ptr.p_complex[i-m].y-t*buf.ptr.p_double[2*i+1];
2837 }
2838 }
2839 else
2840 {
2841
2842 /*
2843 * non-circular, just copy
2844 */
2845 ae_vector_set_length(r, m+n-1, _state);
2846 for(i=0; i<=m+n-2; i++)
2847 {
2848 r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
2849 r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
2850 }
2851 }
2852 }
2853 ae_frame_leave(_state);
2854 return;
2855 }
2856
2857 /*
2858 * overlap-add method for
2859 * circular and non-circular convolutions.
2860 *
2861 * First part of code (separate FFTs of input blocks) is the same
2862 * for all types of convolution. Second part (overlapping outputs)
2863 * differs for different types of convolution. We just copy output
2864 * when convolution is non-circular. We wrap it around, if it is
2865 * circular.
2866 */
2867 if( alg==2 )
2868 {
2869 ae_vector_set_length(&buf, 2*(q+n-1), _state);
2870
2871 /*
2872 * prepare R
2873 */
2874 if( circular )
2875 {
2876 ae_vector_set_length(r, m, _state);
2877 for(i=0; i<=m-1; i++)
2878 {
2879 r->ptr.p_complex[i] = ae_complex_from_i(0);
2880 }
2881 }
2882 else
2883 {
2884 ae_vector_set_length(r, m+n-1, _state);
2885 for(i=0; i<=m+n-2; i++)
2886 {
2887 r->ptr.p_complex[i] = ae_complex_from_i(0);
2888 }
2889 }
2890
2891 /*
2892 * pre-calculated FFT(B)
2893 */
2894 ae_vector_set_length(&bbuf, q+n-1, _state);
2895 ae_v_cmove(&bbuf.ptr.p_complex[0], 1, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
2896 for(j=n; j<=q+n-2; j++)
2897 {
2898 bbuf.ptr.p_complex[j] = ae_complex_from_i(0);
2899 }
2900 fftc1d(&bbuf, q+n-1, _state);
2901
2902 /*
2903 * prepare FFT plan for chunks of A
2904 */
2905 ftcomplexfftplan(q+n-1, 1, &plan, _state);
2906
2907 /*
2908 * main overlap-add cycle
2909 */
2910 i = 0;
2911 while(i<=m-1)
2912 {
2913 p = ae_minint(q, m-i, _state);
2914 for(j=0; j<=p-1; j++)
2915 {
2916 buf.ptr.p_double[2*j+0] = a->ptr.p_complex[i+j].x;
2917 buf.ptr.p_double[2*j+1] = a->ptr.p_complex[i+j].y;
2918 }
2919 for(j=p; j<=q+n-2; j++)
2920 {
2921 buf.ptr.p_double[2*j+0] = (double)(0);
2922 buf.ptr.p_double[2*j+1] = (double)(0);
2923 }
2924 ftapplyplan(&plan, &buf, 0, 1, _state);
2925 for(j=0; j<=q+n-2; j++)
2926 {
2927 ax = buf.ptr.p_double[2*j+0];
2928 ay = buf.ptr.p_double[2*j+1];
2929 bx = bbuf.ptr.p_complex[j].x;
2930 by = bbuf.ptr.p_complex[j].y;
2931 tx = ax*bx-ay*by;
2932 ty = ax*by+ay*bx;
2933 buf.ptr.p_double[2*j+0] = tx;
2934 buf.ptr.p_double[2*j+1] = -ty;
2935 }
2936 ftapplyplan(&plan, &buf, 0, 1, _state);
2937 t = (double)1/(double)(q+n-1);
2938 if( circular )
2939 {
2940 j1 = ae_minint(i+p+n-2, m-1, _state)-i;
2941 j2 = j1+1;
2942 }
2943 else
2944 {
2945 j1 = p+n-2;
2946 j2 = j1+1;
2947 }
2948 for(j=0; j<=j1; j++)
2949 {
2950 r->ptr.p_complex[i+j].x = r->ptr.p_complex[i+j].x+buf.ptr.p_double[2*j+0]*t;
2951 r->ptr.p_complex[i+j].y = r->ptr.p_complex[i+j].y-buf.ptr.p_double[2*j+1]*t;
2952 }
2953 for(j=j2; j<=p+n-2; j++)
2954 {
2955 r->ptr.p_complex[j-j2].x = r->ptr.p_complex[j-j2].x+buf.ptr.p_double[2*j+0]*t;
2956 r->ptr.p_complex[j-j2].y = r->ptr.p_complex[j-j2].y-buf.ptr.p_double[2*j+1]*t;
2957 }
2958 i = i+p;
2959 }
2960 ae_frame_leave(_state);
2961 return;
2962 }
2963 ae_frame_leave(_state);
2964 }
2965
2966
2967 /*************************************************************************
2968 1-dimensional real convolution.
2969
2970 Extended subroutine which allows to choose convolution algorithm.
2971 Intended for internal use, ALGLIB users should call ConvR1D().
2972
2973 INPUT PARAMETERS
2974 A - array[0..M-1] - complex function to be transformed
2975 M - problem size
2976 B - array[0..N-1] - complex function to be transformed
2977 N - problem size, N<=M
2978 Alg - algorithm type:
2979 *-2 auto-select Q for overlap-add
2980 *-1 auto-select algorithm and parameters
2981 * 0 straightforward formula for small N's
2982 * 1 general FFT-based code
2983 * 2 overlap-add with length Q
2984 Q - length for overlap-add
2985
2986 OUTPUT PARAMETERS
2987 R - convolution: A*B. array[0..N+M-1].
2988
2989 -- ALGLIB --
2990 Copyright 21.07.2009 by Bochkanov Sergey
2991 *************************************************************************/
convr1dx(ae_vector * a,ae_int_t m,ae_vector * b,ae_int_t n,ae_bool circular,ae_int_t alg,ae_int_t q,ae_vector * r,ae_state * _state)2992 void convr1dx(/* Real */ ae_vector* a,
2993 ae_int_t m,
2994 /* Real */ ae_vector* b,
2995 ae_int_t n,
2996 ae_bool circular,
2997 ae_int_t alg,
2998 ae_int_t q,
2999 /* Real */ ae_vector* r,
3000 ae_state *_state)
3001 {
3002 ae_frame _frame_block;
3003 double v;
3004 ae_int_t i;
3005 ae_int_t j;
3006 ae_int_t p;
3007 ae_int_t ptotal;
3008 ae_int_t i1;
3009 ae_int_t i2;
3010 ae_int_t j1;
3011 ae_int_t j2;
3012 double ax;
3013 double ay;
3014 double bx;
3015 double by;
3016 double tx;
3017 double ty;
3018 double flopcand;
3019 double flopbest;
3020 ae_int_t algbest;
3021 fasttransformplan plan;
3022 ae_vector buf;
3023 ae_vector buf2;
3024 ae_vector buf3;
3025
3026 ae_frame_make(_state, &_frame_block);
3027 memset(&plan, 0, sizeof(plan));
3028 memset(&buf, 0, sizeof(buf));
3029 memset(&buf2, 0, sizeof(buf2));
3030 memset(&buf3, 0, sizeof(buf3));
3031 ae_vector_clear(r);
3032 _fasttransformplan_init(&plan, _state, ae_true);
3033 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
3034 ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
3035 ae_vector_init(&buf3, 0, DT_REAL, _state, ae_true);
3036
3037 ae_assert(n>0&&m>0, "ConvC1DX: incorrect N or M!", _state);
3038 ae_assert(n<=m, "ConvC1DX: N<M assumption is false!", _state);
3039
3040 /*
3041 * handle special cases
3042 */
3043 if( ae_minint(m, n, _state)<=2 )
3044 {
3045 alg = 0;
3046 }
3047
3048 /*
3049 * Auto-select
3050 */
3051 if( alg<0 )
3052 {
3053
3054 /*
3055 * Initial candidate: straightforward implementation.
3056 *
3057 * If we want to use auto-fitted overlap-add,
3058 * flop count is initialized by large real number - to force
3059 * another algorithm selection
3060 */
3061 algbest = 0;
3062 if( alg==-1 )
3063 {
3064 flopbest = 0.15*m*n;
3065 }
3066 else
3067 {
3068 flopbest = ae_maxrealnumber;
3069 }
3070
3071 /*
3072 * Another candidate - generic FFT code
3073 */
3074 if( alg==-1 )
3075 {
3076 if( (circular&&ftbaseissmooth(m, _state))&&m%2==0 )
3077 {
3078
3079 /*
3080 * special code for circular convolution of a sequence with a smooth length
3081 */
3082 flopcand = 3*ftbasegetflopestimate(m/2, _state)+(double)(6*m)/(double)2;
3083 if( ae_fp_less(flopcand,flopbest) )
3084 {
3085 algbest = 1;
3086 flopbest = flopcand;
3087 }
3088 }
3089 else
3090 {
3091
3092 /*
3093 * general cyclic/non-cyclic convolution
3094 */
3095 p = ftbasefindsmootheven(m+n-1, _state);
3096 flopcand = 3*ftbasegetflopestimate(p/2, _state)+(double)(6*p)/(double)2;
3097 if( ae_fp_less(flopcand,flopbest) )
3098 {
3099 algbest = 1;
3100 flopbest = flopcand;
3101 }
3102 }
3103 }
3104
3105 /*
3106 * Another candidate - overlap-add
3107 */
3108 q = 1;
3109 ptotal = 1;
3110 while(ptotal<n)
3111 {
3112 ptotal = ptotal*2;
3113 }
3114 while(ptotal<=m+n-1)
3115 {
3116 p = ptotal-n+1;
3117 flopcand = ae_iceil((double)m/(double)p, _state)*(2*ftbasegetflopestimate(ptotal/2, _state)+1*(ptotal/2));
3118 if( ae_fp_less(flopcand,flopbest) )
3119 {
3120 flopbest = flopcand;
3121 algbest = 2;
3122 q = p;
3123 }
3124 ptotal = ptotal*2;
3125 }
3126 alg = algbest;
3127 convr1dx(a, m, b, n, circular, alg, q, r, _state);
3128 ae_frame_leave(_state);
3129 return;
3130 }
3131
3132 /*
3133 * straightforward formula for
3134 * circular and non-circular convolutions.
3135 *
3136 * Very simple code, no further comments needed.
3137 */
3138 if( alg==0 )
3139 {
3140
3141 /*
3142 * Special case: N=1
3143 */
3144 if( n==1 )
3145 {
3146 ae_vector_set_length(r, m, _state);
3147 v = b->ptr.p_double[0];
3148 ae_v_moved(&r->ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1), v);
3149 ae_frame_leave(_state);
3150 return;
3151 }
3152
3153 /*
3154 * use straightforward formula
3155 */
3156 if( circular )
3157 {
3158
3159 /*
3160 * circular convolution
3161 */
3162 ae_vector_set_length(r, m, _state);
3163 v = b->ptr.p_double[0];
3164 ae_v_moved(&r->ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1), v);
3165 for(i=1; i<=n-1; i++)
3166 {
3167 v = b->ptr.p_double[i];
3168 i1 = 0;
3169 i2 = i-1;
3170 j1 = m-i;
3171 j2 = m-1;
3172 ae_v_addd(&r->ptr.p_double[i1], 1, &a->ptr.p_double[j1], 1, ae_v_len(i1,i2), v);
3173 i1 = i;
3174 i2 = m-1;
3175 j1 = 0;
3176 j2 = m-i-1;
3177 ae_v_addd(&r->ptr.p_double[i1], 1, &a->ptr.p_double[j1], 1, ae_v_len(i1,i2), v);
3178 }
3179 }
3180 else
3181 {
3182
3183 /*
3184 * non-circular convolution
3185 */
3186 ae_vector_set_length(r, m+n-1, _state);
3187 for(i=0; i<=m+n-2; i++)
3188 {
3189 r->ptr.p_double[i] = (double)(0);
3190 }
3191 for(i=0; i<=n-1; i++)
3192 {
3193 v = b->ptr.p_double[i];
3194 ae_v_addd(&r->ptr.p_double[i], 1, &a->ptr.p_double[0], 1, ae_v_len(i,i+m-1), v);
3195 }
3196 }
3197 ae_frame_leave(_state);
3198 return;
3199 }
3200
3201 /*
3202 * general FFT-based code for
3203 * circular and non-circular convolutions.
3204 *
3205 * First, if convolution is circular, we test whether M is smooth or not.
3206 * If it is smooth, we just use M-length FFT to calculate convolution.
3207 * If it is not, we calculate non-circular convolution and wrap it arount.
3208 *
3209 * If convolution is non-circular, we use zero-padding + FFT.
3210 *
3211 * We assume that M+N-1>2 - we should call small case code otherwise
3212 */
3213 if( alg==1 )
3214 {
3215 ae_assert(m+n-1>2, "ConvR1DX: internal error!", _state);
3216 if( (circular&&ftbaseissmooth(m, _state))&&m%2==0 )
3217 {
3218
3219 /*
3220 * special code for circular convolution with smooth even M
3221 */
3222 ae_vector_set_length(&buf, m, _state);
3223 ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1));
3224 ae_vector_set_length(&buf2, m, _state);
3225 ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
3226 for(i=n; i<=m-1; i++)
3227 {
3228 buf2.ptr.p_double[i] = (double)(0);
3229 }
3230 ae_vector_set_length(&buf3, m, _state);
3231 ftcomplexfftplan(m/2, 1, &plan, _state);
3232 fftr1dinternaleven(&buf, m, &buf3, &plan, _state);
3233 fftr1dinternaleven(&buf2, m, &buf3, &plan, _state);
3234 buf.ptr.p_double[0] = buf.ptr.p_double[0]*buf2.ptr.p_double[0];
3235 buf.ptr.p_double[1] = buf.ptr.p_double[1]*buf2.ptr.p_double[1];
3236 for(i=1; i<=m/2-1; i++)
3237 {
3238 ax = buf.ptr.p_double[2*i+0];
3239 ay = buf.ptr.p_double[2*i+1];
3240 bx = buf2.ptr.p_double[2*i+0];
3241 by = buf2.ptr.p_double[2*i+1];
3242 tx = ax*bx-ay*by;
3243 ty = ax*by+ay*bx;
3244 buf.ptr.p_double[2*i+0] = tx;
3245 buf.ptr.p_double[2*i+1] = ty;
3246 }
3247 fftr1dinvinternaleven(&buf, m, &buf3, &plan, _state);
3248 ae_vector_set_length(r, m, _state);
3249 ae_v_move(&r->ptr.p_double[0], 1, &buf.ptr.p_double[0], 1, ae_v_len(0,m-1));
3250 }
3251 else
3252 {
3253
3254 /*
3255 * M is non-smooth or non-even, general code (circular/non-circular):
3256 * * first part is the same for circular and non-circular
3257 * convolutions. zero padding, FFTs, inverse FFTs
3258 * * second part differs:
3259 * * for non-circular convolution we just copy array
3260 * * for circular convolution we add array tail to its head
3261 */
3262 p = ftbasefindsmootheven(m+n-1, _state);
3263 ae_vector_set_length(&buf, p, _state);
3264 ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1));
3265 for(i=m; i<=p-1; i++)
3266 {
3267 buf.ptr.p_double[i] = (double)(0);
3268 }
3269 ae_vector_set_length(&buf2, p, _state);
3270 ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
3271 for(i=n; i<=p-1; i++)
3272 {
3273 buf2.ptr.p_double[i] = (double)(0);
3274 }
3275 ae_vector_set_length(&buf3, p, _state);
3276 ftcomplexfftplan(p/2, 1, &plan, _state);
3277 fftr1dinternaleven(&buf, p, &buf3, &plan, _state);
3278 fftr1dinternaleven(&buf2, p, &buf3, &plan, _state);
3279 buf.ptr.p_double[0] = buf.ptr.p_double[0]*buf2.ptr.p_double[0];
3280 buf.ptr.p_double[1] = buf.ptr.p_double[1]*buf2.ptr.p_double[1];
3281 for(i=1; i<=p/2-1; i++)
3282 {
3283 ax = buf.ptr.p_double[2*i+0];
3284 ay = buf.ptr.p_double[2*i+1];
3285 bx = buf2.ptr.p_double[2*i+0];
3286 by = buf2.ptr.p_double[2*i+1];
3287 tx = ax*bx-ay*by;
3288 ty = ax*by+ay*bx;
3289 buf.ptr.p_double[2*i+0] = tx;
3290 buf.ptr.p_double[2*i+1] = ty;
3291 }
3292 fftr1dinvinternaleven(&buf, p, &buf3, &plan, _state);
3293 if( circular )
3294 {
3295
3296 /*
3297 * circular, add tail to head
3298 */
3299 ae_vector_set_length(r, m, _state);
3300 ae_v_move(&r->ptr.p_double[0], 1, &buf.ptr.p_double[0], 1, ae_v_len(0,m-1));
3301 if( n>=2 )
3302 {
3303 ae_v_add(&r->ptr.p_double[0], 1, &buf.ptr.p_double[m], 1, ae_v_len(0,n-2));
3304 }
3305 }
3306 else
3307 {
3308
3309 /*
3310 * non-circular, just copy
3311 */
3312 ae_vector_set_length(r, m+n-1, _state);
3313 ae_v_move(&r->ptr.p_double[0], 1, &buf.ptr.p_double[0], 1, ae_v_len(0,m+n-2));
3314 }
3315 }
3316 ae_frame_leave(_state);
3317 return;
3318 }
3319
3320 /*
3321 * overlap-add method
3322 */
3323 if( alg==2 )
3324 {
3325 ae_assert((q+n-1)%2==0, "ConvR1DX: internal error!", _state);
3326 ae_vector_set_length(&buf, q+n-1, _state);
3327 ae_vector_set_length(&buf2, q+n-1, _state);
3328 ae_vector_set_length(&buf3, q+n-1, _state);
3329 ftcomplexfftplan((q+n-1)/2, 1, &plan, _state);
3330
3331 /*
3332 * prepare R
3333 */
3334 if( circular )
3335 {
3336 ae_vector_set_length(r, m, _state);
3337 for(i=0; i<=m-1; i++)
3338 {
3339 r->ptr.p_double[i] = (double)(0);
3340 }
3341 }
3342 else
3343 {
3344 ae_vector_set_length(r, m+n-1, _state);
3345 for(i=0; i<=m+n-2; i++)
3346 {
3347 r->ptr.p_double[i] = (double)(0);
3348 }
3349 }
3350
3351 /*
3352 * pre-calculated FFT(B)
3353 */
3354 ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
3355 for(j=n; j<=q+n-2; j++)
3356 {
3357 buf2.ptr.p_double[j] = (double)(0);
3358 }
3359 fftr1dinternaleven(&buf2, q+n-1, &buf3, &plan, _state);
3360
3361 /*
3362 * main overlap-add cycle
3363 */
3364 i = 0;
3365 while(i<=m-1)
3366 {
3367 p = ae_minint(q, m-i, _state);
3368 ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[i], 1, ae_v_len(0,p-1));
3369 for(j=p; j<=q+n-2; j++)
3370 {
3371 buf.ptr.p_double[j] = (double)(0);
3372 }
3373 fftr1dinternaleven(&buf, q+n-1, &buf3, &plan, _state);
3374 buf.ptr.p_double[0] = buf.ptr.p_double[0]*buf2.ptr.p_double[0];
3375 buf.ptr.p_double[1] = buf.ptr.p_double[1]*buf2.ptr.p_double[1];
3376 for(j=1; j<=(q+n-1)/2-1; j++)
3377 {
3378 ax = buf.ptr.p_double[2*j+0];
3379 ay = buf.ptr.p_double[2*j+1];
3380 bx = buf2.ptr.p_double[2*j+0];
3381 by = buf2.ptr.p_double[2*j+1];
3382 tx = ax*bx-ay*by;
3383 ty = ax*by+ay*bx;
3384 buf.ptr.p_double[2*j+0] = tx;
3385 buf.ptr.p_double[2*j+1] = ty;
3386 }
3387 fftr1dinvinternaleven(&buf, q+n-1, &buf3, &plan, _state);
3388 if( circular )
3389 {
3390 j1 = ae_minint(i+p+n-2, m-1, _state)-i;
3391 j2 = j1+1;
3392 }
3393 else
3394 {
3395 j1 = p+n-2;
3396 j2 = j1+1;
3397 }
3398 ae_v_add(&r->ptr.p_double[i], 1, &buf.ptr.p_double[0], 1, ae_v_len(i,i+j1));
3399 if( p+n-2>=j2 )
3400 {
3401 ae_v_add(&r->ptr.p_double[0], 1, &buf.ptr.p_double[j2], 1, ae_v_len(0,p+n-2-j2));
3402 }
3403 i = i+p;
3404 }
3405 ae_frame_leave(_state);
3406 return;
3407 }
3408 ae_frame_leave(_state);
3409 }
3410
3411
3412 #endif
3413 #if defined(AE_COMPILE_CORR) || !defined(AE_PARTIAL_BUILD)
3414
3415
3416 /*************************************************************************
3417 1-dimensional complex cross-correlation.
3418
3419 For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).
3420
3421 Correlation is calculated using reduction to convolution. Algorithm with
3422 max(N,N)*log(max(N,N)) complexity is used (see ConvC1D() for more info
3423 about performance).
3424
3425 IMPORTANT:
3426 for historical reasons subroutine accepts its parameters in reversed
3427 order: CorrC1D(Signal, Pattern) = Pattern x Signal (using traditional
3428 definition of cross-correlation, denoting cross-correlation as "x").
3429
3430 INPUT PARAMETERS
3431 Signal - array[0..N-1] - complex function to be transformed,
3432 signal containing pattern
3433 N - problem size
3434 Pattern - array[0..M-1] - complex function to be transformed,
3435 pattern to search withing signal
3436 M - problem size
3437
3438 OUTPUT PARAMETERS
3439 R - cross-correlation, array[0..N+M-2]:
3440 * positive lags are stored in R[0..N-1],
3441 R[i] = sum(conj(pattern[j])*signal[i+j]
3442 * negative lags are stored in R[N..N+M-2],
3443 R[N+M-1-i] = sum(conj(pattern[j])*signal[-i+j]
3444
3445 NOTE:
3446 It is assumed that pattern domain is [0..M-1]. If Pattern is non-zero
3447 on [-K..M-1], you can still use this subroutine, just shift result by K.
3448
3449 -- ALGLIB --
3450 Copyright 21.07.2009 by Bochkanov Sergey
3451 *************************************************************************/
corrc1d(ae_vector * signal,ae_int_t n,ae_vector * pattern,ae_int_t m,ae_vector * r,ae_state * _state)3452 void corrc1d(/* Complex */ ae_vector* signal,
3453 ae_int_t n,
3454 /* Complex */ ae_vector* pattern,
3455 ae_int_t m,
3456 /* Complex */ ae_vector* r,
3457 ae_state *_state)
3458 {
3459 ae_frame _frame_block;
3460 ae_vector p;
3461 ae_vector b;
3462 ae_int_t i;
3463
3464 ae_frame_make(_state, &_frame_block);
3465 memset(&p, 0, sizeof(p));
3466 memset(&b, 0, sizeof(b));
3467 ae_vector_clear(r);
3468 ae_vector_init(&p, 0, DT_COMPLEX, _state, ae_true);
3469 ae_vector_init(&b, 0, DT_COMPLEX, _state, ae_true);
3470
3471 ae_assert(n>0&&m>0, "CorrC1D: incorrect N or M!", _state);
3472 ae_vector_set_length(&p, m, _state);
3473 for(i=0; i<=m-1; i++)
3474 {
3475 p.ptr.p_complex[m-1-i] = ae_c_conj(pattern->ptr.p_complex[i], _state);
3476 }
3477 convc1d(&p, m, signal, n, &b, _state);
3478 ae_vector_set_length(r, m+n-1, _state);
3479 ae_v_cmove(&r->ptr.p_complex[0], 1, &b.ptr.p_complex[m-1], 1, "N", ae_v_len(0,n-1));
3480 if( m+n-2>=n )
3481 {
3482 ae_v_cmove(&r->ptr.p_complex[n], 1, &b.ptr.p_complex[0], 1, "N", ae_v_len(n,m+n-2));
3483 }
3484 ae_frame_leave(_state);
3485 }
3486
3487
3488 /*************************************************************************
3489 1-dimensional circular complex cross-correlation.
3490
3491 For given Pattern/Signal returns corr(Pattern,Signal) (circular).
3492 Algorithm has linearithmic complexity for any M/N.
3493
3494 IMPORTANT:
3495 for historical reasons subroutine accepts its parameters in reversed
3496 order: CorrC1DCircular(Signal, Pattern) = Pattern x Signal (using
3497 traditional definition of cross-correlation, denoting cross-correlation
3498 as "x").
3499
3500 INPUT PARAMETERS
3501 Signal - array[0..N-1] - complex function to be transformed,
3502 periodic signal containing pattern
3503 N - problem size
3504 Pattern - array[0..M-1] - complex function to be transformed,
3505 non-periodic pattern to search withing signal
3506 M - problem size
3507
3508 OUTPUT PARAMETERS
3509 R - convolution: A*B. array[0..M-1].
3510
3511
3512 -- ALGLIB --
3513 Copyright 21.07.2009 by Bochkanov Sergey
3514 *************************************************************************/
corrc1dcircular(ae_vector * signal,ae_int_t m,ae_vector * pattern,ae_int_t n,ae_vector * c,ae_state * _state)3515 void corrc1dcircular(/* Complex */ ae_vector* signal,
3516 ae_int_t m,
3517 /* Complex */ ae_vector* pattern,
3518 ae_int_t n,
3519 /* Complex */ ae_vector* c,
3520 ae_state *_state)
3521 {
3522 ae_frame _frame_block;
3523 ae_vector p;
3524 ae_vector b;
3525 ae_int_t i1;
3526 ae_int_t i2;
3527 ae_int_t i;
3528 ae_int_t j2;
3529
3530 ae_frame_make(_state, &_frame_block);
3531 memset(&p, 0, sizeof(p));
3532 memset(&b, 0, sizeof(b));
3533 ae_vector_clear(c);
3534 ae_vector_init(&p, 0, DT_COMPLEX, _state, ae_true);
3535 ae_vector_init(&b, 0, DT_COMPLEX, _state, ae_true);
3536
3537 ae_assert(n>0&&m>0, "ConvC1DCircular: incorrect N or M!", _state);
3538
3539 /*
3540 * normalize task: make M>=N,
3541 * so A will be longer (at least - not shorter) that B.
3542 */
3543 if( m<n )
3544 {
3545 ae_vector_set_length(&b, m, _state);
3546 for(i1=0; i1<=m-1; i1++)
3547 {
3548 b.ptr.p_complex[i1] = ae_complex_from_i(0);
3549 }
3550 i1 = 0;
3551 while(i1<n)
3552 {
3553 i2 = ae_minint(i1+m-1, n-1, _state);
3554 j2 = i2-i1;
3555 ae_v_cadd(&b.ptr.p_complex[0], 1, &pattern->ptr.p_complex[i1], 1, "N", ae_v_len(0,j2));
3556 i1 = i1+m;
3557 }
3558 corrc1dcircular(signal, m, &b, m, c, _state);
3559 ae_frame_leave(_state);
3560 return;
3561 }
3562
3563 /*
3564 * Task is normalized
3565 */
3566 ae_vector_set_length(&p, n, _state);
3567 for(i=0; i<=n-1; i++)
3568 {
3569 p.ptr.p_complex[n-1-i] = ae_c_conj(pattern->ptr.p_complex[i], _state);
3570 }
3571 convc1dcircular(signal, m, &p, n, &b, _state);
3572 ae_vector_set_length(c, m, _state);
3573 ae_v_cmove(&c->ptr.p_complex[0], 1, &b.ptr.p_complex[n-1], 1, "N", ae_v_len(0,m-n));
3574 if( m-n+1<=m-1 )
3575 {
3576 ae_v_cmove(&c->ptr.p_complex[m-n+1], 1, &b.ptr.p_complex[0], 1, "N", ae_v_len(m-n+1,m-1));
3577 }
3578 ae_frame_leave(_state);
3579 }
3580
3581
3582 /*************************************************************************
3583 1-dimensional real cross-correlation.
3584
3585 For given Pattern/Signal returns corr(Pattern,Signal) (non-circular).
3586
3587 Correlation is calculated using reduction to convolution. Algorithm with
3588 max(N,N)*log(max(N,N)) complexity is used (see ConvC1D() for more info
3589 about performance).
3590
3591 IMPORTANT:
3592 for historical reasons subroutine accepts its parameters in reversed
3593 order: CorrR1D(Signal, Pattern) = Pattern x Signal (using traditional
3594 definition of cross-correlation, denoting cross-correlation as "x").
3595
3596 INPUT PARAMETERS
3597 Signal - array[0..N-1] - real function to be transformed,
3598 signal containing pattern
3599 N - problem size
3600 Pattern - array[0..M-1] - real function to be transformed,
3601 pattern to search withing signal
3602 M - problem size
3603
3604 OUTPUT PARAMETERS
3605 R - cross-correlation, array[0..N+M-2]:
3606 * positive lags are stored in R[0..N-1],
3607 R[i] = sum(pattern[j]*signal[i+j]
3608 * negative lags are stored in R[N..N+M-2],
3609 R[N+M-1-i] = sum(pattern[j]*signal[-i+j]
3610
3611 NOTE:
3612 It is assumed that pattern domain is [0..M-1]. If Pattern is non-zero
3613 on [-K..M-1], you can still use this subroutine, just shift result by K.
3614
3615 -- ALGLIB --
3616 Copyright 21.07.2009 by Bochkanov Sergey
3617 *************************************************************************/
corrr1d(ae_vector * signal,ae_int_t n,ae_vector * pattern,ae_int_t m,ae_vector * r,ae_state * _state)3618 void corrr1d(/* Real */ ae_vector* signal,
3619 ae_int_t n,
3620 /* Real */ ae_vector* pattern,
3621 ae_int_t m,
3622 /* Real */ ae_vector* r,
3623 ae_state *_state)
3624 {
3625 ae_frame _frame_block;
3626 ae_vector p;
3627 ae_vector b;
3628 ae_int_t i;
3629
3630 ae_frame_make(_state, &_frame_block);
3631 memset(&p, 0, sizeof(p));
3632 memset(&b, 0, sizeof(b));
3633 ae_vector_clear(r);
3634 ae_vector_init(&p, 0, DT_REAL, _state, ae_true);
3635 ae_vector_init(&b, 0, DT_REAL, _state, ae_true);
3636
3637 ae_assert(n>0&&m>0, "CorrR1D: incorrect N or M!", _state);
3638 ae_vector_set_length(&p, m, _state);
3639 for(i=0; i<=m-1; i++)
3640 {
3641 p.ptr.p_double[m-1-i] = pattern->ptr.p_double[i];
3642 }
3643 convr1d(&p, m, signal, n, &b, _state);
3644 ae_vector_set_length(r, m+n-1, _state);
3645 ae_v_move(&r->ptr.p_double[0], 1, &b.ptr.p_double[m-1], 1, ae_v_len(0,n-1));
3646 if( m+n-2>=n )
3647 {
3648 ae_v_move(&r->ptr.p_double[n], 1, &b.ptr.p_double[0], 1, ae_v_len(n,m+n-2));
3649 }
3650 ae_frame_leave(_state);
3651 }
3652
3653
3654 /*************************************************************************
3655 1-dimensional circular real cross-correlation.
3656
3657 For given Pattern/Signal returns corr(Pattern,Signal) (circular).
3658 Algorithm has linearithmic complexity for any M/N.
3659
3660 IMPORTANT:
3661 for historical reasons subroutine accepts its parameters in reversed
3662 order: CorrR1DCircular(Signal, Pattern) = Pattern x Signal (using
3663 traditional definition of cross-correlation, denoting cross-correlation
3664 as "x").
3665
3666 INPUT PARAMETERS
3667 Signal - array[0..N-1] - real function to be transformed,
3668 periodic signal containing pattern
3669 N - problem size
3670 Pattern - array[0..M-1] - real function to be transformed,
3671 non-periodic pattern to search withing signal
3672 M - problem size
3673
3674 OUTPUT PARAMETERS
3675 R - convolution: A*B. array[0..M-1].
3676
3677
3678 -- ALGLIB --
3679 Copyright 21.07.2009 by Bochkanov Sergey
3680 *************************************************************************/
corrr1dcircular(ae_vector * signal,ae_int_t m,ae_vector * pattern,ae_int_t n,ae_vector * c,ae_state * _state)3681 void corrr1dcircular(/* Real */ ae_vector* signal,
3682 ae_int_t m,
3683 /* Real */ ae_vector* pattern,
3684 ae_int_t n,
3685 /* Real */ ae_vector* c,
3686 ae_state *_state)
3687 {
3688 ae_frame _frame_block;
3689 ae_vector p;
3690 ae_vector b;
3691 ae_int_t i1;
3692 ae_int_t i2;
3693 ae_int_t i;
3694 ae_int_t j2;
3695
3696 ae_frame_make(_state, &_frame_block);
3697 memset(&p, 0, sizeof(p));
3698 memset(&b, 0, sizeof(b));
3699 ae_vector_clear(c);
3700 ae_vector_init(&p, 0, DT_REAL, _state, ae_true);
3701 ae_vector_init(&b, 0, DT_REAL, _state, ae_true);
3702
3703 ae_assert(n>0&&m>0, "ConvC1DCircular: incorrect N or M!", _state);
3704
3705 /*
3706 * normalize task: make M>=N,
3707 * so A will be longer (at least - not shorter) that B.
3708 */
3709 if( m<n )
3710 {
3711 ae_vector_set_length(&b, m, _state);
3712 for(i1=0; i1<=m-1; i1++)
3713 {
3714 b.ptr.p_double[i1] = (double)(0);
3715 }
3716 i1 = 0;
3717 while(i1<n)
3718 {
3719 i2 = ae_minint(i1+m-1, n-1, _state);
3720 j2 = i2-i1;
3721 ae_v_add(&b.ptr.p_double[0], 1, &pattern->ptr.p_double[i1], 1, ae_v_len(0,j2));
3722 i1 = i1+m;
3723 }
3724 corrr1dcircular(signal, m, &b, m, c, _state);
3725 ae_frame_leave(_state);
3726 return;
3727 }
3728
3729 /*
3730 * Task is normalized
3731 */
3732 ae_vector_set_length(&p, n, _state);
3733 for(i=0; i<=n-1; i++)
3734 {
3735 p.ptr.p_double[n-1-i] = pattern->ptr.p_double[i];
3736 }
3737 convr1dcircular(signal, m, &p, n, &b, _state);
3738 ae_vector_set_length(c, m, _state);
3739 ae_v_move(&c->ptr.p_double[0], 1, &b.ptr.p_double[n-1], 1, ae_v_len(0,m-n));
3740 if( m-n+1<=m-1 )
3741 {
3742 ae_v_move(&c->ptr.p_double[m-n+1], 1, &b.ptr.p_double[0], 1, ae_v_len(m-n+1,m-1));
3743 }
3744 ae_frame_leave(_state);
3745 }
3746
3747
3748 #endif
3749
3750 }
3751
3752