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

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

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 # выходной порт коэффициента передачи
);

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!("Дальность, м")

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

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

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)), м")

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")

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

print("Оценка скорости: $(round(tgtspeed;sigdigits=4)) м/с")