using CommonLib.math;
using System;
using System.Collections.Generic;
using System.Text;

namespace CommonLib.math
{
    /// <summary>
    /// 快速傅里叶变换
    /// </summary>
    public class FFT
    {
        /// <summary>
        /// 获取FFT数列
        /// </summary>
        /// <param name="chdata">数字信号数组</param>
        /// <param name="result_lenth">获取结果的长度，数字信号的数组长度需要不小于该值的两倍，否则将无法正常计算</param>
        /// <returns></returns>
        public double[] GetFftValueFromChData(short[] chdata, int result_lenth)
        {
            double[] result = new double[0];
            double[] allvalue = new double[0];//快速计算结果为双倍长度
            int n = result_lenth * 2;// 470;// this.allChanData.Length;
            short[] x = chdata;//被计算的数组 
            complex[] y = new complex[n];//接收复数结果的数组 
            if (chdata.Length >= n)
            {
                result = new Double[n];//接收幅值结果的数组 
                y = airthm.dft(x, n);//重点耗时项
                allvalue = airthm.amplitude(y, n);

                result = new double[allvalue.Length / 2];
                for (int i = 0; i < result.Length; i++)
                { 
                    result[i] = allvalue[i]; 
                }
              
            }
            else { }
            return result;
        }
        //double [,]X_sn;
        double pi = System.Math.PI;
        public double[,] Wcreat(int N, int FFT_IFFT_elect)
        {
            double[,] Wp = new double[2, N / 2];
            if (FFT_IFFT_elect == 0)
            {
                for (int i = 0; i < N / 2; i++)
                {
                    Wp[0, i] = System.Math.Cos(2 * pi / N * i);
                    Wp[1, i] = System.Math.Sin(-2 * pi / N * i);
                }
            }
            else
            {
                for (int i = 0; i < N / 2; i++)
                {
                    Wp[0, i] = System.Math.Cos(2 * pi / N * i);
                    Wp[1, i] = System.Math.Sin(2 * pi / N * i);
                }
            }
            return Wp;
        }
        /// <summary>
        /// 进行傅里叶变换
        /// </summary>
        /// <param name="X_sn">数列/原始信号</param>
        /// <param name="N">数列的长度</param>
        /// <param name="FFT_IFFT_elect"></param>
        /// <returns></returns>
        public double[,] FFT_T(double[,] X_sn, int N, int FFT_IFFT_elect)
        {
            double[,] Wp = new double[2, N / 2];
            FFT Xn = new FFT();
            Wp = Xn.Wcreat(N, FFT_IFFT_elect);
            //测试


            //
            double tem = System.Math.Log(N, 2);
            int M = (int)tem;
            for (int L = 1; L <= M; L++)
            {
                double M_Ld = System.Math.Pow(2, M - L);//计算2的M-L次方
                int M_L = (int)M_Ld;
                double L_1d = System.Math.Pow(2, L - 1);
                int L_1 = (int)L_1d;
                for (int j = 0; j < M_L; j++)
                {
                    int J = j * (int)System.Math.Pow(2, L);
                    for (int k = 0; k < L_1; k++)
                    {
                        double[] T = new double[2];
                        double p_k = k * (int)System.Math.Pow(2, (M - L));
                        int P = (int)p_k;
                        T = Xn.C_multi(X_sn[0, J + L_1 + k], X_sn[1, J + L_1 + k], Wp[0, P], Wp[1, P]);
                        X_sn[0, J + L_1 + k] = X_sn[0, J + k] - T[0];
                        X_sn[1, J + L_1 + k] = X_sn[1, J + k] - T[1];
                        X_sn[0, J + k] = X_sn[0, J + k] + T[0];
                        X_sn[1, J + k] = X_sn[1, J + k] + T[1];
                    }
                }
            }
            return X_sn;
        }
        public double[] C_multi(double a, double b, double c, double d)
        {
            double[] R_value = new double[2];
            R_value[0] = a * c - b * d;
            R_value[1] = a * d + b * c;
            return R_value;
        }
        public int R_value(int num, int M)//将一个数取反
        {
            double R_num = 0;
            for (int i = M; i > 0; i--)
            {
                double tem = System.Math.Pow(2, i);
                int t = (int)tem;
                int j = (num % t) / (t / 2);
                R_num = R_num + j * System.Math.Pow(2, M - i);
            }
            return (int)R_num;
        }
        public double[,] R_Xn(double[,] xn, int N)
        {
            double tem = System.Math.Log(N, 2);//求得序列点数对2为底的对数 例：N=4则tem=2，N=8则tem=3；
            int M = (int)tem;//将该对数取整
            FFT R_fft = new FFT();
            for (int i = 0; i < N / 2; i++)
            {
                int tem1 = R_fft.R_value(i, M);//把索引各种倒换，得出一个换位置的索引，还没看懂
                double[] tem2 = new double[2];
                if (tem1 != i)//防止计算出的索引值与i重合造成计算的浪费
                {
                    tem2[0] = xn[0, i];//将索引位置的值互换
                    tem2[1] = xn[1, i];
                    xn[0, i] = xn[0, tem1];
                    xn[1, i] = xn[1, tem1];
                    xn[0, tem1] = tem2[0];
                    xn[1, tem1] = tem2[1];
                }
            }
            return xn;
        }
    }
}