Frequency adaptation when exposed to interference

This example shows how to model radar systems with frequency adaptation when exposed to 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 = "Power, dBW",
        xticks=0:50:250,yticks=(-1:0.2:1),
    )
    Plots.xlabel!("Time, iss")
    Plots.ylabel!("Doppler frequency, MHz")
    Plots.title!("Signal spectrogram")

    return fig
end;

Introduction

Electronically scanned active phased array antennas can help solve a variety of tasks using the same antenna array equipment. They can be used in radars, electronic warfare and communication systems. However, the environment in which these types of systems operate is complex and may introduce unwanted interference. For example, a jammer on a repeater may repeat the received radar signal and retransmit it to confuse the radar. Frequency adaptation can be an effective method of countering signals generated by interference sources and contributing to the efficient operation of these systems.

In this example, a scenario with a stationary monostatic radar and a moving aircraft target is modeled (a similar scenario is considered in the example Center frequency tuning under impact помех). The aircraft will then generate a replacement signal that confuses the radar. As soon as the radar detects the source of interference, center frequency tuning techniques can be used to allow the radar to bypass interference.

1. Simulation of the system in a non-interference environment

Let's assume that an X-band monostatic radar is located at the origin.

fc = 10e9; # X-band radar carrier frequency, Hz
fs = 2e6; # sampling rate of the model, Hz
c = 3e8; # signal propagation speed, m/s
lambda = c/fc; # wavelength, m

radar_pos = [0;0;0]; # Radar position
radar_vel = [0;0;0]; # Radar speed

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

antenna = EngeePhased.URA(
    Element=EngeePhased.CosineAntennaElement(), # the cosine element of AR
    Size=[8 8], # antenna array size
    ElementSpacing=[lambda/2 lambda/2],  # the distance between the elements
)

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 characterize the directivity of the AR, we will form a 2-dimensional Taylor window function.:

taper = TaylorSpectrumWindow(8); # weighted fuction Taylor window
taperURA = taper.*taper' # формирование 2д взвешанной функции
antenna.Taper = taperURA # weighing the AP

# Visualization of the window function
surface(taperURA,title="2-dimensional Taylor window fuction")

The radiation pattern of the antenna array can be constructed using the function pattern:

pattern(antenna,fc,Type="powerdb")
plot!(title="The bottom of the rectangular AR",colorbar_title="Power, dBW")

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

release!(antenna) 
# Formation of system objects of radar elements

bw = 1e5 # the frequency domain band of the signal
pw = 1e-4 # Long-lasting pulse

wav = EngeePhased.LinearFMWaveform( # LCHM generator
    SampleRate=fs,  # sampling rate of the model
    PulseWidth= pw, # pulse duration
    SweepBandwidth= bw, # the frequency domain band of the signal
    PRF=4e3, # pulse repetition rate
    FrequencyOffsetSource="Input port" # turning on the input to adjust the center frequency
)
tx = EngeePhased.Transmitter( # The transmitter
    Gain=20, # gain, dB
    PeakPower=500 # peak power, W
)
txArray = EngeePhased.WidebandRadiator( # Broadband transmission system
    SampleRate=fs,  # sampling rate of the model
    Sensor=antenna, # defining the AR geometry
    CarrierFrequency=fc # carrier frequency
)
rxArray = EngeePhased.WidebandCollector(
    SampleRate=fs,  # sampling rate of the model
    Sensor=antenna, # defining the AR geometry
    CarrierFrequency=fc # carrier frequency
)
rxPreamp = EngeePhased.ReceiverPreamp( # The preamp
    Gain=10, # gain, dB
    NoiseFigure=5, # noise factor, dB
    SampleRate=fs,  # sampling rate of the model
);

Next, we'll create the target and distribution environment objects.:

target = EngeePhased.RadarTarget(
    MeanRCS=100, # the average EPR value of the goal
    OperatingFrequency=fc # carrier frequency
)

target_pos = [8000;1000;1000] # the initial position of the goal
target_vel = [100;0;0] # The target's velocity vector

envout = EngeePhased.WidebandFreeSpace(
    TwoWayPropagation=false, # Accounting for bidirectional propagation
    SampleRate = fs, # sampling rate
    OperatingFrequency=fc, # carrier frequency
    PropagationSpeed=c # The speed of propagation
)
envin = deepcopy(envout); # duplication of the reverse distribution channel

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

tgtRng, tgtAng = rangeangle(target_pos, radar_pos)           
x = deepcopy(wav(0)) # The generator
xt = tx(x) # The transmitter
xtarray = txArray(xt,tgtAng) # transmitting antenna array
yp = envout(xtarray,radar_pos,target_pos,radar_vel,target_vel) # direct distribution channel
yr = target(yp) # purpose
ye = envin(yr,target_pos,radar_pos,target_vel,radar_vel) # distribution environment
yt = rxArray(ye,tgtAng) # receiving antenna array
yt = rxPreamp(yt);  # The receiver

Let's calculate the direction and implement the ray formation using BeamscanEstimator2D and SubbandPhaseShiftBeamformer

estimator = EngeePhased.BeamscanEstimator2D(
    SensorArray=antenna, # AR geometry
    DOAOutputPort=true, # Turning on the output of the corners
    OperatingFrequency=fc, # carrier frequency
    NumSignals=1, # number of signals
    AzimuthScanAngles=-40:40,  # azimuth scanning range
    ElevationScanAngles=-60:60 # scanning range by location angle
)

doa = estimator(yt)[2] # calculation of the bearing estimate using CO BeamscanEstimator2D

beamformer = EngeePhased.SubbandPhaseShiftBeamformer(
    SensorArray=antenna, # AR geometry
    OperatingFrequency=fc, # carrier frequency
    DirectionSource="Input port", # the method of setting viewing angles
    SampleRate=fs, # sampling rate
    WeightsOutputPort=true # weight factors
)

ybf = beamformer(yt,doa)[1]; # applications of the CO-generation algorithm SubbandPhaseShiftBeamformer

Let's perform consistent filtering. To do this, we calculate the filter coefficients by means of the function getMatchedFilter

mfcoeff1 = getMatchedFilter(wav) # Coefficients of the matched filter
mf1 = EngeePhased.MatchedFilter(Coefficients=mfcoeff1);

Let's apply consistent filtering to the generated beam and visualize the result of the radar operation.:

y1 = mf1(ybf) # applying a consistent filter

nSamples = wav.SampleRate/wav.PRF # number of counts per pulse period
t = ((0:nSamples-1).-(length(mfcoeff1)-1))./fs # time delay grid
r = t*c/2 # conversion from delay to range

# The result of the radar operation
plot(r/1000,abs.(y1),xlabel="Range, km",ylabel="The amplitude",lab="",xticks=-1e6:5:1e6,title="SF output")

2. Adding an interference component

Now let's simulate the interference signal by simulating a useful signal and a time delay.:

 # Simulating an interfering signal by duplicating the signal generator
jwav = deepcopy(wav)

xj = jwav(0) # formation of the interference component
Npad = ceil(Int64,3500/(c/fs)) # number of delay elements (1 count - 150 meters)
xj = circshift(xj,Npad) # signal shift to propagation
txjam = EngeePhased.Transmitter(
    Gain=10, # amplification of interference
    PeakPower=5 # Interference power
)

# Calculation of reflected signal + interference
xj = txjam(xj) # interference radiation
ye = envin(yr.+xj,target_pos,radar_pos,target_vel,radar_vel) # direct distribution channel
yt = rxArray(ye,tgtAng) # AR reception
yt = rxPreamp(yt) # The preamp
ybfj = beamformer(yt,doa)[1] # the beam shaper
y1j = mf1(ybfj); # filtering the total signal

Let's see how the output graph of the model has been transformed.:

plot(r/1000,abs.(y1j),xlabel="Range, km",ylabel="The amplitude",
    lab="",xticks=-1e6:5:1e6,title="SF output when exposed to interference")

We can see that after adding interference, 2 correlation peaks formed at the output of the matched filter, with 1 peak from interference prevailing and obscuring the true position of the target.

3. Tuning the central frequency

To reduce the influence of the interference component, you can use the frequency tuning approach, which consists in changing the center frequency of the signal to the value of several signal bands.:

deltaf = 2.5 * bw # center frequency offset by 250 kHz
xh = deepcopy(wav(deltaf)); # a signal with an offset center frequency

Let's construct a spectrogram of the sum of 2 signals - an unbiased and an offset FM signal.:

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

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

Let's repeat the simulation of the radar operation scenario for an offset signal:

xth = tx(xh) # The transmitter
xtharray = txArray(xth, tgtAng) # transmitting headlight
yph = envout(xtharray,radar_pos,target_pos,radar_vel,target_vel) # direct distribution channel
yrh = target(yph) # purpose

yeh = envin(yrh.+xj,target_pos,radar_pos,target_vel,radar_vel) # reverse distribution channel
yth = rxArray(yeh,tgtAng) # receiving lights
yth = rxPreamp(yth) # The preamp
ybfh = beamformer(yth,doa)[1]; # beam formation

To isolate the useful component and remove the interfering component, a band-pass filter is required, tuned to an offset frequency with a ZS band.

Let's set the coefficients of the low-pass Butterworth filter of the 9th order and shift the center band by 250 kHz:

include("$(@__DIR__)/data_filter_coeff.jl") # reading coefficients from a jl file
# center frequency offset by 250 kHz
bf2 = buttercoef.*exp.(1im*2*pi*deltaf*Vector(0:length(buttercoef)-1)/fs);

We implement filtering in the SW band:

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

# bandpass filtering and matched filtering
yb2 = mf2(DSP.filt(bf2,1,ybfh))

plot(r./1000,abs.(yb2),lab="")
xlabel!("Range, km")
ylabel!("The amplitude")
title!("SF output after perestroika")

From the graph of the signal at the SF output, it can be seen that the correlation peak of the interference component has been successfully suppressed.
Thus, the use of center frequency tuning allowed

Conclusion

In the example, we considered the simulation of a frequency-tunable radar using the system objects of the EngeePhased library. This approach made it possible to reduce the impact of the interference component and successfully isolate the useful signal.