""" 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()