Улучшена устойчивость расчёта биномиальных коэффициентов

Добавлен метод GetBinomialCoefficients(int N) для безопасного вычисления биномиальных коэффициентов в формате double без переполнения. Переведены фильтры Butterworth и Chebyshev на использование этого метода при формировании коэффициентов B, что повысило численную устойчивость. В ChebyshevLowPass.cs заменён LINQ на явный цикл. В тестах снят [Ignore] с класса, но добавлен [Ignore] к тесту импульсной характеристики для высоких порядков (требуется SOS-реализация).

Author: Infarh

MathCore.DSP/Filters/AnalogBasedFilter.cs

@@ -145,6 +145,23 @@ public static (Complex[] Poles, Complex[] Zeros) ToBandStop(
         }
     }
 
+    /// <summary>Получить биномиальные коэффициенты C(N,k) в вещественном формате без целочисленного переполнения</summary>
+    /// <param name="N">Порядок бинома</param>
+    /// <returns>Массив коэффициентов C(N,0..N)</returns>
+    /// <exception cref="ArgumentOutOfRangeException">Если N меньше 0</exception>
+    protected static double[] GetBinomialCoefficients(int N)
+    {
+        if (N < 0) throw new ArgumentOutOfRangeException(nameof(N), N, "Порядок бинома не может быть отрицательным");
+
+        var coefficients = new double[N + 1];
+        coefficients[0] = 1;
+
+        for (var k = 1; k <= N; k++)
+            coefficients[k] = coefficients[k - 1] * (N - (k - 1d)) / k;
+
+        return coefficients;
+    }
+
     /// <summary>Спецификация фильтра</summary>
     public readonly ref struct Specification
     {

MathCore.DSP/Filters/ButterworthHighPass.cs

@@ -26,8 +26,9 @@ public static (double[] A, double[] B) GetPolynoms(Specification Spec)
         var g_norm = z_poles.Multiply(z => (1 + z) / 2).Abs;
 
         var B = new double[N + 1];
+        var binomial = GetBinomialCoefficients(N);
         for (var i = 0; i < B.Length; i++)
-            B[i] = BinomialCoefficient(N, i) * (i % 2 == 0 ? g_norm : -g_norm);
+            B[i] = binomial[i] * (i % 2 == 0 ? g_norm : -g_norm);
         var A = GetCoefficientsInverted(z_poles).ToRe();
 
         return (A, B);

MathCore.DSP/Filters/ButterworthLowPass.cs

@@ -163,8 +163,9 @@ public static (double[] A, double[] B) GetPolynoms(Specification Spec)
         var g_norm = z_poles.Multiply(z => (1 - z) / 2).Re;
 
         var B = new double[N + 1];
+        var binomial = GetBinomialCoefficients(N);
         for (var i = 0; i < B.Length; i++)
-            B[i] = g_norm * BinomialCoefficient(N, i);
+            B[i] = g_norm * binomial[i];
 
         var A = GetCoefficientsInverted(z_poles).ToRe();
 

MathCore.DSP/Filters/ChebyshevHighPass.cs

@@ -27,8 +27,9 @@ private static (double[] A, double[] B) InitializeI(Specification Spec)
             : 1 * k_poles;
 
         var B = new double[N + 1];
+        var binomial = GetBinomialCoefficients(N);
         for (var i = 0; i < B.Length; i++)
-            B[i] = SpecialFunctions.BinomialCoefficient(N, i) * (i % 2 == 0 ? g_norm : -g_norm);
+            B[i] = binomial[i] * (i % 2 == 0 ? g_norm : -g_norm);
         var A = Polynom.Array.GetCoefficientsInverted(z_poles).ToRe();
 
         return (A, B);

MathCore.DSP/Filters/ChebyshevLowPass.cs

@@ -113,9 +113,10 @@ private static (double[] A, double[] B) InitializeI(Specification Spec)
 
         var g_norm = (N.IsOdd() ? 1 : Spec.Gp) / (2.Pow(N) / z_poles.Multiply(z => 1 - z).Re);
 
-        var B = Enumerable
-           .Range(0, N + 1)
-           .ToArray(i => g_norm * BinomialCoefficient(N, i));
+        var binomial = GetBinomialCoefficients(N);
+        var B = new double[N + 1];
+        for (var i = 0; i < B.Length; i++)
+            B[i] = g_norm * binomial[i];
 
         return (A, B);
     }

Tests/MathCore.DSP.Tests/Filters/IirCoefficientsStabilityTests.cs

@@ -3,7 +3,6 @@
 namespace MathCore.DSP.Tests.Filters;
 
 [TestClass]
-[Ignore("Диагностические стресс-тесты устойчивости; включить после внесения исправлений численной устойчивости")]
 public class IirCoefficientsStabilityTests
 {
     private sealed record StressCase(string Name, double Dt, Func<DSP.Filters.IIR> Factory);
@@ -17,6 +16,7 @@ public void FrequencyResponse_ShouldRemainFinite_OnStressSpecifications() =>
         RunStressValidation(AssertFrequencyResponseIsFinite);
 
     [TestMethod]
+    [Ignore("Диагностический тест: для высоких порядков требуется секционная реализация (SOS) вместо прямой формы")]
     public void ImpulseResponse_ShouldRemainFinite_OnStressSpecifications() =>
         RunStressValidation((case_name, filter, dt) => AssertImpulseResponseIsFinite(case_name, filter));