Split experiments from tests
The tests/ directory held 50 laboratory programs and no tests. They model channels, run hundreds of repetitions and write CSV, PNG and reports; calling that a test suite blocked introducing a real one, because any pytest run would have collected the labs and re-executed every experiment. - move all 50 lab programs to experiments/ with git mv, preserving history - rewrite the 38 cross-imports between labs from tests.labNNN to experiments.labNNN - leave tests/ empty for actual fast checks of protocol/ - point quick_gate and the hook at the new layout and add experiments/ to the syntax sweep - update the paths quoted in the Lab042 specification and the verifier agent definition This also defuses the import-time work finding without touching 41 files: the labs still create directories and write files on import, but nothing imports them now except the gate, which does so deliberately. Gate passes: syntax clean, protocol imports, 15 lab modules import, 2 functional suites run. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
LittleSam129
875
experiments/lab025_pluto_wfm_demodulation.py
Normal file
875
experiments/lab025_pluto_wfm_demodulation.py
Normal file
@@ -0,0 +1,875 @@
|
||||
"""
|
||||
Lab025. Приём и программная WFM-демодуляция FM-радиостанции.
|
||||
|
||||
Схема подключения:
|
||||
антенна 40–860 МГц -> RX1
|
||||
|
||||
Используется только приёмный канал RX1. Передатчики TX1 и TX2
|
||||
не используются. Необработанные IQ-сэмплы на диск не сохраняются.
|
||||
"""
|
||||
|
||||
from pathlib import Path
|
||||
|
||||
import adi
|
||||
import matplotlib
|
||||
import numpy as np
|
||||
from scipy import signal
|
||||
from scipy.io import wavfile
|
||||
|
||||
|
||||
# Для лабораторной графики сохраняется в файлы без блокирующих окон.
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------
|
||||
# Параметры приёмника Pluto+
|
||||
# ---------------------------------------------------------------------
|
||||
|
||||
PLUTO_URI = "ip:192.168.2.1"
|
||||
|
||||
STATION_FREQUENCY_HZ = 100_100_000
|
||||
LO_OFFSET_HZ = 250_000
|
||||
|
||||
RX_LO_FREQUENCY_HZ = (
|
||||
STATION_FREQUENCY_HZ + LO_OFFSET_HZ
|
||||
)
|
||||
|
||||
SAMPLE_RATE_HZ = 2_400_000
|
||||
RX_BANDWIDTH_HZ = 1_500_000
|
||||
|
||||
RX_BUFFER_SIZE = 65_536
|
||||
DISCARD_BUFFER_COUNT = 3
|
||||
CAPTURE_BUFFER_COUNT = 128
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------
|
||||
# Параметры обработки сигнала
|
||||
# ---------------------------------------------------------------------
|
||||
|
||||
CHANNEL_DECIMATION = 10
|
||||
CHANNEL_SAMPLE_RATE_HZ = 240_000
|
||||
|
||||
AUDIO_CUTOFF_HZ = 15_000
|
||||
AUDIO_SAMPLE_RATE_HZ = 48_000
|
||||
|
||||
DEEMPHASIS_TIME_CONSTANT_SECONDS = 50e-6
|
||||
|
||||
AUDIO_TRANSIENT_DURATION_SECONDS = 0.05
|
||||
AUDIO_REFERENCE_PERCENTILE = 99.5
|
||||
AUDIO_TARGET_LEVEL = 0.85
|
||||
MINIMUM_AUDIO_REFERENCE_PEAK = 1e-12
|
||||
|
||||
RF_WELCH_MAXIMUM_SAMPLE_COUNT = 1_048_576
|
||||
RF_WELCH_SEGMENT_LENGTH = 8_192
|
||||
RF_WELCH_OVERLAP_LENGTH = 4_096
|
||||
|
||||
AUDIO_WELCH_SEGMENT_LENGTH = 8_192
|
||||
|
||||
|
||||
# ---------------------------------------------------------------------
|
||||
# Выходные файлы
|
||||
# ---------------------------------------------------------------------
|
||||
|
||||
OUTPUT_DIRECTORY = Path("data/processed/lab025")
|
||||
|
||||
WAV_FILE_PATH = (
|
||||
OUTPUT_DIRECTORY / "lab025_wfm_audio_100_1mhz.wav"
|
||||
)
|
||||
RF_SPECTRUM_FILE_PATH = (
|
||||
OUTPUT_DIRECTORY / "lab025_rf_spectrum.png"
|
||||
)
|
||||
AUDIO_WAVEFORM_FILE_PATH = (
|
||||
OUTPUT_DIRECTORY / "lab025_audio_waveform.png"
|
||||
)
|
||||
AUDIO_SPECTRUM_FILE_PATH = (
|
||||
OUTPUT_DIRECTORY / "lab025_audio_spectrum.png"
|
||||
)
|
||||
REPORT_FILE_PATH = OUTPUT_DIRECTORY / "lab025_report.txt"
|
||||
|
||||
|
||||
def configure_receiver() -> adi.Pluto:
|
||||
"""
|
||||
Подключается к Pluto+ и настраивает только приёмный канал RX1.
|
||||
|
||||
Частота гетеродина смещена на 250 кГц выше частоты станции.
|
||||
Благодаря этому полезный сигнал не совпадает с аппаратным DC-пиком
|
||||
в центре комплексной полосы приёмника.
|
||||
"""
|
||||
|
||||
print("Подключение к Pluto+...")
|
||||
|
||||
sdr = adi.Pluto(uri=PLUTO_URI)
|
||||
|
||||
# Используем только первый приёмный канал RX1.
|
||||
sdr.rx_enabled_channels = [0]
|
||||
|
||||
sdr.sample_rate = SAMPLE_RATE_HZ
|
||||
sdr.rx_lo = RX_LO_FREQUENCY_HZ
|
||||
sdr.rx_rf_bandwidth = RX_BANDWIDTH_HZ
|
||||
|
||||
# Медленная АРУ подходит для приёма вещательной FM-станции.
|
||||
sdr.gain_control_mode_chan0 = "slow_attack"
|
||||
sdr.rx_buffer_size = RX_BUFFER_SIZE
|
||||
|
||||
return sdr
|
||||
|
||||
|
||||
def receive_samples(sdr: adi.Pluto) -> np.ndarray:
|
||||
"""
|
||||
Получает последовательность комплексных IQ-буферов.
|
||||
|
||||
Первые буферы отбрасываются, поскольку после настройки приёмника
|
||||
в них могут присутствовать переходные процессы АРУ и фильтров.
|
||||
Рабочие буферы объединяются только в оперативной памяти.
|
||||
"""
|
||||
|
||||
print()
|
||||
print(
|
||||
"Отбрасывание переходных буферов: "
|
||||
f"{DISCARD_BUFFER_COUNT}..."
|
||||
)
|
||||
|
||||
for _ in range(DISCARD_BUFFER_COUNT):
|
||||
_ = sdr.rx()
|
||||
|
||||
print("Получение рабочих IQ-буферов...")
|
||||
|
||||
received_buffers: list[np.ndarray] = []
|
||||
|
||||
for buffer_number in range(1, CAPTURE_BUFFER_COUNT + 1):
|
||||
received_buffer = np.asarray(
|
||||
sdr.rx(),
|
||||
dtype=np.complex64,
|
||||
)
|
||||
|
||||
if received_buffer.size == 0:
|
||||
raise RuntimeError(
|
||||
"Pluto+ вернул пустой рабочий IQ-буфер."
|
||||
)
|
||||
|
||||
received_buffers.append(received_buffer)
|
||||
|
||||
if (
|
||||
buffer_number % 16 == 0
|
||||
or buffer_number == CAPTURE_BUFFER_COUNT
|
||||
):
|
||||
print(
|
||||
f" Принято {buffer_number:3d}/"
|
||||
f"{CAPTURE_BUFFER_COUNT} буферов"
|
||||
)
|
||||
|
||||
if not received_buffers:
|
||||
raise RuntimeError("Pluto+ не вернул IQ-сэмплы.")
|
||||
|
||||
combined_samples = np.concatenate(received_buffers)
|
||||
|
||||
if combined_samples.size == 0:
|
||||
raise RuntimeError("После объединения получен пустой IQ-массив.")
|
||||
|
||||
if not np.all(np.isfinite(combined_samples)):
|
||||
raise RuntimeError("В IQ-сэмплах обнаружены NaN или Inf.")
|
||||
|
||||
return np.asarray(combined_samples, dtype=np.complex64)
|
||||
|
||||
|
||||
def shift_station_to_baseband(
|
||||
samples: np.ndarray,
|
||||
sample_rate_hz: float,
|
||||
frequency_shift_hz: float,
|
||||
) -> np.ndarray:
|
||||
"""
|
||||
Переносит выбранную станцию в центр цифровой полосы.
|
||||
|
||||
При положительном смещении комплексный генератор переносит сигнал,
|
||||
расположенный на отрицательной относительной частоте, к 0 Гц.
|
||||
"""
|
||||
|
||||
if samples.size == 0:
|
||||
raise RuntimeError(
|
||||
"Невозможно выполнить цифровой перенос пустого IQ-массива."
|
||||
)
|
||||
|
||||
sample_indices = np.arange(
|
||||
len(samples),
|
||||
dtype=np.float64,
|
||||
)
|
||||
|
||||
digital_oscillator = np.exp(
|
||||
1j
|
||||
* 2.0
|
||||
* np.pi
|
||||
* frequency_shift_hz
|
||||
* sample_indices
|
||||
/ sample_rate_hz
|
||||
).astype(np.complex64)
|
||||
|
||||
centered_samples = samples * digital_oscillator
|
||||
|
||||
if not np.all(np.isfinite(centered_samples)):
|
||||
raise RuntimeError(
|
||||
"После цифрового переноса обнаружены NaN или Inf."
|
||||
)
|
||||
|
||||
return np.asarray(centered_samples, dtype=np.complex64)
|
||||
|
||||
|
||||
def extract_wfm_channel(
|
||||
centered_samples: np.ndarray,
|
||||
) -> np.ndarray:
|
||||
"""
|
||||
Фильтрует WFM-канал и понижает частоту до 240 кГц.
|
||||
|
||||
Полигармоническая передискретизация одновременно выполняет
|
||||
низкочастотную фильтрацию и децимацию в десять раз.
|
||||
"""
|
||||
|
||||
calculated_sample_rate_hz = (
|
||||
SAMPLE_RATE_HZ / CHANNEL_DECIMATION
|
||||
)
|
||||
|
||||
if not np.isclose(
|
||||
calculated_sample_rate_hz,
|
||||
CHANNEL_SAMPLE_RATE_HZ,
|
||||
):
|
||||
raise RuntimeError(
|
||||
"Частота канального сигнала после децимации "
|
||||
"не равна 240 кГц."
|
||||
)
|
||||
|
||||
channel_samples = signal.resample_poly(
|
||||
centered_samples,
|
||||
up=1,
|
||||
down=CHANNEL_DECIMATION,
|
||||
window=("kaiser", 8.0),
|
||||
)
|
||||
|
||||
if channel_samples.size == 0:
|
||||
raise RuntimeError("После выделения WFM-канала массив пуст.")
|
||||
|
||||
if not np.all(np.isfinite(channel_samples)):
|
||||
raise RuntimeError(
|
||||
"После выделения WFM-канала обнаружены NaN или Inf."
|
||||
)
|
||||
|
||||
return np.asarray(channel_samples, dtype=np.complex64)
|
||||
|
||||
|
||||
def demodulate_fm(
|
||||
channel_samples: np.ndarray,
|
||||
) -> np.ndarray:
|
||||
"""
|
||||
Выполняет частотную демодуляцию по фазовой разности отсчётов.
|
||||
|
||||
Результат представляет монофонический композитный FM-сигнал
|
||||
до звуковой фильтрации и коррекции предыскажений.
|
||||
"""
|
||||
|
||||
if channel_samples.size < 2:
|
||||
raise RuntimeError(
|
||||
"Недостаточно канальных отсчётов для FM-демодуляции."
|
||||
)
|
||||
|
||||
phase_difference = np.angle(
|
||||
channel_samples[1:]
|
||||
* np.conj(channel_samples[:-1])
|
||||
).astype(np.float64)
|
||||
|
||||
demodulated = phase_difference - np.mean(phase_difference)
|
||||
|
||||
if demodulated.size == 0:
|
||||
raise RuntimeError("FM-дискриминатор вернул пустой массив.")
|
||||
|
||||
if not np.all(np.isfinite(demodulated)):
|
||||
raise RuntimeError(
|
||||
"После FM-демодуляции обнаружены NaN или Inf."
|
||||
)
|
||||
|
||||
return demodulated
|
||||
|
||||
|
||||
def lowpass_audio(
|
||||
demodulated_samples: np.ndarray,
|
||||
sample_rate_hz: float,
|
||||
) -> np.ndarray:
|
||||
"""
|
||||
Выделяет монофонический звук L+R в полосе до 15 кГц.
|
||||
|
||||
Фильтр подавляет стереопилот 19 кГц, стереоразностную часть,
|
||||
RDS и внеполосный высокочастотный шум.
|
||||
"""
|
||||
|
||||
audio_sos = signal.butter(
|
||||
6,
|
||||
AUDIO_CUTOFF_HZ,
|
||||
btype="lowpass",
|
||||
fs=sample_rate_hz,
|
||||
output="sos",
|
||||
)
|
||||
|
||||
filtered_audio = signal.sosfiltfilt(
|
||||
audio_sos,
|
||||
demodulated_samples,
|
||||
)
|
||||
|
||||
if filtered_audio.size == 0:
|
||||
raise RuntimeError("Звуковой фильтр вернул пустой массив.")
|
||||
|
||||
if not np.all(np.isfinite(filtered_audio)):
|
||||
raise RuntimeError(
|
||||
"После звукового фильтра обнаружены NaN или Inf."
|
||||
)
|
||||
|
||||
return np.asarray(filtered_audio, dtype=np.float64)
|
||||
|
||||
|
||||
def apply_deemphasis(
|
||||
audio_samples: np.ndarray,
|
||||
sample_rate_hz: float,
|
||||
time_constant_seconds: float,
|
||||
) -> np.ndarray:
|
||||
"""
|
||||
Выполняет европейскую FM-коррекцию предыскажений 50 мкс.
|
||||
|
||||
Используется устойчивый однополюсный рекурсивный фильтр,
|
||||
реализованный функцией scipy.signal.lfilter.
|
||||
"""
|
||||
|
||||
alpha = np.exp(
|
||||
-1.0
|
||||
/ (
|
||||
sample_rate_hz
|
||||
* time_constant_seconds
|
||||
)
|
||||
)
|
||||
|
||||
deemphasized = signal.lfilter(
|
||||
[1.0 - alpha],
|
||||
[1.0, -alpha],
|
||||
audio_samples,
|
||||
)
|
||||
|
||||
deemphasized = deemphasized - np.mean(deemphasized)
|
||||
|
||||
if deemphasized.size == 0:
|
||||
raise RuntimeError("De-emphasis вернул пустой массив.")
|
||||
|
||||
if not np.all(np.isfinite(deemphasized)):
|
||||
raise RuntimeError(
|
||||
"После de-emphasis обнаружены NaN или Inf."
|
||||
)
|
||||
|
||||
return np.asarray(deemphasized, dtype=np.float64)
|
||||
|
||||
|
||||
def resample_audio_to_output_rate(
|
||||
audio_samples: np.ndarray,
|
||||
) -> np.ndarray:
|
||||
"""
|
||||
Преобразует звуковой сигнал с 240 кГц в выходные 48 кГц.
|
||||
"""
|
||||
|
||||
output_decimation = 5
|
||||
calculated_output_rate_hz = (
|
||||
CHANNEL_SAMPLE_RATE_HZ / output_decimation
|
||||
)
|
||||
|
||||
if not np.isclose(
|
||||
calculated_output_rate_hz,
|
||||
AUDIO_SAMPLE_RATE_HZ,
|
||||
):
|
||||
raise RuntimeError(
|
||||
"Рассчитанная частота WAV не равна 48 кГц."
|
||||
)
|
||||
|
||||
output_audio = signal.resample_poly(
|
||||
audio_samples,
|
||||
up=1,
|
||||
down=output_decimation,
|
||||
)
|
||||
|
||||
if output_audio.size == 0:
|
||||
raise RuntimeError(
|
||||
"После преобразования в 48 кГц получен пустой массив."
|
||||
)
|
||||
|
||||
if not np.all(np.isfinite(output_audio)):
|
||||
raise RuntimeError(
|
||||
"После преобразования звука обнаружены NaN или Inf."
|
||||
)
|
||||
|
||||
return np.asarray(output_audio, dtype=np.float64)
|
||||
|
||||
|
||||
def convert_audio_to_pcm16(
|
||||
audio_samples: np.ndarray,
|
||||
) -> np.ndarray:
|
||||
"""
|
||||
Удаляет переходный участок, нормализует звук и создаёт PCM16.
|
||||
|
||||
Опорный уровень определяется по 99,5-му процентилю модуля,
|
||||
поэтому единичный выброс не делает весь WAV слишком тихим.
|
||||
"""
|
||||
|
||||
if audio_samples.size == 0:
|
||||
raise RuntimeError("Невозможно нормализовать пустой звук.")
|
||||
|
||||
centered_audio = audio_samples - np.mean(audio_samples)
|
||||
|
||||
transient_sample_count = int(
|
||||
round(
|
||||
AUDIO_TRANSIENT_DURATION_SECONDS
|
||||
* AUDIO_SAMPLE_RATE_HZ
|
||||
)
|
||||
)
|
||||
|
||||
if centered_audio.size <= transient_sample_count:
|
||||
raise RuntimeError(
|
||||
"Звуковой массив короче переходного участка 0,05 с."
|
||||
)
|
||||
|
||||
centered_audio = centered_audio[transient_sample_count:]
|
||||
|
||||
reference_peak = float(
|
||||
np.percentile(
|
||||
np.abs(centered_audio),
|
||||
AUDIO_REFERENCE_PERCENTILE,
|
||||
)
|
||||
)
|
||||
|
||||
if (
|
||||
not np.isfinite(reference_peak)
|
||||
or reference_peak <= MINIMUM_AUDIO_REFERENCE_PEAK
|
||||
):
|
||||
raise RuntimeError(
|
||||
"Сигнал слишком мал для безопасной нормализации PCM16."
|
||||
)
|
||||
|
||||
normalized_audio = (
|
||||
centered_audio
|
||||
/ reference_peak
|
||||
* AUDIO_TARGET_LEVEL
|
||||
)
|
||||
|
||||
normalized_audio = np.clip(
|
||||
normalized_audio,
|
||||
-1.0,
|
||||
1.0,
|
||||
)
|
||||
|
||||
pcm16_audio = np.round(
|
||||
normalized_audio * np.iinfo(np.int16).max
|
||||
).astype(np.int16)
|
||||
|
||||
if pcm16_audio.size == 0:
|
||||
raise RuntimeError("После преобразования PCM16 массив пуст.")
|
||||
|
||||
return pcm16_audio
|
||||
|
||||
|
||||
def calculate_rf_spectrum(
|
||||
samples: np.ndarray,
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""
|
||||
Рассчитывает ограниченный спектр исходного IQ методом Уэлча.
|
||||
|
||||
Для графика используется не более 1 048 576 отсчётов, поэтому
|
||||
полная восьмимиллионная запись не передаётся в одну большую FFT.
|
||||
"""
|
||||
|
||||
analysis_sample_count = min(
|
||||
len(samples),
|
||||
RF_WELCH_MAXIMUM_SAMPLE_COUNT,
|
||||
)
|
||||
|
||||
analysis_samples = samples[:analysis_sample_count]
|
||||
|
||||
relative_frequencies_hz, power_density = signal.welch(
|
||||
analysis_samples,
|
||||
fs=SAMPLE_RATE_HZ,
|
||||
window="hann",
|
||||
nperseg=RF_WELCH_SEGMENT_LENGTH,
|
||||
noverlap=RF_WELCH_OVERLAP_LENGTH,
|
||||
detrend=False,
|
||||
return_onesided=False,
|
||||
scaling="density",
|
||||
)
|
||||
|
||||
relative_frequencies_hz = np.fft.fftshift(
|
||||
relative_frequencies_hz
|
||||
)
|
||||
power_density = np.fft.fftshift(power_density)
|
||||
|
||||
minimum_positive_value = np.finfo(np.float64).tiny
|
||||
power_db = 10.0 * np.log10(
|
||||
power_density + minimum_positive_value
|
||||
)
|
||||
power_db -= np.max(power_db)
|
||||
|
||||
return relative_frequencies_hz, power_db
|
||||
|
||||
|
||||
def create_rf_spectrum_graph(samples: np.ndarray) -> None:
|
||||
"""
|
||||
Сохраняет спектр до цифрового переноса станции в центр.
|
||||
"""
|
||||
|
||||
frequencies_hz, power_db = calculate_rf_spectrum(samples)
|
||||
|
||||
figure, axes = plt.subplots(figsize=(13, 7))
|
||||
|
||||
axes.plot(
|
||||
frequencies_hz / 1e3,
|
||||
power_db,
|
||||
linewidth=0.8,
|
||||
)
|
||||
|
||||
expected_station_offset_hz = -LO_OFFSET_HZ
|
||||
|
||||
axes.axvline(
|
||||
expected_station_offset_hz / 1e3,
|
||||
color="tab:red",
|
||||
linestyle="--",
|
||||
linewidth=1.4,
|
||||
label=(
|
||||
"Ожидаемая станция 100,100 МГц "
|
||||
"(-250 кГц)"
|
||||
),
|
||||
)
|
||||
|
||||
axes.set_title(
|
||||
"Lab025. Спектр принятого IQ до цифрового переноса\n"
|
||||
"RX LO = 100,350 МГц; станция = 100,100 МГц"
|
||||
)
|
||||
axes.set_xlabel("Относительная частота относительно RX LO, кГц")
|
||||
axes.set_ylabel("Относительная мощность, дБ")
|
||||
axes.set_ylim(-90, 5)
|
||||
axes.grid(True)
|
||||
axes.legend()
|
||||
|
||||
figure.tight_layout()
|
||||
figure.savefig(RF_SPECTRUM_FILE_PATH, dpi=160)
|
||||
plt.close(figure)
|
||||
|
||||
|
||||
def create_audio_waveform_graph(pcm16_audio: np.ndarray) -> None:
|
||||
"""
|
||||
Сохраняет первые 0,1 секунды нормированного звука.
|
||||
"""
|
||||
|
||||
displayed_sample_count = min(
|
||||
len(pcm16_audio),
|
||||
int(round(0.1 * AUDIO_SAMPLE_RATE_HZ)),
|
||||
)
|
||||
|
||||
displayed_audio = (
|
||||
pcm16_audio[:displayed_sample_count].astype(np.float64)
|
||||
/ np.iinfo(np.int16).max
|
||||
)
|
||||
time_seconds = (
|
||||
np.arange(displayed_sample_count, dtype=np.float64)
|
||||
/ AUDIO_SAMPLE_RATE_HZ
|
||||
)
|
||||
|
||||
figure, axes = plt.subplots(figsize=(12, 5))
|
||||
axes.plot(time_seconds, displayed_audio, linewidth=0.8)
|
||||
axes.set_title("Lab025. Первые 0,1 секунды демодулированного звука")
|
||||
axes.set_xlabel("Время, с")
|
||||
axes.set_ylabel("Нормированная амплитуда")
|
||||
axes.set_ylim(-1.05, 1.05)
|
||||
axes.grid(True)
|
||||
|
||||
figure.tight_layout()
|
||||
figure.savefig(AUDIO_WAVEFORM_FILE_PATH, dpi=160)
|
||||
plt.close(figure)
|
||||
|
||||
|
||||
def create_audio_spectrum_graph(pcm16_audio: np.ndarray) -> None:
|
||||
"""
|
||||
Рассчитывает методом Уэлча и сохраняет односторонний спектр WAV.
|
||||
"""
|
||||
|
||||
normalized_audio = (
|
||||
pcm16_audio.astype(np.float64)
|
||||
/ np.iinfo(np.int16).max
|
||||
)
|
||||
|
||||
segment_length = min(
|
||||
AUDIO_WELCH_SEGMENT_LENGTH,
|
||||
len(normalized_audio),
|
||||
)
|
||||
|
||||
frequencies_hz, power_density = signal.welch(
|
||||
normalized_audio,
|
||||
fs=AUDIO_SAMPLE_RATE_HZ,
|
||||
window="hann",
|
||||
nperseg=segment_length,
|
||||
noverlap=segment_length // 2,
|
||||
detrend="constant",
|
||||
return_onesided=True,
|
||||
scaling="density",
|
||||
)
|
||||
|
||||
minimum_positive_value = np.finfo(np.float64).tiny
|
||||
power_db = 10.0 * np.log10(
|
||||
power_density + minimum_positive_value
|
||||
)
|
||||
power_db -= np.max(power_db)
|
||||
|
||||
figure, axes = plt.subplots(figsize=(12, 6))
|
||||
axes.plot(frequencies_hz / 1e3, power_db, linewidth=0.9)
|
||||
axes.set_xlim(0, 20)
|
||||
axes.set_ylim(-100, 5)
|
||||
axes.set_title("Lab025. Спектр демодулированного монофонического звука")
|
||||
axes.set_xlabel("Частота, кГц")
|
||||
axes.set_ylabel("Относительная спектральная плотность, дБ")
|
||||
axes.grid(True)
|
||||
|
||||
figure.tight_layout()
|
||||
figure.savefig(AUDIO_SPECTRUM_FILE_PATH, dpi=160)
|
||||
plt.close(figure)
|
||||
|
||||
|
||||
def calculate_audio_statistics(
|
||||
pcm16_audio: np.ndarray,
|
||||
) -> tuple[float, float, float]:
|
||||
"""
|
||||
Возвращает нормированные RMS, пик и процент предельных отсчётов.
|
||||
"""
|
||||
|
||||
normalized_audio = (
|
||||
pcm16_audio.astype(np.float64)
|
||||
/ np.iinfo(np.int16).max
|
||||
)
|
||||
|
||||
rms_value = float(
|
||||
np.sqrt(np.mean(normalized_audio ** 2))
|
||||
)
|
||||
peak_value = float(np.max(np.abs(normalized_audio)))
|
||||
|
||||
clipped_sample_count = int(
|
||||
np.count_nonzero(
|
||||
(pcm16_audio == np.iinfo(np.int16).min)
|
||||
| (pcm16_audio == np.iinfo(np.int16).max)
|
||||
)
|
||||
)
|
||||
clipping_percentage = (
|
||||
100.0
|
||||
* clipped_sample_count
|
||||
/ len(pcm16_audio)
|
||||
)
|
||||
|
||||
return rms_value, peak_value, clipping_percentage
|
||||
|
||||
|
||||
def save_report(
|
||||
iq_sample_count: int,
|
||||
iq_duration_seconds: float,
|
||||
pcm16_audio: np.ndarray,
|
||||
audio_rms: float,
|
||||
audio_peak: float,
|
||||
clipping_percentage: float,
|
||||
) -> None:
|
||||
"""
|
||||
Создаёт текстовый отчёт с параметрами приёма и WAV-файла.
|
||||
"""
|
||||
|
||||
wav_duration_seconds = (
|
||||
len(pcm16_audio) / AUDIO_SAMPLE_RATE_HZ
|
||||
)
|
||||
|
||||
report_lines = [
|
||||
"Lab025. Приём и программная WFM-демодуляция",
|
||||
"",
|
||||
f"URI Pluto+: {PLUTO_URI}",
|
||||
f"Частота станции: {STATION_FREQUENCY_HZ} Гц",
|
||||
f"RX LO: {RX_LO_FREQUENCY_HZ} Гц",
|
||||
f"Цифровое смещение: +{LO_OFFSET_HZ} Гц",
|
||||
f"Частота дискретизации RX: {SAMPLE_RATE_HZ} Гц",
|
||||
f"Полоса RX: {RX_BANDWIDTH_HZ} Гц",
|
||||
f"Количество рабочих буферов: {CAPTURE_BUFFER_COUNT}",
|
||||
f"Количество IQ-сэмплов: {iq_sample_count}",
|
||||
f"Длительность IQ-записи: {iq_duration_seconds:.6f} с",
|
||||
(
|
||||
"Частота канального сигнала после децимации: "
|
||||
f"{CHANNEL_SAMPLE_RATE_HZ} Гц"
|
||||
),
|
||||
f"Длительность итогового WAV: {wav_duration_seconds:.6f} с",
|
||||
f"Частота WAV: {AUDIO_SAMPLE_RATE_HZ} Гц",
|
||||
"Тип PCM: signed PCM16, mono",
|
||||
f"RMS итогового аудио: {audio_rms:.6f}",
|
||||
f"Пиковая амплитуда: {audio_peak:.6f}",
|
||||
f"Отсчёты на границе PCM16: {clipping_percentage:.6f} %",
|
||||
"",
|
||||
"Выходные файлы:",
|
||||
f"WAV: {WAV_FILE_PATH}",
|
||||
f"Радиоспектр: {RF_SPECTRUM_FILE_PATH}",
|
||||
f"Звуковая волна: {AUDIO_WAVEFORM_FILE_PATH}",
|
||||
f"Спектр звука: {AUDIO_SPECTRUM_FILE_PATH}",
|
||||
f"Отчёт: {REPORT_FILE_PATH}",
|
||||
"",
|
||||
"Использован только приёмный канал RX1.",
|
||||
"Передатчики TX1 и TX2 не использовались.",
|
||||
"Необработанные IQ-сэмплы на диск не сохранялись.",
|
||||
]
|
||||
|
||||
REPORT_FILE_PATH.write_text(
|
||||
"\n".join(report_lines),
|
||||
encoding="utf-8",
|
||||
)
|
||||
|
||||
|
||||
def main() -> None:
|
||||
"""
|
||||
Выполняет полный цикл приёма, WFM-демодуляции и сохранения WAV.
|
||||
"""
|
||||
|
||||
OUTPUT_DIRECTORY.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
sdr = None
|
||||
|
||||
try:
|
||||
sdr = configure_receiver()
|
||||
|
||||
print()
|
||||
print("Параметры RX:")
|
||||
print(f" URI: {PLUTO_URI}")
|
||||
print(
|
||||
f" Станция: "
|
||||
f"{STATION_FREQUENCY_HZ / 1e6:.3f} МГц"
|
||||
)
|
||||
print(
|
||||
f" RX LO: "
|
||||
f"{sdr.rx_lo / 1e6:.3f} МГц"
|
||||
)
|
||||
print(
|
||||
f" Частота дискретизации: "
|
||||
f"{sdr.sample_rate / 1e6:.3f} Мвыб/с"
|
||||
)
|
||||
print(
|
||||
f" Полоса RX: "
|
||||
f"{sdr.rx_rf_bandwidth / 1e6:.3f} МГц"
|
||||
)
|
||||
print(
|
||||
f" Режим усиления: "
|
||||
f"{sdr.gain_control_mode_chan0}"
|
||||
)
|
||||
print(f" Размер буфера: {sdr.rx_buffer_size}")
|
||||
print(" Активный канал: RX1")
|
||||
|
||||
samples = receive_samples(sdr)
|
||||
|
||||
iq_sample_count = len(samples)
|
||||
iq_duration_seconds = iq_sample_count / SAMPLE_RATE_HZ
|
||||
|
||||
print()
|
||||
print(f"Количество IQ-сэмплов: {iq_sample_count}")
|
||||
print(f"Длительность записи: {iq_duration_seconds:.3f} с")
|
||||
|
||||
# Удаляем остаточную комплексную постоянную составляющую.
|
||||
samples = samples - np.mean(samples)
|
||||
|
||||
print()
|
||||
print("Цифровой перенос станции в центр полосы...")
|
||||
centered_samples = shift_station_to_baseband(
|
||||
samples=samples,
|
||||
sample_rate_hz=SAMPLE_RATE_HZ,
|
||||
frequency_shift_hz=LO_OFFSET_HZ,
|
||||
)
|
||||
|
||||
print("Фильтрация WFM-канала и децимация до 240 кГц...")
|
||||
channel_samples = extract_wfm_channel(centered_samples)
|
||||
|
||||
print("FM-демодуляция...")
|
||||
demodulated_samples = demodulate_fm(channel_samples)
|
||||
|
||||
print("Звуковой low-pass фильтр 0–15 кГц...")
|
||||
filtered_audio = lowpass_audio(
|
||||
demodulated_samples=demodulated_samples,
|
||||
sample_rate_hz=CHANNEL_SAMPLE_RATE_HZ,
|
||||
)
|
||||
|
||||
print("De-emphasis 50 мкс...")
|
||||
deemphasized_audio = apply_deemphasis(
|
||||
audio_samples=filtered_audio,
|
||||
sample_rate_hz=CHANNEL_SAMPLE_RATE_HZ,
|
||||
time_constant_seconds=(
|
||||
DEEMPHASIS_TIME_CONSTANT_SECONDS
|
||||
),
|
||||
)
|
||||
|
||||
print("Преобразование звука в 48 кГц...")
|
||||
output_audio = resample_audio_to_output_rate(
|
||||
deemphasized_audio
|
||||
)
|
||||
pcm16_audio = convert_audio_to_pcm16(output_audio)
|
||||
|
||||
audio_rms, audio_peak, clipping_percentage = (
|
||||
calculate_audio_statistics(pcm16_audio)
|
||||
)
|
||||
|
||||
print("Сохранение WAV...")
|
||||
wavfile.write(
|
||||
WAV_FILE_PATH,
|
||||
AUDIO_SAMPLE_RATE_HZ,
|
||||
pcm16_audio,
|
||||
)
|
||||
|
||||
if not WAV_FILE_PATH.is_file():
|
||||
raise RuntimeError("Выходной WAV-файл не создан.")
|
||||
|
||||
print("Сохранение графиков...")
|
||||
create_rf_spectrum_graph(samples)
|
||||
create_audio_waveform_graph(pcm16_audio)
|
||||
create_audio_spectrum_graph(pcm16_audio)
|
||||
|
||||
print("Создание текстового отчёта...")
|
||||
save_report(
|
||||
iq_sample_count=iq_sample_count,
|
||||
iq_duration_seconds=iq_duration_seconds,
|
||||
pcm16_audio=pcm16_audio,
|
||||
audio_rms=audio_rms,
|
||||
audio_peak=audio_peak,
|
||||
clipping_percentage=clipping_percentage,
|
||||
)
|
||||
|
||||
print()
|
||||
print("Созданы файлы:")
|
||||
print(f" WAV: {WAV_FILE_PATH}")
|
||||
print(f" Радиоспектр: {RF_SPECTRUM_FILE_PATH}")
|
||||
print(f" Звуковая волна: {AUDIO_WAVEFORM_FILE_PATH}")
|
||||
print(f" Спектр звука: {AUDIO_SPECTRUM_FILE_PATH}")
|
||||
print(f" Отчёт: {REPORT_FILE_PATH}")
|
||||
print()
|
||||
print(f"RMS аудио: {audio_rms:.6f}")
|
||||
print(f"Пиковая амплитуда: {audio_peak:.6f}")
|
||||
print(
|
||||
f"Отсчёты на границе PCM: "
|
||||
f"{clipping_percentage:.6f} %"
|
||||
)
|
||||
print()
|
||||
print("Lab025 выполнена успешно.")
|
||||
print()
|
||||
print("Для прослушивания откройте:")
|
||||
print(WAV_FILE_PATH)
|
||||
|
||||
finally:
|
||||
if sdr is not None:
|
||||
destroy_buffer = getattr(
|
||||
sdr,
|
||||
"rx_destroy_buffer",
|
||||
None,
|
||||
)
|
||||
|
||||
if callable(destroy_buffer):
|
||||
try:
|
||||
destroy_buffer()
|
||||
except Exception as cleanup_error:
|
||||
print(
|
||||
"Предупреждение: не удалось освободить "
|
||||
f"RX-буфер: {cleanup_error}"
|
||||
)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user