source: branches/HeuristicLab.ALGLIB-2.5.0/ALGLIB-2.5.0/fft.cs @ 5192

/*************************************************************************
Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/

using System;

namespace alglib
{
    public class fft
    {
        /*************************************************************************
        1-dimensional complex FFT.

        Array size N may be arbitrary number (composite or prime).  Composite  N's
        are handled with cache-oblivious variation of  a  Cooley-Tukey  algorithm.
        Small prime-factors are transformed using hard coded  codelets (similar to
        FFTW codelets, but without low-level  optimization),  large  prime-factors
        are handled with Bluestein's algorithm.

        Fastests transforms are for smooth N's (prime factors are 2, 3,  5  only),
        most fast for powers of 2. When N have prime factors  larger  than  these,
        but orders of magnitude smaller than N, computations will be about 4 times
        slower than for nearby highly composite N's. When N itself is prime, speed
        will be 6 times lower.

        Algorithm has O(N*logN) complexity for any N (composite or prime).

        INPUT PARAMETERS
            A   -   array[0..N-1] - complex function to be transformed
            N   -   problem size
            
        OUTPUT PARAMETERS
            A   -   DFT of a input array, array[0..N-1]
                    A_out[j] = SUM(A_in[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)


          -- ALGLIB --
             Copyright 29.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftc1d(ref AP.Complex[] a,
            int n)
        {
            ftbase.ftplan plan = new ftbase.ftplan();
            int i = 0;
            double[] buf = new double[0];

            System.Diagnostics.Debug.Assert(n>0, "FFTC1D: incorrect N!");
            
            //
            // Special case: N=1, FFT is just identity transform.
            // After this block we assume that N is strictly greater than 1.
            //
            if( n==1 )
            {
                return;
            }
            
            //
            // convert input array to the more convinient format
            //
            buf = new double[2*n];
            for(i=0; i<=n-1; i++)
            {
                buf[2*i+0] = a[i].x;
                buf[2*i+1] = a[i].y;
            }
            
            //
            // Generate plan and execute it.
            //
            // Plan is a combination of a successive factorizations of N and
            // precomputed data. It is much like a FFTW plan, but is not stored
            // between subroutine calls and is much simpler.
            //
            ftbase.ftbasegeneratecomplexfftplan(n, ref plan);
            ftbase.ftbaseexecuteplan(ref buf, 0, n, ref plan);
            
            //
            // result
            //
            for(i=0; i<=n-1; i++)
            {
                a[i].x = buf[2*i+0];
                a[i].y = buf[2*i+1];
            }
        }


        /*************************************************************************
        1-dimensional complex inverse FFT.

        Array size N may be arbitrary number (composite or prime).  Algorithm  has
        O(N*logN) complexity for any N (composite or prime).

        See FFTC1D() description for more information about algorithm performance.

        INPUT PARAMETERS
            A   -   array[0..N-1] - complex array to be transformed
            N   -   problem size

        OUTPUT PARAMETERS
            A   -   inverse DFT of a input array, array[0..N-1]
                    A_out[j] = SUM(A_in[k]/N*exp(+2*pi*sqrt(-1)*j*k/N), k = 0..N-1)


          -- ALGLIB --
             Copyright 29.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftc1dinv(ref AP.Complex[] a,
            int n)
        {
            int i = 0;

            System.Diagnostics.Debug.Assert(n>0, "FFTC1DInv: incorrect N!");
            
            //
            // Inverse DFT can be expressed in terms of the DFT as
            //
            //     invfft(x) = fft(x')'/N
            //
            // here x' means conj(x).
            //
            for(i=0; i<=n-1; i++)
            {
                a[i].y = -a[i].y;
            }
            fftc1d(ref a, n);
            for(i=0; i<=n-1; i++)
            {
                a[i].x = a[i].x/n;
                a[i].y = -(a[i].y/n);
            }
        }


        /*************************************************************************
        1-dimensional real FFT.

        Algorithm has O(N*logN) complexity for any N (composite or prime).

        INPUT PARAMETERS
            A   -   array[0..N-1] - real function to be transformed
            N   -   problem size

        OUTPUT PARAMETERS
            F   -   DFT of a input array, array[0..N-1]
                    F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)

        NOTE:
            F[] satisfies symmetry property F[k] = conj(F[N-k]),  so just one half
        of  array  is  usually needed. But for convinience subroutine returns full
        complex array (with frequencies above N/2), so its result may be  used  by
        other FFT-related subroutines.


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftr1d(ref double[] a,
            int n,
            ref AP.Complex[] f)
        {
            int i = 0;
            int n2 = 0;
            int idx = 0;
            AP.Complex hn = 0;
            AP.Complex hmnc = 0;
            AP.Complex v = 0;
            double[] buf = new double[0];
            ftbase.ftplan plan = new ftbase.ftplan();
            int i_ = 0;

            System.Diagnostics.Debug.Assert(n>0, "FFTR1D: incorrect N!");
            
            //
            // Special cases:
            // * N=1, FFT is just identity transform.
            // * N=2, FFT is simple too
            //
            // After this block we assume that N is strictly greater than 2
            //
            if( n==1 )
            {
                f = new AP.Complex[1];
                f[0] = a[0];
                return;
            }
            if( n==2 )
            {
                f = new AP.Complex[2];
                f[0].x = a[0]+a[1];
                f[0].y = 0;
                f[1].x = a[0]-a[1];
                f[1].y = 0;
                return;
            }
            
            //
            // Choose between odd-size and even-size FFTs
            //
            if( n%2==0 )
            {
                
                //
                // even-size real FFT, use reduction to the complex task
                //
                n2 = n/2;
                buf = new double[n];
                for(i_=0; i_<=n-1;i_++)
                {
                    buf[i_] = a[i_];
                }
                ftbase.ftbasegeneratecomplexfftplan(n2, ref plan);
                ftbase.ftbaseexecuteplan(ref buf, 0, n2, ref plan);
                f = new AP.Complex[n];
                for(i=0; i<=n2; i++)
                {
                    idx = 2*(i%n2);
                    hn.x = buf[idx+0];
                    hn.y = buf[idx+1];
                    idx = 2*((n2-i)%n2);
                    hmnc.x = buf[idx+0];
                    hmnc.y = -buf[idx+1];
                    v.x = -Math.Sin(-(2*Math.PI*i/n));
                    v.y = Math.Cos(-(2*Math.PI*i/n));
                    f[i] = hn+hmnc-v*(hn-hmnc);
                    f[i].x = 0.5*f[i].x;
                    f[i].y = 0.5*f[i].y;
                }
                for(i=n2+1; i<=n-1; i++)
                {
                    f[i] = AP.Math.Conj(f[n-i]);
                }
                return;
            }
            else
            {
                
                //
                // use complex FFT
                //
                f = new AP.Complex[n];
                for(i=0; i<=n-1; i++)
                {
                    f[i] = a[i];
                }
                fftc1d(ref f, n);
                return;
            }
        }


        /*************************************************************************
        1-dimensional real inverse FFT.

        Algorithm has O(N*logN) complexity for any N (composite or prime).

        INPUT PARAMETERS
            F   -   array[0..floor(N/2)] - frequencies from forward real FFT
            N   -   problem size

        OUTPUT PARAMETERS
            A   -   inverse DFT of a input array, array[0..N-1]

        NOTE:
            F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just  one
        half of frequencies array is needed - elements from 0 to floor(N/2).  F[0]
        is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd,  then
        F[floor(N/2)] has no special properties.

        Relying on properties noted above, FFTR1DInv subroutine uses only elements
        from 0th to floor(N/2)-th. It ignores imaginary part of F[0],  and in case
        N is even it ignores imaginary part of F[floor(N/2)] too.  So you can pass
        either frequencies array with N elements or reduced array with roughly N/2
        elements - subroutine will successfully transform both.


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftr1dinv(ref AP.Complex[] f,
            int n,
            ref double[] a)
        {
            int i = 0;
            double[] h = new double[0];
            AP.Complex[] fh = new AP.Complex[0];

            System.Diagnostics.Debug.Assert(n>0, "FFTR1DInv: incorrect N!");
            
            //
            // Special case: N=1, FFT is just identity transform.
            // After this block we assume that N is strictly greater than 1.
            //
            if( n==1 )
            {
                a = new double[1];
                a[0] = f[0].x;
                return;
            }
            
            //
            // inverse real FFT is reduced to the inverse real FHT,
            // which is reduced to the forward real FHT,
            // which is reduced to the forward real FFT.
            //
            // Don't worry, it is really compact and efficient reduction :)
            //
            h = new double[n];
            a = new double[n];
            h[0] = f[0].x;
            for(i=1; i<=(int)Math.Floor((double)(n)/(double)(2))-1; i++)
            {
                h[i] = f[i].x-f[i].y;
                h[n-i] = f[i].x+f[i].y;
            }
            if( n%2==0 )
            {
                h[(int)Math.Floor((double)(n)/(double)(2))] = f[(int)Math.Floor((double)(n)/(double)(2))].x;
            }
            else
            {
                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;
                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;
            }
            fftr1d(ref h, n, ref fh);
            for(i=0; i<=n-1; i++)
            {
                a[i] = (fh[i].x-fh[i].y)/n;
            }
        }


        /*************************************************************************
        Internal subroutine. Never call it directly!


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftr1dinternaleven(ref double[] a,
            int n,
            ref double[] buf,
            ref ftbase.ftplan plan)
        {
            double x = 0;
            double y = 0;
            int i = 0;
            int n2 = 0;
            int idx = 0;
            AP.Complex hn = 0;
            AP.Complex hmnc = 0;
            AP.Complex v = 0;
            int i_ = 0;

            System.Diagnostics.Debug.Assert(n>0 & n%2==0, "FFTR1DEvenInplace: incorrect N!");
            
            //
            // Special cases:
            // * N=2
            //
            // After this block we assume that N is strictly greater than 2
            //
            if( n==2 )
            {
                x = a[0]+a[1];
                y = a[0]-a[1];
                a[0] = x;
                a[1] = y;
                return;
            }
            
            //
            // even-size real FFT, use reduction to the complex task
            //
            n2 = n/2;
            for(i_=0; i_<=n-1;i_++)
            {
                buf[i_] = a[i_];
            }
            ftbase.ftbaseexecuteplan(ref buf, 0, n2, ref plan);
            a[0] = buf[0]+buf[1];
            for(i=1; i<=n2-1; i++)
            {
                idx = 2*(i%n2);
                hn.x = buf[idx+0];
                hn.y = buf[idx+1];
                idx = 2*(n2-i);
                hmnc.x = buf[idx+0];
                hmnc.y = -buf[idx+1];
                v.x = -Math.Sin(-(2*Math.PI*i/n));
                v.y = Math.Cos(-(2*Math.PI*i/n));
                v = hn+hmnc-v*(hn-hmnc);
                a[2*i+0] = 0.5*v.x;
                a[2*i+1] = 0.5*v.y;
            }
            a[1] = buf[0]-buf[1];
        }


        /*************************************************************************
        Internal subroutine. Never call it directly!


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void fftr1dinvinternaleven(ref double[] a,
            int n,
            ref double[] buf,
            ref ftbase.ftplan plan)
        {
            double x = 0;
            double y = 0;
            double t = 0;
            int i = 0;
            int n2 = 0;

            System.Diagnostics.Debug.Assert(n>0 & n%2==0, "FFTR1DInvInternalEven: incorrect N!");
            
            //
            // Special cases:
            // * N=2
            //
            // After this block we assume that N is strictly greater than 2
            //
            if( n==2 )
            {
                x = 0.5*(a[0]+a[1]);
                y = 0.5*(a[0]-a[1]);
                a[0] = x;
                a[1] = y;
                return;
            }
            
            //
            // inverse real FFT is reduced to the inverse real FHT,
            // which is reduced to the forward real FHT,
            // which is reduced to the forward real FFT.
            //
            // Don't worry, it is really compact and efficient reduction :)
            //
            n2 = n/2;
            buf[0] = a[0];
            for(i=1; i<=n2-1; i++)
            {
                x = a[2*i+0];
                y = a[2*i+1];
                buf[i] = x-y;
                buf[n-i] = x+y;
            }
            buf[n2] = a[1];
            fftr1dinternaleven(ref buf, n, ref a, ref plan);
            a[0] = buf[0]/n;
            t = (double)(1)/(double)(n);
            for(i=1; i<=n2-1; i++)
            {
                x = buf[2*i+0];
                y = buf[2*i+1];
                a[i] = t*(x-y);
                a[n-i] = t*(x+y);
            }
            a[n2] = buf[1]/n;
        }
    }
}