Engee documentation
Notebook

Range and Doppler estimation using pulse compression

This example demonstrates the effect of pulse compression when the transmitted pulse is correlated with the received signal. Radar and sonar systems use pulse compression to improve signal-to-noise ratio (SNR) and range resolution by reducing the duration of reflected signals. This example also demonstrates Doppler processing, in which the radial velocity of a target is determined from the Doppler shift that occurs as the target moves.

Determining Target Range and Velocity

The radar system shown below transmits a pulsed sequence of linear frequency modulated (LFM) waves to a moving target and receives echo signals. By applying matched filtering and Doppler processing, the radar system can effectively determine the range and velocity of the target.

Let's specify the requirements for the radar system. In this example, the carrier frequency fc = 3 GHz and the sampling frequency fs = 1 mHz are used.

In [ ]: 
using Pkg; Pkg.add("DSP")
using FFTW,DSP
fc = 3e9; # несущая частота, Гц
fs = 1e6; # частота дискретизации, Гц
c = 299_792_458; # м/с

Create system objects to model the radar system. The system is a monostatic system. The transmitter is located at (0,0,0,0) m and is stationary and the target is at (5000,5000,0) m with a velocity of (25,25,0) m/s and has an effective scattering area (EDA) of 1 square metre.

In [ ]: 
antenna = EngeePhased.IsotropicAntennaElement(FrequencyRange=[1e8 10e9]); # антенный элемент
transmitter = EngeePhased.Transmitter(Gain=20,InUseOutputPort=true); # передатчик
txloc = [0;0;0]; # положение передатчика
tgtloc = [5000;5000;0]; # координаты цели, дальность ~= 7071 м
tgtvel = [25;25;0]; # скорость цели,  относительная скорость ~= 35.4 м/с
target = EngeePhased.RadarTarget(Model="Nonfluctuating",MeanRCS=1,OperatingFrequency=fc);
antennaplatform = EngeePhased.Platform(InitialPosition=txloc); # объект передатчика в пространстве
targetplatform = EngeePhased.Platform(InitialPosition=tgtloc,Velocity=tgtvel); # объект цели в пространстве
radiator = EngeePhased.Radiator( # излучатель
   PropagationSpeed=c, # скорость распространения
   OperatingFrequency=fc, # несущая частота cигнала
   Sensor=antenna # антенный элемент
);
channel = EngeePhased.FreeSpace( # канал распространения
   PropagationSpeed=c, # скорость распространения
   OperatingFrequency=fc, # несущая частота
   TwoWayPropagation=false # двухнаправленное направление
);
collector = EngeePhased.Collector(
   PropagationSpeed=c, # скорость распространения
   OperatingFrequency=fc, # несущая частота cигнала
   Sensor=antenna # антенный элемент
);
receiver = EngeePhased.ReceiverPreamp( # передатчик
   NoiseFigure=0, # Фактор шума
   EnableInputPort=true, # входной порт синхронизации с передатчиком
   SeedSource="Property", # источник начального значения генератора ПСП для фазового шума
   Seed=2024 # значение генератора ПСП для фазового шума
);

Generate an LFM signal with a pulse width of 10 µs, a pulse repetition rate of 10 kHz and a sweep bandwidth of 100 kHz. Matched filter coefficients are generated from this signal.

In [ ]: 
waveform = EngeePhased.LinearFMWaveform( # генератор ЛЧМ сигнала
    PulseWidth=10e-6, # длительность импульса
    PRF=10e3, # частота следования импульсов
    OutputFormat="Pulses", # тип выхода
    NumPulses=1, # количество импульсов
    SweepBandwidth=1e5 # время наростания пилы
);
wav = waveform(); # формирование ЛЧМ-сигнала 
maxrange = c/(2*waveform.PRF); # Максимальная одназначная дальность обнаружения цели
SNR = npwgnthresh(1e-6,1,"noncoherent"); # расчет ОСШ на основе порога обнаружения
lambda = c/fc; # длина волны
tau = waveform.PulseWidth; # длительность импульса
Ts = 290; # Шумовая температура приемника

# Согласованный фильтр
Mfilter = EngeePhased.MatchedFilter( 
   Coefficients=getMatchedFilter(waveform), # коэффициенты фильтра
   GainOutputPort=true # выходной порт коэффициента передачи
);

The system emits 64 pulses to improve the Doppler resolution. The signal data will be stored in the variable rxsig, with the Data Matrix storing fast time samples (time within each pulse) in each column and slow time samples (time between pulses) in each row.

In [ ]: 
numPulses = 64; # количество импульсов
rxsig = zeros(ComplexF64,length(wav),numPulses); # выделение памяти под выходной сигнал

for n = 1:numPulses
    tgtloc,tgtvel = targetplatform(1/waveform.PRF); # параметры движения цели
    tgtrng,tgtang = rangeangle(tgtloc,txloc); # параметры движения радара
    
    txsig, txstatus = transmitter(wav); # формирование передающего сигнала 
    txsig = radiator(txsig,tgtang); # излучение ЭМ-волны в пространство
    txsig = channel(txsig,txloc,tgtloc,[0;0;0],tgtvel); # канал: передатчик-цель
    txsig = target(txsig); # Отражение сигнала с учетом ЭПР цели
    txsig = channel(txsig,tgtloc,txloc,tgtvel,[0;0;0]); # канал: цель-приемник
    txsig = collector(txsig,tgtang); # формирование сигнала в приемной антенне
    rxsig[:,n] = receiver(txsig,xor.(txstatus,1)); # предусиление принятого сигнала
end

prf = waveform.PRF; # частота следования импульсов
fs = waveform.SampleRate; # частота дискретизации
response = EngeePhased.RangeDopplerResponse(
    DopplerFFTLengthSource="Property", # задание источника длины БПФ
    DopplerFFTLength=2048, # Длина БПФ по доплеру
    SampleRate=fs, # частота дискретизации
    DopplerOutput="Speed", # тип выхода 
    OperatingFrequency=fc, # несущая частота
    PRFSource="Property", # источник задания ЧСИ
    PRF=prf # частота следования импульсов
);
filt = getMatchedFilter(waveform); # коэффициенты согласованного фильтра
resp,rng_grid,dop_grid = response(rxsig,filt); # отклик по дальности и скорости
rng_grid = range(extrema(rng_grid)...,length = length(rng_grid)) # cетка дальности
dop_grid = range(extrema(dop_grid)...,length = length(dop_grid)) # cетка скорости

# построение долероского-дальностного портрета
heatmap(dop_grid,rng_grid,DSP.pow2db.((abs.(resp)).^2),ylims=(0,12e3), 
    fontfamily=font(14,"Computer Modern"), xticks=(-200:50:200),
    colorbar_title = "Мощность,дБВт",
    title=("Дальностно-доплероский портрет"),
    margin=5Plots.mm
)
xlabel!("Скорость, м/с")
ylabel!("Дальность, м")
Out[0]:

The response indicates that the target is travelling at about -40 m/s and that the target is moving away from the transmitter as the velocity is negative.

Let's calculate the range estimate corresponding to the signal propagation speed. The plot of the first slow time sample shows the largest peak around 7000 m, which corresponds to the range speed response plot.

In [ ]: 
fasttime = unigrid(0,1/fs,1/prf,"[)"); # сетка быстрого времени
rangebins = (c*fasttime/2); # разрешение по дальности

# Построение зависимости амплитуда от дальности
plot(rangebins.*1e-3,
    abs.(rxsig[:,1].*1e6),
    lab="",xlabel="Дальность, км",
    ylabel="Амплитуда, мкВ"
)
Out[0]:

Range determination

Let's set the detection threshold at which the probability of false alarm is less than $10^{-6}$. Perform incoherent accumulation of 64 pulses, assuming that the signal is in white Gaussian noise. By the largest peak exceeding the threshold value, determine the target range estimate.

In [ ]: 
pfa = 1e-6; # вероятность ложной тревоги
NoiseBandwidth = 5e6/2; # шумовая полоса

# Расчет мощности шума
npower = noisepow(
    NoiseBandwidth, # Шумовая полоса 
    receiver.NoiseFigure, # фактор шума
    receiver.ReferenceTemperature # шумовая температура
);

thresh = npwgnthresh(pfa,numPulses,"noncoherent"); # порог обнаружение
thresh = npower*DSP.db2pow(thresh); # чувствительность приемника

pks,range_detect = findpeaks(pulsint(rxsig,"noncoherent"),"MinPeakHeight",thresh); # обнаружение пиков
range_estimate = rangebins[range_detect[1]]; # определение оценки дальности
print("Дальность: $(round(range_estimate;sigdigits=6)), м")

Дальность: 7344.92, м

Determination of velocity

We select a range sample that includes the detected target and plot the power spectral density estimation using the periodagram function.

In [ ]: 
out = DSP.periodogram(rxsig[range_detect[1],:][:]; onesided=false, 
    nfft=256, fs=prf, window=nothing)
spec_power = fftshift(DSP.pow2db.(out.power))
freq = fftshift(fftfreq(256,prf))
plot(freq.*1e-3,spec_power,lab="",gridalpha=0.15)
xlabel!("Частота,кГц", fontfamily="Computer Modern")
ylabel!("Мощность,дБВт", fontfamily="Computer Modern")
title!("Спектральная плотность мощности", fontfamily="Computer Modern")
Out[0]:

The peak frequency corresponds to the Doppler shift divided by 2, which can be converted to target velocity. A positive velocity means that the target is approaching the transmitter and a negative velocity means that the target is moving away from the transmitter.

In [ ]: 
Y,I = maximum(spec_power),argmax(spec_power);
lambda = c/fc;
tgtspeed = dop2speed(freq[I]/2,lambda);
print("Оценка доплеровского сдвига: $(freq[I]/2) Гц")

Оценка доплеровского сдвига: -351.5625 Гц
In [ ]: 
print("Оценка скорости: $(round(tgtspeed;sigdigits=4)) м/с")

Оценка скорости: -35.13 м/с