Skip to content

Commit 6f715fc

Browse files
committed
Оптимизация метода PhaseDemodulation
Изменения в методе `PhaseDemodulation` в классе `SampleSI16Ex.cs` включают оптимизацию алгоритма фазовой демодуляции, что позволяет уменьшить количество вызовов функции `atan2` и избежать создания промежуточного массива фаз. Вместо этого используется прямое вычисление разности фаз через комплексное произведение, что улучшает производительность и экономит память. Также добавлена функция `UnwrapSinglePhaseDiff` для упрощенного разворачивания разностей фаз с накоплением коррекции. Обновлена документация в файле `PhaseDemodulation.md`, чтобы отразить изменения в алгоритме и его оптимизации, включая описание ключевых улучшений и математического обоснования. Добавлены новые тесты в `SampleSI16ExTests.cs`, которые проверяют корректность работы оптимизированного метода фазовой демодуляции, включая тесты на пустые массивы, постоянные сигналы, синусоидальные сигналы и производительность. Тесты также сравнивают точность оптимизированной реализации с наивной версией.
1 parent 5c04cdb commit 6f715fc

3 files changed

Lines changed: 363 additions & 64 deletions

File tree

MathCore.DSP/Samples/Extensions/SampleSI16Ex.cs

Lines changed: 48 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -16,65 +16,65 @@ public static float[] PhaseDemodulation(this Span<SampleSI16> samples, double f0
1616

1717
var result = new float[samples.Length];
1818
result[0] = 0f; // Первый отсчёт всегда 0, так как нет предыдущего для вычисления производной
19-
20-
// Вычисляем массив фаз
21-
var phases = new float[samples.Length];
22-
for (var i = 0; i < samples.Length; i++)
23-
{
24-
phases[i] = GetPhase(samples[i]);
25-
}
26-
27-
// Разворачиваем фазы (phase unwrapping)
28-
UnwrapPhases(phases);
29-
30-
// Вычисляем мгновенную частоту как производную фазы
31-
var dt = 1.0 / fd;
19+
3220
var scale_factor = (float)(fd / (2.0 * Math.PI));
33-
21+
var f0_offset = (float)f0;
22+
var unwrapped_correction = 0f;
23+
3424
for (var i = 1; i < samples.Length; i++)
3525
{
36-
// Производная фазы дает мгновенную частоту в рад/с
37-
// Делим на 2π для перевода в Гц
38-
var instantaneous_frequency = (phases[i] - phases[i - 1]) * scale_factor;
39-
40-
// Вычитаем центральную частоту f0
41-
result[i] = instantaneous_frequency - (float)f0;
26+
// Вычисляем разность фаз напрямую через комплексное произведение
27+
// Δφ = arg(z_curr * z_prev*) = atan2(Im(z_curr * z_prev*), Re(z_curr * z_prev*))
28+
var curr = samples[i];
29+
var prev = samples[i - 1];
30+
31+
// z_curr * z_prev* = (I1 + jQ1) * (I2 - jQ2) = (I1*I2 + Q1*Q2) + j(Q1*I2 - I1*Q2)
32+
var real_part = curr.I * prev.I + curr.Q * prev.Q;
33+
var imag_part = curr.Q * prev.I - curr.I * prev.Q;
34+
35+
// Разность фаз с одним вызовом atan2
36+
#if NET8_0_OR_GREATER
37+
var phase_diff = MathF.Atan2(imag_part, real_part);
38+
#else
39+
var phase_diff = (float)Math.Atan2(imag_part, real_part);
40+
#endif
41+
42+
// Unwrapping - приводим к диапазону (-π, π] и накапливаем коррекцию
43+
phase_diff = UnwrapSinglePhaseDiff(phase_diff, ref unwrapped_correction);
44+
45+
// Мгновенная частота = производная фазы / (2π)
46+
var instantaneous_frequency = phase_diff * scale_factor;
47+
48+
// Вычитаем центральную частоту
49+
result[i] = instantaneous_frequency - f0_offset;
4250
}
43-
51+
4452
return result;
4553
}
4654

47-
/// <summary>Быстрое вычисление фазы с использованием статических функций SampleSI16</summary>
55+
/// <summary>Unwrapping одной разности фаз с накоплением коррекции</summary>
4856
[MethodImpl(MethodImplOptions.AggressiveInlining)]
49-
private static float GetPhase(SampleSI16 sample)
50-
{
51-
// Используем уже оптимизированную функцию GetArg из SampleSI16
52-
return SampleSI16.GetArg(sample);
53-
}
54-
55-
/// <summary>Разворачивание фаз - устранение скачков ±2π</summary>
56-
private static void UnwrapPhases(Span<float> phases)
57+
private static float UnwrapSinglePhaseDiff(float phase_diff, ref float accumulated_correction)
5758
{
58-
if (phases.Length <= 1) return;
59-
60-
const float two_pi = (float)(2.0 * Math.PI);
59+
#if NET8_0_OR_GREATER
60+
const float pi = MathF.PI;
61+
#else
6162
const float pi = (float)Math.PI;
62-
63-
for (var i = 1; i < phases.Length; i++)
63+
#endif
64+
const float two_pi = 2f * pi;
65+
66+
// Приводим разность к диапазону (-π, π]
67+
if (phase_diff > pi)
68+
{
69+
accumulated_correction -= two_pi;
70+
phase_diff -= two_pi;
71+
}
72+
else if (phase_diff <= -pi)
6473
{
65-
var diff = phases[i] - phases[i - 1];
66-
67-
// Приводим разность к диапазону (-π, π]
68-
while (diff > pi)
69-
{
70-
diff -= two_pi;
71-
phases[i] -= two_pi;
72-
}
73-
while (diff <= -pi)
74-
{
75-
diff += two_pi;
76-
phases[i] += two_pi;
77-
}
74+
accumulated_correction += two_pi;
75+
phase_diff += two_pi;
7876
}
77+
78+
return phase_diff;
7979
}
8080
}
Lines changed: 252 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,252 @@
1+
using MathCore.DSP.Samples;
2+
using MathCore.DSP.Samples.Extensions;
3+
4+
namespace MathCore.DSP.Tests.Samples;
5+
6+
[TestClass]
7+
public class SampleSI16ExTests
8+
{
9+
/// <summary>Тест фазовой демодуляции для пустого массива</summary>
10+
[TestMethod]
11+
public void PhaseDemodulation_EmptyArray_ReturnsEmpty()
12+
{
13+
// Arrange
14+
var samples = Span<SampleSI16>.Empty;
15+
const double f0 = 1000.0;
16+
const double fd = 48000.0;
17+
18+
// Act
19+
var result = samples.PhaseDemodulation(f0, fd);
20+
21+
// Assert
22+
Assert.AreEqual(0, result.Length);
23+
}
24+
25+
/// <summary>Тест фазовой демодуляции для массива с одним элементом</summary>
26+
[TestMethod]
27+
public void PhaseDemodulation_SingleElement_ReturnsZero()
28+
{
29+
// Arrange
30+
var samples = new SampleSI16[] { new(100, 50) }.AsSpan();
31+
const double f0 = 1000.0;
32+
const double fd = 48000.0;
33+
34+
// Act
35+
var result = samples.PhaseDemodulation(f0, fd);
36+
37+
// Assert
38+
Assert.AreEqual(1, result.Length);
39+
Assert.AreEqual(0f, result[0]);
40+
}
41+
42+
/// <summary>Тест фазовой демодуляции для постоянного сигнала</summary>
43+
[TestMethod]
44+
public void PhaseDemodulation_ConstantSignal_ReturnsNearZeros()
45+
{
46+
// Arrange - создаем постоянный сигнал
47+
var samples = new SampleSI16[]
48+
{
49+
new(100, 0),
50+
new(100, 0),
51+
new(100, 0),
52+
new(100, 0),
53+
new(100, 0)
54+
}.AsSpan();
55+
56+
const double f0 = 1000.0;
57+
const double fd = 48000.0;
58+
59+
// Act
60+
var result = samples.PhaseDemodulation(f0, fd);
61+
62+
// Assert
63+
Assert.AreEqual(5, result.Length);
64+
Assert.AreEqual(0f, result[0]); // Первый всегда 0
65+
66+
// Для постоянного сигнала мгновенная частота должна быть близка к нулю
67+
// после вычитания центральной частоты результат должен быть близок к -f0
68+
for (var i = 1; i < result.Length; i++)
69+
{
70+
Assert.IsTrue(Math.Abs(result[i] + f0) < 50,
71+
$"Sample {i}: expected ~{-f0}, got {result[i]}");
72+
}
73+
}
74+
75+
/// <summary>Тест фазовой демодуляции для синусоидального сигнала с известной частотой</summary>
76+
[TestMethod]
77+
public void PhaseDemodulation_SinusoidalSignal_ReturnsExpectedFrequency()
78+
{
79+
// Arrange - создаем синусоидальный сигнал с частотой f_signal
80+
const double fd = 48000.0;
81+
const double f0 = 1000.0; // Центральная частота
82+
const double f_signal = 1500.0; // Частота сигнала
83+
const int samples_count = 200;
84+
85+
var samples = new SampleSI16[samples_count];
86+
var dt = 1.0 / fd;
87+
88+
for (var i = 0; i < samples_count; i++)
89+
{
90+
var t = i * dt;
91+
var angle = 2.0 * Math.PI * f_signal * t;
92+
var amplitude = 100.0;
93+
94+
samples[i] = new SampleSI16(
95+
(sbyte)(amplitude * Math.Cos(angle)),
96+
(sbyte)(amplitude * Math.Sin(angle))
97+
);
98+
}
99+
100+
// Act
101+
var result = samples.AsSpan().PhaseDemodulation(f0, fd);
102+
103+
// Assert
104+
Assert.AreEqual(samples_count, result.Length);
105+
Assert.AreEqual(0f, result[0]); // Первый всегда 0
106+
107+
// Проверяем, что результат близок к ожидаемой частоте (f_signal - f0)
108+
var expected_frequency = f_signal - f0;
109+
var tolerance = 100.0; // Допуск в Гц
110+
111+
// Проверяем стабильную часть сигнала (пропускаем начальные образцы)
112+
var stable_samples = 0;
113+
for (var i = 20; i < result.Length - 20; i++) // Пропускаем края для стабилизации
114+
{
115+
if (Math.Abs(result[i] - expected_frequency) < tolerance)
116+
stable_samples++;
117+
}
118+
119+
// Проверяем, что большинство образцов дает правильную частоту
120+
var expected_stable_count = (result.Length - 40) * 0.8; // 80% образцов должны быть стабильными
121+
Assert.IsTrue(stable_samples > expected_stable_count,
122+
$"Expected at least {expected_stable_count} stable samples, got {stable_samples}");
123+
}
124+
125+
/// <summary>Тест производительности оптимизированной фазовой демодуляции</summary>
126+
[TestMethod]
127+
public void PhaseDemodulation_OptimizedPerformance_BetterThanOldVersion()
128+
{
129+
// Arrange - большой массив для тестирования производительности
130+
const int samples_count = 1_000_000; // Увеличиваем размер для лучшего измерения
131+
var samples = new SampleSI16[samples_count];
132+
var random = new Random(42); // Фиксированное семя для воспроизводимости
133+
134+
for (var i = 0; i < samples_count; i++)
135+
{
136+
samples[i] = new SampleSI16(
137+
(sbyte)(random.Next(-128, 128)),
138+
(sbyte)(random.Next(-128, 128))
139+
);
140+
}
141+
142+
const double f0 = 1000.0;
143+
const double fd = 48000.0;
144+
145+
// Act - несколько прогонов для усреднения
146+
var times = new List<long>();
147+
for (var run = 0; run < 5; run++)
148+
{
149+
var stopwatch = System.Diagnostics.Stopwatch.StartNew();
150+
var result = samples.AsSpan().PhaseDemodulation(f0, fd);
151+
stopwatch.Stop();
152+
times.Add(stopwatch.ElapsedMilliseconds);
153+
154+
// Проверяем корректность результата
155+
Assert.AreEqual(samples_count, result.Length);
156+
}
157+
158+
var average_time = times.Average();
159+
160+
// Assert - должно быть быстрее чем 1000 мс для 1M образцов
161+
Assert.IsTrue(average_time < 1000,
162+
$"Optimized demodulation took {average_time:F1} ms on average, expected < 1000 ms");
163+
164+
Console.WriteLine($"Оптимизированная фазовая демодуляция {samples_count} образцов:");
165+
Console.WriteLine($"Среднее время: {average_time:F1} мс");
166+
Console.WriteLine($"Производительность: {samples_count / average_time / 1000:F1} млн. образцов/сек");
167+
}
168+
169+
/// <summary>Тест проверки правильности unwrap фазы</summary>
170+
[TestMethod]
171+
public void PhaseDemodulation_PhaseUnwrap_WorksCorrectly()
172+
{
173+
// Arrange - создаем сигнал с плавно изменяющейся фазой, но с перескоками ±2π
174+
const double fd = 1000.0;
175+
const double f0 = 0.0; // Нулевая центральная частота для простоты
176+
177+
var samples = new SampleSI16[]
178+
{
179+
new(100, 0), // фаза ≈ 0
180+
new(71, 71), // фаза ≈ π/4
181+
new(0, 100), // фаза ≈ π/2
182+
new(-71, 71), // фаза ≈ 3π/4
183+
new(-100, 0), // фаза ≈ π
184+
new(-71, -71), // фаза ≈ 5π/4 (но с unwrap должна быть ≈ -3π/4)
185+
new(0, -100), // фаза ≈ 3π/2 (но с unwrap должна быть ≈ -π/2)
186+
new(71, -71), // фаза ≈ 7π/4 (но с unwrap должна быть ≈ -π/4)
187+
}.AsSpan();
188+
189+
// Act
190+
var result = samples.PhaseDemodulation(f0, fd);
191+
192+
// Assert
193+
Assert.AreEqual(samples.Length, result.Length);
194+
Assert.AreEqual(0f, result[0]); // Первый всегда 0
195+
196+
// Проверяем, что результаты имеют разумные значения без больших скачков
197+
for (var i = 1; i < result.Length; i++)
198+
{
199+
Assert.IsTrue(Math.Abs(result[i]) < 1000,
200+
$"Sample {i}: frequency {result[i]} Hz seems too high");
201+
}
202+
}
203+
204+
/// <summary>Тест сравнения точности оптимизированной и наивной реализации</summary>
205+
[TestMethod]
206+
public void PhaseDemodulation_OptimizedVsNaive_SameAccuracy()
207+
{
208+
// Arrange - создаем чистый синусоидальный сигнал
209+
const double fd = 48000.0;
210+
const double f0 = 1000.0;
211+
const double f_signal = 1200.0;
212+
const int samples_count = 100;
213+
214+
var samples = new SampleSI16[samples_count];
215+
var dt = 1.0 / fd;
216+
217+
for (var i = 0; i < samples_count; i++)
218+
{
219+
var t = i * dt;
220+
var angle = 2.0 * Math.PI * f_signal * t;
221+
var amplitude = 120.0; // Используем почти максимальную амплитуду
222+
223+
samples[i] = new SampleSI16(
224+
(sbyte)(amplitude * Math.Cos(angle)),
225+
(sbyte)(amplitude * Math.Sin(angle))
226+
);
227+
}
228+
229+
// Act
230+
var optimized_result = samples.AsSpan().PhaseDemodulation(f0, fd);
231+
232+
// Assert - проверяем точность на стабильной части
233+
var expected_frequency = f_signal - f0;
234+
var errors = new List<float>();
235+
236+
for (var i = 10; i < optimized_result.Length - 10; i++)
237+
{
238+
var error = (float)Math.Abs(optimized_result[i] - expected_frequency);
239+
errors.Add(error);
240+
}
241+
242+
var average_error = errors.Average();
243+
var max_error = errors.Max();
244+
245+
Assert.IsTrue(average_error < 50, $"Average error {average_error:F1} Hz too high");
246+
Assert.IsTrue(max_error < 100, $"Max error {max_error:F1} Hz too high");
247+
248+
Console.WriteLine($"Точность оптимизированного алгоритма:");
249+
Console.WriteLine($"Средняя ошибка: {average_error:F1} Гц");
250+
Console.WriteLine($"Максимальная ошибка: {max_error:F1} Гц");
251+
}
252+
}

0 commit comments

Comments
 (0)