// <copyright file="Fourier.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2018 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>

using System;
using IStation.Numerics.LinearAlgebra;
using IStation.Numerics.LinearAlgebra.Storage;
using IStation.Numerics.Providers.FourierTransform;
using Complex = System.Numerics.Complex;

namespace IStation.Numerics.IntegralTransforms
{
    /// <summary>
    /// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
    /// </summary>
    public static partial class Fourier
    {
        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
        public static void Forward(Complex32[] samples)
        {
            FourierTransformControl.Provider.Forward(samples, FourierTransformScaling.SymmetricScaling);
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
        public static void Forward(Complex[] samples)
        {
            FourierTransformControl.Provider.Forward(samples, FourierTransformScaling.SymmetricScaling);
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Forward(Complex32[] samples, FourierOptions options)
        {
            switch (options)
            {
                case FourierOptions.NoScaling:
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.Forward(samples, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.InverseExponent:
                    FourierTransformControl.Provider.Backward(samples, FourierTransformScaling.SymmetricScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.NoScaling:
                case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.Backward(samples, FourierTransformScaling.NoScaling);
                    break;
                default:
                    FourierTransformControl.Provider.Forward(samples, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Forward(Complex[] samples, FourierOptions options)
        {
            switch (options)
            {
                case FourierOptions.NoScaling:
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.Forward(samples, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.InverseExponent:
                    FourierTransformControl.Provider.Backward(samples, FourierTransformScaling.SymmetricScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.NoScaling:
                case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.Backward(samples, FourierTransformScaling.NoScaling);
                    break;
                default:
                    FourierTransformControl.Provider.Forward(samples, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="real">Real part of the sample vector, where the FFT is evaluated in place.</param>
        /// <param name="imaginary">Imaginary part of the sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Forward(float[] real, float[] imaginary, FourierOptions options = FourierOptions.Default)
        {
            if (real.Length != imaginary.Length)
            {
                throw new ArgumentException("The array arguments must have the same length.");
            }

            // TODO: consider to support this natively by the provider, without the need for copying
            // TODO: otherwise, consider ArrayPool

            Complex32[] data = new Complex32[real.Length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = new Complex32(real[i], imaginary[i]);
            }

            Forward(data, options);

            for (int i = 0; i < data.Length; i++)
            {
                real[i] = data[i].Real;
                imaginary[i] = data[i].Imaginary;
            }
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="real">Real part of the sample vector, where the FFT is evaluated in place.</param>
        /// <param name="imaginary">Imaginary part of the sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Forward(double[] real, double[] imaginary, FourierOptions options = FourierOptions.Default)
        {
            if (real.Length != imaginary.Length)
            {
                throw new ArgumentException("The array arguments must have the same length.");
            }

            // TODO: consider to support this natively by the provider, without the need for copying
            // TODO: otherwise, consider ArrayPool

            Complex[] data = new Complex[real.Length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = new Complex(real[i], imaginary[i]);
            }

            Forward(data, options);

            for (int i = 0; i < data.Length; i++)
            {
                real[i] = data[i].Real;
                imaginary[i] = data[i].Imaginary;
            }
        }

        /// <summary>
        /// Packed Real-Complex forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
        /// Since for real-valued time samples the complex spectrum is conjugate-even (symmetry),
        /// the spectrum can be fully reconstructed from the positive frequencies only (first half).
        /// The data array needs to be N+2 (if N is even) or N+1 (if N is odd) long in order to support such a packed spectrum.
        /// </summary>
        /// <param name="data">Data array of length N+2 (if N is even) or N+1 (if N is odd).</param>
        /// <param name="n">The number of samples.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void ForwardReal(float[] data, int n, FourierOptions options = FourierOptions.Default)
        {
            int length = n.IsEven() ? n + 2 : n + 1;
            if (data.Length < length)
            {
                throw new ArgumentException($"The given array is too small. It must be at least {length} long.");
            }

            if ((options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent)
            {
                throw new NotSupportedException();
            }

            switch (options)
            {
                case FourierOptions.NoScaling:
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.ForwardReal(data, n, FourierTransformScaling.NoScaling);
                    break;
                default:
                    FourierTransformControl.Provider.ForwardReal(data, n, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Packed Real-Complex forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
        /// Since for real-valued time samples the complex spectrum is conjugate-even (symmetry),
        /// the spectrum can be fully reconstructed form the positive frequencies only (first half).
        /// The data array needs to be N+2 (if N is even) or N+1 (if N is odd) long in order to support such a packed spectrum.
        /// </summary>
        /// <param name="data">Data array of length N+2 (if N is even) or N+1 (if N is odd).</param>
        /// <param name="n">The number of samples.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void ForwardReal(double[] data, int n, FourierOptions options = FourierOptions.Default)
        {
            int length = n.IsEven() ? n + 2 : n + 1;
            if (data.Length < length)
            {
                throw new ArgumentException($"The given array is too small. It must be at least {length} long.");
            }

            if ((options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent)
            {
                throw new NotSupportedException();
            }

            switch (options)
            {
                case FourierOptions.NoScaling:
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.ForwardReal(data, n, FourierTransformScaling.NoScaling);
                    break;
                default:
                    FourierTransformControl.Provider.ForwardReal(data, n, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to multiple dimensional sample data.
        /// </summary>
        /// <param name="samples">Sample data, where the FFT is evaluated in place.</param>
        /// <param name="dimensions">
        /// The data size per dimension. The first dimension is the major one.
        /// For example, with two dimensions "rows" and "columns" the samples are assumed to be organized row by row.
        /// </param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void ForwardMultiDim(Complex32[] samples, int[] dimensions, FourierOptions options = FourierOptions.Default)
        {
            switch (options)
            {
                case FourierOptions.NoScaling:
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.ForwardMultidim(samples, dimensions, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.InverseExponent:
                    FourierTransformControl.Provider.BackwardMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.NoScaling:
                case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.BackwardMultidim(samples, dimensions, FourierTransformScaling.NoScaling);
                    break;
                default:
                    FourierTransformControl.Provider.ForwardMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to multiple dimensional sample data.
        /// </summary>
        /// <param name="samples">Sample data, where the FFT is evaluated in place.</param>
        /// <param name="dimensions">
        /// The data size per dimension. The first dimension is the major one.
        /// For example, with two dimensions "rows" and "columns" the samples are assumed to be organized row by row.
        /// </param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void ForwardMultiDim(Complex[] samples, int[] dimensions, FourierOptions options = FourierOptions.Default)
        {
            switch (options)
            {
                case FourierOptions.NoScaling:
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.ForwardMultidim(samples, dimensions, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.InverseExponent:
                    FourierTransformControl.Provider.BackwardMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.NoScaling:
                case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.BackwardMultidim(samples, dimensions, FourierTransformScaling.NoScaling);
                    break;
                default:
                    FourierTransformControl.Provider.ForwardMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to two dimensional sample data.
        /// </summary>
        /// <param name="samplesRowWise">Sample data, organized row by row, where the FFT is evaluated in place</param>
        /// <param name="rows">The number of rows.</param>
        /// <param name="columns">The number of columns.</param>
        /// <remarks>Data available organized column by column instead of row by row can be processed directly by swapping the rows and columns arguments.</remarks>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Forward2D(Complex32[] samplesRowWise, int rows, int columns, FourierOptions options = FourierOptions.Default)
        {
            ForwardMultiDim(samplesRowWise, new[] { rows, columns }, options);
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to two dimensional sample data.
        /// </summary>
        /// <param name="samplesRowWise">Sample data, organized row by row, where the FFT is evaluated in place</param>
        /// <param name="rows">The number of rows.</param>
        /// <param name="columns">The number of columns.</param>
        /// <remarks>Data available organized column by column instead of row by row can be processed directly by swapping the rows and columns arguments.</remarks>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Forward2D(Complex[] samplesRowWise, int rows, int columns, FourierOptions options = FourierOptions.Default)
        {
            ForwardMultiDim(samplesRowWise, new[] { rows, columns }, options);
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to a two dimensional data in form of a matrix.
        /// </summary>
        /// <param name="samples">Sample matrix, where the FFT is evaluated in place</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Forward2D(Matrix<Complex32> samples, FourierOptions options = FourierOptions.Default)
        {
            var rowMajorArray = samples.AsRowMajorArray();
            if (rowMajorArray != null)
            {
                ForwardMultiDim(rowMajorArray, new[] { samples.RowCount, samples.ColumnCount }, options);
                return;
            }

            var columnMajorArray = samples.AsColumnMajorArray();
            if (columnMajorArray != null)
            {
                ForwardMultiDim(columnMajorArray, new[] { samples.ColumnCount, samples.RowCount }, options);
                return;
            }

            // Fall Back
            columnMajorArray = samples.ToColumnMajorArray();
            ForwardMultiDim(columnMajorArray, new[] { samples.ColumnCount, samples.RowCount }, options);
            var denseStorage = new DenseColumnMajorMatrixStorage<Complex32>(samples.RowCount, samples.ColumnCount, columnMajorArray);
            denseStorage.CopyToUnchecked(samples.Storage, ExistingData.Clear);
        }

        /// <summary>
        /// Applies the forward Fast Fourier Transform (FFT) to a two dimensional data in form of a matrix.
        /// </summary>
        /// <param name="samples">Sample matrix, where the FFT is evaluated in place</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Forward2D(Matrix<Complex> samples, FourierOptions options = FourierOptions.Default)
        {
            var rowMajorArray = samples.AsRowMajorArray();
            if (rowMajorArray != null)
            {
                ForwardMultiDim(rowMajorArray, new[] { samples.RowCount, samples.ColumnCount }, options);
                return;
            }

            var columnMajorArray = samples.AsColumnMajorArray();
            if (columnMajorArray != null)
            {
                ForwardMultiDim(columnMajorArray, new[] { samples.ColumnCount, samples.RowCount }, options);
                return;
            }

            // Fall Back
            columnMajorArray = samples.ToColumnMajorArray();
            ForwardMultiDim(columnMajorArray, new[] { samples.ColumnCount, samples.RowCount }, options);
            var denseStorage = new DenseColumnMajorMatrixStorage<Complex>(samples.RowCount, samples.ColumnCount, columnMajorArray);
            denseStorage.CopyToUnchecked(samples.Storage, ExistingData.Clear);
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
        public static void Inverse(Complex32[] spectrum)
        {
            FourierTransformControl.Provider.Backward(spectrum, FourierTransformScaling.SymmetricScaling);
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
        public static void Inverse(Complex[] spectrum)
        {
            FourierTransformControl.Provider.Backward(spectrum, FourierTransformScaling.SymmetricScaling);
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Inverse(Complex32[] spectrum, FourierOptions options)
        {
            switch (options)
            {
                case FourierOptions.NoScaling:
                    FourierTransformControl.Provider.Backward(spectrum, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.Backward(spectrum, FourierTransformScaling.BackwardScaling);
                    break;
                case FourierOptions.InverseExponent:
                    FourierTransformControl.Provider.Forward(spectrum, FourierTransformScaling.SymmetricScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.NoScaling:
                    FourierTransformControl.Provider.Forward(spectrum, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.Forward(spectrum, FourierTransformScaling.ForwardScaling);
                    break;
                default:
                    FourierTransformControl.Provider.Backward(spectrum, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Inverse(Complex[] spectrum, FourierOptions options)
        {
            switch (options)
            {
                case FourierOptions.NoScaling:
                    FourierTransformControl.Provider.Backward(spectrum, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.Backward(spectrum, FourierTransformScaling.BackwardScaling);
                    break;
                case FourierOptions.InverseExponent:
                    FourierTransformControl.Provider.Forward(spectrum, FourierTransformScaling.SymmetricScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.NoScaling:
                    FourierTransformControl.Provider.Forward(spectrum, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.Forward(spectrum, FourierTransformScaling.ForwardScaling);
                    break;
                default:
                    FourierTransformControl.Provider.Backward(spectrum, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="real">Real part of the sample vector, where the iFFT is evaluated in place.</param>
        /// <param name="imaginary">Imaginary part of the sample vector, where the iFFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Inverse(float[] real, float[] imaginary, FourierOptions options = FourierOptions.Default)
        {
            if (real.Length != imaginary.Length)
            {
                throw new ArgumentException("The array arguments must have the same length.");
            }

            // TODO: consider to support this natively by the provider, without the need for copying
            // TODO: otherwise, consider ArrayPool

            Complex32[] data = new Complex32[real.Length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = new Complex32(real[i], imaginary[i]);
            }

            Inverse(data, options);

            for (int i = 0; i < data.Length; i++)
            {
                real[i] = data[i].Real;
                imaginary[i] = data[i].Imaginary;
            }
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
        /// </summary>
        /// <param name="real">Real part of the sample vector, where the iFFT is evaluated in place.</param>
        /// <param name="imaginary">Imaginary part of the sample vector, where the iFFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Inverse(double[] real, double[] imaginary, FourierOptions options = FourierOptions.Default)
        {
            if (real.Length != imaginary.Length)
            {
                throw new ArgumentException("The array arguments must have the same length.");
            }

            // TODO: consider to support this natively by the provider, without the need for copying
            // TODO: otherwise, consider ArrayPool

            Complex[] data = new Complex[real.Length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = new Complex(real[i], imaginary[i]);
            }

            Inverse(data, options);

            for (int i = 0; i < data.Length; i++)
            {
                real[i] = data[i].Real;
                imaginary[i] = data[i].Imaginary;
            }
        }

        /// <summary>
        /// Packed Real-Complex inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
        /// Since for real-valued time samples the complex spectrum is conjugate-even (symmetry),
        /// the spectrum can be fully reconstructed form the positive frequencies only (first half).
        /// The data array needs to be N+2 (if N is even) or N+1 (if N is odd) long in order to support such a packed spectrum.
        /// </summary>
        /// <param name="data">Data array of length N+2 (if N is even) or N+1 (if N is odd).</param>
        /// <param name="n">The number of samples.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void InverseReal(float[] data, int n, FourierOptions options = FourierOptions.Default)
        {
            int length = n.IsEven() ? n + 2 : n + 1;
            if (data.Length < length)
            {
                throw new ArgumentException($"The given array is too small. It must be at least {length} long.");
            }

            if ((options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent)
            {
                throw new NotSupportedException();
            }

            switch (options)
            {
                case FourierOptions.NoScaling:
                    FourierTransformControl.Provider.BackwardReal(data, n, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.BackwardReal(data, n, FourierTransformScaling.BackwardScaling);
                    break;
                default:
                    FourierTransformControl.Provider.BackwardReal(data, n, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Packed Real-Complex inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
        /// Since for real-valued time samples the complex spectrum is conjugate-even (symmetry),
        /// the spectrum can be fully reconstructed form the positive frequencies only (first half).
        /// The data array needs to be N+2 (if N is even) or N+1 (if N is odd) long in order to support such a packed spectrum.
        /// </summary>
        /// <param name="data">Data array of length N+2 (if N is even) or N+1 (if N is odd).</param>
        /// <param name="n">The number of samples.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void InverseReal(double[] data, int n, FourierOptions options = FourierOptions.Default)
        {
            int length = n.IsEven() ? n + 2 : n + 1;
            if (data.Length < length)
            {
                throw new ArgumentException($"The given array is too small. It must be at least {length} long.");
            }

            if ((options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent)
            {
                throw new NotSupportedException();
            }

            switch (options)
            {
                case FourierOptions.NoScaling:
                    FourierTransformControl.Provider.BackwardReal(data, n, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.BackwardReal(data, n, FourierTransformScaling.BackwardScaling);
                    break;
                default:
                    FourierTransformControl.Provider.BackwardReal(data, n, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to multiple dimensional sample data.
        /// </summary>
        /// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
        /// <param name="dimensions">
        /// The data size per dimension. The first dimension is the major one.
        /// For example, with two dimensions "rows" and "columns" the samples are assumed to be organized row by row.
        /// </param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void InverseMultiDim(Complex32[] spectrum, int[] dimensions, FourierOptions options = FourierOptions.Default)
        {
            switch (options)
            {
                case FourierOptions.NoScaling:
                    FourierTransformControl.Provider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.BackwardScaling);
                    break;
                case FourierOptions.InverseExponent:
                    FourierTransformControl.Provider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.SymmetricScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.NoScaling:
                    FourierTransformControl.Provider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.ForwardScaling);
                    break;
                default:
                    FourierTransformControl.Provider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to multiple dimensional sample data.
        /// </summary>
        /// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
        /// <param name="dimensions">
        /// The data size per dimension. The first dimension is the major one.
        /// For example, with two dimensions "rows" and "columns" the samples are assumed to be organized row by row.
        /// </param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void InverseMultiDim(Complex[] spectrum, int[] dimensions, FourierOptions options = FourierOptions.Default)
        {
            switch (options)
            {
                case FourierOptions.NoScaling:
                    FourierTransformControl.Provider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.BackwardScaling);
                    break;
                case FourierOptions.InverseExponent:
                    FourierTransformControl.Provider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.SymmetricScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.NoScaling:
                    FourierTransformControl.Provider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.NoScaling);
                    break;
                case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
                    FourierTransformControl.Provider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.ForwardScaling);
                    break;
                default:
                    FourierTransformControl.Provider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.SymmetricScaling);
                    break;
            }
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to two dimensional sample data.
        /// </summary>
        /// <param name="spectrumRowWise">Sample data, organized row by row, where the iFFT is evaluated in place</param>
        /// <param name="rows">The number of rows.</param>
        /// <param name="columns">The number of columns.</param>
        /// <remarks>Data available organized column by column instead of row by row can be processed directly by swapping the rows and columns arguments.</remarks>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Inverse2D(Complex32[] spectrumRowWise, int rows, int columns, FourierOptions options = FourierOptions.Default)
        {
            InverseMultiDim(spectrumRowWise, new[] { rows, columns }, options);
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to two dimensional sample data.
        /// </summary>
        /// <param name="spectrumRowWise">Sample data, organized row by row, where the iFFT is evaluated in place</param>
        /// <param name="rows">The number of rows.</param>
        /// <param name="columns">The number of columns.</param>
        /// <remarks>Data available organized column by column instead of row by row can be processed directly by swapping the rows and columns arguments.</remarks>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Inverse2D(Complex[] spectrumRowWise, int rows, int columns, FourierOptions options = FourierOptions.Default)
        {
            InverseMultiDim(spectrumRowWise, new[] { rows, columns }, options);
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to a two dimensional data in form of a matrix.
        /// </summary>
        /// <param name="spectrum">Sample matrix, where the iFFT is evaluated in place</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Inverse2D(Matrix<Complex32> spectrum, FourierOptions options = FourierOptions.Default)
        {
            var rowMajorArray = spectrum.AsRowMajorArray();
            if (rowMajorArray != null)
            {
                InverseMultiDim(rowMajorArray, new[] { spectrum.RowCount, spectrum.ColumnCount }, options);
                return;
            }

            var columnMajorArray = spectrum.AsColumnMajorArray();
            if (columnMajorArray != null)
            {
                InverseMultiDim(columnMajorArray, new[] { spectrum.ColumnCount, spectrum.RowCount }, options);
                return;
            }

            // Fall Back
            columnMajorArray = spectrum.ToColumnMajorArray();
            InverseMultiDim(columnMajorArray, new[] { spectrum.ColumnCount, spectrum.RowCount }, options);
            var denseStorage = new DenseColumnMajorMatrixStorage<Complex32>(spectrum.RowCount, spectrum.ColumnCount, columnMajorArray);
            denseStorage.CopyToUnchecked(spectrum.Storage, ExistingData.Clear);
        }

        /// <summary>
        /// Applies the inverse Fast Fourier Transform (iFFT) to a two dimensional data in form of a matrix.
        /// </summary>
        /// <param name="spectrum">Sample matrix, where the iFFT is evaluated in place</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        public static void Inverse2D(Matrix<Complex> spectrum, FourierOptions options = FourierOptions.Default)
        {
            var rowMajorArray = spectrum.AsRowMajorArray();
            if (rowMajorArray != null)
            {
                InverseMultiDim(rowMajorArray, new[] { spectrum.RowCount, spectrum.ColumnCount }, options);
                return;
            }

            var columnMajorArray = spectrum.AsColumnMajorArray();
            if (columnMajorArray != null)
            {
                InverseMultiDim(columnMajorArray, new[] { spectrum.ColumnCount, spectrum.RowCount }, options);
                return;
            }

            // Fall Back
            columnMajorArray = spectrum.ToColumnMajorArray();
            InverseMultiDim(columnMajorArray, new[] { spectrum.ColumnCount, spectrum.RowCount }, options);
            var denseStorage = new DenseColumnMajorMatrixStorage<Complex>(spectrum.RowCount, spectrum.ColumnCount, columnMajorArray);
            denseStorage.CopyToUnchecked(spectrum.Storage, ExistingData.Clear);
        }

        /// <summary>
        /// Naive forward DFT, useful e.g. to verify faster algorithms.
        /// </summary>
        /// <param name="samples">Time-space sample vector.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        /// <returns>Corresponding frequency-space vector.</returns>
        [Obsolete("Use Forward instead. Will be dropped in version 5.0 and behave like Forward until then.")]
        public static Complex32[] NaiveForward(Complex32[] samples, FourierOptions options = FourierOptions.Default)
        {
            var result = new Complex32[samples.Length];
            samples.Copy(result);
            Forward(result, options);
            return result;
        }

        /// <summary>
        /// Naive forward DFT, useful e.g. to verify faster algorithms.
        /// </summary>
        /// <param name="samples">Time-space sample vector.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        /// <returns>Corresponding frequency-space vector.</returns>
        [Obsolete("Use Forward instead. Will be dropped in version 5.0 and behave like Forward until then.")]
        public static Complex[] NaiveForward(Complex[] samples, FourierOptions options = FourierOptions.Default)
        {
            var result = new Complex[samples.Length];
            samples.Copy(result);
            Forward(result, options);
            return result;
        }

        /// <summary>
        /// Naive inverse DFT, useful e.g. to verify faster algorithms.
        /// </summary>
        /// <param name="spectrum">Frequency-space sample vector.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        /// <returns>Corresponding time-space vector.</returns>
        [Obsolete("Use Inverse instead. Will be dropped in version 5.0 and behave like Inverse until then.")]
        public static Complex32[] NaiveInverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default)
        {
            var result = new Complex32[spectrum.Length];
            spectrum.Copy(result);
            Inverse(result, options);
            return result;
        }

        /// <summary>
        /// Naive inverse DFT, useful e.g. to verify faster algorithms.
        /// </summary>
        /// <param name="spectrum">Frequency-space sample vector.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        /// <returns>Corresponding time-space vector.</returns>
        [Obsolete("Use Inverse instead. Will be dropped in version 5.0 and behave like Inverse until then.")]
        public static Complex[] NaiveInverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default)
        {
            var result = new Complex[spectrum.Length];
            spectrum.Copy(result);
            Inverse(result, options);
            return result;
        }

        /// <summary>
        /// Radix-2 forward FFT for power-of-two sized sample vectors.
        /// </summary>
        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        /// <exception cref="ArgumentException"/>
        [Obsolete("Use Forward instead. Will be dropped in version 5.0 and behave like Forward until then.")]
        public static void Radix2Forward(Complex32[] samples, FourierOptions options = FourierOptions.Default)
        {
            Forward(samples, options);
        }

        /// <summary>
        /// Radix-2 forward FFT for power-of-two sized sample vectors.
        /// </summary>
        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        /// <exception cref="ArgumentException"/>
        [Obsolete("Use Forward instead. Will be dropped in version 5.0 and behave like Forward until then.")]
        public static void Radix2Forward(Complex[] samples, FourierOptions options = FourierOptions.Default)
        {
            Forward(samples, options);
        }

        /// <summary>
        /// Radix-2 inverse FFT for power-of-two sized sample vectors.
        /// </summary>
        /// <param name="spectrum">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        /// <exception cref="ArgumentException"/>
        [Obsolete("Use Inverse instead. Will be dropped in version 5.0 and behave like Inverse until then.")]
        public static void Radix2Inverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default)
        {
            Inverse(spectrum, options);
        }

        /// <summary>
        /// Radix-2 inverse FFT for power-of-two sized sample vectors.
        /// </summary>
        /// <param name="spectrum">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        /// <exception cref="ArgumentException"/>
        [Obsolete("Use Inverse instead. Will be dropped in version 5.0 and behave like Inverse until then.")]
        public static void Radix2Inverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default)
        {
            Inverse(spectrum, options);
        }

        /// <summary>
        /// Bluestein forward FFT for arbitrary sized sample vectors.
        /// </summary>
        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        [Obsolete("Use Forward instead. Will be dropped in version 5.0 and behave like Forward until then.")]
        public static void BluesteinForward(Complex32[] samples, FourierOptions options = FourierOptions.Default)
        {
            Forward(samples, options);
        }

        /// <summary>
        /// Bluestein forward FFT for arbitrary sized sample vectors.
        /// </summary>
        /// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        [Obsolete("Use Forward instead. Will be dropped in version 5.0 and behave like Forward until then.")]
        public static void BluesteinForward(Complex[] samples, FourierOptions options = FourierOptions.Default)
        {
            Forward(samples, options);
        }

        /// <summary>
        /// Bluestein inverse FFT for arbitrary sized sample vectors.
        /// </summary>
        /// <param name="spectrum">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        [Obsolete("Use Inverse instead. Will be dropped in version 5.0 and behave like Inverse until then.")]
        public static void BluesteinInverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default)
        {
            Inverse(spectrum, options);
        }

        /// <summary>
        /// Bluestein inverse FFT for arbitrary sized sample vectors.
        /// </summary>
        /// <param name="spectrum">Sample vector, where the FFT is evaluated in place.</param>
        /// <param name="options">Fourier Transform Convention Options.</param>
        [Obsolete("Use Inverse instead. Will be dropped in version 5.0 and behave like Inverse until then.")]
        public static void BluesteinInverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default)
        {
            Inverse(spectrum, options);
        }

        /// <summary>
        /// Generate the frequencies corresponding to each index in frequency space.
        /// The frequency space has a resolution of sampleRate/N.
        /// Index 0 corresponds to the DC part, the following indices correspond to
        /// the positive frequencies up to the Nyquist frequency (sampleRate/2),
        /// followed by the negative frequencies wrapped around.
        /// </summary>
        /// <param name="length">Number of samples.</param>
        /// <param name="sampleRate">The sampling rate of the time-space data.</param>
        public static double[] FrequencyScale(int length, double sampleRate)
        {
            double[] scale = new double[length];
            double f = 0, step = sampleRate / length;
            int secondHalf = (length >> 1) + 1;
            for (int i = 0; i < secondHalf; i++)
            {
                scale[i] = f;
                f += step;
            }

            f = -step * (secondHalf - 2);
            for (int i = secondHalf; i < length; i++)
            {
                scale[i] = f;
                f += step;
            }

            return scale;
        }
    }
}