Line data Source code
1 : /*************************************************************************
2 : ALGLIB 3.17.0 (source code generated 2020-12-27)
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 : *************************************************************************/
89 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
94 0 : if( setjmp(_break_jump) )
95 : {
96 : #if !defined(AE_NO_EXCEPTIONS)
97 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
104 0 : if( _xparams.flags!=0x0 )
105 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
106 0 : alglib_impl::fftc1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
107 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
108 0 : 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)
141 0 : 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 0 : n = a.length();
148 0 : alglib_impl::ae_state_init(&_alglib_env_state);
149 0 : if( setjmp(_break_jump) )
150 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
151 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
152 0 : if( _xparams.flags!=0x0 )
153 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
154 0 : alglib_impl::fftc1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
155 :
156 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
157 0 : 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 : *************************************************************************/
181 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
186 0 : if( setjmp(_break_jump) )
187 : {
188 : #if !defined(AE_NO_EXCEPTIONS)
189 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
196 0 : if( _xparams.flags!=0x0 )
197 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
198 0 : alglib_impl::fftc1dinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
199 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
200 0 : 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)
224 0 : 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 0 : n = a.length();
231 0 : alglib_impl::ae_state_init(&_alglib_env_state);
232 0 : if( setjmp(_break_jump) )
233 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
234 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
235 0 : if( _xparams.flags!=0x0 )
236 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
237 0 : alglib_impl::fftc1dinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
238 :
239 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
240 0 : 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 : *************************************************************************/
267 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
272 0 : if( setjmp(_break_jump) )
273 : {
274 : #if !defined(AE_NO_EXCEPTIONS)
275 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
282 0 : if( _xparams.flags!=0x0 )
283 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
284 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
286 0 : 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)
313 0 : 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 0 : n = a.length();
320 0 : alglib_impl::ae_state_init(&_alglib_env_state);
321 0 : if( setjmp(_break_jump) )
322 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
323 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
324 0 : if( _xparams.flags!=0x0 )
325 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
326 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
329 0 : 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 : *************************************************************************/
368 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
373 0 : if( setjmp(_break_jump) )
374 : {
375 : #if !defined(AE_NO_EXCEPTIONS)
376 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
383 0 : if( _xparams.flags!=0x0 )
384 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
385 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
387 0 : 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)
426 0 : 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 0 : n = f.length();
433 0 : alglib_impl::ae_state_init(&_alglib_env_state);
434 0 : if( setjmp(_break_jump) )
435 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
436 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
437 0 : if( _xparams.flags!=0x0 )
438 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
439 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
442 0 : 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 : *************************************************************************/
465 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
470 0 : if( setjmp(_break_jump) )
471 : {
472 : #if !defined(AE_NO_EXCEPTIONS)
473 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
480 0 : if( _xparams.flags!=0x0 )
481 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
482 0 : alglib_impl::fhtr1d(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
483 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
484 0 : 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 : *************************************************************************/
503 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
508 0 : if( setjmp(_break_jump) )
509 : {
510 : #if !defined(AE_NO_EXCEPTIONS)
511 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
518 0 : if( _xparams.flags!=0x0 )
519 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
520 0 : alglib_impl::fhtr1dinv(const_cast<alglib_impl::ae_vector*>(a.c_ptr()), n, &_alglib_env_state);
521 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
522 0 : 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 : *************************************************************************/
556 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
561 0 : if( setjmp(_break_jump) )
562 : {
563 : #if !defined(AE_NO_EXCEPTIONS)
564 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
571 0 : if( _xparams.flags!=0x0 )
572 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
573 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
575 0 : 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 : *************************************************************************/
604 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
609 0 : if( setjmp(_break_jump) )
610 : {
611 : #if !defined(AE_NO_EXCEPTIONS)
612 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
619 0 : if( _xparams.flags!=0x0 )
620 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
621 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
623 0 : 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 : *************************************************************************/
654 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
659 0 : if( setjmp(_break_jump) )
660 : {
661 : #if !defined(AE_NO_EXCEPTIONS)
662 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
669 0 : if( _xparams.flags!=0x0 )
670 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
671 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
673 0 : 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 : *************************************************************************/
702 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
707 0 : if( setjmp(_break_jump) )
708 : {
709 : #if !defined(AE_NO_EXCEPTIONS)
710 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
717 0 : if( _xparams.flags!=0x0 )
718 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
719 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
721 0 : 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 : *************************************************************************/
746 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
751 0 : if( setjmp(_break_jump) )
752 : {
753 : #if !defined(AE_NO_EXCEPTIONS)
754 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
761 0 : if( _xparams.flags!=0x0 )
762 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
763 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
765 0 : 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 : *************************************************************************/
794 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
799 0 : if( setjmp(_break_jump) )
800 : {
801 : #if !defined(AE_NO_EXCEPTIONS)
802 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
809 0 : if( _xparams.flags!=0x0 )
810 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
811 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
813 0 : 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 : *************************************************************************/
838 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
843 0 : if( setjmp(_break_jump) )
844 : {
845 : #if !defined(AE_NO_EXCEPTIONS)
846 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
853 0 : if( _xparams.flags!=0x0 )
854 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
855 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
857 0 : 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 : *************************************************************************/
886 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
891 0 : if( setjmp(_break_jump) )
892 : {
893 : #if !defined(AE_NO_EXCEPTIONS)
894 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
901 0 : if( _xparams.flags!=0x0 )
902 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
903 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
905 0 : 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 : *************************************************************************/
946 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
951 0 : if( setjmp(_break_jump) )
952 : {
953 : #if !defined(AE_NO_EXCEPTIONS)
954 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
961 0 : if( _xparams.flags!=0x0 )
962 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
963 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
965 0 : 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 : *************************************************************************/
995 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
1000 0 : if( setjmp(_break_jump) )
1001 : {
1002 : #if !defined(AE_NO_EXCEPTIONS)
1003 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1010 0 : if( _xparams.flags!=0x0 )
1011 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1012 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1014 0 : 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 : *************************************************************************/
1053 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
1058 0 : if( setjmp(_break_jump) )
1059 : {
1060 : #if !defined(AE_NO_EXCEPTIONS)
1061 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1068 0 : if( _xparams.flags!=0x0 )
1069 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1070 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1072 0 : 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 : *************************************************************************/
1102 0 : 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 0 : alglib_impl::ae_state_init(&_alglib_env_state);
1107 0 : if( setjmp(_break_jump) )
1108 : {
1109 : #if !defined(AE_NO_EXCEPTIONS)
1110 0 : _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 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1117 0 : if( _xparams.flags!=0x0 )
1118 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1119 0 : 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 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1121 0 : 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 : *************************************************************************/
1182 0 : 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 0 : ae_frame_make(_state, &_frame_block);
1190 0 : memset(&plan, 0, sizeof(plan));
1191 0 : memset(&buf, 0, sizeof(buf));
1192 0 : _fasttransformplan_init(&plan, _state, ae_true);
1193 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
1194 :
1195 0 : ae_assert(n>0, "FFTC1D: incorrect N!", _state);
1196 0 : ae_assert(a->cnt>=n, "FFTC1D: Length(A)<N!", _state);
1197 0 : 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 0 : if( n==1 )
1204 : {
1205 0 : ae_frame_leave(_state);
1206 0 : return;
1207 : }
1208 :
1209 : /*
1210 : * convert input array to the more convinient format
1211 : */
1212 0 : ae_vector_set_length(&buf, 2*n, _state);
1213 0 : for(i=0; i<=n-1; i++)
1214 : {
1215 0 : buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
1216 0 : 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 0 : ftcomplexfftplan(n, 1, &plan, _state);
1227 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
1228 :
1229 : /*
1230 : * result
1231 : */
1232 0 : for(i=0; i<=n-1; i++)
1233 : {
1234 0 : a->ptr.p_complex[i].x = buf.ptr.p_double[2*i+0];
1235 0 : a->ptr.p_complex[i].y = buf.ptr.p_double[2*i+1];
1236 : }
1237 0 : 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 : *************************************************************************/
1261 0 : void fftc1dinv(/* Complex */ ae_vector* a, ae_int_t n, ae_state *_state)
1262 : {
1263 : ae_int_t i;
1264 :
1265 :
1266 0 : ae_assert(n>0, "FFTC1DInv: incorrect N!", _state);
1267 0 : ae_assert(a->cnt>=n, "FFTC1DInv: Length(A)<N!", _state);
1268 0 : 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 0 : for(i=0; i<=n-1; i++)
1278 : {
1279 0 : a->ptr.p_complex[i].y = -a->ptr.p_complex[i].y;
1280 : }
1281 0 : fftc1d(a, n, _state);
1282 0 : for(i=0; i<=n-1; i++)
1283 : {
1284 0 : a->ptr.p_complex[i].x = a->ptr.p_complex[i].x/n;
1285 0 : a->ptr.p_complex[i].y = -a->ptr.p_complex[i].y/n;
1286 : }
1287 0 : }
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 : *************************************************************************/
1313 0 : 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 0 : ae_frame_make(_state, &_frame_block);
1329 0 : memset(&buf, 0, sizeof(buf));
1330 0 : memset(&plan, 0, sizeof(plan));
1331 0 : ae_vector_clear(f);
1332 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
1333 0 : _fasttransformplan_init(&plan, _state, ae_true);
1334 :
1335 0 : ae_assert(n>0, "FFTR1D: incorrect N!", _state);
1336 0 : ae_assert(a->cnt>=n, "FFTR1D: Length(A)<N!", _state);
1337 0 : 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 0 : if( n==1 )
1347 : {
1348 0 : ae_vector_set_length(f, 1, _state);
1349 0 : f->ptr.p_complex[0] = ae_complex_from_d(a->ptr.p_double[0]);
1350 0 : ae_frame_leave(_state);
1351 0 : return;
1352 : }
1353 0 : if( n==2 )
1354 : {
1355 0 : ae_vector_set_length(f, 2, _state);
1356 0 : f->ptr.p_complex[0].x = a->ptr.p_double[0]+a->ptr.p_double[1];
1357 0 : f->ptr.p_complex[0].y = (double)(0);
1358 0 : f->ptr.p_complex[1].x = a->ptr.p_double[0]-a->ptr.p_double[1];
1359 0 : f->ptr.p_complex[1].y = (double)(0);
1360 0 : ae_frame_leave(_state);
1361 0 : return;
1362 : }
1363 :
1364 : /*
1365 : * Choose between odd-size and even-size FFTs
1366 : */
1367 0 : if( n%2==0 )
1368 : {
1369 :
1370 : /*
1371 : * even-size real FFT, use reduction to the complex task
1372 : */
1373 0 : n2 = n/2;
1374 0 : ae_vector_set_length(&buf, n, _state);
1375 0 : ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,n-1));
1376 0 : ftcomplexfftplan(n2, 1, &plan, _state);
1377 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
1378 0 : ae_vector_set_length(f, n, _state);
1379 0 : for(i=0; i<=n2; i++)
1380 : {
1381 0 : idx = 2*(i%n2);
1382 0 : hn.x = buf.ptr.p_double[idx+0];
1383 0 : hn.y = buf.ptr.p_double[idx+1];
1384 0 : idx = 2*((n2-i)%n2);
1385 0 : hmnc.x = buf.ptr.p_double[idx+0];
1386 0 : hmnc.y = -buf.ptr.p_double[idx+1];
1387 0 : v.x = -ae_sin(-2*ae_pi*i/n, _state);
1388 0 : v.y = ae_cos(-2*ae_pi*i/n, _state);
1389 0 : f->ptr.p_complex[i] = ae_c_sub(ae_c_add(hn,hmnc),ae_c_mul(v,ae_c_sub(hn,hmnc)));
1390 0 : f->ptr.p_complex[i].x = 0.5*f->ptr.p_complex[i].x;
1391 0 : f->ptr.p_complex[i].y = 0.5*f->ptr.p_complex[i].y;
1392 : }
1393 0 : for(i=n2+1; i<=n-1; i++)
1394 : {
1395 0 : 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 0 : ae_vector_set_length(f, n, _state);
1405 0 : for(i=0; i<=n-1; i++)
1406 : {
1407 0 : f->ptr.p_complex[i] = ae_complex_from_d(a->ptr.p_double[i]);
1408 : }
1409 0 : fftc1d(f, n, _state);
1410 : }
1411 0 : 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 : *************************************************************************/
1450 0 : 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 0 : ae_frame_make(_state, &_frame_block);
1461 0 : memset(&h, 0, sizeof(h));
1462 0 : memset(&fh, 0, sizeof(fh));
1463 0 : ae_vector_clear(a);
1464 0 : ae_vector_init(&h, 0, DT_REAL, _state, ae_true);
1465 0 : ae_vector_init(&fh, 0, DT_COMPLEX, _state, ae_true);
1466 :
1467 0 : ae_assert(n>0, "FFTR1DInv: incorrect N!", _state);
1468 0 : ae_assert(f->cnt>=ae_ifloor((double)n/(double)2, _state)+1, "FFTR1DInv: Length(F)<Floor(N/2)+1!", _state);
1469 0 : ae_assert(ae_isfinite(f->ptr.p_complex[0].x, _state), "FFTR1DInv: F contains infinite or NAN values!", _state);
1470 0 : for(i=1; i<=ae_ifloor((double)n/(double)2, _state)-1; i++)
1471 : {
1472 0 : 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 0 : 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 0 : if( n%2!=0 )
1476 : {
1477 0 : 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 0 : if( n==1 )
1485 : {
1486 0 : ae_vector_set_length(a, 1, _state);
1487 0 : a->ptr.p_double[0] = f->ptr.p_complex[0].x;
1488 0 : ae_frame_leave(_state);
1489 0 : 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 0 : ae_vector_set_length(&h, n, _state);
1500 0 : ae_vector_set_length(a, n, _state);
1501 0 : h.ptr.p_double[0] = f->ptr.p_complex[0].x;
1502 0 : for(i=1; i<=ae_ifloor((double)n/(double)2, _state)-1; i++)
1503 : {
1504 0 : h.ptr.p_double[i] = f->ptr.p_complex[i].x-f->ptr.p_complex[i].y;
1505 0 : h.ptr.p_double[n-i] = f->ptr.p_complex[i].x+f->ptr.p_complex[i].y;
1506 : }
1507 0 : if( n%2==0 )
1508 : {
1509 0 : 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 0 : 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 0 : 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 0 : fftr1d(&h, n, &fh, _state);
1517 0 : for(i=0; i<=n-1; i++)
1518 : {
1519 0 : a->ptr.p_double[i] = (fh.ptr.p_complex[i].x-fh.ptr.p_complex[i].y)/n;
1520 : }
1521 0 : 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 : *************************************************************************/
1532 0 : 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 0 : 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 0 : if( n==2 )
1557 : {
1558 0 : x = a->ptr.p_double[0]+a->ptr.p_double[1];
1559 0 : y = a->ptr.p_double[0]-a->ptr.p_double[1];
1560 0 : a->ptr.p_double[0] = x;
1561 0 : a->ptr.p_double[1] = y;
1562 0 : return;
1563 : }
1564 :
1565 : /*
1566 : * even-size real FFT, use reduction to the complex task
1567 : */
1568 0 : n2 = n/2;
1569 0 : ae_v_move(&buf->ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,n-1));
1570 0 : ftapplyplan(plan, buf, 0, 1, _state);
1571 0 : a->ptr.p_double[0] = buf->ptr.p_double[0]+buf->ptr.p_double[1];
1572 0 : for(i=1; i<=n2-1; i++)
1573 : {
1574 0 : idx = 2*(i%n2);
1575 0 : hn.x = buf->ptr.p_double[idx+0];
1576 0 : hn.y = buf->ptr.p_double[idx+1];
1577 0 : idx = 2*(n2-i);
1578 0 : hmnc.x = buf->ptr.p_double[idx+0];
1579 0 : hmnc.y = -buf->ptr.p_double[idx+1];
1580 0 : v.x = -ae_sin(-2*ae_pi*i/n, _state);
1581 0 : v.y = ae_cos(-2*ae_pi*i/n, _state);
1582 0 : v = ae_c_sub(ae_c_add(hn,hmnc),ae_c_mul(v,ae_c_sub(hn,hmnc)));
1583 0 : a->ptr.p_double[2*i+0] = 0.5*v.x;
1584 0 : a->ptr.p_double[2*i+1] = 0.5*v.y;
1585 : }
1586 0 : 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 : *************************************************************************/
1597 0 : 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 0 : 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 0 : if( n==2 )
1619 : {
1620 0 : x = 0.5*(a->ptr.p_double[0]+a->ptr.p_double[1]);
1621 0 : y = 0.5*(a->ptr.p_double[0]-a->ptr.p_double[1]);
1622 0 : a->ptr.p_double[0] = x;
1623 0 : a->ptr.p_double[1] = y;
1624 0 : 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 0 : n2 = n/2;
1635 0 : buf->ptr.p_double[0] = a->ptr.p_double[0];
1636 0 : for(i=1; i<=n2-1; i++)
1637 : {
1638 0 : x = a->ptr.p_double[2*i+0];
1639 0 : y = a->ptr.p_double[2*i+1];
1640 0 : buf->ptr.p_double[i] = x-y;
1641 0 : buf->ptr.p_double[n-i] = x+y;
1642 : }
1643 0 : buf->ptr.p_double[n2] = a->ptr.p_double[1];
1644 0 : fftr1dinternaleven(buf, n, a, plan, _state);
1645 0 : a->ptr.p_double[0] = buf->ptr.p_double[0]/n;
1646 0 : t = (double)1/(double)n;
1647 0 : for(i=1; i<=n2-1; i++)
1648 : {
1649 0 : x = buf->ptr.p_double[2*i+0];
1650 0 : y = buf->ptr.p_double[2*i+1];
1651 0 : a->ptr.p_double[i] = t*(x-y);
1652 0 : a->ptr.p_double[n-i] = t*(x+y);
1653 : }
1654 0 : 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 : *************************************************************************/
1679 0 : 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 0 : ae_frame_make(_state, &_frame_block);
1686 0 : memset(&fa, 0, sizeof(fa));
1687 0 : ae_vector_init(&fa, 0, DT_COMPLEX, _state, ae_true);
1688 :
1689 0 : 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 0 : if( n==1 )
1696 : {
1697 0 : ae_frame_leave(_state);
1698 0 : return;
1699 : }
1700 :
1701 : /*
1702 : * Reduce FHt to real FFT
1703 : */
1704 0 : fftr1d(a, n, &fa, _state);
1705 0 : for(i=0; i<=n-1; i++)
1706 : {
1707 0 : a->ptr.p_double[i] = fa.ptr.p_complex[i].x-fa.ptr.p_complex[i].y;
1708 : }
1709 0 : 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 : *************************************************************************/
1729 0 : void fhtr1dinv(/* Real */ ae_vector* a, ae_int_t n, ae_state *_state)
1730 : {
1731 : ae_int_t i;
1732 :
1733 :
1734 0 : 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 0 : if( n==1 )
1741 : {
1742 0 : 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 0 : fhtr1d(a, n, _state);
1751 0 : for(i=0; i<=n-1; i++)
1752 : {
1753 0 : 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 : *************************************************************************/
1791 0 : 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 0 : ae_vector_clear(r);
1800 :
1801 0 : 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 0 : if( m<n )
1808 : {
1809 0 : convc1d(b, n, a, m, r, _state);
1810 0 : return;
1811 : }
1812 0 : 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 : *************************************************************************/
1842 0 : 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 0 : ae_frame_make(_state, &_frame_block);
1861 0 : memset(&buf, 0, sizeof(buf));
1862 0 : memset(&buf2, 0, sizeof(buf2));
1863 0 : memset(&plan, 0, sizeof(plan));
1864 0 : ae_vector_clear(r);
1865 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
1866 0 : ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
1867 0 : _fasttransformplan_init(&plan, _state, ae_true);
1868 :
1869 0 : ae_assert((n>0&&m>0)&&n<=m, "ConvC1DInv: incorrect N or M!", _state);
1870 0 : p = ftbasefindsmooth(m, _state);
1871 0 : ftcomplexfftplan(p, 1, &plan, _state);
1872 0 : ae_vector_set_length(&buf, 2*p, _state);
1873 0 : for(i=0; i<=m-1; i++)
1874 : {
1875 0 : buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
1876 0 : buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
1877 : }
1878 0 : for(i=m; i<=p-1; i++)
1879 : {
1880 0 : buf.ptr.p_double[2*i+0] = (double)(0);
1881 0 : buf.ptr.p_double[2*i+1] = (double)(0);
1882 : }
1883 0 : ae_vector_set_length(&buf2, 2*p, _state);
1884 0 : for(i=0; i<=n-1; i++)
1885 : {
1886 0 : buf2.ptr.p_double[2*i+0] = b->ptr.p_complex[i].x;
1887 0 : buf2.ptr.p_double[2*i+1] = b->ptr.p_complex[i].y;
1888 : }
1889 0 : for(i=n; i<=p-1; i++)
1890 : {
1891 0 : buf2.ptr.p_double[2*i+0] = (double)(0);
1892 0 : buf2.ptr.p_double[2*i+1] = (double)(0);
1893 : }
1894 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
1895 0 : ftapplyplan(&plan, &buf2, 0, 1, _state);
1896 0 : for(i=0; i<=p-1; i++)
1897 : {
1898 0 : c1.x = buf.ptr.p_double[2*i+0];
1899 0 : c1.y = buf.ptr.p_double[2*i+1];
1900 0 : c2.x = buf2.ptr.p_double[2*i+0];
1901 0 : c2.y = buf2.ptr.p_double[2*i+1];
1902 0 : c3 = ae_c_div(c1,c2);
1903 0 : buf.ptr.p_double[2*i+0] = c3.x;
1904 0 : buf.ptr.p_double[2*i+1] = -c3.y;
1905 : }
1906 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
1907 0 : t = (double)1/(double)p;
1908 0 : ae_vector_set_length(r, m-n+1, _state);
1909 0 : for(i=0; i<=m-n; i++)
1910 : {
1911 0 : r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
1912 0 : r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
1913 : }
1914 0 : ae_frame_leave(_state);
1915 0 : }
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 : *************************************************************************/
1946 0 : 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 0 : ae_frame_make(_state, &_frame_block);
1960 0 : memset(&buf, 0, sizeof(buf));
1961 0 : ae_vector_clear(c);
1962 0 : ae_vector_init(&buf, 0, DT_COMPLEX, _state, ae_true);
1963 :
1964 0 : 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 0 : if( m<n )
1971 : {
1972 0 : ae_vector_set_length(&buf, m, _state);
1973 0 : for(i1=0; i1<=m-1; i1++)
1974 : {
1975 0 : buf.ptr.p_complex[i1] = ae_complex_from_i(0);
1976 : }
1977 0 : i1 = 0;
1978 0 : while(i1<n)
1979 : {
1980 0 : i2 = ae_minint(i1+m-1, n-1, _state);
1981 0 : j2 = i2-i1;
1982 0 : ae_v_cadd(&buf.ptr.p_complex[0], 1, &r->ptr.p_complex[i1], 1, "N", ae_v_len(0,j2));
1983 0 : i1 = i1+m;
1984 : }
1985 0 : convc1dcircular(s, m, &buf, m, c, _state);
1986 0 : ae_frame_leave(_state);
1987 0 : return;
1988 : }
1989 0 : convc1dx(s, m, r, n, ae_true, -1, 0, c, _state);
1990 0 : 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 : *************************************************************************/
2020 0 : 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 0 : ae_frame_make(_state, &_frame_block);
2042 0 : memset(&buf, 0, sizeof(buf));
2043 0 : memset(&buf2, 0, sizeof(buf2));
2044 0 : memset(&cbuf, 0, sizeof(cbuf));
2045 0 : memset(&plan, 0, sizeof(plan));
2046 0 : ae_vector_clear(r);
2047 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2048 0 : ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
2049 0 : ae_vector_init(&cbuf, 0, DT_COMPLEX, _state, ae_true);
2050 0 : _fasttransformplan_init(&plan, _state, ae_true);
2051 :
2052 0 : 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 0 : if( m<n )
2059 : {
2060 0 : ae_vector_set_length(&cbuf, m, _state);
2061 0 : for(i=0; i<=m-1; i++)
2062 : {
2063 0 : cbuf.ptr.p_complex[i] = ae_complex_from_i(0);
2064 : }
2065 0 : i1 = 0;
2066 0 : while(i1<n)
2067 : {
2068 0 : i2 = ae_minint(i1+m-1, n-1, _state);
2069 0 : j2 = i2-i1;
2070 0 : ae_v_cadd(&cbuf.ptr.p_complex[0], 1, &b->ptr.p_complex[i1], 1, "N", ae_v_len(0,j2));
2071 0 : i1 = i1+m;
2072 : }
2073 0 : convc1dcircularinv(a, m, &cbuf, m, r, _state);
2074 0 : ae_frame_leave(_state);
2075 0 : return;
2076 : }
2077 :
2078 : /*
2079 : * Task is normalized
2080 : */
2081 0 : ftcomplexfftplan(m, 1, &plan, _state);
2082 0 : ae_vector_set_length(&buf, 2*m, _state);
2083 0 : for(i=0; i<=m-1; i++)
2084 : {
2085 0 : buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
2086 0 : buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
2087 : }
2088 0 : ae_vector_set_length(&buf2, 2*m, _state);
2089 0 : for(i=0; i<=n-1; i++)
2090 : {
2091 0 : buf2.ptr.p_double[2*i+0] = b->ptr.p_complex[i].x;
2092 0 : buf2.ptr.p_double[2*i+1] = b->ptr.p_complex[i].y;
2093 : }
2094 0 : for(i=n; i<=m-1; i++)
2095 : {
2096 0 : buf2.ptr.p_double[2*i+0] = (double)(0);
2097 0 : buf2.ptr.p_double[2*i+1] = (double)(0);
2098 : }
2099 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
2100 0 : ftapplyplan(&plan, &buf2, 0, 1, _state);
2101 0 : for(i=0; i<=m-1; i++)
2102 : {
2103 0 : c1.x = buf.ptr.p_double[2*i+0];
2104 0 : c1.y = buf.ptr.p_double[2*i+1];
2105 0 : c2.x = buf2.ptr.p_double[2*i+0];
2106 0 : c2.y = buf2.ptr.p_double[2*i+1];
2107 0 : c3 = ae_c_div(c1,c2);
2108 0 : buf.ptr.p_double[2*i+0] = c3.x;
2109 0 : buf.ptr.p_double[2*i+1] = -c3.y;
2110 : }
2111 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
2112 0 : t = (double)1/(double)m;
2113 0 : ae_vector_set_length(r, m, _state);
2114 0 : for(i=0; i<=m-1; i++)
2115 : {
2116 0 : r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
2117 0 : r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
2118 : }
2119 0 : 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 : *************************************************************************/
2145 0 : 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 0 : ae_vector_clear(r);
2154 :
2155 0 : 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 0 : if( m<n )
2162 : {
2163 0 : convr1d(b, n, a, m, r, _state);
2164 0 : return;
2165 : }
2166 0 : 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 : *************************************************************************/
2196 0 : 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 0 : ae_frame_make(_state, &_frame_block);
2215 0 : memset(&buf, 0, sizeof(buf));
2216 0 : memset(&buf2, 0, sizeof(buf2));
2217 0 : memset(&buf3, 0, sizeof(buf3));
2218 0 : memset(&plan, 0, sizeof(plan));
2219 0 : ae_vector_clear(r);
2220 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2221 0 : ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
2222 0 : ae_vector_init(&buf3, 0, DT_REAL, _state, ae_true);
2223 0 : _fasttransformplan_init(&plan, _state, ae_true);
2224 :
2225 0 : ae_assert((n>0&&m>0)&&n<=m, "ConvR1DInv: incorrect N or M!", _state);
2226 0 : p = ftbasefindsmootheven(m, _state);
2227 0 : ae_vector_set_length(&buf, p, _state);
2228 0 : ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1));
2229 0 : for(i=m; i<=p-1; i++)
2230 : {
2231 0 : buf.ptr.p_double[i] = (double)(0);
2232 : }
2233 0 : ae_vector_set_length(&buf2, p, _state);
2234 0 : ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2235 0 : for(i=n; i<=p-1; i++)
2236 : {
2237 0 : buf2.ptr.p_double[i] = (double)(0);
2238 : }
2239 0 : ae_vector_set_length(&buf3, p, _state);
2240 0 : ftcomplexfftplan(p/2, 1, &plan, _state);
2241 0 : fftr1dinternaleven(&buf, p, &buf3, &plan, _state);
2242 0 : fftr1dinternaleven(&buf2, p, &buf3, &plan, _state);
2243 0 : buf.ptr.p_double[0] = buf.ptr.p_double[0]/buf2.ptr.p_double[0];
2244 0 : buf.ptr.p_double[1] = buf.ptr.p_double[1]/buf2.ptr.p_double[1];
2245 0 : for(i=1; i<=p/2-1; i++)
2246 : {
2247 0 : c1.x = buf.ptr.p_double[2*i+0];
2248 0 : c1.y = buf.ptr.p_double[2*i+1];
2249 0 : c2.x = buf2.ptr.p_double[2*i+0];
2250 0 : c2.y = buf2.ptr.p_double[2*i+1];
2251 0 : c3 = ae_c_div(c1,c2);
2252 0 : buf.ptr.p_double[2*i+0] = c3.x;
2253 0 : buf.ptr.p_double[2*i+1] = c3.y;
2254 : }
2255 0 : fftr1dinvinternaleven(&buf, p, &buf3, &plan, _state);
2256 0 : ae_vector_set_length(r, m-n+1, _state);
2257 0 : ae_v_move(&r->ptr.p_double[0], 1, &buf.ptr.p_double[0], 1, ae_v_len(0,m-n));
2258 0 : ae_frame_leave(_state);
2259 0 : }
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 : *************************************************************************/
2284 0 : 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 0 : ae_frame_make(_state, &_frame_block);
2298 0 : memset(&buf, 0, sizeof(buf));
2299 0 : ae_vector_clear(c);
2300 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2301 :
2302 0 : 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 0 : if( m<n )
2309 : {
2310 0 : ae_vector_set_length(&buf, m, _state);
2311 0 : for(i1=0; i1<=m-1; i1++)
2312 : {
2313 0 : buf.ptr.p_double[i1] = (double)(0);
2314 : }
2315 0 : i1 = 0;
2316 0 : while(i1<n)
2317 : {
2318 0 : i2 = ae_minint(i1+m-1, n-1, _state);
2319 0 : j2 = i2-i1;
2320 0 : ae_v_add(&buf.ptr.p_double[0], 1, &r->ptr.p_double[i1], 1, ae_v_len(0,j2));
2321 0 : i1 = i1+m;
2322 : }
2323 0 : convr1dcircular(s, m, &buf, m, c, _state);
2324 0 : ae_frame_leave(_state);
2325 0 : return;
2326 : }
2327 :
2328 : /*
2329 : * reduce to usual convolution
2330 : */
2331 0 : convr1dx(s, m, r, n, ae_true, -1, 0, c, _state);
2332 0 : 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 : *************************************************************************/
2362 0 : 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 0 : ae_frame_make(_state, &_frame_block);
2385 0 : memset(&buf, 0, sizeof(buf));
2386 0 : memset(&buf2, 0, sizeof(buf2));
2387 0 : memset(&buf3, 0, sizeof(buf3));
2388 0 : memset(&cbuf, 0, sizeof(cbuf));
2389 0 : memset(&cbuf2, 0, sizeof(cbuf2));
2390 0 : memset(&plan, 0, sizeof(plan));
2391 0 : ae_vector_clear(r);
2392 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2393 0 : ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
2394 0 : ae_vector_init(&buf3, 0, DT_REAL, _state, ae_true);
2395 0 : ae_vector_init(&cbuf, 0, DT_COMPLEX, _state, ae_true);
2396 0 : ae_vector_init(&cbuf2, 0, DT_COMPLEX, _state, ae_true);
2397 0 : _fasttransformplan_init(&plan, _state, ae_true);
2398 :
2399 0 : 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 0 : if( m<n )
2406 : {
2407 0 : ae_vector_set_length(&buf, m, _state);
2408 0 : for(i=0; i<=m-1; i++)
2409 : {
2410 0 : buf.ptr.p_double[i] = (double)(0);
2411 : }
2412 0 : i1 = 0;
2413 0 : while(i1<n)
2414 : {
2415 0 : i2 = ae_minint(i1+m-1, n-1, _state);
2416 0 : j2 = i2-i1;
2417 0 : ae_v_add(&buf.ptr.p_double[0], 1, &b->ptr.p_double[i1], 1, ae_v_len(0,j2));
2418 0 : i1 = i1+m;
2419 : }
2420 0 : convr1dcircularinv(a, m, &buf, m, r, _state);
2421 0 : ae_frame_leave(_state);
2422 0 : return;
2423 : }
2424 :
2425 : /*
2426 : * Task is normalized
2427 : */
2428 0 : if( m%2==0 )
2429 : {
2430 :
2431 : /*
2432 : * size is even, use fast even-size FFT
2433 : */
2434 0 : ae_vector_set_length(&buf, m, _state);
2435 0 : ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1));
2436 0 : ae_vector_set_length(&buf2, m, _state);
2437 0 : ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2438 0 : for(i=n; i<=m-1; i++)
2439 : {
2440 0 : buf2.ptr.p_double[i] = (double)(0);
2441 : }
2442 0 : ae_vector_set_length(&buf3, m, _state);
2443 0 : ftcomplexfftplan(m/2, 1, &plan, _state);
2444 0 : fftr1dinternaleven(&buf, m, &buf3, &plan, _state);
2445 0 : fftr1dinternaleven(&buf2, m, &buf3, &plan, _state);
2446 0 : buf.ptr.p_double[0] = buf.ptr.p_double[0]/buf2.ptr.p_double[0];
2447 0 : buf.ptr.p_double[1] = buf.ptr.p_double[1]/buf2.ptr.p_double[1];
2448 0 : for(i=1; i<=m/2-1; i++)
2449 : {
2450 0 : c1.x = buf.ptr.p_double[2*i+0];
2451 0 : c1.y = buf.ptr.p_double[2*i+1];
2452 0 : c2.x = buf2.ptr.p_double[2*i+0];
2453 0 : c2.y = buf2.ptr.p_double[2*i+1];
2454 0 : c3 = ae_c_div(c1,c2);
2455 0 : buf.ptr.p_double[2*i+0] = c3.x;
2456 0 : buf.ptr.p_double[2*i+1] = c3.y;
2457 : }
2458 0 : fftr1dinvinternaleven(&buf, m, &buf3, &plan, _state);
2459 0 : ae_vector_set_length(r, m, _state);
2460 0 : 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 0 : fftr1d(a, m, &cbuf, _state);
2469 0 : ae_vector_set_length(&buf2, m, _state);
2470 0 : ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2471 0 : for(i=n; i<=m-1; i++)
2472 : {
2473 0 : buf2.ptr.p_double[i] = (double)(0);
2474 : }
2475 0 : fftr1d(&buf2, m, &cbuf2, _state);
2476 0 : for(i=0; i<=ae_ifloor((double)m/(double)2, _state); i++)
2477 : {
2478 0 : cbuf.ptr.p_complex[i] = ae_c_div(cbuf.ptr.p_complex[i],cbuf2.ptr.p_complex[i]);
2479 : }
2480 0 : fftr1dinv(&cbuf, m, r, _state);
2481 : }
2482 0 : 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 : *************************************************************************/
2511 0 : 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 0 : ae_frame_make(_state, &_frame_block);
2547 0 : memset(&bbuf, 0, sizeof(bbuf));
2548 0 : memset(&plan, 0, sizeof(plan));
2549 0 : memset(&buf, 0, sizeof(buf));
2550 0 : memset(&buf2, 0, sizeof(buf2));
2551 0 : ae_vector_clear(r);
2552 0 : ae_vector_init(&bbuf, 0, DT_COMPLEX, _state, ae_true);
2553 0 : _fasttransformplan_init(&plan, _state, ae_true);
2554 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
2555 0 : ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
2556 :
2557 0 : ae_assert(n>0&&m>0, "ConvC1DX: incorrect N or M!", _state);
2558 0 : ae_assert(n<=m, "ConvC1DX: N<M assumption is false!", _state);
2559 :
2560 : /*
2561 : * Auto-select
2562 : */
2563 0 : 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 0 : algbest = 0;
2574 0 : if( alg==-1 )
2575 : {
2576 0 : flopbest = (double)(2*m*n);
2577 : }
2578 : else
2579 : {
2580 0 : flopbest = ae_maxrealnumber;
2581 : }
2582 :
2583 : /*
2584 : * Another candidate - generic FFT code
2585 : */
2586 0 : if( alg==-1 )
2587 : {
2588 0 : if( circular&&ftbaseissmooth(m, _state) )
2589 : {
2590 :
2591 : /*
2592 : * special code for circular convolution of a sequence with a smooth length
2593 : */
2594 0 : flopcand = 3*ftbasegetflopestimate(m, _state)+6*m;
2595 0 : if( ae_fp_less(flopcand,flopbest) )
2596 : {
2597 0 : algbest = 1;
2598 0 : flopbest = flopcand;
2599 : }
2600 : }
2601 : else
2602 : {
2603 :
2604 : /*
2605 : * general cyclic/non-cyclic convolution
2606 : */
2607 0 : p = ftbasefindsmooth(m+n-1, _state);
2608 0 : flopcand = 3*ftbasegetflopestimate(p, _state)+6*p;
2609 0 : if( ae_fp_less(flopcand,flopbest) )
2610 : {
2611 0 : algbest = 1;
2612 0 : flopbest = flopcand;
2613 : }
2614 : }
2615 : }
2616 :
2617 : /*
2618 : * Another candidate - overlap-add
2619 : */
2620 0 : q = 1;
2621 0 : ptotal = 1;
2622 0 : while(ptotal<n)
2623 : {
2624 0 : ptotal = ptotal*2;
2625 : }
2626 0 : while(ptotal<=m+n-1)
2627 : {
2628 0 : p = ptotal-n+1;
2629 0 : flopcand = ae_iceil((double)m/(double)p, _state)*(2*ftbasegetflopestimate(ptotal, _state)+8*ptotal);
2630 0 : if( ae_fp_less(flopcand,flopbest) )
2631 : {
2632 0 : flopbest = flopcand;
2633 0 : algbest = 2;
2634 0 : q = p;
2635 : }
2636 0 : ptotal = ptotal*2;
2637 : }
2638 0 : alg = algbest;
2639 0 : convc1dx(a, m, b, n, circular, alg, q, r, _state);
2640 0 : ae_frame_leave(_state);
2641 0 : 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 0 : if( alg==0 )
2651 : {
2652 :
2653 : /*
2654 : * Special case: N=1
2655 : */
2656 0 : if( n==1 )
2657 : {
2658 0 : ae_vector_set_length(r, m, _state);
2659 0 : v = b->ptr.p_complex[0];
2660 0 : ae_v_cmovec(&r->ptr.p_complex[0], 1, &a->ptr.p_complex[0], 1, "N", ae_v_len(0,m-1), v);
2661 0 : ae_frame_leave(_state);
2662 0 : return;
2663 : }
2664 :
2665 : /*
2666 : * use straightforward formula
2667 : */
2668 0 : if( circular )
2669 : {
2670 :
2671 : /*
2672 : * circular convolution
2673 : */
2674 0 : ae_vector_set_length(r, m, _state);
2675 0 : v = b->ptr.p_complex[0];
2676 0 : ae_v_cmovec(&r->ptr.p_complex[0], 1, &a->ptr.p_complex[0], 1, "N", ae_v_len(0,m-1), v);
2677 0 : for(i=1; i<=n-1; i++)
2678 : {
2679 0 : v = b->ptr.p_complex[i];
2680 0 : i1 = 0;
2681 0 : i2 = i-1;
2682 0 : j1 = m-i;
2683 0 : j2 = m-1;
2684 0 : ae_v_caddc(&r->ptr.p_complex[i1], 1, &a->ptr.p_complex[j1], 1, "N", ae_v_len(i1,i2), v);
2685 0 : i1 = i;
2686 0 : i2 = m-1;
2687 0 : j1 = 0;
2688 0 : j2 = m-i-1;
2689 0 : 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 0 : ae_vector_set_length(r, m+n-1, _state);
2699 0 : for(i=0; i<=m+n-2; i++)
2700 : {
2701 0 : r->ptr.p_complex[i] = ae_complex_from_i(0);
2702 : }
2703 0 : for(i=0; i<=n-1; i++)
2704 : {
2705 0 : v = b->ptr.p_complex[i];
2706 0 : 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 0 : ae_frame_leave(_state);
2710 0 : 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 0 : if( alg==1 )
2724 : {
2725 0 : if( circular&&ftbaseissmooth(m, _state) )
2726 : {
2727 :
2728 : /*
2729 : * special code for circular convolution with smooth M
2730 : */
2731 0 : ftcomplexfftplan(m, 1, &plan, _state);
2732 0 : ae_vector_set_length(&buf, 2*m, _state);
2733 0 : for(i=0; i<=m-1; i++)
2734 : {
2735 0 : buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
2736 0 : buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
2737 : }
2738 0 : ae_vector_set_length(&buf2, 2*m, _state);
2739 0 : for(i=0; i<=n-1; i++)
2740 : {
2741 0 : buf2.ptr.p_double[2*i+0] = b->ptr.p_complex[i].x;
2742 0 : buf2.ptr.p_double[2*i+1] = b->ptr.p_complex[i].y;
2743 : }
2744 0 : for(i=n; i<=m-1; i++)
2745 : {
2746 0 : buf2.ptr.p_double[2*i+0] = (double)(0);
2747 0 : buf2.ptr.p_double[2*i+1] = (double)(0);
2748 : }
2749 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
2750 0 : ftapplyplan(&plan, &buf2, 0, 1, _state);
2751 0 : for(i=0; i<=m-1; i++)
2752 : {
2753 0 : ax = buf.ptr.p_double[2*i+0];
2754 0 : ay = buf.ptr.p_double[2*i+1];
2755 0 : bx = buf2.ptr.p_double[2*i+0];
2756 0 : by = buf2.ptr.p_double[2*i+1];
2757 0 : tx = ax*bx-ay*by;
2758 0 : ty = ax*by+ay*bx;
2759 0 : buf.ptr.p_double[2*i+0] = tx;
2760 0 : buf.ptr.p_double[2*i+1] = -ty;
2761 : }
2762 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
2763 0 : t = (double)1/(double)m;
2764 0 : ae_vector_set_length(r, m, _state);
2765 0 : for(i=0; i<=m-1; i++)
2766 : {
2767 0 : r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
2768 0 : 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 0 : p = ftbasefindsmooth(m+n-1, _state);
2783 0 : ftcomplexfftplan(p, 1, &plan, _state);
2784 0 : ae_vector_set_length(&buf, 2*p, _state);
2785 0 : for(i=0; i<=m-1; i++)
2786 : {
2787 0 : buf.ptr.p_double[2*i+0] = a->ptr.p_complex[i].x;
2788 0 : buf.ptr.p_double[2*i+1] = a->ptr.p_complex[i].y;
2789 : }
2790 0 : for(i=m; i<=p-1; i++)
2791 : {
2792 0 : buf.ptr.p_double[2*i+0] = (double)(0);
2793 0 : buf.ptr.p_double[2*i+1] = (double)(0);
2794 : }
2795 0 : ae_vector_set_length(&buf2, 2*p, _state);
2796 0 : for(i=0; i<=n-1; i++)
2797 : {
2798 0 : buf2.ptr.p_double[2*i+0] = b->ptr.p_complex[i].x;
2799 0 : buf2.ptr.p_double[2*i+1] = b->ptr.p_complex[i].y;
2800 : }
2801 0 : for(i=n; i<=p-1; i++)
2802 : {
2803 0 : buf2.ptr.p_double[2*i+0] = (double)(0);
2804 0 : buf2.ptr.p_double[2*i+1] = (double)(0);
2805 : }
2806 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
2807 0 : ftapplyplan(&plan, &buf2, 0, 1, _state);
2808 0 : for(i=0; i<=p-1; i++)
2809 : {
2810 0 : ax = buf.ptr.p_double[2*i+0];
2811 0 : ay = buf.ptr.p_double[2*i+1];
2812 0 : bx = buf2.ptr.p_double[2*i+0];
2813 0 : by = buf2.ptr.p_double[2*i+1];
2814 0 : tx = ax*bx-ay*by;
2815 0 : ty = ax*by+ay*bx;
2816 0 : buf.ptr.p_double[2*i+0] = tx;
2817 0 : buf.ptr.p_double[2*i+1] = -ty;
2818 : }
2819 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
2820 0 : t = (double)1/(double)p;
2821 0 : if( circular )
2822 : {
2823 :
2824 : /*
2825 : * circular, add tail to head
2826 : */
2827 0 : ae_vector_set_length(r, m, _state);
2828 0 : for(i=0; i<=m-1; i++)
2829 : {
2830 0 : r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
2831 0 : r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
2832 : }
2833 0 : for(i=m; i<=m+n-2; i++)
2834 : {
2835 0 : r->ptr.p_complex[i-m].x = r->ptr.p_complex[i-m].x+t*buf.ptr.p_double[2*i+0];
2836 0 : 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 0 : ae_vector_set_length(r, m+n-1, _state);
2846 0 : for(i=0; i<=m+n-2; i++)
2847 : {
2848 0 : r->ptr.p_complex[i].x = t*buf.ptr.p_double[2*i+0];
2849 0 : r->ptr.p_complex[i].y = -t*buf.ptr.p_double[2*i+1];
2850 : }
2851 : }
2852 : }
2853 0 : ae_frame_leave(_state);
2854 0 : 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 0 : if( alg==2 )
2868 : {
2869 0 : ae_vector_set_length(&buf, 2*(q+n-1), _state);
2870 :
2871 : /*
2872 : * prepare R
2873 : */
2874 0 : if( circular )
2875 : {
2876 0 : ae_vector_set_length(r, m, _state);
2877 0 : for(i=0; i<=m-1; i++)
2878 : {
2879 0 : r->ptr.p_complex[i] = ae_complex_from_i(0);
2880 : }
2881 : }
2882 : else
2883 : {
2884 0 : ae_vector_set_length(r, m+n-1, _state);
2885 0 : for(i=0; i<=m+n-2; i++)
2886 : {
2887 0 : r->ptr.p_complex[i] = ae_complex_from_i(0);
2888 : }
2889 : }
2890 :
2891 : /*
2892 : * pre-calculated FFT(B)
2893 : */
2894 0 : ae_vector_set_length(&bbuf, q+n-1, _state);
2895 0 : ae_v_cmove(&bbuf.ptr.p_complex[0], 1, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
2896 0 : for(j=n; j<=q+n-2; j++)
2897 : {
2898 0 : bbuf.ptr.p_complex[j] = ae_complex_from_i(0);
2899 : }
2900 0 : fftc1d(&bbuf, q+n-1, _state);
2901 :
2902 : /*
2903 : * prepare FFT plan for chunks of A
2904 : */
2905 0 : ftcomplexfftplan(q+n-1, 1, &plan, _state);
2906 :
2907 : /*
2908 : * main overlap-add cycle
2909 : */
2910 0 : i = 0;
2911 0 : while(i<=m-1)
2912 : {
2913 0 : p = ae_minint(q, m-i, _state);
2914 0 : for(j=0; j<=p-1; j++)
2915 : {
2916 0 : buf.ptr.p_double[2*j+0] = a->ptr.p_complex[i+j].x;
2917 0 : buf.ptr.p_double[2*j+1] = a->ptr.p_complex[i+j].y;
2918 : }
2919 0 : for(j=p; j<=q+n-2; j++)
2920 : {
2921 0 : buf.ptr.p_double[2*j+0] = (double)(0);
2922 0 : buf.ptr.p_double[2*j+1] = (double)(0);
2923 : }
2924 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
2925 0 : for(j=0; j<=q+n-2; j++)
2926 : {
2927 0 : ax = buf.ptr.p_double[2*j+0];
2928 0 : ay = buf.ptr.p_double[2*j+1];
2929 0 : bx = bbuf.ptr.p_complex[j].x;
2930 0 : by = bbuf.ptr.p_complex[j].y;
2931 0 : tx = ax*bx-ay*by;
2932 0 : ty = ax*by+ay*bx;
2933 0 : buf.ptr.p_double[2*j+0] = tx;
2934 0 : buf.ptr.p_double[2*j+1] = -ty;
2935 : }
2936 0 : ftapplyplan(&plan, &buf, 0, 1, _state);
2937 0 : t = (double)1/(double)(q+n-1);
2938 0 : if( circular )
2939 : {
2940 0 : j1 = ae_minint(i+p+n-2, m-1, _state)-i;
2941 0 : j2 = j1+1;
2942 : }
2943 : else
2944 : {
2945 0 : j1 = p+n-2;
2946 0 : j2 = j1+1;
2947 : }
2948 0 : for(j=0; j<=j1; j++)
2949 : {
2950 0 : r->ptr.p_complex[i+j].x = r->ptr.p_complex[i+j].x+buf.ptr.p_double[2*j+0]*t;
2951 0 : r->ptr.p_complex[i+j].y = r->ptr.p_complex[i+j].y-buf.ptr.p_double[2*j+1]*t;
2952 : }
2953 0 : for(j=j2; j<=p+n-2; j++)
2954 : {
2955 0 : r->ptr.p_complex[j-j2].x = r->ptr.p_complex[j-j2].x+buf.ptr.p_double[2*j+0]*t;
2956 0 : r->ptr.p_complex[j-j2].y = r->ptr.p_complex[j-j2].y-buf.ptr.p_double[2*j+1]*t;
2957 : }
2958 0 : i = i+p;
2959 : }
2960 0 : ae_frame_leave(_state);
2961 0 : return;
2962 : }
2963 0 : 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 : *************************************************************************/
2992 0 : 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 0 : ae_frame_make(_state, &_frame_block);
3027 0 : memset(&plan, 0, sizeof(plan));
3028 0 : memset(&buf, 0, sizeof(buf));
3029 0 : memset(&buf2, 0, sizeof(buf2));
3030 0 : memset(&buf3, 0, sizeof(buf3));
3031 0 : ae_vector_clear(r);
3032 0 : _fasttransformplan_init(&plan, _state, ae_true);
3033 0 : ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
3034 0 : ae_vector_init(&buf2, 0, DT_REAL, _state, ae_true);
3035 0 : ae_vector_init(&buf3, 0, DT_REAL, _state, ae_true);
3036 :
3037 0 : ae_assert(n>0&&m>0, "ConvC1DX: incorrect N or M!", _state);
3038 0 : ae_assert(n<=m, "ConvC1DX: N<M assumption is false!", _state);
3039 :
3040 : /*
3041 : * handle special cases
3042 : */
3043 0 : if( ae_minint(m, n, _state)<=2 )
3044 : {
3045 0 : alg = 0;
3046 : }
3047 :
3048 : /*
3049 : * Auto-select
3050 : */
3051 0 : 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 0 : algbest = 0;
3062 0 : if( alg==-1 )
3063 : {
3064 0 : flopbest = 0.15*m*n;
3065 : }
3066 : else
3067 : {
3068 0 : flopbest = ae_maxrealnumber;
3069 : }
3070 :
3071 : /*
3072 : * Another candidate - generic FFT code
3073 : */
3074 0 : if( alg==-1 )
3075 : {
3076 0 : 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 0 : flopcand = 3*ftbasegetflopestimate(m/2, _state)+(double)(6*m)/(double)2;
3083 0 : if( ae_fp_less(flopcand,flopbest) )
3084 : {
3085 0 : algbest = 1;
3086 0 : flopbest = flopcand;
3087 : }
3088 : }
3089 : else
3090 : {
3091 :
3092 : /*
3093 : * general cyclic/non-cyclic convolution
3094 : */
3095 0 : p = ftbasefindsmootheven(m+n-1, _state);
3096 0 : flopcand = 3*ftbasegetflopestimate(p/2, _state)+(double)(6*p)/(double)2;
3097 0 : if( ae_fp_less(flopcand,flopbest) )
3098 : {
3099 0 : algbest = 1;
3100 0 : flopbest = flopcand;
3101 : }
3102 : }
3103 : }
3104 :
3105 : /*
3106 : * Another candidate - overlap-add
3107 : */
3108 0 : q = 1;
3109 0 : ptotal = 1;
3110 0 : while(ptotal<n)
3111 : {
3112 0 : ptotal = ptotal*2;
3113 : }
3114 0 : while(ptotal<=m+n-1)
3115 : {
3116 0 : p = ptotal-n+1;
3117 0 : flopcand = ae_iceil((double)m/(double)p, _state)*(2*ftbasegetflopestimate(ptotal/2, _state)+1*(ptotal/2));
3118 0 : if( ae_fp_less(flopcand,flopbest) )
3119 : {
3120 0 : flopbest = flopcand;
3121 0 : algbest = 2;
3122 0 : q = p;
3123 : }
3124 0 : ptotal = ptotal*2;
3125 : }
3126 0 : alg = algbest;
3127 0 : convr1dx(a, m, b, n, circular, alg, q, r, _state);
3128 0 : ae_frame_leave(_state);
3129 0 : 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 0 : if( alg==0 )
3139 : {
3140 :
3141 : /*
3142 : * Special case: N=1
3143 : */
3144 0 : if( n==1 )
3145 : {
3146 0 : ae_vector_set_length(r, m, _state);
3147 0 : v = b->ptr.p_double[0];
3148 0 : ae_v_moved(&r->ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1), v);
3149 0 : ae_frame_leave(_state);
3150 0 : return;
3151 : }
3152 :
3153 : /*
3154 : * use straightforward formula
3155 : */
3156 0 : if( circular )
3157 : {
3158 :
3159 : /*
3160 : * circular convolution
3161 : */
3162 0 : ae_vector_set_length(r, m, _state);
3163 0 : v = b->ptr.p_double[0];
3164 0 : ae_v_moved(&r->ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1), v);
3165 0 : for(i=1; i<=n-1; i++)
3166 : {
3167 0 : v = b->ptr.p_double[i];
3168 0 : i1 = 0;
3169 0 : i2 = i-1;
3170 0 : j1 = m-i;
3171 0 : j2 = m-1;
3172 0 : ae_v_addd(&r->ptr.p_double[i1], 1, &a->ptr.p_double[j1], 1, ae_v_len(i1,i2), v);
3173 0 : i1 = i;
3174 0 : i2 = m-1;
3175 0 : j1 = 0;
3176 0 : j2 = m-i-1;
3177 0 : 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 0 : ae_vector_set_length(r, m+n-1, _state);
3187 0 : for(i=0; i<=m+n-2; i++)
3188 : {
3189 0 : r->ptr.p_double[i] = (double)(0);
3190 : }
3191 0 : for(i=0; i<=n-1; i++)
3192 : {
3193 0 : v = b->ptr.p_double[i];
3194 0 : 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 0 : ae_frame_leave(_state);
3198 0 : 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 0 : if( alg==1 )
3214 : {
3215 0 : ae_assert(m+n-1>2, "ConvR1DX: internal error!", _state);
3216 0 : if( (circular&&ftbaseissmooth(m, _state))&&m%2==0 )
3217 : {
3218 :
3219 : /*
3220 : * special code for circular convolution with smooth even M
3221 : */
3222 0 : ae_vector_set_length(&buf, m, _state);
3223 0 : ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1));
3224 0 : ae_vector_set_length(&buf2, m, _state);
3225 0 : ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
3226 0 : for(i=n; i<=m-1; i++)
3227 : {
3228 0 : buf2.ptr.p_double[i] = (double)(0);
3229 : }
3230 0 : ae_vector_set_length(&buf3, m, _state);
3231 0 : ftcomplexfftplan(m/2, 1, &plan, _state);
3232 0 : fftr1dinternaleven(&buf, m, &buf3, &plan, _state);
3233 0 : fftr1dinternaleven(&buf2, m, &buf3, &plan, _state);
3234 0 : buf.ptr.p_double[0] = buf.ptr.p_double[0]*buf2.ptr.p_double[0];
3235 0 : buf.ptr.p_double[1] = buf.ptr.p_double[1]*buf2.ptr.p_double[1];
3236 0 : for(i=1; i<=m/2-1; i++)
3237 : {
3238 0 : ax = buf.ptr.p_double[2*i+0];
3239 0 : ay = buf.ptr.p_double[2*i+1];
3240 0 : bx = buf2.ptr.p_double[2*i+0];
3241 0 : by = buf2.ptr.p_double[2*i+1];
3242 0 : tx = ax*bx-ay*by;
3243 0 : ty = ax*by+ay*bx;
3244 0 : buf.ptr.p_double[2*i+0] = tx;
3245 0 : buf.ptr.p_double[2*i+1] = ty;
3246 : }
3247 0 : fftr1dinvinternaleven(&buf, m, &buf3, &plan, _state);
3248 0 : ae_vector_set_length(r, m, _state);
3249 0 : 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 0 : p = ftbasefindsmootheven(m+n-1, _state);
3263 0 : ae_vector_set_length(&buf, p, _state);
3264 0 : ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[0], 1, ae_v_len(0,m-1));
3265 0 : for(i=m; i<=p-1; i++)
3266 : {
3267 0 : buf.ptr.p_double[i] = (double)(0);
3268 : }
3269 0 : ae_vector_set_length(&buf2, p, _state);
3270 0 : ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
3271 0 : for(i=n; i<=p-1; i++)
3272 : {
3273 0 : buf2.ptr.p_double[i] = (double)(0);
3274 : }
3275 0 : ae_vector_set_length(&buf3, p, _state);
3276 0 : ftcomplexfftplan(p/2, 1, &plan, _state);
3277 0 : fftr1dinternaleven(&buf, p, &buf3, &plan, _state);
3278 0 : fftr1dinternaleven(&buf2, p, &buf3, &plan, _state);
3279 0 : buf.ptr.p_double[0] = buf.ptr.p_double[0]*buf2.ptr.p_double[0];
3280 0 : buf.ptr.p_double[1] = buf.ptr.p_double[1]*buf2.ptr.p_double[1];
3281 0 : for(i=1; i<=p/2-1; i++)
3282 : {
3283 0 : ax = buf.ptr.p_double[2*i+0];
3284 0 : ay = buf.ptr.p_double[2*i+1];
3285 0 : bx = buf2.ptr.p_double[2*i+0];
3286 0 : by = buf2.ptr.p_double[2*i+1];
3287 0 : tx = ax*bx-ay*by;
3288 0 : ty = ax*by+ay*bx;
3289 0 : buf.ptr.p_double[2*i+0] = tx;
3290 0 : buf.ptr.p_double[2*i+1] = ty;
3291 : }
3292 0 : fftr1dinvinternaleven(&buf, p, &buf3, &plan, _state);
3293 0 : if( circular )
3294 : {
3295 :
3296 : /*
3297 : * circular, add tail to head
3298 : */
3299 0 : ae_vector_set_length(r, m, _state);
3300 0 : ae_v_move(&r->ptr.p_double[0], 1, &buf.ptr.p_double[0], 1, ae_v_len(0,m-1));
3301 0 : if( n>=2 )
3302 : {
3303 0 : 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 0 : ae_vector_set_length(r, m+n-1, _state);
3313 0 : 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 0 : ae_frame_leave(_state);
3317 0 : return;
3318 : }
3319 :
3320 : /*
3321 : * overlap-add method
3322 : */
3323 0 : if( alg==2 )
3324 : {
3325 0 : ae_assert((q+n-1)%2==0, "ConvR1DX: internal error!", _state);
3326 0 : ae_vector_set_length(&buf, q+n-1, _state);
3327 0 : ae_vector_set_length(&buf2, q+n-1, _state);
3328 0 : ae_vector_set_length(&buf3, q+n-1, _state);
3329 0 : ftcomplexfftplan((q+n-1)/2, 1, &plan, _state);
3330 :
3331 : /*
3332 : * prepare R
3333 : */
3334 0 : if( circular )
3335 : {
3336 0 : ae_vector_set_length(r, m, _state);
3337 0 : for(i=0; i<=m-1; i++)
3338 : {
3339 0 : r->ptr.p_double[i] = (double)(0);
3340 : }
3341 : }
3342 : else
3343 : {
3344 0 : ae_vector_set_length(r, m+n-1, _state);
3345 0 : for(i=0; i<=m+n-2; i++)
3346 : {
3347 0 : r->ptr.p_double[i] = (double)(0);
3348 : }
3349 : }
3350 :
3351 : /*
3352 : * pre-calculated FFT(B)
3353 : */
3354 0 : ae_v_move(&buf2.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
3355 0 : for(j=n; j<=q+n-2; j++)
3356 : {
3357 0 : buf2.ptr.p_double[j] = (double)(0);
3358 : }
3359 0 : fftr1dinternaleven(&buf2, q+n-1, &buf3, &plan, _state);
3360 :
3361 : /*
3362 : * main overlap-add cycle
3363 : */
3364 0 : i = 0;
3365 0 : while(i<=m-1)
3366 : {
3367 0 : p = ae_minint(q, m-i, _state);
3368 0 : ae_v_move(&buf.ptr.p_double[0], 1, &a->ptr.p_double[i], 1, ae_v_len(0,p-1));
3369 0 : for(j=p; j<=q+n-2; j++)
3370 : {
3371 0 : buf.ptr.p_double[j] = (double)(0);
3372 : }
3373 0 : fftr1dinternaleven(&buf, q+n-1, &buf3, &plan, _state);
3374 0 : buf.ptr.p_double[0] = buf.ptr.p_double[0]*buf2.ptr.p_double[0];
3375 0 : buf.ptr.p_double[1] = buf.ptr.p_double[1]*buf2.ptr.p_double[1];
3376 0 : for(j=1; j<=(q+n-1)/2-1; j++)
3377 : {
3378 0 : ax = buf.ptr.p_double[2*j+0];
3379 0 : ay = buf.ptr.p_double[2*j+1];
3380 0 : bx = buf2.ptr.p_double[2*j+0];
3381 0 : by = buf2.ptr.p_double[2*j+1];
3382 0 : tx = ax*bx-ay*by;
3383 0 : ty = ax*by+ay*bx;
3384 0 : buf.ptr.p_double[2*j+0] = tx;
3385 0 : buf.ptr.p_double[2*j+1] = ty;
3386 : }
3387 0 : fftr1dinvinternaleven(&buf, q+n-1, &buf3, &plan, _state);
3388 0 : if( circular )
3389 : {
3390 0 : j1 = ae_minint(i+p+n-2, m-1, _state)-i;
3391 0 : j2 = j1+1;
3392 : }
3393 : else
3394 : {
3395 0 : j1 = p+n-2;
3396 0 : j2 = j1+1;
3397 : }
3398 0 : ae_v_add(&r->ptr.p_double[i], 1, &buf.ptr.p_double[0], 1, ae_v_len(i,i+j1));
3399 0 : if( p+n-2>=j2 )
3400 : {
3401 0 : 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 0 : i = i+p;
3404 : }
3405 0 : ae_frame_leave(_state);
3406 0 : return;
3407 : }
3408 0 : 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 : *************************************************************************/
3452 0 : 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 0 : ae_frame_make(_state, &_frame_block);
3465 0 : memset(&p, 0, sizeof(p));
3466 0 : memset(&b, 0, sizeof(b));
3467 0 : ae_vector_clear(r);
3468 0 : ae_vector_init(&p, 0, DT_COMPLEX, _state, ae_true);
3469 0 : ae_vector_init(&b, 0, DT_COMPLEX, _state, ae_true);
3470 :
3471 0 : ae_assert(n>0&&m>0, "CorrC1D: incorrect N or M!", _state);
3472 0 : ae_vector_set_length(&p, m, _state);
3473 0 : for(i=0; i<=m-1; i++)
3474 : {
3475 0 : p.ptr.p_complex[m-1-i] = ae_c_conj(pattern->ptr.p_complex[i], _state);
3476 : }
3477 0 : convc1d(&p, m, signal, n, &b, _state);
3478 0 : ae_vector_set_length(r, m+n-1, _state);
3479 0 : ae_v_cmove(&r->ptr.p_complex[0], 1, &b.ptr.p_complex[m-1], 1, "N", ae_v_len(0,n-1));
3480 0 : if( m+n-2>=n )
3481 : {
3482 0 : 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 0 : ae_frame_leave(_state);
3485 0 : }
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 : *************************************************************************/
3515 0 : 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 0 : ae_frame_make(_state, &_frame_block);
3531 0 : memset(&p, 0, sizeof(p));
3532 0 : memset(&b, 0, sizeof(b));
3533 0 : ae_vector_clear(c);
3534 0 : ae_vector_init(&p, 0, DT_COMPLEX, _state, ae_true);
3535 0 : ae_vector_init(&b, 0, DT_COMPLEX, _state, ae_true);
3536 :
3537 0 : 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 0 : if( m<n )
3544 : {
3545 0 : ae_vector_set_length(&b, m, _state);
3546 0 : for(i1=0; i1<=m-1; i1++)
3547 : {
3548 0 : b.ptr.p_complex[i1] = ae_complex_from_i(0);
3549 : }
3550 0 : i1 = 0;
3551 0 : while(i1<n)
3552 : {
3553 0 : i2 = ae_minint(i1+m-1, n-1, _state);
3554 0 : j2 = i2-i1;
3555 0 : ae_v_cadd(&b.ptr.p_complex[0], 1, &pattern->ptr.p_complex[i1], 1, "N", ae_v_len(0,j2));
3556 0 : i1 = i1+m;
3557 : }
3558 0 : corrc1dcircular(signal, m, &b, m, c, _state);
3559 0 : ae_frame_leave(_state);
3560 0 : return;
3561 : }
3562 :
3563 : /*
3564 : * Task is normalized
3565 : */
3566 0 : ae_vector_set_length(&p, n, _state);
3567 0 : for(i=0; i<=n-1; i++)
3568 : {
3569 0 : p.ptr.p_complex[n-1-i] = ae_c_conj(pattern->ptr.p_complex[i], _state);
3570 : }
3571 0 : convc1dcircular(signal, m, &p, n, &b, _state);
3572 0 : ae_vector_set_length(c, m, _state);
3573 0 : ae_v_cmove(&c->ptr.p_complex[0], 1, &b.ptr.p_complex[n-1], 1, "N", ae_v_len(0,m-n));
3574 0 : if( m-n+1<=m-1 )
3575 : {
3576 0 : 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 0 : 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 : *************************************************************************/
3618 0 : 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 0 : ae_frame_make(_state, &_frame_block);
3631 0 : memset(&p, 0, sizeof(p));
3632 0 : memset(&b, 0, sizeof(b));
3633 0 : ae_vector_clear(r);
3634 0 : ae_vector_init(&p, 0, DT_REAL, _state, ae_true);
3635 0 : ae_vector_init(&b, 0, DT_REAL, _state, ae_true);
3636 :
3637 0 : ae_assert(n>0&&m>0, "CorrR1D: incorrect N or M!", _state);
3638 0 : ae_vector_set_length(&p, m, _state);
3639 0 : for(i=0; i<=m-1; i++)
3640 : {
3641 0 : p.ptr.p_double[m-1-i] = pattern->ptr.p_double[i];
3642 : }
3643 0 : convr1d(&p, m, signal, n, &b, _state);
3644 0 : ae_vector_set_length(r, m+n-1, _state);
3645 0 : ae_v_move(&r->ptr.p_double[0], 1, &b.ptr.p_double[m-1], 1, ae_v_len(0,n-1));
3646 0 : if( m+n-2>=n )
3647 : {
3648 0 : ae_v_move(&r->ptr.p_double[n], 1, &b.ptr.p_double[0], 1, ae_v_len(n,m+n-2));
3649 : }
3650 0 : ae_frame_leave(_state);
3651 0 : }
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 : *************************************************************************/
3681 0 : 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 0 : ae_frame_make(_state, &_frame_block);
3697 0 : memset(&p, 0, sizeof(p));
3698 0 : memset(&b, 0, sizeof(b));
3699 0 : ae_vector_clear(c);
3700 0 : ae_vector_init(&p, 0, DT_REAL, _state, ae_true);
3701 0 : ae_vector_init(&b, 0, DT_REAL, _state, ae_true);
3702 :
3703 0 : 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 0 : if( m<n )
3710 : {
3711 0 : ae_vector_set_length(&b, m, _state);
3712 0 : for(i1=0; i1<=m-1; i1++)
3713 : {
3714 0 : b.ptr.p_double[i1] = (double)(0);
3715 : }
3716 0 : i1 = 0;
3717 0 : while(i1<n)
3718 : {
3719 0 : i2 = ae_minint(i1+m-1, n-1, _state);
3720 0 : j2 = i2-i1;
3721 0 : ae_v_add(&b.ptr.p_double[0], 1, &pattern->ptr.p_double[i1], 1, ae_v_len(0,j2));
3722 0 : i1 = i1+m;
3723 : }
3724 0 : corrr1dcircular(signal, m, &b, m, c, _state);
3725 0 : ae_frame_leave(_state);
3726 0 : return;
3727 : }
3728 :
3729 : /*
3730 : * Task is normalized
3731 : */
3732 0 : ae_vector_set_length(&p, n, _state);
3733 0 : for(i=0; i<=n-1; i++)
3734 : {
3735 0 : p.ptr.p_double[n-1-i] = pattern->ptr.p_double[i];
3736 : }
3737 0 : convr1dcircular(signal, m, &p, n, &b, _state);
3738 0 : ae_vector_set_length(c, m, _state);
3739 0 : ae_v_move(&c->ptr.p_double[0], 1, &b.ptr.p_double[n-1], 1, ae_v_len(0,m-n));
3740 0 : if( m-n+1<=m-1 )
3741 : {
3742 0 : 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 0 : ae_frame_leave(_state);
3745 : }
3746 :
3747 :
3748 : #endif
3749 :
3750 : }
3751 :
|
