Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/fft.cs @ 3839

Last change on this file since 3839 was 3839, checked in by mkommend, 14 years ago

implemented first version of LR (ticket #1012)

File size: 16.3 KB
Line 
1/*************************************************************************
2Copyright (c) 2009, 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
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class fft
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 AP.Complex[] a,
57            int n)
58        {
59            ftbase.ftplan plan = new ftbase.ftplan();
60            int i = 0;
61            double[] buf = new double[0];
62
63            System.Diagnostics.Debug.Assert(n>0, "FFTC1D: incorrect N!");
64           
65            //
66            // Special case: N=1, FFT is just identity transform.
67            // After this block we assume that N is strictly greater than 1.
68            //
69            if( n==1 )
70            {
71                return;
72            }
73           
74            //
75            // convert input array to the more convinient format
76            //
77            buf = new double[2*n];
78            for(i=0; i<=n-1; i++)
79            {
80                buf[2*i+0] = a[i].x;
81                buf[2*i+1] = a[i].y;
82            }
83           
84            //
85            // Generate plan and execute it.
86            //
87            // Plan is a combination of a successive factorizations of N and
88            // precomputed data. It is much like a FFTW plan, but is not stored
89            // between subroutine calls and is much simpler.
90            //
91            ftbase.ftbasegeneratecomplexfftplan(n, ref plan);
92            ftbase.ftbaseexecuteplan(ref buf, 0, n, ref plan);
93           
94            //
95            // result
96            //
97            for(i=0; i<=n-1; i++)
98            {
99                a[i].x = buf[2*i+0];
100                a[i].y = buf[2*i+1];
101            }
102        }
103
104
105        /*************************************************************************
106        1-dimensional complex inverse FFT.
107
108        Array size N may be arbitrary number (composite or prime).  Algorithm  has
109        O(N*logN) complexity for any N (composite or prime).
110
111        See FFTC1D() description for more information about algorithm performance.
112
113        INPUT PARAMETERS
114            A   -   array[0..N-1] - complex array to be transformed
115            N   -   problem size
116
117        OUTPUT PARAMETERS
118            A   -   inverse DFT of a input array, array[0..N-1]
119                    A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
120
121
122          -- ALGLIB --
123             Copyright 29.05.2009 by Bochkanov Sergey
124        *************************************************************************/
125        public static void fftc1dinv(ref AP.Complex[] a,
126            int n)
127        {
128            int i = 0;
129
130            System.Diagnostics.Debug.Assert(n>0, "FFTC1DInv: incorrect N!");
131           
132            //
133            // Inverse DFT can be expressed in terms of the DFT as
134            //
135            //     invfft(x) = fft(x')'/N
136            //
137            // here x' means conj(x).
138            //
139            for(i=0; i<=n-1; i++)
140            {
141                a[i].y = -a[i].y;
142            }
143            fftc1d(ref a, n);
144            for(i=0; i<=n-1; i++)
145            {
146                a[i].x = a[i].x/n;
147                a[i].y = -(a[i].y/n);
148            }
149        }
150
151
152        /*************************************************************************
153        1-dimensional real FFT.
154
155        Algorithm has O(N*logN) complexity for any N (composite or prime).
156
157        INPUT PARAMETERS
158            A   -   array[0..N-1] - real function to be transformed
159            N   -   problem size
160
161        OUTPUT PARAMETERS
162            F   -   DFT of a input array, array[0..N-1]
163                    F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)
164
165        NOTE:
166            F[] satisfies symmetry property F[k] = conj(F[N-k]),  so just one half
167        of  array  is  usually needed. But for convinience subroutine returns full
168        complex array (with frequencies above N/2), so its result may be  used  by
169        other FFT-related subroutines.
170
171
172          -- ALGLIB --
173             Copyright 01.06.2009 by Bochkanov Sergey
174        *************************************************************************/
175        public static void fftr1d(ref double[] a,
176            int n,
177            ref AP.Complex[] f)
178        {
179            int i = 0;
180            int n2 = 0;
181            int idx = 0;
182            AP.Complex hn = 0;
183            AP.Complex hmnc = 0;
184            AP.Complex v = 0;
185            double[] buf = new double[0];
186            ftbase.ftplan plan = new ftbase.ftplan();
187            int i_ = 0;
188
189            System.Diagnostics.Debug.Assert(n>0, "FFTR1D: incorrect N!");
190           
191            //
192            // Special cases:
193            // * N=1, FFT is just identity transform.
194            // * N=2, FFT is simple too
195            //
196            // After this block we assume that N is strictly greater than 2
197            //
198            if( n==1 )
199            {
200                f = new AP.Complex[1];
201                f[0] = a[0];
202                return;
203            }
204            if( n==2 )
205            {
206                f = new AP.Complex[2];
207                f[0].x = a[0]+a[1];
208                f[0].y = 0;
209                f[1].x = a[0]-a[1];
210                f[1].y = 0;
211                return;
212            }
213           
214            //
215            // Choose between odd-size and even-size FFTs
216            //
217            if( n%2==0 )
218            {
219               
220                //
221                // even-size real FFT, use reduction to the complex task
222                //
223                n2 = n/2;
224                buf = new double[n];
225                for(i_=0; i_<=n-1;i_++)
226                {
227                    buf[i_] = a[i_];
228                }
229                ftbase.ftbasegeneratecomplexfftplan(n2, ref plan);
230                ftbase.ftbaseexecuteplan(ref buf, 0, n2, ref plan);
231                f = new AP.Complex[n];
232                for(i=0; i<=n2; i++)
233                {
234                    idx = 2*(i%n2);
235                    hn.x = buf[idx+0];
236                    hn.y = buf[idx+1];
237                    idx = 2*((n2-i)%n2);
238                    hmnc.x = buf[idx+0];
239                    hmnc.y = -buf[idx+1];
240                    v.x = -Math.Sin(-(2*Math.PI*i/n));
241                    v.y = Math.Cos(-(2*Math.PI*i/n));
242                    f[i] = hn+hmnc-v*(hn-hmnc);
243                    f[i].x = 0.5*f[i].x;
244                    f[i].y = 0.5*f[i].y;
245                }
246                for(i=n2+1; i<=n-1; i++)
247                {
248                    f[i] = AP.Math.Conj(f[n-i]);
249                }
250                return;
251            }
252            else
253            {
254               
255                //
256                // use complex FFT
257                //
258                f = new AP.Complex[n];
259                for(i=0; i<=n-1; i++)
260                {
261                    f[i] = a[i];
262                }
263                fftc1d(ref f, n);
264                return;
265            }
266        }
267
268
269        /*************************************************************************
270        1-dimensional real inverse FFT.
271
272        Algorithm has O(N*logN) complexity for any N (composite or prime).
273
274        INPUT PARAMETERS
275            F   -   array[0..floor(N/2)] - frequencies from forward real FFT
276            N   -   problem size
277
278        OUTPUT PARAMETERS
279            A   -   inverse DFT of a input array, array[0..N-1]
280
281        NOTE:
282            F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just  one
283        half of frequencies array is needed - elements from 0 to floor(N/2).  F[0]
284        is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd,  then
285        F[floor(N/2)] has no special properties.
286
287        Relying on properties noted above, FFTR1DInv subroutine uses only elements
288        from 0th to floor(N/2)-th. It ignores imaginary part of F[0],  and in case
289        N is even it ignores imaginary part of F[floor(N/2)] too.  So you can pass
290        either frequencies array with N elements or reduced array with roughly N/2
291        elements - subroutine will successfully transform both.
292
293
294          -- ALGLIB --
295             Copyright 01.06.2009 by Bochkanov Sergey
296        *************************************************************************/
297        public static void fftr1dinv(ref AP.Complex[] f,
298            int n,
299            ref double[] a)
300        {
301            int i = 0;
302            double[] h = new double[0];
303            AP.Complex[] fh = new AP.Complex[0];
304
305            System.Diagnostics.Debug.Assert(n>0, "FFTR1DInv: incorrect N!");
306           
307            //
308            // Special case: N=1, FFT is just identity transform.
309            // After this block we assume that N is strictly greater than 1.
310            //
311            if( n==1 )
312            {
313                a = new double[1];
314                a[0] = f[0].x;
315                return;
316            }
317           
318            //
319            // inverse real FFT is reduced to the inverse real FHT,
320            // which is reduced to the forward real FHT,
321            // which is reduced to the forward real FFT.
322            //
323            // Don't worry, it is really compact and efficient reduction :)
324            //
325            h = new double[n];
326            a = new double[n];
327            h[0] = f[0].x;
328            for(i=1; i<=(int)Math.Floor((double)(n)/(double)(2))-1; i++)
329            {
330                h[i] = f[i].x-f[i].y;
331                h[n-i] = f[i].x+f[i].y;
332            }
333            if( n%2==0 )
334            {
335                h[(int)Math.Floor((double)(n)/(double)(2))] = f[(int)Math.Floor((double)(n)/(double)(2))].x;
336            }
337            else
338            {
339                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;
340                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;
341            }
342            fftr1d(ref h, n, ref fh);
343            for(i=0; i<=n-1; i++)
344            {
345                a[i] = (fh[i].x-fh[i].y)/n;
346            }
347        }
348
349
350        /*************************************************************************
351        Internal subroutine. Never call it directly!
352
353
354          -- ALGLIB --
355             Copyright 01.06.2009 by Bochkanov Sergey
356        *************************************************************************/
357        public static void fftr1dinternaleven(ref double[] a,
358            int n,
359            ref double[] buf,
360            ref ftbase.ftplan plan)
361        {
362            double x = 0;
363            double y = 0;
364            int i = 0;
365            int n2 = 0;
366            int idx = 0;
367            AP.Complex hn = 0;
368            AP.Complex hmnc = 0;
369            AP.Complex v = 0;
370            int i_ = 0;
371
372            System.Diagnostics.Debug.Assert(n>0 & n%2==0, "FFTR1DEvenInplace: incorrect N!");
373           
374            //
375            // Special cases:
376            // * N=2
377            //
378            // After this block we assume that N is strictly greater than 2
379            //
380            if( n==2 )
381            {
382                x = a[0]+a[1];
383                y = a[0]-a[1];
384                a[0] = x;
385                a[1] = y;
386                return;
387            }
388           
389            //
390            // even-size real FFT, use reduction to the complex task
391            //
392            n2 = n/2;
393            for(i_=0; i_<=n-1;i_++)
394            {
395                buf[i_] = a[i_];
396            }
397            ftbase.ftbaseexecuteplan(ref buf, 0, n2, ref plan);
398            a[0] = buf[0]+buf[1];
399            for(i=1; i<=n2-1; i++)
400            {
401                idx = 2*(i%n2);
402                hn.x = buf[idx+0];
403                hn.y = buf[idx+1];
404                idx = 2*(n2-i);
405                hmnc.x = buf[idx+0];
406                hmnc.y = -buf[idx+1];
407                v.x = -Math.Sin(-(2*Math.PI*i/n));
408                v.y = Math.Cos(-(2*Math.PI*i/n));
409                v = hn+hmnc-v*(hn-hmnc);
410                a[2*i+0] = 0.5*v.x;
411                a[2*i+1] = 0.5*v.y;
412            }
413            a[1] = buf[0]-buf[1];
414        }
415
416
417        /*************************************************************************
418        Internal subroutine. Never call it directly!
419
420
421          -- ALGLIB --
422             Copyright 01.06.2009 by Bochkanov Sergey
423        *************************************************************************/
424        public static void fftr1dinvinternaleven(ref double[] a,
425            int n,
426            ref double[] buf,
427            ref ftbase.ftplan plan)
428        {
429            double x = 0;
430            double y = 0;
431            double t = 0;
432            int i = 0;
433            int n2 = 0;
434
435            System.Diagnostics.Debug.Assert(n>0 & n%2==0, "FFTR1DInvInternalEven: incorrect N!");
436           
437            //
438            // Special cases:
439            // * N=2
440            //
441            // After this block we assume that N is strictly greater than 2
442            //
443            if( n==2 )
444            {
445                x = 0.5*(a[0]+a[1]);
446                y = 0.5*(a[0]-a[1]);
447                a[0] = x;
448                a[1] = y;
449                return;
450            }
451           
452            //
453            // inverse real FFT is reduced to the inverse real FHT,
454            // which is reduced to the forward real FHT,
455            // which is reduced to the forward real FFT.
456            //
457            // Don't worry, it is really compact and efficient reduction :)
458            //
459            n2 = n/2;
460            buf[0] = a[0];
461            for(i=1; i<=n2-1; i++)
462            {
463                x = a[2*i+0];
464                y = a[2*i+1];
465                buf[i] = x-y;
466                buf[n-i] = x+y;
467            }
468            buf[n2] = a[1];
469            fftr1dinternaleven(ref buf, n, ref a, ref plan);
470            a[0] = buf[0]/n;
471            t = (double)(1)/(double)(n);
472            for(i=1; i<=n2-1; i++)
473            {
474                x = buf[2*i+0];
475                y = buf[2*i+1];
476                a[i] = t*(x-y);
477                a[n-i] = t*(x+y);
478            }
479            a[n2] = buf[1]/n;
480        }
481    }
482}
Note: See TracBrowser for help on using the repository browser.