Engee documentation
Notebook

Estimation of range and Doppler using pulse compression

This example demonstrates the pulse compression effect when the transmitted pulse correlates 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 the target is determined by the Doppler shift that occurs when the target moves.

Determining the range and speed of targets

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

Let's set 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 simulate the radar system. The system is monostatic. The transmitter is located at (0,0,0) m and is stationary, and the target is located at (5000,500,0) m at a speed of (25,25,0) m/s and has an effective scattering area (ESR) of 1 sq. m.

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 # значение генератора ПСП для фазового шума
);

Create an FM signal with a pulse duration of 10 microseconds, a pulse repetition rate of 10 kHz, and a scan bandwidth of 100 kHz. Matched filter coefficients are generated based on 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 # выходной порт коэффициента передачи
);

To improve the Doppler resolution, the system emits 64 pulses. The signal data will be stored in the variable rxsig, and the Data Matrix stores 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 moving at a speed of about -40 m/s, and that the target is moving away from the transmitter because the speed is negative.

We will calculate the range estimate corresponding to the propagation velocity of the signal. The graph of the first slow time sample shows the highest peak in the region of 7000 m, which corresponds to the graph of the response to the range speed.

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 detection

Let's set a detection threshold at which the probability of a false alarm will be lower. . Let's perform an incoherent accumulation of 64 pulses, assuming that the signal is in white Gaussian noise. Based on the largest peak exceeding the threshold value, we will 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, м

Speed detection

Let's select a range sample that includes the detected target, and plot the power spectral density estimate 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 the target speed. A positive velocity means that the target is approaching the radiator, and a negative velocity means that the target is moving away from the radiator.

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 м/с