Frequency adaptation under interference

This example shows how to model radar systems with frequency adaptation under interference.

Connecting libraries and auxiliary functions

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

using DSP

function TaylorSpectrumWindow(N, nbar = 4, sll = 30)
    function calcFm(m, sp2, A, NBAR)
        n = Array(1:NBAR-1)
        p = [Array(1:m-1); Array((m+1):(NBAR-1))]
        Num = prod((1 .- (m^2/sp2)./(A^2 .+(n .-0.5).^2)));
        Den = prod((1 .- m.^2 ./p.^2));
    
        return ((-1)^(m+1).*Num)./(2 .*Den)
    end

    A = acosh((10^(sll/20)))/π
    sp2 = nbar^2 / (A^2 + (nbar - 0.5)^2)
    Fm = zeros(nbar - 1)
    sum = zeros(N)
    k = Array(0:N-1)
    xi = @. (k-0.5*N+0.5)/N
    for m = 1:(nbar-1)
        Fm[m] = calcFm(m,sp2,A,nbar);
        @. sum = Fm[m]*cos(2*pi*m*xi)+sum;
    end
    @. sum = 2sum + 1
    return sum
end

function calc_spectrogramm(x,t_step;down_lim = -40)
    plotlyjs()
    out,f,t = EngeePhased.Functions.spectrogram(x;
        window = DSP.kaiser(128,10),
        noverlap=120,nfft=256,fs=1/t_step,
        centered=true,freqloc = "yaxis",out = :data
    )
    out_sig = map(x -> x < down_lim ? down_lim : x,DSP.pow2db.(abs.(out)))
    fig = Plots.heatmap(
        t[:].*1e6,f*1e-6,out_sig,color= :jet,
        gridalpha=0.5,colorbar_title = "Мощность, дБВт",
        xticks=0:50:250,yticks=(-1:0.2:1),
    )
    Plots.xlabel!("Время, мкс")
    Plots.ylabel!("Частота Доплера, МГц")
    Plots.title!("Cпектрограмма сигнала")

    return fig
end;

Introduction

Active phased electronically scanned antenna arrays can help solve many applications using the same antenna array hardware. They can be used in radar, REB and communication systems. However, the environment in which these types of systems operate is complex and can introduce unwanted interference. For example, a jamming device on a repeater may repeat a received radar signal and retransmit it to confuse the radar. Frequency adaptation can be an effective method of counteracting signals generated by interference sources and contribute to the efficient operation of these systems.

This example simulates a scenario with a stationary monostatic radar and a moving aircraft target (a similar scenario is discussed in the example Tuning the centre frequency under interference). The aircraft will then generate a spoofing signal that confuses the radar. Once the radar detects the interference source, centre frequency hopping techniques can be used to allow the radar to bypass the interference.

1. Modelling the system in an interference-free environment

Assume that there is a monostatic X-band radar at the origin.

fc = 10e9; # несущая частота работы РЛС X-диапазона, Гц 
fs = 2e6; # частота дискретизации модели, Гц
c = 3e8; # скорость распросранения сигнала, м/с
lambda = c/fc; # длина волны, м

radar_pos = [0;0;0]; # положение РЛС
radar_vel = [0;0;0]; # скорость РЛС

The receiving antenna consists of 64 elements (8x8) of an equidistant rectangular antenna array with a half wavelength spacing between the elements:

antenna = EngeePhased.URA(
    Element=EngeePhased.CosineAntennaElement(), # косинусный элемент АР 
    Size=[8 8], # размер антенной решетки
    ElementSpacing=[lambda/2 lambda/2],  # расстояние между элементами
)

URA:
    Element=CosineAntennaElement: FrequencyRange=[0.0 1.0e20] CosinePower=[1.5 1.5]
    Size=[8 8]
    ElementSpacing=[0.015 0.015]
    Lattice=Rectangular
    ArrayNormal=x
    Taper=1

To better characterise the directivity of the AR, we form a 2-dimensional Taylor window function:

taper = TaylorSpectrumWindow(8); # взвешанная фукция окно Тэйлора
taperURA = taper.*taper' # формирование 2д взвешанной функции
antenna.Taper = taperURA # взвешивание АР

# Визуализация оконной функции
surface(taperURA,title="2-мерная оконная фукция Тейлора")

The antenna array directivity diagram can be plotted using the function pattern:

pattern(antenna,fc,Type="powerdb")
plot!(title="ДН прямоугольной АР",colorbar_title="Мощность, дБВт")

In the radar system under consideration, the signal source will be a linear frequency modulated (LFM) signal:

release!(antenna) 
# Формирование системных объектов элементов РЛС 

bw = 1e5 # полоса сигнала в частотной области
pw = 1e-4 # длительной импульса

wav = EngeePhased.LinearFMWaveform( # генератор ЛЧМ
    SampleRate=fs,  # частота дискретизации модели
    PulseWidth= pw, # длительность импульса
    SweepBandwidth= bw, # полоса сигнала в частотной области
    PRF=4e3, # частота следования импульсов
    FrequencyOffsetSource="Input port" # включение входа для настройки центральной частоты
)
tx = EngeePhased.Transmitter( # передатчик
    Gain=20, # усиление, дБ
    PeakPower=500 # пиковая мощность, Вт
)
txArray = EngeePhased.WidebandRadiator( # Широкополосная передающая АР
    SampleRate=fs,  # частота дискретизации модели
    Sensor=antenna, # задание геометрии АР
    CarrierFrequency=fc # несущая частота
)
rxArray = EngeePhased.WidebandCollector(
    SampleRate=fs,  # частота дискретизации модели
    Sensor=antenna, # задание геометрии АР
    CarrierFrequency=fc # несущая частота
)
rxPreamp = EngeePhased.ReceiverPreamp( # предусилитель
    Gain=10, # усиление, дБ
    NoiseFigure=5, # фактор шума, дБ 
    SampleRate=fs,  # частота дискретизации модели
);

Next, let's create the target and distribution environment objects:

target = EngeePhased.RadarTarget(
    MeanRCS=100, # среднее значение ЭПР цели
    OperatingFrequency=fc # несущая частота 
)

target_pos = [8000;1000;1000] # начальное положение цели
target_vel = [100;0;0] # вектор скорости цели

envout = EngeePhased.WidebandFreeSpace(
    TwoWayPropagation=false, # Учет двунаправленного распространения
    SampleRate = fs, # частота дискретизации
    OperatingFrequency=fc, # несущая частота
    PropagationSpeed=c # скорость распространения
)
envin = deepcopy(envout); # дублирование канала обратного распространения

Using the previously defined system objects of the system, we calculate the reflected signal from the target in the absence of interference:

tgtRng, tgtAng = rangeangle(target_pos, radar_pos)           
x = deepcopy(wav(0)) # Генератор
xt = tx(x) # передатчик
xtarray = txArray(xt,tgtAng) # передающая антенная решетка
yp = envout(xtarray,radar_pos,target_pos,radar_vel,target_vel) # канал прямого распространения
yr = target(yp) # цель
ye = envin(yr,target_pos,radar_pos,target_vel,radar_vel) # среда распространения
yt = rxArray(ye,tgtAng) # приемная антенная решетка
yt = rxPreamp(yt);  # приемник

Let's calculate the direction and realise the beam forming with the help of BeamscanEstimator2D and SubbandPhaseShiftBeamformer

estimator = EngeePhased.BeamscanEstimator2D(
    SensorArray=antenna, # геометрия АР 
    DOAOutputPort=true, # Включение выхода уголов
    OperatingFrequency=fc, # несущая частота
    NumSignals=1, # количество сигналов
    AzimuthScanAngles=-40:40,  # диапазон сканирования по азимуту
    ElevationScanAngles=-60:60 # диапазон сканирования по углу места
)

doa = estimator(yt)[2] # расчет оценки пеленга с помощью СО BeamscanEstimator2D

beamformer = EngeePhased.SubbandPhaseShiftBeamformer(
    SensorArray=antenna, # геометрия АР 
    OperatingFrequency=fc, # несущая частота
    DirectionSource="Input port", # способ задания углов визирования
    SampleRate=fs, # частота дискретизации 
    WeightsOutputPort=true # весовые коэффициенты 
)

ybf = beamformer(yt,doa)[1]; # применения алгоритма формирования СО SubbandPhaseShiftBeamformer

Let's perform matched filtering, for this purpose let's calculate filter coefficients by means of the function getMatchedFilter

mfcoeff1 = getMatchedFilter(wav) # Коэффициентов согласованнго фильтра
mf1 = EngeePhased.MatchedFilter(Coefficients=mfcoeff1);

Apply matched filtering to the generated beam and visualise the result of radar operation:

y1 = mf1(ybf) # применение согласованного фильтра

nSamples = wav.SampleRate/wav.PRF # количество отсчетов в периоде импульса
t = ((0:nSamples-1).-(length(mfcoeff1)-1))./fs # временная сетка задержек
r = t*c/2 # пересчет из задержки в дальность

# Результат работы РЛС
plot(r/1000,abs.(y1),xlabel="Дальность, км",ylabel="Амплитуда",lab="",xticks=-1e6:5:1e6,title="Выход СФ")

2. Adding an interference component

Let us now simulate the interfering signal by simulating the useful signal and time delay:

 # имитация помехового сигнала путемдублирование генератор сигнала
jwav = deepcopy(wav)

xj = jwav(0) # формирование помеховой составляющей
Npad = ceil(Int64,3500/(c/fs)) # количество элементов задержки (1 отсчет - 150 метров)
xj = circshift(xj,Npad) # сдвиг сигнала на распространения
txjam = EngeePhased.Transmitter(
    Gain=10, # усиление помехи
    PeakPower=5 # мощность помехи
)

# Расчет отраженного сигнала + помеха
xj = txjam(xj) # излучение помехи
ye = envin(yr.+xj,target_pos,radar_pos,target_vel,radar_vel) # канал прямого распространения 
yt = rxArray(ye,tgtAng) # приемная АР                               
yt = rxPreamp(yt) # предусилитель                      
ybfj = beamformer(yt,doa)[1] # формирователь луча
y1j = mf1(ybfj); # фильтрация суммарного сигнала 

Let's see how the output graph of the model is transformed:

plot(r/1000,abs.(y1j),xlabel="Дальность, км",ylabel="Амплитуда",
    lab="",xticks=-1e6:5:1e6,title="Выход СФ при воздействии помехи")

We can notice that after adding the interference, 2 correlation peaks are formed at the output of the matched filter, with one peak from the interference dominating and shading the true position of the target.

3. Tuning the centre frequency

To reduce the influence of the interference component, a frequency tuning approach can be used, which consists in changing the centre frequency of the signal by the value of several signal bands:

deltaf = 2.5 * bw # смещение центральной частоты на 250 кГц
xh = deepcopy(wav(deltaf)); # сигнал со смещенной центральной частотой

Let's construct a spectrogram of the sum of 2 signals - unbiased and biased LFM signal:

calc_spectrogramm(x.+xh,1/fs;down_lim = -40)

The spectrogram shows that the spectra of the signals are shifted relative to each other by about 250 kHz.

Let us repeat the modelling of the radar operation scenario for the shifted signal:

xth = tx(xh) # передатчик                         
xtharray = txArray(xth, tgtAng) # передающая ФАР
yph = envout(xtharray,radar_pos,target_pos,radar_vel,target_vel) # канал прямого распространения
yrh = target(yph) # цель

yeh = envin(yrh.+xj,target_pos,radar_pos,target_vel,radar_vel) # канал обратного распространения
yth = rxArray(yeh,tgtAng) # приемная ФАР                              
yth = rxPreamp(yth) # предусилитель                            
ybfh = beamformer(yth,doa)[1]; # формирование луча

A bandpass filter tuned to an offset frequency with a bandwidth of the WL is required to extract the useful component and remove the interfering component.

Set the coefficients of the low-pass Butterworth filter of order 9 and shift the centre bandwidth to 250 kHz:

include("$(@__DIR__)/data_filter_coeff.jl") # считывание коэффициентов из jl файла
# смещение центральной частоты на 250 кГц
bf2 = buttercoef.*exp.(1im*2*pi*deltaf*Vector(0:length(buttercoef)-1)/fs); 

Let's realise the filtering in the WF band:

mfcoeff2 = getMatchedFilter(wav,deltaf)
mf2 = EngeePhased.MatchedFilter(Coefficients=mfcoeff2)

# выделение полосым фильтром и согласованная фильтрация
yb2 = mf2(DSP.filt(bf2,1,ybfh))

plot(r./1000,abs.(yb2),lab="")
xlabel!("Дальность, км")
ylabel!("Амплитуда")
title!("Выход СФ после перестройки")

From the graph of the signal at the output of the NF, it can be seen that the correlation peak of the noise component is successfully suppressed. Thus, the application of the centre frequency tuning allowed

Conclusion

In the case study, the modelling of a frequency tunable radar was examined using the EngeePhased library system objects. This approach reduced the impact of the interfering component and successfully isolated the useful signal.