FM-приёмник c ФАПЧ

using Main: engee, plot, plot!, histogram, hline!, savefig, Animation, frame, gif
Kd = 0.5                 # коэффициент фазового детектора, 1/рад
τ = 0.02                 # постоянная времени каждого ФНЧ, с
Kp, Ki = 15.0, 100.0     # коэффициенты ПИ-регулятора
ω0 = 62.8319             # несущая частота, рад/с
Am, Ω = 1.5, π           # амплитуда отклонения частоты и частота модуляции, рад/с

pll_F(s) = 1 / (τ * s + 1)^2         # два последовательных ФНЧ
pll_C(s) = Kp + Ki / s               # ПИ-регулятор
pll_L(s) = Kd * pll_F(s) * pll_C(s) / s      # разомкнутая петля: ГУН интегрирует частоту в фазу
pll_H(s) = pll_L(s) / (1 + pll_L(s))         # от отклонения частоты входа к сигналу управления

β = Am / Ω
ripple = Kp * Kd * abs(pll_F(2im * ω0))   # амплитуда составляющей 2ω0 после фильтров и Kp
Hm = pll_H(im * Ω)

println("Индекс модуляции β = ", round(β, digits = 3))
println("Ослабление составляющей 2ω0 двумя ФНЧ: ", round(abs(pll_F(2im * ω0)), digits = 3))
println("Оценка пульсаций 2ω0 в сигнале управления: ", round(ripple, digits = 2), " рад/с")
println("|H| на частоте модуляции: ", round(abs(Hm), digits = 3),
        ", фаза ", round(rad2deg(angle(Hm)), digits = 1), "°")

using Statistics

model_name = "PLL_FM"
model_file = joinpath(@__DIR__, "PLL_FM.engee")

model = engee.load(model_file; force=true)

# Запуск моделирования; сигналы To Workspace появятся в рабочей области
engee.run(model)

# To Workspace хранит сигнал в виде таблицы: time, value.
function signal_data(signal)
    table = collect(signal)
    return table[!, 1], table[!, 2]
end

t_mod,  x_mod  = signal_data(Core.eval(Main, :(mod_sig)))
t_fm,   x_fm   = signal_data(Core.eval(Main, :(fm_in)))
t_vco,  x_vco  = signal_data(Core.eval(Main, :(pll_vco_sig)))
t_dem,  x_dem  = signal_data(Core.eval(Main, :(demod_sig)))
t_ctrl, x_ctrl = signal_data(Core.eval(Main, :(ctrl_sig)))
t_freq, x_freq = signal_data(Core.eval(Main, :(vco_freq)))

omega0 = 62.8319
freq_deviation = x_freq .- omega0

println("Число отсчётов: ", length(t_mod))
println("Длительность: ", last(t_mod), " с")
println("RMS mod_sig: ", round(sqrt(mean(x_mod .^ 2)), digits=4))
println("RMS demod_sig: ", round(sqrt(mean(x_dem .^ 2)), digits=4))

engee.close(model_name; force = true)

# 1. Модулирующий и демодулированный сигналы.
p1 = plot(
    t_mod, x_mod;
    label="Модуляция",
    linewidth=2,
    ylabel="Амплитуда",
    title="Модулирующий и демодулированный сигналы",
    legend=:topright, legendfontsize=8, xlabel="Время, с",
    grid=true
)
plot!(
    p1, t_dem, x_dem;
    label="Выход демодулятора",
    linewidth=1.5,
    linestyle=:dash
)

# 2. FM-сигнал и сигнал VCO на коротком участке.
p2 = plot(
    t_fm, x_fm;
    label="Входной FM-сигнал",
    linewidth=1.2,
    xlims=(2.0, 2.4),
    ylabel="Амплитуда",
    title="FM-сигнал и VCO",
    legend=:topright, legendfontsize=8, xlabel="Время, с",
    grid=true
)
plot!(
    p2, t_vco, x_vco;
    label="ГУН",
    linewidth=1.2,
    linestyle=:dash
)

# 3. Управляющий сигнал и отклонение частоты VCO.
p3 = plot(
    t_ctrl, x_ctrl;
    label="Сигнал управления",
    linewidth=2,
    ylabel="рад/с",
    title="Управление VCO",
    legend=:topright, legendfontsize=8, xlabel="Время, с",
    grid=true
)
plot!(
    p3, t_freq, freq_deviation;
    label="Отклонение частоты",
    linewidth=1.2,
    linestyle=:dash
)

# 4. Увеличенный фрагмент реакции PLL.
p4 = plot(
    t_mod, x_mod;
    label="Модуляция",
    linewidth=2,
    xlims=(2.0, 4.0),
    ylabel="Амплитуда",
    title="Реакция PLL",
    legend=:topright, legendfontsize=8, xlabel="Время, с",
    grid=true
)
plot!(
    p4, t_dem, x_dem;
    label="Выход демодулятора",
    linewidth=1.5,
    linestyle=:dash
)

plot(
    p1, p2, p3, p4;
    layout=(2, 2),
    size=(1400, 1050)
)

# Амплитуда гармоники частоты f (Гц) в сигнале x(t) на участке t ≥ t_start
function harmonic_amplitude(t, x, f; t_start = 2.0)
    mask = t .>= t_start
    tt, xx = t[mask], x[mask]
    c = sum(xx .* cos.(2π * f .* tt))
    s = sum(xx .* sin.(2π * f .* tt))
    return 2 * hypot(c, s) / length(tt)
end

a_mod = harmonic_amplitude(t_mod, x_mod, 0.5)
a_dem = harmonic_amplitude(t_dem, x_dem, 0.5)
a_2f0 = harmonic_amplitude(t_dem, x_dem, 20.0)

println("mod_sig, 0.5 Гц:   ", round(a_mod, digits = 3), " рад/с")
println("demod_sig, 0.5 Гц: ", round(a_dem, digits = 3), " рад/с; отношение ", round(a_dem / a_mod, digits = 3),
        " (теория |H| = ", round(abs(Hm), digits = 3), ")")
println("demod_sig, 20 Гц:  ", round(a_2f0, digits = 3), " рад/с")
@assert length(t_mod) > 100 && all(isfinite, x_dem)
string("PLL FM passed; samples=", length(t_mod), "; modulation gain=", round(a_dem/a_mod; digits=3))