Free cookie consent management tool by TermsFeed Policy Generator

source: stable/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.7.0/ALGLIB-3.7.0/fasttransforms.cs @ 16204

Last change on this file since 16204 was 9268, checked in by mkommend, 12 years ago

#2007: Added AlgLib 3.7.0 to the external libraries.

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