klhztrader/KLHZ.Trader.Core.Math/Declisions/Utils/FFT.cs

using KLHZ.Trader.Core.Math.Declisions.Dtos;
using MathNet.Numerics;
using MathNet.Numerics.IntegralTransforms;

namespace KLHZ.Trader.Core.Math.Declisions.Utils
{
    public static class FFT
    {
        public static void Test()
        {
            var da = new List<float>();
            for (int i = 0; i < 1000; i++)
            {
                da.Add((float)System.Math.Sin(0.01 * i) + (float)System.Math.Cos(0.01 * i));
            }

            var start = da.ToArray();
            var arrv = da.Select(d => new Complex32(d, 0)).ToArray();
            Fourier.Forward(arrv);

            Fourier.Inverse(arrv);
            var res = arrv.Select(a => a.Real).ToArray();

            for (int i = 0; i < 1000; i++)
            {
                var d = res[i] - start[i];
            }
        }

        public static TimeSpan GetMainHarmonictPeriod(double[] data, TimeSpan period)
        {
            var arrv = data.Select(d => new Complex32((float)d, 0)).ToArray();
            Fourier.Forward(arrv);
            var res = new List<Harmonic>();
            var times = new TimeSpan[arrv.Length];
            var mags = arrv.Select(a => a.Magnitude).ToArray();
            var phases = arrv.Select(a => a.Phase).ToArray();
            var max = 0f;
            var kmax = 1;
            var sum = 0f;
            for (int i=0;i< arrv.Length / 2; i++)
            {
                if (i == 0)
                {
                    res.Add(new Harmonic()
                    {
                        Magnitude = arrv[i].Magnitude,
                        Period = TimeSpan.MaxValue,
                        Phase = arrv[i].Phase,
                    });
                }
                else
                {
                    times[i] = CaclHarmonycPeriod(period, data.Length, i);
                    res.Add(new Harmonic()
                    {
                        Magnitude = arrv[i].Magnitude,
                        Period = times[i],
                        Phase = arrv[i].Phase,
                    });
                    var mag = arrv[i].Magnitude;
                    sum += mag;
                    if (mag > max)
                    {
                        max = mag;
                        kmax = i;
                    }
                }
            }
            return CaclHarmonycPeriod(period, data.Length, kmax);
        }

        public static Harmonic[] GetHarmonics(double[] data, TimeSpan period, int count)
        {
            var arrv = data.Select(d => new Complex32((float)d, 0)).ToArray();
            Fourier.Forward(arrv);
            var res = new List<Harmonic>();

            for (int i = 1; i < count; i++)
            {
                res.Add(new Harmonic()
                {
                    Imaginary = arrv[i].Imaginary,
                    Real = arrv[i].Real,
                    Magnitude = arrv[i].Magnitude,
                    Period = CaclHarmonycPeriod(period, data.Length, i),
                    Phase = arrv[i].Phase,
                });
            }

            return res.ToArray();
        }

        public static double Calc(Harmonic[] harmonics, DateTime startTime, DateTime endTime, DateTime currentTime)
        {
            var sum = 0d;
            var timeSpan = currentTime - startTime;
            for (int i = 0; i < harmonics.Length; i++)
            {
                var currentPhase = System.Math.PI * 2 * timeSpan.TotalSeconds * (endTime - startTime).TotalSeconds / harmonics[i].Period.TotalSeconds / harmonics[i].Period.TotalSeconds + harmonics[i].Phase;
                var value = harmonics[i].Real * System.Math.Cos(currentPhase) + harmonics[i].Imaginary * System.Math.Sign(currentPhase);
                sum += value;
            }
            return sum;
        }

        public static TimeSpan CaclHarmonycPeriod(TimeSpan signalLength, int signalLengthItems, int harmonyNumber)
        {
            var fdiscretisation = signalLengthItems / signalLength.TotalSeconds;
            var fharm = harmonyNumber * fdiscretisation / signalLengthItems;
            return TimeSpan.FromSeconds(1 / fharm);
        }

        public static double CalcCurrentPhase(Complex32[] spectr, int harmonyNumber)
        {
            var item = spectr[harmonyNumber];
            return System.Math.Atan(item.Imaginary / item.Real);
        }

        public static (double min, double max) CalcPhaseRangeFoxMax(double aSin, double aCos, double initPhase, double level)
        {
            return CalcPhaseRange(aSin, aCos,initPhase, level, CheckMaxValue);
        }

        public static (double min, double max) CalcPhaseRangeFoxMin(double aSin, double aCos, double initPhase, double level)
        {
            return CalcPhaseRange(aSin, aCos, initPhase, level, CheckMinValue);
        }

        internal static (double min, double max) CalcPhaseRange(double aSin, double aCos, double initPhase, double level, Func<double, double, double, double,bool> comparer)
        {
            var x = new List<double>();
            var xIndexes = new List<int>();
            var y = new List<double>();
            var max = double.MinValue;
            var min = double.MaxValue;
            var _2pi = 2 * System.Math.PI;
            for (double i = 0; i <= 10 * System.Math.PI; i += 0.01)
            {
                var df = (((i) / _2pi) % 1) * _2pi;
                var val = aSin * System.Math.Sin(i + initPhase) + aCos * System.Math.Cos(i + initPhase);
                if (val > max)
                {
                    max = val;
                }
                if (val < min)
                {
                    min = val;
                }
                x.Add(df);
                y.Add(val);
            }

            int start = -2;
            int prevI = -2;
            int end = -2;

            var drange = -2;

            var minPhaseTmp = 0d;
            var maxPhaseTmp = 0d;

            var minPhase = 0d;
            var maxPhase = 0d;
            for (int i = 0; i < y.Count; i++)
            {
                if (comparer(y[i],min,max,level))
                {
                    if (start < 0)
                    {
                        minPhaseTmp = x[i];
                        start = i;
                    }
                }
                else if (start >= 0 && i - prevI == 1)
                {
                    end = prevI;
                    maxPhaseTmp = x[end];
                    var drangeTmp = end - start;
                    if (drangeTmp > drange)
                    {
                        drange = drangeTmp;
                        minPhase = minPhaseTmp;
                        maxPhase = maxPhaseTmp;
                    }
                    start = -2;
                }
                prevI = i;
            }


            return (minPhase, maxPhase);
        }

        internal static bool CheckMaxValue(double val, double min, double max, double relativeValue)
        {
            return (val > (1 - relativeValue) * max);
        }

        internal static bool CheckMinValue(double val, double min, double max, double relativeValue)
        {
            return (val < (1 - relativeValue) * min);
        }
    }
}