''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'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 >>>
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'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 Sub FFTC1D(ByRef A() As Complex, ByVal N As Long)
    Dim Plan As FTPlan
    Dim I As Long
    Dim Buf() As Double

    
    '
    ' Special case: N=1, FFT is just identity transform.
    ' After this block we assume that N is strictly greater than 1.
    '
    If N=1# then
        Exit Sub
    End If
    
    '
    ' convert input array to the more convinient format
    '
    ReDim Buf(0 To 2#*N-1)
    For I=0# To N-1# Step 1
        Buf(2#*I+0#) = A(I).X
        Buf(2#*I+1#) = A(I).Y
    Next I
    
    '
    ' 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.
    '
    Call FTBaseGenerateComplexFFTPlan(N, Plan)
    Call FTBaseExecutePlan(Buf, 0#, N, Plan)
    
    '
    ' result
    '
    For I=0# To N-1# Step 1
        A(I).X = Buf(2#*I+0#)
        A(I).Y = Buf(2#*I+1#)
    Next I
End Sub


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'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 Sub FFTC1DInv(ByRef A() As Complex, ByVal N As Long)
    Dim I As Long

    
    '
    ' Inverse DFT can be expressed in terms of the DFT as
    '
    '     invfft(x) = fft(x')'/N
    '
    ' here x' means conj(x).
    '
    For I=0# To N-1# Step 1
        A(I).Y = -A(I).Y
    Next I
    Call FFTC1D(A, N)
    For I=0# To N-1# Step 1
        A(I).X = A(I).X/N
        A(I).Y = -(A(I).Y/N)
    Next I
End Sub


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'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 Sub FFTR1D(ByRef A() As Double, ByVal N As Long, ByRef F() As Complex)
    Dim I As Long
    Dim N2 As Long
    Dim Idx As Long
    Dim Hn As Complex
    Dim HmnC As Complex
    Dim V As Complex
    Dim Buf() As Double
    Dim Plan As FTPlan
    Dim i_ As Long

    
    '
    ' 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# then
        ReDim F(0 To 1#-1)
        F(0#) = C_Complex(A(0#))
        Exit Sub
    End If
    If N=2# then
        ReDim F(0 To 2#-1)
        F(0#).X = A(0#)+A(1#)
        F(0#).Y = 0#
        F(1#).X = A(0#)-A(1#)
        F(1#).Y = 0#
        Exit Sub
    End If
    
    '
    ' Choose between odd-size and even-size FFTs
    '
    If N Mod 2#=0# then
        
        '
        ' even-size real FFT, use reduction to the complex task
        '
        N2 = N\2#
        ReDim Buf(0 To N-1)
        For i_ = 0# To N-1# Step 1
            Buf(i_) = A(i_)
        Next i_
        Call FTBaseGenerateComplexFFTPlan(N2, Plan)
        Call FTBaseExecutePlan(Buf, 0#, N2, Plan)
        ReDim F(0 To N-1)
        For I=0# To N2 Step 1
            Idx = 2#*(I Mod N2)
            Hn.X = Buf(Idx+0#)
            Hn.Y = Buf(Idx+1#)
            Idx = 2#*((N2-I) Mod N2)
            HmnC.X = Buf(Idx+0#)
            HmnC.Y = -Buf(Idx+1#)
            V.X = -Sin(-(2#*Pi()*I/N))
            V.Y = Cos(-(2#*Pi()*I/N))
            F(I) = C_Sub(C_Add(Hn,HmnC),C_Mul(V,C_Sub(Hn,HmnC)))
            F(I).X = 0.5*F(I).X
            F(I).Y = 0.5*F(I).Y
        Next I
        For I=N2+1# To N-1# Step 1
            F(I) = Conj(F(N-I))
        Next I
        Exit Sub
    Else
        
        '
        ' use complex FFT
        '
        ReDim F(0 To N-1)
        For I=0# To N-1# Step 1
            F(I) = C_Complex(A(I))
        Next I
        Call FFTC1D(F, N)
        Exit Sub
    End If
End Sub


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'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 Sub FFTR1DInv(ByRef F() As Complex, _
         ByVal N As Long, _
         ByRef A() As Double)
    Dim I As Long
    Dim H() As Double
    Dim FH() As Complex

    
    '
    ' Special case: N=1, FFT is just identity transform.
    ' After this block we assume that N is strictly greater than 1.
    '
    If N=1# then
        ReDim A(0 To 1#-1)
        A(0#) = F(0#).X
        Exit Sub
    End If
    
    '
    ' 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 :)
    '
    ReDim H(0 To N-1)
    ReDim A(0 To N-1)
    H(0#) = F(0#).X
    For I=1# To Int(N/2#)-1# Step 1
        H(I) = F(I).X-F(I).Y
        H(N-I) = F(I).X+F(I).Y
    Next I
    If N Mod 2#=0# then
        H(Int(N/2#)) = F(Int(N/2#)).X
    Else
        H(Int(N/2#)) = F(Int(N/2#)).X-F(Int(N/2#)).Y
        H(Int(N/2#)+1#) = F(Int(N/2#)).X+F(Int(N/2#)).Y
    End If
    Call FFTR1D(H, N, FH)
    For I=0# To N-1# Step 1
        A(I) = (FH(I).X-FH(I).Y)/N
    Next I
End Sub


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Internal subroutine. Never call it directly!
'
'
'  -- ALGLIB --
'     Copyright 01.06.2009 by Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Sub FFTR1DInternalEven(ByRef A() As Double, _
         ByVal N As Long, _
         ByRef Buf() As Double, _
         ByRef Plan As FTPlan)
    Dim X As Double
    Dim Y As Double
    Dim I As Long
    Dim N2 As Long
    Dim Idx As Long
    Dim Hn As Complex
    Dim HmnC As Complex
    Dim V As Complex
    Dim i_ As Long

    
    '
    ' Special cases:
    ' * N=2
    '
    ' After this block we assume that N is strictly greater than 2
    '
    If N=2# then
        X = A(0#)+A(1#)
        Y = A(0#)-A(1#)
        A(0#) = X
        A(1#) = Y
        Exit Sub
    End If
    
    '
    ' even-size real FFT, use reduction to the complex task
    '
    N2 = N\2#
    For i_ = 0# To N-1# Step 1
        Buf(i_) = A(i_)
    Next i_
    Call FTBaseExecutePlan(Buf, 0#, N2, Plan)
    A(0#) = Buf(0#)+Buf(1#)
    For I=1# To N2-1# Step 1
        Idx = 2#*(I Mod 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 = -Sin(-(2#*Pi()*I/N))
        V.Y = Cos(-(2#*Pi()*I/N))
        V = C_Sub(C_Add(Hn,HmnC),C_Mul(V,C_Sub(Hn,HmnC)))
        A(2#*I+0#) = 0.5*V.X
        A(2#*I+1#) = 0.5*V.Y
    Next I
    A(1#) = Buf(0#)-Buf(1#)
End Sub


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Internal subroutine. Never call it directly!
'
'
'  -- ALGLIB --
'     Copyright 01.06.2009 by Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Sub FFTR1DInvInternalEven(ByRef A() As Double, _
         ByVal N As Long, _
         ByRef Buf() As Double, _
         ByRef Plan As FTPlan)
    Dim X As Double
    Dim Y As Double
    Dim T As Double
    Dim I As Long
    Dim N2 As Long

    
    '
    ' Special cases:
    ' * N=2
    '
    ' After this block we assume that N is strictly greater than 2
    '
    If N=2# then
        X = 0.5*(A(0#)+A(1#))
        Y = 0.5*(A(0#)-A(1#))
        A(0#) = X
        A(1#) = Y
        Exit Sub
    End If
    
    '
    ' 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# To N2-1# Step 1
        X = A(2#*I+0#)
        Y = A(2#*I+1#)
        Buf(I) = X-Y
        Buf(N-I) = X+Y
    Next I
    Buf(N2) = A(1#)
    Call FFTR1DInternalEven(Buf, N, A, Plan)
    A(0#) = Buf(0#)/N
    T = 1#/N
    For I=1# To N2-1# Step 1
        X = Buf(2#*I+0#)
        Y = Buf(2#*I+1#)
        A(I) = T*(X-Y)
        A(N-I) = T*(X+Y)
    Next I
    A(N2) = Buf(1#)/N
End Sub