function run_model( name_model)
    Path = (@__DIR__) * "/" * name_model * ".engee"
    if name_model in [m.name for m in engee.get_all_models()] 
        model = engee.open( name_model ) 
        model_output = engee.run( model, verbose=true );
    else
        model = engee.load( Path, force=true ) 
        model_output = engee.run( model, verbose=true );
        engee.close( name_model, force=true );
    end
    sleep(0.1)
    return model_output
end
run_model("PUCCH_model") 

using EngeeDSP
using Statistics, LinearAlgebra
using FFTW  
using StatsBase 

function unwrap(phase::Vector{Float64})
    unwrapped = copy(phase)
    for i in 2:length(phase)
        diff = phase[i] - phase[i-1]
        if diff > π
            unwrapped[i] -= 2π
        elseif diff < -π
            unwrapped[i] += 2π
        end
    end
    return unwrapped
end

function compute_papr(signal::AbstractArray{<:Complex})
    power = abs2.(signal)
    papr_db = 10 * log10(maximum(power) / mean(power))
    return round(papr_db, digits=2)
end

function compute_spectrum(signal::AbstractArray{<:Complex}, fs=1.0)
    N = length(signal)
    fft_result = fft(signal)
    freqs = fftfreq(N, fs)
    pwr = abs2.(fft_result) / N
    return freqs, pwr
end

function analyze_constellation(signal::AbstractArray{<:Complex}, title)
    p = scatter(real.(signal), imag.(signal),
                title=title, xlabel="In-phase (I)", ylabel="Quadrature (Q)",
                legend=false, markersize=4, markerstrokewidth=0, alpha=0.7,
                aspect_ratio=:equal)
    return p
end

function compute_signal_stats(signal::AbstractArray{<:Complex}, name)
    println("\n" * "═"^50)
    println("statistics: $name")
    println("═"^50)
    println("   Data type:          $(eltype(signal))")
    println("   Size:              $(size(signal))")
    println("   PAPR: $(compute_papr(signal)) dB")
    println("   Energy:             $(round(sum(abs2.(signal)), digits=2))")
    println("   Average amplitude: $(round(mean(abs.(signal)), digits=4))")
    println("   Average phase: $(round(mean(angle.(signal)), digits=4)) glad")
end

y = collect(simout["PUCCH_model/PUCCH1/# Parse and validate inputs.y"]).value[end]
oSeq1 = collect(simout["PUCCH_model/PUCCH1/Spreading Format-1.1"]).value[end]
final = collect(simout["PUCCH_model/PUCCH1/# Get the PUCCH format 1 sequence.Seq"]).value[end]

compute_signal_stats(y, "Signal y (after nrPUCCH1)")
compute_signal_stats(oSeq1, "Orthogonal sequence 1")
compute_signal_stats(final, "The final sequence")

p1 = analyze_constellation(y, "PUCCH1: The constellation of the 'y' signal")
p2 = analyze_constellation(final, "PUCCH1: The Final Constellation")
p3 = plot(abs.(final), title="PUCCH1: The amplitude envelope",
            xlabel="The counts", ylabel="The amplitude", legend=false, grid=true)
display(plot(p1, p2, p3, layout=(1,3), size=(1200, 400)))
savefig("pucch1_analysis.png")

prbs_bits = collect(simout["PUCCH_model/PUCCH2/PUCCH PRBS-1.1"]).value[end]
modulated = collect(simout["PUCCH_model/PUCCH2/Symbol Modulate-1.1"]).value[end]
println("Scrambled Bits (PRBS):")
println("   Number of bits: $(length(prbs_bits))")
println("   Balance (1/0): $(sum(prbs_bits))/$(sum(.!prbs_bits))")

compute_signal_stats(modulated, "QPSK modulated characters")
unique_symbols = unique(round.(modulated, digits=6))
println("\Punic points of the QPSK constellation:")
for sym in unique_symbols
    println("   I=$(round(real(sym), digits=4)), Q=$(round(imag(sym), digits=4))")
end

phase = angle.(modulated)
instant_freq = diff(unwrap(phase)) 
p1 = histogram(instant_freq, bins=30, title="PUCCH2: Instantaneous frequency",
                xlabel="Frequency", ylabel="Frequency of occurrence", legend=false)
p2 = analyze_constellation(modulated, "PUCCH2: The QPSK Constellation")

display(plot(p1, p2, layout=(1,2), size=(1200, 400)))
savefig("pucch2_analysis.png")

converted = collect(simout["PUCCH_model/PUCCH3/Convert in Int.1"]).value[end]
final = collect(simout["PUCCH_model/PUCCH3/Transform Precoding-1.1"]).value[end]
println("\Converted integer data:")
println("   Range: $(minimum(converted)) - $(maximum(converted))")
println("   Average: $(round(mean(converted), digits=2))")
compute_signal_stats(final, "The signal after Transform Precoding")

N = length(final)
fs = 1.0 
fft_result = fft(final)
freqs = fftfreq(N, fs)
pwr = abs2.(fft_result) / N
p1 = plot(freqs, pwr, title="PUCCH3: Power spectral Density",
            xlabel="Frequency", ylabel="Power (dB)", legend=false, grid=true,
            xlim=(-0.5, 0.5))
p2 = analyze_constellation(final, "PUCCH3: Constellation after DFT")

max_pwr = maximum(pwr)
half_pwr = max_pwr / 2 
indices = findall(pwr .> half_pwr)
if !isempty(indices)
    min_freq = minimum(freqs[indices])
    max_freq = maximum(freqs[indices])
    bandwidth = max_freq - min_freq
    println("\Lane width mark-up:")
    println("   Bandwidth at 3 dB level: $(round(bandwidth, digits=4))")
else
    println("We have not been able to estimate the bandwidth at the 3 dB level")
end
display(plot(p1, p2, layout=(1,2), size=(900, 400)))
savefig("pucch3_analysis.png")

prbs = collect(simout["PUCCH_model/PUCCH4/PUCCH PRBS.1"]).value[end]
modulated = collect(simout["PUCCH_model/PUCCH4/Symbol Modulate.1"]).value[end]
after_spreading = collect(simout["PUCCH_model/PUCCH4/Spreading.Out1"]).value[end]
final = collect(simout["PUCCH_model/PUCCH4/Transform Precoding.1"]).value[end]
compute_signal_stats(modulated, "Before Block-wise Spreading")
compute_signal_stats(after_spreading, "After Spreading")
compute_signal_stats(final, "After Transform Precoding")
gain_factor = mean(abs2.(after_spreading)) / mean(abs2.(modulated))
println("Efficiency of Block-wise Spreading:")
println("   Power gain: $(round(gain_factor, digits=2))")

plots = []
push!(plots, analyze_constellation(modulated, "Before Spreading"))
push!(plots, analyze_constellation(after_spreading, "After Spreading"))
push!(plots, analyze_constellation(final, "The final signal"))
p4 = plot(abs.(modulated), label="Before Spreading", linewidth=1.5, alpha=0.8)
plot!(p4, abs.(after_spreading), label="After Spreading", linewidth=1, alpha=0.7)
plot!(p4, abs.(final), label="After Precoding", linewidth=1, alpha=0.7,
        title="PUCCH4: Comparing amplitudes", xlabel="The counts", ylabel="The amplitude", grid=true)
push!(plots, p4)
display(plot(plots..., layout=(2,2), size=(1000, 800)))
savefig("pucch4_analysis.png")