/*
 *	Karinto Library Project
 *
 *	This software is distributed under a zlib-style license.
 *	See license.txt for more information.
 */

using System;
using System.Collections.Generic;
using System.Text;
using System.Runtime.InteropServices;
using System.Diagnostics;
using Karinto.Numeric;

/*
 *  Fft class uses FFTSS which is an open source library for 
 *  computing the Fast Fourier Transform.
 *  FFTSS is distributed under a BSD-style license.
 *  http://www.ssisc.org/fftss/
 */

namespace Karinto
{
#if WITH_FFTSS

    using Plan = IntPtr;

    public unsafe class Fft
    {
        #region type definitions
        
        [Flags]
        private enum Flags : uint
        {
            Measure = 0,
            DestroyInput = 1,
            Unaligned = 2,
            ConserveMemory = 4,
            Exhaustive = 8,
            PreserveInput = 16,
            Patient = 32,
            Estimate = 64,
        }

        public enum Direction : int
        {
            Forward = -1,
            Backward = 1,
        }

        delegate void* _Malloc(int size);
        delegate void* _Free(void* ptr);
        delegate Plan _PlanDft1d(Int32 n, double* input, double* output,
                                                Direction sign, Flags flags);
        delegate Plan _PlanDft2d(Int32 nx, Int32 ny, Int32 py,
                double* input, double* output, Direction sign, Flags flags);
        delegate void _Execute(Plan p);
        delegate void _Destroy(Plan p);

        #endregion

        #region interop
        // Int32 の箇所は C/C++ での long 
        private class Fftss
        {

            [DllImport(@"libfftss.dll",
                EntryPoint = "fftss_malloc",
                ExactSpelling = true)]
            public static extern void* Malloc(int size);

            [DllImport(@"libfftss.dll",
                EntryPoint = "fftss_free",
                ExactSpelling = true)]
            public static extern void* Free(void* ptr);

            [DllImport(@"libfftss.dll",
                 EntryPoint = "fftss_plan_dft_1d",
                 ExactSpelling = true)]
            public static extern Plan PlanDft1d(
                                        Int32 n, double* input, double* output, 
                                        Direction sign, Flags flags);

            [DllImport(@"libfftss.dll",
                 EntryPoint = "fftss_plan_dft_2d",
                 ExactSpelling = true)]
            public static extern Plan PlanDft2d(
                Int32 nx, Int32 ny, Int32 py,
                double* input, double* output, Direction sign, Flags flags);

            [DllImport(@"libfftss.dll",
                EntryPoint = "fftss_execute",
                ExactSpelling = true)]
            public static extern void Execute(Plan p);

            [DllImport(@"libfftss.dll",
                EntryPoint = "fftss_destroy_plan",
                ExactSpelling = true)]
            public static extern void Destroy(Plan p);
        }
        
        // for x64
        private class FftssX64 : Fftss
        {
            [DllImport(@"x64\libfftss.dll",
                EntryPoint = "fftss_malloc",
                ExactSpelling = true)]
            public static new extern void* Malloc(int size);

            [DllImport(@"x64\libfftss.dll",
                EntryPoint = "fftss_free",
                ExactSpelling = true)]
            public static new extern void* Free(void* ptr);

            [DllImport(@"x64\libfftss.dll",
                 EntryPoint = "fftss_plan_dft_1d",
                 ExactSpelling = true)]
            public static new extern Plan PlanDft1d(
                                        Int32 n, double* input, double* output, 
                                        Direction sign, Flags flags);

            [DllImport(@"x64\libfftss.dll",
                 EntryPoint = "fftss_plan_dft_2d",
                 ExactSpelling = true)]
            public static new extern Plan PlanDft2d(
                Int32 nx, Int32 ny, Int32 py,
                double* input, double* output, Direction sign, Flags flags);

            [DllImport(@"x64\libfftss.dll",
                EntryPoint = "fftss_execute",
                ExactSpelling = true)]
            public static new extern void Execute(Plan p);

            [DllImport(@"x64\libfftss.dll",
                EntryPoint = "fftss_destroy_plan",
                ExactSpelling = true)]
            public static new extern void Destroy(Plan p);
        }
        #endregion

        #region fields
        private static readonly _Malloc Malloc;
        private static readonly _Free Free;
        private static readonly _PlanDft1d PlanDft1d;
        private static readonly _PlanDft2d PlanDft2d;
        private static readonly _Execute Execute;
        private static readonly _Destroy Destroy;

        private const int sizeOfComplex = 16;
        private double* alignedSrc;
        private double* alignedDest;
        int size;
        Plan plan;
        Direction direction;

        #endregion

        #region constructors
        static Fft()
        {
            try
            {
                Fftss.Free(Fftss.Malloc(1));

                Malloc = Fftss.Malloc;
                Free = Fftss.Free;
                PlanDft1d = Fftss.PlanDft1d;
                PlanDft2d = Fftss.PlanDft2d;
                Execute = Fftss.Execute;
                Destroy = Fftss.Destroy;
            }
            catch (Exception ex)
            {
                if (ex is DllNotFoundException || 
                    ex is BadImageFormatException)
                {
                    FftssX64.Free(FftssX64.Malloc(1));
                    Malloc = FftssX64.Malloc;
                    Free = FftssX64.Free;
                    PlanDft1d = FftssX64.PlanDft1d;
                    PlanDft2d = FftssX64.PlanDft2d;
                    Execute = FftssX64.Execute;
                    Destroy = FftssX64.Destroy;
                }
                else
                {
                    throw ex;
                }
            }
        }

        public Fft(int size, Direction direction)
        {
            this.size = size;
            this.direction = direction;
            alignedSrc = (double*)Malloc(size * sizeOfComplex);
            alignedDest = (double*)Malloc(size * sizeOfComplex);
            plan = PlanDft1d(
                    size, alignedSrc, alignedDest, direction, Flags.Measure);
        }

        #endregion

        #region destructors

        ~Fft()
        {
            if (plan != IntPtr.Zero)
            {
                Destroy(plan);
            }
            Free(alignedSrc);
            Free(alignedDest);
        }

        #endregion

        #region public methods
        public Complex[] Transform(Complex[] input)
        {
            Complex[] output = new Complex[size];
            fixed (void* pInput = input)
            {
                CopyDoubleArray((double*)pInput, alignedSrc, size * 2);
            }

            Execute(plan);

            fixed (void* pOutput = output)
            {
                CopyDoubleArray(alignedDest, (double*)pOutput, size * 2);
            }
            return output;
        }

        public void Transform(double[] input, out Complex[] output)
        {
            output = new Complex[size];
            fixed (void* pInput = input)
            {
                double* src = (double*)pInput;
                double* dest = alignedSrc;
                for (int i = 0; i < size; ++i)
                {
                    *dest++ = *src++;
                    *dest++ = 0.0;
                }
            }
            Execute(plan);

            fixed (void* pOutput = output)
            {
                CopyDoubleArray(alignedDest, (double*)pOutput, size * 2);
            }
        }

        public void Transform(Complex[] input, out double[] output)
        {
            output = new double[size];
            fixed (void* pInput = input)
            {
                CopyDoubleArray((double*)pInput, alignedSrc, size * 2);
            }
            Execute(plan);

            fixed (void* pOutput = output)
            {
                double* dest = (double*)pOutput;
                double* src = alignedDest;
                for (int i = 0; i < size; ++i)
                {
                    *dest = *src++;
                    *dest += 2;
                };
            }
        }

        public void Scale(Complex[] target)
        { 
            Complex scale = direction == Direction.Forward ? 1.0 / size: size;
            for (int i = 0; i < size; ++i)
            {
                target[i] *= scale;
            }
        }

        public static Complex[] TransformForward(Complex[] input)
        {
            return Transform(input, Direction.Forward);
        }

        public static Complex[] TransformBackward(Complex[] input)
        {
            return Transform(input, Direction.Backward);
        }

        #endregion

        #region private methods

        private static Complex[] Transform(Complex[] input, Direction direction)
        {
            int n = input.Length;

            Complex[] output = new Complex[n];
            double* alignedSrc = (double*)Malloc(n * sizeOfComplex);
            double* alignedDest = (double*)Malloc(n * sizeOfComplex);

            fixed (void* pInput = input)
            {
                CopyDoubleArray((double*)pInput, alignedSrc, n * 2);
            }

            Plan plan = PlanDft1d(n, alignedSrc, alignedDest,
                    direction, Flags.PreserveInput | Flags.Estimate);
            Execute(plan);
            Destroy(plan);

            fixed (void* pOutput = output)
            {
                CopyDoubleArray(alignedDest, (double*)pOutput, n * 2);
            }
            Free(alignedSrc);
            Free(alignedDest);

            return output;
        }

        private static void CopyDoubleArray(double* src, double* dest, int size)
        {
            for (int i = 0; i < size; ++i)
            {
                *dest++ = *src++;
            }
        }

        /*
        private static void DebugWriteComplexArray(double* src, int size)
        {
            for (int i = 0; i < size; ++i)
            {
                double re = *src++;
                double im = *src++;

                Debug.WriteLine( re.ToString()+ "\t" + im.ToString());
            }
            
        }
        */

        #endregion
    }
#endif
}