Генератор сигналов GPS

В примере рассматривается работа блока gpsWaveformGenerator, а именно генерируемые с его помощью сигналы различных типов (P-, L1C-, L2C- и L5-коды).

Введение

GPS — это космическая радионавигационная инфраструктура, состоящая из группировки спутников, передающих навигационные сигналы, и сети наземных станций и станций управления спутниками для мониторинга и контроля этих спутников.

Спутники GPS используют три разные частоты для гражданских применений: L1 на 1575,42 МГц, L2 на 1227,60 МГц и L5 на 1176,45 МГц.

GPS использует пять различных гражданских сигналов: код грубого захвата (C/A) и точный код (P) в унаследованном диапазоне L1, модернизированный гражданский код L1 (L1C) в диапазоне L1, гражданский средний код (CM) и гражданский длинный код (CL) в диапазоне L2C, а также синфазный код (I5-код) и квадратурный код (Q5-код) в диапазоне L5.

P-код

P-код — дальномерный код, использует метод модуляции BPSK (Binary Phase Shift Keying) для передачи на частотах L1 и L2 с чиповой скоростью 10,23 Мчип/с. Более высокая чиповая скорость P-кода способствует повышению точности по сравнению с C/A-кодом.

P-код является общедоступным, однако его шифрование с помощью Y-кода для формирования P(Y)-кода ограничивает его использование военными целями.

Ниже представлен скрипт для генерации и отображения спектрограммы Р-кода.

SampleRate = 4 * 10.23e6 задает частоту дискретизации 40,92 МГц, то есть 4 отсчета на один чип P‑кода. Это удовлетворяет теореме Котельникова, поскольку ширина главного лепестка спектра P‑кода составляет 20,46 МГц (от –10,23 до +10,23 МГц), и позволяет увидеть несколько боковых лепестков.

navdata — вектор из двух случайных битов (0 или 1), имитирующий навигационное сообщение GPS (в реальности скорость передачи данных 50 бит/с, P‑код модулируется этими битами через BPSK).

p = welch_pgram(waveform_vec, 1024, 0; 

fs=gpswaveobj.SampleRate,
window = hann,
onesided=false)

Метод Уэлча с окном Ханна дает сглаженную оценку, усредняя периодограммы отдельных сегментов.

Подключение вспомогательных файлов и библиотек

Для визуализации результатов работы блока необходимо подключить сторонние библиотеки.

let
    installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
    list_packages = ["DSP"]
    for pack in list_packages
        pack in installed_packages || Pkg.add(pack)
    end
end

import DSP
using EngeeDSP, EngeeSatellites

# Создание и настройка генератора
gpswaveobj = gpsWaveformGenerator()
gpswaveobj.EnablePCode = true
gpswaveobj.SampleRate = 4 * 10.23e6

# Генерация сигнала
navdata = rand([0,1], 2, 1)
waveform = gpswaveobj(navdata)

# Преобразуем матрицу сигнала в одномерный вектор
waveform_vec = vec(waveform)

# nfft=1024 — размер БПФ (чем больше, тем выше разрешение)
# Задаем длину окна n=1024 и перекрытие noverlap=0 (0%)
p = DSP.welch_pgram(
     waveform_vec, 1024, 0; 
     fs=gpswaveobj.SampleRate,
     window = DSP.hann,
     onesided=false
)

# Извлекаем частоты и спектральную плотность мощности
freqs = DSP.freq(p)      # вектор частот в Гц
psd = DSP.power(p)       # вектор PSD в Вт/Гц

# Центрируем спектр
freqs_shifted = EngeeDSP.Functions.fftshift(freqs)
psd_shifted = EngeeDSP.Functions.fftshift(psd)

# Построение спектра в dBW/Гц
plot(freqs_shifted/1e6, 10*log10.(psd_shifted .+ eps()), 
     xlabel="Частота (МГц)", 
     ylabel="Мощность (dBW/Гц)", 
     title="Спектр GPS P-кода",
     legend=false)

Сигнал L1C

L1C — Этот сигнал открыт для гражданского использования и передается на частоте L1. Он использует схему модуляции с двоичным смещением несущей (BOC)

Глобальное отличие от предыдущего скрипта - наличие двух спутников и соответственно, двух сигналов. В основе лежат те же основные принципы, что и в генерации P-кода.

prn = [4, 70]                 # идентификаторы спутников
numsat = length(prn)          # 2 спутника
numbits = 100                 # длина сообщения в битах
navdata = rand([0, 1], numbits, numsat)
fs = 25e6                     # частота дискретизации 25 МГц

gpswaveobj = gpsWaveformGenerator(; SignalType="l1c", PRNID=prn, SampleRate=fs)
waveform = gpswaveobj(navdata)  # Это матрица размером N×2, где N - количесство отсчетов

# Настройка параметров для Welch (сглаживание спектра)
# Если сигнал короче 2048 точек, используем его полную длину
seg_len = min(2048, size(waveform, 1))

# Построение графика
# Создаем пустую область с подписями осей и легендой
plot_obj = plot(; xlabel="Частота (МГц)", ylabel="Мощность (dBW/Гц)", 
                title="Спектр GPS L1C (PRN 4 и 70)", legend=:topleft)

# Цикл по двум каналам (спутникам)
for i in 1:numsat
    # Берем только i-й столбец (канал)
    channel = waveform[:, i]
    
    # Считаем спектр методом Уэлча
    p_i = DSP.welch_pgram(channel, seg_len, 0; fs=fs, onesided=false)
    
    # Извлекаем и центрируем частоты
    freqs_i = DSP.freq(p_i)
    psd_i = DSP.power(p_i)
    freqs_shifted = EngeeDSP.Functions.fftshift(freqs_i)
    psd_shifted = EngeeDSP.Functions.fftshift(psd_i)
    
    # Добавляем график
    plot!(plot_obj, freqs_shifted / 1e6, 10 * log10.(psd_shifted .+ eps()), 
          label="PRN $(prn[i])")
end

# Отображаем график
display(plot_obj)

Сигнал L2C

SignalType="l2c" – основой является гражданский сигнал L2C. Он передаётся на частоте 1227,60 МГц и содержит два кода, мультиплексированных по времени:

CM (Civil Moderate) – код длиной 10 230 чипов (период 20 мс), модулирован данными CNAV (25 бит/с).
CL (Civil Long) – код длиной 767 250 чипов (период 1,5 с), без данных (пилотная компонента).
Оба кода имеют чиповую скорость 511,5 кчип/с, но мультиплексируются таким образом, что суммарная чиповая скорость составляет 1,023 Мчип/с (как у C/A).

waveform = gpswaveobj(lnavdata, cnavdata)

В отличие от предыдущего примера L1C, здесь данные передаются кортежем из двух матриц, потому что L2C использует два разных навигационных потока: LNAV и CNAV. Генератор использует их для модуляции соответствующих компонентов.

Число спутников увеличено до четырёх, но в остальном принцип работы остаётся тем же.

# Задаем параметры
prn = [7, 11, 20, 28]
numsat = length(prn)
numbits = 10

# Генерация случайных данных LNAV и CNAV (матрицы 10x4)
lnavdata = rand([0, 1], numbits, numsat)
cnavdata = rand([0, 1], numbits, numsat)

# Создаем генератор с настройками
gpswaveobj = gpsWaveformGenerator(;
    SignalType="l2c",
    PRNID=prn,
    SampleRate=15e6,
    EnablePCode=true
)

# Генерация сигнала
# ВНИМАНИЕ: Передаем данные в виде кортежа (аналог cell-массива MATLAB)
waveform = gpswaveobj(lnavdata, cnavdata)

# Проверка размера
println("Размер сгенерированного сигнала: ", size(waveform))
# Ожидаемый вывод: (3000000, 4)

# Настройка БПФ для сглаживания
seg_len = min(2048, size(waveform, 1))

# Построение графика
plot_obj = plot(; xlabel="Частота (МГц)", ylabel="Мощность (dBW/Гц)", 
                title="Спектр GPS L2C (PRN 7, 11, 20, 28)", legend=:topleft)

for i in 1:numsat
    channel = waveform[:, i] # Берем канал i-го спутника
    
    # Вычисляем спектр
    p_i = DSP.welch_pgram(channel, seg_len, 0; fs=15e6, onesided=false, nfft=seg_len)
    
    freqs_i = DSP.freq(p_i)
    psd_i = DSP.power(p_i)
    freqs_shifted = EngeeDSP.Functions.fftshift(freqs_i)
    psd_shifted = EngeeDSP.Functions.fftshift(psd_i)
    
    # Рисуем график для этого спутника
    plot!(plot_obj, freqs_shifted / 1e6, 10 * log10.(psd_shifted .+ eps()), 
          label="PRN $(prn[i])")
end

display(plot_obj) # Отображаем график

Размер сгенерированного сигнала: (3000000, 4)

Сигнал L5

Сигнал GPS L5 — это один из современных гражданских сигналов системы GPS, передаваемый на частоте 1176,45 МГц (диапазон L5). Он был введён в рамках программы модернизации GPS и предназначен в первую очередь для обеспечения высокой точности и целостности, необходимых для авиационных и других ответственных приложений (Safety-of-Life).

Основное тело скрипта остаётся схожим с предыдущими, однако к сигналу добавляются искажения, вызванные доплеровским смещением, возникающим из-за движения спутника и приемника. Введение разных сдвигов (+3 кГц и –4,6 кГц) разносит их спектры, что приводит к появлению двух близких максимумов внутри общей огибающей – интересный и наглядный эффект.

Внесение искажений: доплеровский сдвиг частоты

# Внесение искажений 

# Частотный сдвиг
t = (0:size(waveform, 1) - 1) / fs                 # Вектор времени
freq_offsets = [3e3, -4.6e3]                       # Сдвиги для двух каналов

# Создаем комплексную экспоненту и умножаем на сигнал
fwave = waveform .* exp.(im * 2 * pi * t .* freq_offsets')

# Внесение временной задержки
delay_val = 5.23e-4
delay_samples = round(Int, delay_val * fs)         # 13075 отсчетов

# Физическая задержка: добавляем нули в начало и отбрасываем конец
dwave = vcat(zeros(ComplexF64, delay_samples, numsat), fwave[1:end-delay_samples, :])

# Суммирование двух каналов
bbwave = sum(dwave, dims=2)                        # Размер (10_000_000 x 1)
bbwave_vec = vec(bbwave)                           # Превращаем в одномерный вектор

# Генерация сигнала GPS L5
prn = [185, 189]
numsat = length(prn)
numbits = 40
msg = rand([0, 1], numbits, numsat)

fs = 25e6
gpswaveobj = gpsWaveformGenerator(;
    SignalType="l5",
    PRNID=prn,
    InitialTime=42123,
    SampleRate=fs
)

# Генерируем сигнал. Получаем матрицу (10_000_000 x 2)
waveform = gpswaveobj(msg)
println("Размер исходного сигнала: ", size(waveform))


# Внесение искажений 
# Частотный сдвиг
t = (0:size(waveform, 1) - 1) / fs                 # Вектор времени
freq_offsets = [3e3, -4.6e3]                       # Сдвиги для двух каналов

# Создаем комплексную экспоненту и умножаем на сигнал
fwave = waveform .* exp.(im * 2 * pi * t .* freq_offsets')

# Внесение временной задержки
delay_val = 5.23e-4
delay_samples = round(Int, delay_val * fs)         # 13075 отсчетов

# Физическая задержка: добавляем нули в начало и отбрасываем конец
dwave = vcat(zeros(ComplexF64, delay_samples, numsat), fwave[1:end-delay_samples, :])

# Суммирование двух каналов
bbwave = sum(dwave, dims=2)                        # Размер (10_000_000 x 1)
bbwave_vec = vec(bbwave)                           # Превращаем в одномерный вектор
println("Размер результирующего сигнала (BB): ", size(bbwave_vec))


# Визуализация спектра
seg_len = min(2048, length(bbwave_vec))

# Вычисляем спектр
p = DSP.welch_pgram(bbwave_vec, seg_len, 0; fs=fs, onesided=false, nfft=seg_len)

freqs = DSP.freq(p)
psd = DSP.power(p)

# Центрируем спектр
freqs_shifted = EngeeDSP.Functions.fftshift(freqs)
psd_shifted = EngeeDSP.Functions.fftshift(psd)

# Рисуем график
plot(freqs_shifted / 1e6, 10 * log10.(psd_shifted .+ eps()),
     xlabel="Частота (МГц)", 
     ylabel="Мощность (dBW/Гц)",
     title="Спектр GPS L5",
     legend=false)

Размер исходного сигнала: (10000000, 2)
Размер результирующего сигнала (BB): (10000000,)

Вывод

Представленные примеры позволяют оценить работу блока gpsWaveformGenerator, его функционал и влияние настроек на результирующий GPS-сигнал на выходе.

Генератор сигналов GPS