using System.Runtime.CompilerServices; namespace QuanTAlib; /// /// AFIRMA: Adaptive FIR Moving Average /// A finite impulse response (FIR) filter that combines windowing functions with sinc-based filtering. /// Provides superior noise reduction while maintaining signal fidelity through adaptive filtering. /// /// /// Implementation: /// Original implementation based on FIR filter design principles /// public class Afirma : AbstractBase { public enum WindowType { Rectangular, Hanning1, Hanning2, Blackman, BlackmanHarris } private readonly int Periods; private readonly int Taps; private readonly WindowType Window; private readonly CircularBuffer _buffer; private readonly double[] _weights; private readonly double _wsum; private readonly double[] _armaBuffer; private readonly int _n; private readonly double _sx2, _sx3, _sx4, _sx5, _sx6, _den; private readonly double _twoPi = 2.0 * Math.PI; private readonly double _fourPi = 4.0 * Math.PI; private readonly double _sixPi = 6.0 * Math.PI; /// The number of periods for the sinc filter calculation. /// The number of filter taps (filter length). Must be odd number. /// The type of window function to apply (Rectangular, Hanning1, Hanning2, Blackman, or BlackmanHarris). /// Thrown when periods or taps is less than 1. public Afirma(int periods, int taps, WindowType window) { if (periods < 1) { throw new ArgumentOutOfRangeException(nameof(periods), "Periods must be greater than or equal to 1."); } if (taps < 1) { throw new ArgumentOutOfRangeException(nameof(taps), "Taps must be greater than or equal to 1."); } Periods = periods; Taps = taps; Window = window; WarmupPeriod = taps; _buffer = new CircularBuffer(taps); _weights = new double[taps]; _wsum = CalculateWeights(); _armaBuffer = new double[taps]; _n = (Taps - 1) / 2; // Precalculate least squares coefficients _sx2 = ((2 * _n) + 1) / 3.0; _sx3 = _n * (_n + 1) / 2.0; _sx4 = _sx2 * ((3 * _n * _n) + (3 * _n) - 1) / 5.0; _sx5 = _sx3 * ((2 * _n * _n) + (2 * _n) - 1) / 3.0; _sx6 = _sx2 * ((3 * Math.Pow(_n, 3) * (_n + 2)) - (3 * _n) + 1) / 7.0; _den = (_sx6 * _sx4 / _sx5) - _sx5; Name = "Afirma"; Init(); } /// The data source object that publishes updates. /// The number of periods for the sinc filter calculation. /// The number of filter taps (filter length). Must be odd number. /// The type of window function to apply (Rectangular, Hanning1, Hanning2, Blackman, or BlackmanHarris). public Afirma(object source, int periods, int taps, WindowType window) : this(periods, taps, window) { var pubEvent = source.GetType().GetEvent("Pub"); pubEvent?.AddEventHandler(source, new ValueSignal(Sub)); } [MethodImpl(MethodImplOptions.AggressiveInlining)] protected override void ManageState(bool isNew) { if (isNew) { _lastValidValue = Input.Value; _index++; } } [MethodImpl(MethodImplOptions.AggressiveInlining)] private static double CalculateSincWeight(double x) { return Math.Abs(x) < 1e-10 ? 1.0 : Math.Sin(x) / x; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private double GetWindowWeight(int k, int tapsMinusOne) { switch (Window) { case WindowType.Rectangular: return 1.0; case WindowType.Hanning1: return 0.50 - (0.50 * Math.Cos(_twoPi * k / tapsMinusOne)); case WindowType.Hanning2: return 0.54 - (0.46 * Math.Cos(_twoPi * k / tapsMinusOne)); case WindowType.Blackman: return 0.42 - (0.50 * Math.Cos(_twoPi * k / tapsMinusOne)) + (0.08 * Math.Cos(_fourPi * k / tapsMinusOne)); case WindowType.BlackmanHarris: return 0.35875 - (0.48829 * Math.Cos(_twoPi * k / tapsMinusOne)) + (0.14128 * Math.Cos(_fourPi * k / tapsMinusOne)) - (0.01168 * Math.Cos(_sixPi * k / tapsMinusOne)); default: return 1.0; } } protected override double Calculation() { ManageState(IsNew); _buffer.Add(Input.Value, Input.IsNew); if (_index >= Taps) { CalculateAdaptiveCoefficients(); } double result = 0.0; for (int k = 0; k < Taps; k++) { result += _buffer[k] * _weights[k]; } IsHot = _index >= WarmupPeriod; return result / _wsum; } [MethodImpl(MethodImplOptions.AggressiveInlining)] private void CalculateAdaptiveCoefficients() { double a0 = _buffer[_n]; double a1 = _buffer[_n] - _buffer[_n + 1]; double sx2y = 0.0; double sx3y = 0.0; for (int i = 0; i <= _n; i++) { double i2 = i * i; sx2y += i2 * _buffer[_n - i]; sx3y += i2 * i * _buffer[_n - i]; } sx2y = 2.0 * sx2y / _n / (_n + 1); sx3y = 2.0 * sx3y / _n / (_n + 1); double p = sx2y - (a0 * _sx2) - (a1 * _sx3); double q = sx3y - (a0 * _sx3) - (a1 * _sx4); double a2 = ((p * _sx6 / _sx5) - q) / _den; double a3 = ((q * _sx4 / _sx5) - p) / _den; for (int k = 0; k <= _n; k++) { double k2 = k * k; _armaBuffer[_n - k] = a0 + (k * a1) + (k2 * a2) + (k2 * k * a3); } } private double CalculateWeights() { double wsum = 0.0; double centerTap = (Taps - 1) / 2.0; int tapsMinusOne = Taps - 1; for (int k = 0; k < Taps; k++) { double windowWeight = GetWindowWeight(k, tapsMinusOne); double x = Math.PI * (k - centerTap) / Periods; double sincWeight = CalculateSincWeight(x); _weights[k] = windowWeight * sincWeight; wsum += _weights[k]; } return wsum; } }