using System;
using System.Collections.Generic;
using System.Drawing;

#if !NETSTANDARD2_0
using System.Drawing.Drawing2D;
#endif
using System.Linq;
using System.Text;

namespace HslCommunication.Algorithms.Fourier
{
    /// <summary>
    /// 离散傅氏变换的快速算法，处理的信号，适合单周期信号数为2的N次方个，支持变换及逆变换
    /// </summary>
    public class FFTHelper
    {


        /// <summary>
        /// 
        /// </summary>
        /// <param name="xreal"></param>
        /// <param name="ximag"></param>
        /// <param name="n"></param>
        private static void bitrp( double[] xreal, double[] ximag, int n )
        {
            // 位反转置换 Bit-reversal Permutation
            int i, j, a, b, p;

            for (i = 1, p = 0; i < n; i *= 2)
            {
                p++;
            }
            for (i = 0; i < n; i++)
            {
                a = i;
                b = 0;
                for (j = 0; j < p; j++)
                {
                    b = b * 2 + a % 2;
                    a = a / 2;
                }
                if (b > i)
                {
                    double t = xreal[i];
                    xreal[i] = xreal[b];
                    xreal[b] = t;

                    t = ximag[i];
                    ximag[i] = ximag[b];
                    ximag[b] = t;
                }
            }
        }


        /// <summary>
        /// 快速傅立叶变换
        /// </summary>
        /// <param name="xreal">实数部分</param>
        /// <returns>变换后的数组值</returns>
        public static double[] FFT( double[] xreal )
        {
            return FFT( xreal, new double[xreal.Length] );
        }


#if !NETSTANDARD2_0

        /// <summary>
        /// 获取FFT变换后的显示图形，需要指定图形的相关参数
        /// </summary>
        /// <param name="xreal">实数部分的值</param>
        /// <param name="width">图形的宽度</param>
        /// <param name="heigh">图形的高度</param>
        /// <param name="lineColor">线条颜色</param>
        /// <returns>等待呈现的图形</returns>
        /// <remarks>
        /// <note type="warning">.net standrard2.0 下不支持。</note>
        /// </remarks>
        public static Bitmap GetFFTImage( double[] xreal,int width,int heigh ,Color lineColor)
        {
            double[] ximag = new double[xreal.Length];                // 构造虚对象
            double[] array = FFT( xreal, ximag );                     // 傅立叶变换

            Bitmap bitmap = new Bitmap( width, heigh );               // 构造图形
            Graphics g = Graphics.FromImage( bitmap );
            g.SmoothingMode = SmoothingMode.HighQuality;
            g.TextRenderingHint = System.Drawing.Text.TextRenderingHint.ClearTypeGridFit;
            g.Clear( Color.White );

            Pen pen_Line = new Pen( Color.DimGray, 1 );               // 定义画笔资源
            Pen pen_Dash = new Pen( Color.LightGray, 1 );
            Pen pen_Fourier = new Pen( lineColor, 1 );
            pen_Dash.DashPattern = new float[2] { 5, 5 };
            pen_Dash.DashStyle = DashStyle.Custom;

            Font Font_Normal = SystemFonts.DefaultFont;               // 定义字体资源

            StringFormat sf_right = new StringFormat( );
            sf_right.Alignment = StringAlignment.Far;
            sf_right.LineAlignment = StringAlignment.Center;

            StringFormat sf_center = new StringFormat( );
            sf_center.LineAlignment = StringAlignment.Center;
            sf_center.Alignment = StringAlignment.Center;

            int padding_top = 20;
            int padding_left = 49;
            int padding_right = 49;
            int padding_down = 30;
            int sections = 9;

            // g.DrawLine( pen_Line, new Point( padding_left, padding_top ), new Point( padding_left, heigh - padding_down ) );
            
            float paint_height = heigh - padding_top - padding_down;
            float paint_width = width - padding_left - padding_right;

            if (array.Length > 1)
            {
                double Max = array.Max( );
                double Min = array.Min( );

                Max = Max - Min > 1 ? Max : Min + 1;
                double Length = Max - Min;
                
                //提取峰值
                List<float> Peaks = new List<float>( );
                if (array.Length >= 2)
                {
                    if (array[0] > array[1])
                    {
                        Peaks.Add( 0 );
                    }

                    for (int i = 1; i < array.Length - 2; i++)
                    {
                        if (array[i - 1] < array[i] && array[i] > array[i + 1])
                        {
                            Peaks.Add( i );
                        }
                    }
                    
                    if (array[array.Length - 1] > array[array.Length - 2])
                    {
                        Peaks.Add( array.Length - 1 );
                    }
                }
                

                //高400
                for (int i = 0; i < sections; i++)
                {
                    RectangleF rec = new RectangleF( -10f, (float)i / (sections - 1) * paint_height, padding_left + 8f, 20f );
                    double n = (sections - 1 - i) * Length / (sections - 1) + Min;
                    g.DrawString( n.ToString( "F1" ), Font_Normal, Brushes.Black, rec, sf_right );
                    g.DrawLine( 
                        pen_Dash, padding_left - 3, paint_height * i / (sections - 1) + padding_top,
                        width - padding_right, paint_height * i / (sections - 1) + padding_top );
                }

                float intervalX = paint_width / array.Length;                        // 横向间隔

                for (int i = 0; i < Peaks.Count; i++)
                {
                    if (array[(int)Peaks[i]] * 200 / Max > 1)
                    {
                        g.DrawLine( pen_Dash, Peaks[i] * intervalX + padding_left + 1, padding_top, Peaks[i] * intervalX + padding_left + 1, heigh - padding_down );
                        RectangleF rec = new RectangleF( Peaks[i] * intervalX + padding_left + 1 - 40, heigh - padding_down + 1, 80f, 20f );

                        g.DrawString( Peaks[i].ToString( ), Font_Normal, Brushes.DeepPink, rec, sf_center );
                    }
                }
                
                for (int i = 0; i < array.Length; i++)
                {
                    PointF point = new PointF( );
                    point.X = i * intervalX + padding_left + 1;
                    point.Y = (float)(paint_height - (array[i] - Min) * paint_height / Length + padding_top);

                    PointF point2 = new PointF( );
                    point2.X = i * intervalX + padding_left + 1;
                    point2.Y = (float)(paint_height - (Min - Min) * paint_height / Length + padding_top);
                    

                    g.DrawLine( Pens.Tomato, point, point2 );
                }
            }
            else
            {
                double Max = 100;
                double Min = 0;
                double Length = Max - Min;
                //高400
                for (int i = 0; i < sections; i++)
                {
                    RectangleF rec = new RectangleF( -10f, (float)i / (sections - 1) * paint_height, padding_left + 8f, 20f );
                    double n = (sections - 1 - i) * Length / (sections - 1) + Min;
                    g.DrawString( n.ToString( "F1" ), Font_Normal, Brushes.Black, rec, sf_right );
                    g.DrawLine(
                        pen_Dash, padding_left - 3, paint_height * i / (sections - 1) + padding_top,
                        width - padding_right, paint_height * i / (sections - 1) + padding_top );
                }
            }


            pen_Dash.Dispose( );
            pen_Line.Dispose( );
            pen_Fourier.Dispose( );
            Font_Normal.Dispose( );
            sf_right.Dispose( );
            sf_center.Dispose( );
            g.Dispose( );
            return bitmap;
        }

#endif


        /// <summary>
        /// 快速傅立叶变换
        /// </summary>
        /// <param name="xreal">实数部分，数组长度最好为2的n次方</param>
        /// <param name="ximag">虚数部分，数组长度最好为2的n次方</param>
        /// <returns>变换后的数组值</returns>
        public static double[] FFT( double[] xreal, double[] ximag )
        {
            //n值为2的N次方
            int n = 2;
            while (n <= xreal.Length)
            {
                n *= 2;
            }
            n /= 2;

            // 快速傅立叶变换，将复数 x 变换后仍保存在 x 中，xreal, ximag 分别是 x 的实部和虚部
            double[] wreal = new double[n / 2];
            double[] wimag = new double[n / 2];
            double treal, timag, ureal, uimag, arg;
            int m, k, j, t, index1, index2;

            bitrp( xreal, ximag, n );

            // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
            arg = (-2 * Math.PI / n);
            treal = Math.Cos( arg );
            timag = Math.Sin( arg );
            wreal[0] = 1.0f;
            wimag[0] = 0.0f;
            for (j = 1; j < n / 2; j++)
            {
                wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
                wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
            }

            for (m = 2; m <= n; m *= 2)
            {
                for (k = 0; k < n; k += m)
                {
                    for (j = 0; j < m / 2; j++)
                    {
                        index1 = k + j;
                        index2 = index1 + m / 2;
                        t = n * j / m;    // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                        treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                        timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                        ureal = xreal[index1];
                        uimag = ximag[index1];
                        xreal[index1] = ureal + treal;
                        ximag[index1] = uimag + timag;
                        xreal[index2] = ureal - treal;
                        ximag[index2] = uimag - timag;
                    }
                }
            }

            double[] result = new double[n];
            for (int i = 0; i < result.Length; i++)
            {
                result[i] = Math.Sqrt( Math.Pow( xreal[i], 2 ) + Math.Pow( ximag[i], 2 ) );
            }

            return result;
        }


        /// <summary>
        /// 快速傅立叶变换的逆变换
        /// </summary>
        /// <param name="xreal">实数部分，数组长度最好为2的n次方</param>
        /// <param name="ximag">虚数部分，数组长度最好为2的n次方</param>
        /// <returns>2的多少次方</returns>
        public static int IFFT( double[] xreal, double[] ximag )
        {
            //n值为2的N次方
            int n = 2;
            while (n <= xreal.Length)
            {
                n *= 2;
            }
            n /= 2;

            // 快速傅立叶逆变换
            double[] wreal = new double[n / 2];
            double[] wimag = new double[n / 2];
            double treal, timag, ureal, uimag, arg;
            int m, k, j, t, index1, index2;

            bitrp( xreal, ximag, n );

            // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
            arg = (2 * Math.PI / n);
            treal = (Math.Cos( arg ));
            timag = (Math.Sin( arg ));
            wreal[0] = 1.0f;
            wimag[0] = 0.0f;
            for (j = 1; j < n / 2; j++)
            {
                wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
                wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
            }

            for (m = 2; m <= n; m *= 2)
            {
                for (k = 0; k < n; k += m)
                {
                    for (j = 0; j < m / 2; j++)
                    {
                        index1 = k + j;
                        index2 = index1 + m / 2;
                        t = n * j / m;    // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                        treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                        timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                        ureal = xreal[index1];
                        uimag = ximag[index1];
                        xreal[index1] = ureal + treal;
                        ximag[index1] = uimag + timag;
                        xreal[index2] = ureal - treal;
                        ximag[index2] = uimag - timag;
                    }
                }
            }

            for (j = 0; j < n; j++)
            {
                xreal[j] /= n;
                ximag[j] /= n;
            }

            return n;
        }
    }
}