Определение ЛАФЧХ линейной системы с помощью PRBS

function interpfreqresp(
    w0::Vector{Float64},
    h0::Vector{ComplexF64},
    w::Vector{Float64})

    h = zeros(ComplexF64, length(w))
    ip0 = findall(w0 .> 0.)
    
    if isempty(ip0)
        w0p1 = 0.
    else
        w0p1 = w0[first(ip0)]
        ip = findall(w .> w0p1)
        h[ip] = _locallinloginterp(log.(w0[ip0]), h0[ip0], log.(w[ip]))
    end

    in0 = findall(w0 .< 0.)
    if isempty(in0)
        w0n1 = 0.
    else
        in0 = in0[end:-1:1]
        w0n1 = w0[first(in0)]
        in = findall(w .< w0n1)
        h[in] = _locallinloginterp(log.(-w0[in0]), h0[in0], log.(-w[in]))
    end

    ix = findall((w .>= w0n1) .& (w .<= w0p1))

    if !isempty(ix)
        ix0 = findall((w0 .>= w0n1) .& (w0 .<= w0p1))
        h[ix] = _locallinloginterp(w0[ix0], h0[ix0], w[ix])
    end

    return h
end

function _locallinloginterp(
    w0::Vector{Float64},
    h0::Vector{ComplexF64},
    w::Vector{Float64})

    h = zeros(ComplexF64, length(w))
    iw = utinterp1(w0, Float64.(collect(1:length(w0))), w)
    iwf = rem.(iw, 1)
    jexact = findall(iwf .== 0.)
    h[jexact] = h0[iw[jexact]]
    jinterp = findall(iwf .> 0.)

    if !isempty(jinterp)
        logh0 = log.(h0)
        logh0 = complex.(real.(logh0), DSP.unwrap(imag.(logh0)))

        for ct in eachindex(jinterp)
            j = jinterp[ct]
            iwj = Int64(floor(iw[j]))
            t = iwf[j]
            aux = (1 - t) * logh0[iwj] + t * logh0[iwj+1]
            h[j] = isfinite(aux) ? exp(aux) : _localrepairsingularity(aux, h0[iwj], h0[iwj+1], t)
        end
    end

    return h
end

function utinterp1(
    x::Vector{Float64},
    y::Vector{Float64},
    xi::Vector{Float64})

    n = length(y)
    nyi = size(xi, 1)
    h = diff(x)

    if any(h .< 0.)
        p = sortperm(x)
        y = y[p]
        x = sort(x)
        h = diff(x)
    end

    duplicatexidx = findall((h .== 0.) .| isnan.(h))

    if !isempty(duplicatexidx)
        deleteat!(x, duplicatexidx)
        deleteat!(y, duplicatexidx)
        deleteat!(h, duplicatexidx)
        n = length(y)
    end

    yi = fill(NaN, nyi)
    inboundsidx = findall((xi .>= first(x)) .& (xi .<= x[n]))
    xi = xi[inboundsidx]

    k = map(y -> findbin(y, x), xi)
    kidx = findall(k .== n)
    isempty(kidx) || (k[idx] .= n - 1)

    for ct in eachindex(xi)
        bin = k[ct]

        if (bin == 1) && (x[bin] == -Inf)
            yi[inboundsidx[ct]] = y[bin + Int64(!isinf(xi[ct]))]
        elseif (bin == n-1) && (x[bin+1] == Inf)
            yi[inboundsidx[ct]] = y[bin + Int64(isinf(xi[ct]))]
        else
            s = (xi[ct] - x[bin]) / h[bin]

            if (s == 0.) || (s == 1.)
                yi[inboundsidx[ct]] = y[bin+Int64(s)]
            else
                yi[inboundsidx[ct]] = (1. - s) * y[bin] + s * y[bin+1]
            end
        end
        end

    return yi
end

function _localrepairsingularity(
    logh::Union{ComplexF64, Float64},
    h1::Vector{ComplexF64},
    h2::Vector{ComplexF64},
    t::Float64)

    h = zeros(size(logh))

    for ct in eachindex(h)
        if isfinite(logh[ct])
            h[ct] = exp(logh[ct])
        else
            y1 =  h1[ct]
            y2 = h2[ct]
            MAP = Int64([y1 == 0.; isinf(y1); ((y1 != 0.) & isfinite(y1))]) *
                Int64([y2 == 0.; isinf(y2); ((y2 != 0.) & isfinite(y2))])'

            if MAP[1,1] == 1 || MAP[2,2] == 1
                h[ct] = y1
            elseif MAP[1,2] == 1 || MAP[2,1] == 1
                h[ct] = NaN
            elseif MAP[1,3] == 1 || MAP[3,1] == 1 
                h[ct] = (1 - t) * y1 + t * y2
            else
                h[ct] = 1. / ((1 - t) / y1 + t / y2)
            end
        end
    end

    return h
end

function findbin(x, edges)
    nbins = length(edges)

    if !isnan(x) && !isempty(edges)
        k = k0 = 0
        k1 = nbins

        if x >= edges[1] && x < edges[nbins]
            k = (k0 + k1) / 2

            while k0 < k1-1
                if x >= edges[Int64(ceil(k))]
                    k0 = k
                else
                    k1 = k
                end

                k = (k0 + k1) / 2
            end

            k = k0
        end

        (x == edges[nbins]) && (k = nbins)

        return Int64(ceil(k))
    else
        return 0
    end
end