Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

fasttransforms.cs @ 17931

Visit:

Last change on this file since 17931 was 17931, checked in by gkronber, 3 years ago
#3117: update alglib to version 3.17
File size: 119.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.17.0/ALGLIB-3.17.0/fasttransforms.cs @ 17931

Download in other formats: