|
| 1 | +using MathCore.DSP.Filters; |
| 2 | + |
| 3 | +namespace MathCore.DSP.Tests.Filters; |
| 4 | + |
| 5 | +[TestClass] |
| 6 | +[Ignore("Диагностические стресс-тесты устойчивости; включить после внесения исправлений численной устойчивости")] |
| 7 | +public class IirCoefficientsStabilityTests |
| 8 | +{ |
| 9 | + private sealed record StressCase(string Name, double Dt, Func<DSP.Filters.IIR> Factory); |
| 10 | + |
| 11 | + [TestMethod] |
| 12 | + public void CalculatedCoefficients_ShouldRemainFinite_OnStressSpecifications() => |
| 13 | + RunStressValidation((case_name, filter, dt) => AssertCoefficientsAreFinite(case_name, filter)); |
| 14 | + |
| 15 | + [TestMethod] |
| 16 | + public void FrequencyResponse_ShouldRemainFinite_OnStressSpecifications() => |
| 17 | + RunStressValidation(AssertFrequencyResponseIsFinite); |
| 18 | + |
| 19 | + [TestMethod] |
| 20 | + public void ImpulseResponse_ShouldRemainFinite_OnStressSpecifications() => |
| 21 | + RunStressValidation((case_name, filter, dt) => AssertImpulseResponseIsFinite(case_name, filter)); |
| 22 | + |
| 23 | + private static void RunStressValidation(Action<string, DSP.Filters.IIR, double> Validator) |
| 24 | + { |
| 25 | + ArgumentNullException.ThrowIfNull(Validator); |
| 26 | + |
| 27 | + var errors = new List<string>(); |
| 28 | + |
| 29 | + foreach (var @case in GetStressCases()) |
| 30 | + { |
| 31 | + try |
| 32 | + { |
| 33 | + var filter = @case.Factory(); |
| 34 | + Validator(@case.Name, filter, @case.Dt); |
| 35 | + } |
| 36 | + catch (Exception error) |
| 37 | + { |
| 38 | + errors.Add($"{@case.Name}: {error.GetType().Name} - {error.Message}"); |
| 39 | + } |
| 40 | + } |
| 41 | + |
| 42 | + if (errors.Count > 0) |
| 43 | + Assert.Fail($"Обнаружены проблемы устойчивости:{Environment.NewLine}{string.Join(Environment.NewLine, errors)}"); |
| 44 | + } |
| 45 | + |
| 46 | + private static IReadOnlyList<StressCase> GetStressCases() |
| 47 | + { |
| 48 | + const double dt = 1d / 48_000; |
| 49 | + const double gp = 0.891250938; |
| 50 | + const double gs = 0.01; |
| 51 | + |
| 52 | + return |
| 53 | + [ |
| 54 | + new("Butterworth-LP", dt, () => new DSP.Filters.ButterworthLowPass(dt, 9_000, 9_300, gp, gs)), |
| 55 | + new("Butterworth-HP", dt, () => new DSP.Filters.ButterworthHighPass(dt, 6_000, 6_300, gp, gs)), |
| 56 | + new("Butterworth-BP", dt, () => new DSP.Filters.ButterworthBandPass(dt, 2_000, 2_300, 5_200, 5_500, gp, gs)), |
| 57 | + new("Butterworth-BS", dt, () => new DSP.Filters.ButterworthBandStop(dt, 1_800, 2_200, 4_800, 5_300, gp, gs)), |
| 58 | + |
| 59 | + new("ChebyshevI-LP", dt, () => new DSP.Filters.ChebyshevLowPass(dt, 9_000, 9_250, gp, gs, ChebyshevType.I)), |
| 60 | + new("ChebyshevII-LP", dt, () => new DSP.Filters.ChebyshevLowPass(dt, 9_000, 9_250, gp, gs, ChebyshevType.II)), |
| 61 | + new("ChebyshevIIC-LP", dt, () => new DSP.Filters.ChebyshevLowPass(dt, 9_000, 9_250, gp, gs, ChebyshevType.IICorrected)), |
| 62 | + |
| 63 | + new("ChebyshevI-HP", dt, () => new DSP.Filters.ChebyshevHighPass(dt, 6_000, 6_200, gp, gs, ChebyshevType.I)), |
| 64 | + new("ChebyshevII-HP", dt, () => new DSP.Filters.ChebyshevHighPass(dt, 6_000, 6_200, gp, gs, ChebyshevType.II)), |
| 65 | + new("ChebyshevIIC-HP", dt, () => new DSP.Filters.ChebyshevHighPass(dt, 6_000, 6_200, gp, gs, ChebyshevType.IICorrected)), |
| 66 | + |
| 67 | + new("ChebyshevI-BP", dt, () => new DSP.Filters.ChebyshevBandPass(dt, 2_000, 2_300, 5_200, 5_500, gp, gs, ChebyshevType.I)), |
| 68 | + new("ChebyshevII-BP", dt, () => new DSP.Filters.ChebyshevBandPass(dt, 2_000, 2_300, 5_200, 5_500, gp, gs, ChebyshevType.II)), |
| 69 | + new("ChebyshevIIC-BP", dt, () => new DSP.Filters.ChebyshevBandPass(dt, 2_000, 2_300, 5_200, 5_500, gp, gs, ChebyshevType.IICorrected)), |
| 70 | + |
| 71 | + new("ChebyshevI-BS", dt, () => new DSP.Filters.ChebyshevBandStop(dt, 1_800, 2_200, 4_800, 5_300, gp, gs, ChebyshevType.I)), |
| 72 | + new("ChebyshevII-BS", dt, () => new DSP.Filters.ChebyshevBandStop(dt, 1_800, 2_200, 4_800, 5_300, gp, gs, ChebyshevType.II)), |
| 73 | + new("ChebyshevIIC-BS", dt, () => new DSP.Filters.ChebyshevBandStop(dt, 1_800, 2_200, 4_800, 5_300, gp, gs, ChebyshevType.IICorrected)), |
| 74 | + |
| 75 | + new("Elliptic-LP", dt, () => new DSP.Filters.EllipticLowPass(dt, 9_000, 9_200, gp, gs)), |
| 76 | + new("Elliptic-HP", dt, () => new DSP.Filters.EllipticHighPass(dt, 6_000, 6_200, gp, gs)), |
| 77 | + new("Elliptic-BP", dt, () => new DSP.Filters.EllipticBandPass(dt, 2_000, 2_300, 5_200, 5_500, gp, gs)), |
| 78 | + new("Elliptic-BS", dt, () => new DSP.Filters.EllipticBandStop(dt, 1_800, 2_200, 4_800, 5_300, gp, gs)) |
| 79 | + ]; |
| 80 | + } |
| 81 | + |
| 82 | + private static void AssertCoefficientsAreFinite(string CaseName, DSP.Filters.IIR Filter) |
| 83 | + { |
| 84 | + for (var i = 0; i < Filter.A.Count; i++) |
| 85 | + Assert.IsTrue(double.IsFinite(Filter.A[i]), $"{CaseName}: A[{i}] не является конечным числом"); |
| 86 | + |
| 87 | + for (var i = 0; i < Filter.B.Count; i++) |
| 88 | + Assert.IsTrue(double.IsFinite(Filter.B[i]), $"{CaseName}: B[{i}] не является конечным числом"); |
| 89 | + |
| 90 | + Assert.IsTrue(double.IsFinite(Filter.A[0]), $"{CaseName}: A[0] не является конечным числом"); |
| 91 | + Assert.IsTrue(Math.Abs(Filter.A[0]) > 1e-30, $"{CaseName}: A[0] слишком близок к нулю"); |
| 92 | + } |
| 93 | + |
| 94 | + private static void AssertFrequencyResponseIsFinite(string CaseName, DSP.Filters.IIR Filter, double Dt) |
| 95 | + { |
| 96 | + var fd = 1d / Dt; |
| 97 | + const int points_count = 256; |
| 98 | + |
| 99 | + for (var i = 0; i <= points_count; i++) |
| 100 | + { |
| 101 | + var f = fd * i / (2d * points_count); |
| 102 | + var h = Filter.FrequencyResponse(f); |
| 103 | + |
| 104 | + Assert.IsTrue(double.IsFinite(h.Re), $"{CaseName}: Re(H) не является конечным числом на частоте {f}"); |
| 105 | + Assert.IsTrue(double.IsFinite(h.Im), $"{CaseName}: Im(H) не является конечным числом на частоте {f}"); |
| 106 | + Assert.IsTrue(double.IsFinite(h.Abs), $"{CaseName}: |H| не является конечным числом на частоте {f}"); |
| 107 | + } |
| 108 | + } |
| 109 | + |
| 110 | + private static void AssertImpulseResponseIsFinite(string CaseName, DSP.Filters.IIR Filter) |
| 111 | + { |
| 112 | + const int samples_count = 2048; |
| 113 | + const double max_abs_limit = 1e12; |
| 114 | + |
| 115 | + Filter.Reset(); |
| 116 | + |
| 117 | + var max_abs = 0d; |
| 118 | + for (var i = 0; i < samples_count; i++) |
| 119 | + { |
| 120 | + var x = i == 0 ? 1d : 0d; |
| 121 | + var y = Filter.Process(x); |
| 122 | + |
| 123 | + Assert.IsTrue(double.IsFinite(y), $"{CaseName}: значение импульсной характеристики не является конечным на отсчёте {i}"); |
| 124 | + |
| 125 | + var abs_y = Math.Abs(y); |
| 126 | + if (abs_y > max_abs) max_abs = abs_y; |
| 127 | + } |
| 128 | + |
| 129 | + Assert.IsTrue(max_abs < max_abs_limit, $"{CaseName}: наблюдается неустойчивый рост импульсной характеристики, max={max_abs}"); |
| 130 | + } |
| 131 | +} |
0 commit comments