using Main: engee, plot, plot!, histogram, hline!, savefig, Animation, frame, gif
Kd = 0.5                 # 鉴相系数，1/rad
τ = 0.02                 # 每个低通滤波器的时间常数，用
Kp, Ki = 15.0, 100.0     # PI控制器系数
ω0 = 62.8319             # 载波频率，rad/s
Am, Ω = 1.5, π           # 频偏幅度和调制频率，rad/s

pll_F(s) = 1 / (τ * s + 1)^2         # 连续两次低通
pll_C(s) = Kp + Ki / s               # PI控制器
pll_L(s) = Kd * pll_F(s) * pll_C(s) / s      # 开环：GONG将频率集成到相位
pll_H(s) = pll_L(s) / (1 + pll_L(s))         # 从输入频率到控制信号的偏差

β = Am / Ω
ripple = Kp * Kd * abs(pll_F(2im * ω0))   # 在滤波器和Kp之后，分量的幅度为2ω0
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)

# 运行模拟；到工作区信号将出现在工作区中
engee.run(model)

# To Workspace将信号存储为表格：时间，值。
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)
)

# 段t中信号x(t)中谐波频率f(Hz)的幅值≥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("解调,0.5赫兹: ", round(a_dem, digits = 3), " рад/с; отношение ", round(a_dem / a_mod, digits = 3),
        " (теория |H| = ", round(abs(Hm), digits = 3), ")")
println("解调,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))