Преобразование сигнала с помощью PCM и PWM

В этом примере сравним импульсно-кодовую модуляцию (PCM) и широтно-импульсную модуляцию (PWM) в модели Engee. Измерим ошибки квантования и кодирования, оценим задержку фильтров и проверим зависимость средней величины PWM от входного сигнала.

Дискретизация и квантование

При дискретизации берём отсчёты . Для восстановления ограниченного по полосе сигнала частота дискретизации должна превышать удвоенную верхнюю частоту спектра.

Квантователь округляет отсчёты с шагом :

Для равномерно распределённой ошибки её среднеквадратическое значение оценивается как .

Кодирование и декодирование

Для битов и диапазона шаг равен . Знаковый код принимает значения от до , а декодер умножает код на . Верхний восстановимый уровень — .

Отношение мощности сигнала к мощности ошибки квантования:

Широтно-импульсная модуляция

PWM-сигнал переключается между 0 и 1. Его среднее значение за период равно коэффициенту заполнения:

При дискретном шаге заполнение изменяется с шагом .

Описание модели

Модель pcm_pwm_analysis показана на рисунке ниже.

Модель содержит ветви PCM и PWM. В PCM последовательно выполняются фильтрация, дискретизация, квантование, кодирование и декодирование.

Сигнал ошибки	Что сравнивается
quant_error	Вход и выход квантователя
pcm_error	Квантованный и декодированный сигналы
total_error	Вход квантователя и декодированный сигнал

Задержку фильтров оцениваем отдельно от ошибок кодирования.

Подготовим функции анализа.

using Main: engee, plot, plot!, histogram, hline!, savefig, Animation, frame, gif
using Statistics

Зададим параметры модели и рассчитаем шаг квантования и ожидаемый SQNR.

# Параметры блоков модели pcm_pwm_analysis
A      = 0.5         # амплитуда Sine Wave
Ω      = 1.0         # частота Sine Wave, рад/с
ωc     = 30.0        # частота среза фильтров Analog Filter Design, рад/с
N      = 8           # порядок фильтров
Ts     = 0.01        # Sample Time блока Zero-Order Hold, с
q_set  = 0.0078125   # Quantization interval блока Quantizer
V      = 1.0         # Peak блоков Uniform Encoder и Uniform Decoder
B      = 8           # Bits блоков Uniform Encoder и Uniform Decoder
T_pwm  = 0.05        # Period блока PWM, с
Ts_pwm = 0.001       # Sample Time блока PWM, с

q = 2V / 2^B
sqnr_theory = 10 * log10((A^2 / 2) / (q^2 / 12))

println("Шаг кодера q = 2V/2^B: ", q, "; совпадает с интервалом Quantizer: ", q == q_set)
println("Предел ошибки квантования q/2: ", q / 2)
println("СКЗ ошибки при равномерном распределении q/√12: ", round(q / sqrt(12), sigdigits = 4))
println("Код уровня 1.0 без насыщения: ", round(Int, 1.0 / q), "; наибольший код: ", 2^(B - 1) - 1)
println("Оценка SQNR для синусоиды амплитуды ", A, ": ", round(sqnr_theory, digits = 1), " дБ")
println("Шаг коэффициента заполнения PWM T_pwm/T: ", Ts_pwm / T_pwm)

Шаг кодера q = 2V/2^B: 0.0078125; совпадает с интервалом Quantizer: true
Предел ошибки квантования q/2: 0.00390625
СКЗ ошибки при равномерном распределении q/√12: 0.002255
Код уровня 1.0 без насыщения: 128; наибольший код: 127
Оценка SQNR для синусоиды амплитуды 0.5: 43.9 дБ
Шаг коэффициента заполнения PWM T_pwm/T: 0.02

Максимальный знаковый код — 127, поэтому уровень 1.0 не представим точно. Для используемой амплитуды сигнала ожидаемый SQNR составляет около 44 дБ.

Передаточную функцию фильтра Баттерворта получим из его полюсов:

Оценим усиление и задержку на частоте сигнала по модулю и фазе .

# Полюсы и частотная характеристика аналогового фильтра Баттерворта
poles = [ωc * exp(im * π * (2k + N - 1) / (2N)) for k in 1:N]
H(s) = prod(-p / (s - p) for p in poles)

H_sig = H(im * Ω)
τ_filter = -angle(H_sig) / Ω

println("Модуль |H(jΩ)| на частоте сигнала: ", round(abs(H_sig), digits = 6))
println("Фазовая задержка одного фильтра: ", round(τ_filter, digits = 4), " с")
println("Групповая задержка на нулевой частоте: ", round(1 / (ωc * sin(π / (2N))), digits = 4), " с")

Модуль |H(jΩ)| на частоте сигнала: 1.0
Фазовая задержка одного фильтра: 0.1709 с
Групповая задержка на нулевой частоте: 0.1709 с

Один фильтр даёт задержку около 0.17 с при усилении, близком к единице. Два последовательных фильтра примерно удваивают задержку.

Загрузим модель из папки примера и запустим расчёт.

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

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

Получим записанные сигналы и закроем модель.

x_sine  = collect(results["Sine Wave.1"])
x_zoh   = collect(results["Zero-Order Hold.1"])
x_quant = collect(results["Quantizer.1"])
x_code  = collect(results["Uniform Encoder.1"])
x_dec   = collect(results["Uniform Decoder.1"])
x_add   = collect(results["Add.1"])
x_duty  = collect(results["Analog Filter Design-1.1"])
x_pwm   = collect(results["PWM.1"])

engee.close(model_name; force = true)

signals = ["Sine Wave" => x_sine, "Zero-Order Hold" => x_zoh, "Quantizer" => x_quant,
           "Uniform Encoder" => x_code, "Uniform Decoder" => x_dec, "Add" => x_add,
           "Analog Filter Design-1" => x_duty, "PWM" => x_pwm]

for (name, s) in signals
    println(rpad(name, 24), length(s.time), " точек, t от ", first(s.time), " до ", last(s.time),
            " с, шаг ≈ ", round(median(diff(s.time)), sigdigits = 3), " с, тип ", eltype(s.value))
end

Sine Wave               20001 точек, t от 0.0 до 20.0 с, шаг ≈ 0.001 с, тип Float64
Zero-Order Hold         2001 точек, t от 0.0 до 20.0 с, шаг ≈ 0.01 с, тип Float64
Quantizer               2001 точек, t от 0.0 до 20.0 с, шаг ≈ 0.01 с, тип Float64
Uniform Encoder         2001 точек, t от 0.0 до 20.0 с, шаг ≈ 0.01 с, тип Int64
Uniform Decoder         2001 точек, t от 0.0 до 20.0 с, шаг ≈ 0.01 с, тип Float64
Add                     2001 точек, t от 0.0 до 20.0 с, шаг ≈ 0.01 с, тип Float64
Analog Filter Design-1  20001 точек, t от 0.0 до 20.0 с, шаг ≈ 0.001 с, тип Float64
PWM                     20001 точек, t от 0.0 до 20.0 с, шаг ≈ 0.001 с, тип Float64

Сигналы записаны на разных временных сетках. Приведём их к общей сетке для сравнения.

Функция hold_at выбирает последний отсчёт, предшествующий заданному моменту времени. Используем её для сопоставления ошибок скрипта с выходами блоков модели.

function hold_at(s, t)
    idx = searchsortedlast.(Ref(s.time), t .+ 1e-6)
    return s.value[clamp.(idx, 1, length(s.time))]
end

rms(x) = sqrt(mean(abs2, x))

t_k = collect(0:Ts:last(x_zoh.time))

u_k   = hold_at(x_zoh, t_k)          # вход квантователя (выход Zero-Order Hold)
y_k   = hold_at(x_quant, t_k)        # выход Quantizer
c_k   = Int.(hold_at(x_code, t_k))   # коды Uniform Encoder
d_k   = hold_at(x_dec, t_k)          # выход Uniform Decoder
add_k = hold_at(x_add, t_k)          # выход Add

quant_error = y_k .- u_k             # Quantizer − Zero-Order Hold
pcm_error   = y_k .- d_k             # Quantizer − Uniform Decoder
total_error = d_k .- u_k             # Uniform Decoder − Zero-Order Hold

println("Моментов дискретизации: ", length(t_k))
println("Наибольшее расхождение pcm_error и выхода Add: ", maximum(abs.(pcm_error .- add_k)))

Моментов дискретизации: 2001
Наибольшее расхождение pcm_error и выхода Add: 0.0

Рассчитанная ошибка PCM совпадает с сигналом соответствующего блока Add.

Ошибка квантования

Проверим границы ошибки квантования и её среднеквадратическое значение.

println("Ошибка квантования Quantizer − Zero-Order Hold")
println("  наибольший модуль: ", round(maximum(abs.(quant_error)), sigdigits = 5), " (предел q/2 = ", q / 2, ")")
println("  СКЗ: ", round(rms(quant_error), sigdigits = 5), " (оценка q/√12 = ", round(q / sqrt(12), sigdigits = 5), ")")
println("  среднее: ", round(mean(quant_error), sigdigits = 3))

Ошибка квантования Quantizer − Zero-Order Hold
  наибольший модуль: 0.003905 (предел q/2 = 0.00390625)
  СКЗ: 0.0022233 (оценка q/√12 = 0.0022553)
  среднее: -4.57e-5

Ошибка не выходит за ±q/2. Сравним её распределение с оценкой q/√12.

histogram(quant_error ./ q;
    bins = range(-0.5, 0.5; length = 21),
    normalize = :pdf,
    label = "quant_error / q",
    xlabel = "Ошибка в долях шага q",
    ylabel = "Плотность",
    title = "Распределение ошибки квантования",
    size = (900, 400))
hline!([1.0]; label = "Равномерное распределение", linestyle = :dash, linewidth = 2)

Гистограмма не обязана быть равномерной: за 20 с моделирования медленный входной сигнал проходит только несколько периодов.

Найдём ошибку кодирования и проверим, что происходит на верхнем уровне 1.0.

nz = findall(abs.(pcm_error) .> q / 2)
share = length(nz) / length(t_k)

println("Ошибка кодека Quantizer − Uniform Decoder (выход Add)")
println("  наибольший модуль: ", maximum(abs.(pcm_error)))
println("  СКЗ: ", round(rms(pcm_error), sigdigits = 5), "; q·√(доля ненулевых): ", round(q * sqrt(share), sigdigits = 5))
println("  ненулевых отсчётов: ", length(nz), " из ", length(t_k), " (", round(100share, digits = 2), " %)")
println("  значения Quantizer в этих отсчётах: ", unique(y_k[nz]))
println("  коды Uniform Encoder в этих отсчётах: ", unique(c_k[nz]))
println("  отсчётов с уровнем Quantizer 1.0: ", count(==(1.0), y_k))
println("  диапазон кодов за расчёт: от ", minimum(c_k), " до ", maximum(c_k))

Ошибка кодека Quantizer − Uniform Decoder (выход Add)
  наибольший модуль: 0.0078125
  СКЗ: 0.0015226; q·√(доля ненулевых): 0.0015226
  ненулевых отсчётов: 76 из 2001 (3.8 %)
  значения Quantizer в этих отсчётах: [1.0]
  коды Uniform Encoder в этих отсчётах: [127]
  отсчётов с уровнем Quantizer 1.0: 76
  диапазон кодов за расчёт: от 0 до 127

Код 127 восстанавливается как 1−q. Поэтому на вершине сигнала появляется ошибка величиной q.

Рассмотрим участок около вершины и отметим границы ошибки.

i0 = isempty(nz) ? argmax(y_k) : first(nz)
win = abs.(t_k .- t_k[i0]) .<= 0.3

plot(t_k[win], u_k[win];
    seriestype = :steppost, label = "Zero-Order Hold", linewidth = 2,
    xlabel = "Время, с", ylabel = "Уровень",
    title = "Отсчёты у вершины синусоиды",
    legend = :outertopright, size = (1000, 450))
plot!(t_k[win], y_k[win]; seriestype = :steppost, label = "Quantizer", linewidth = 2)
plot!(t_k[win], d_k[win]; seriestype = :steppost, label = "Uniform Decoder", linewidth = 2, linestyle = :dash)
hline!([1.0, 1 - q]; label = "Уровни 1.0 и 1 − q", color = :gray, linestyle = :dot)

До уровня 1.0 выходы квантователя и декодера совпадают. На вершине декодер ограничен уровнем 127/128.

Сравним ошибки квантования, кодирования и их суммарный результат.

plot(t_k, quant_error;
    label = "quant_error: Quantizer − Zero-Order Hold", linewidth = 1,
    xlabel = "Время, с", ylabel = "Ошибка",
    title = "Ошибки между точками модели",
    legend = :outertopright, size = (1100, 450))
plot!(t_k, pcm_error; label = "pcm_error: Quantizer − Uniform Decoder", linewidth = 1.5)
plot!(t_k, total_error; label = "total_error: Uniform Decoder − Zero-Order Hold", linewidth = 1, alpha = 0.6)
hline!([q / 2, -q / 2]; label = "±q/2", color = :gray, linestyle = :dash)

Ошибка квантования ограничена ±q/2, а ошибка PCM появляется на вершинах. Поэтому суммарная ошибка может выходить за границы ошибки одного квантователя.

Вычислим SQNR после первых двух секунд, исключив начальный переходный процесс.

steady = t_k .>= 2.0
p_signal = var(u_k[steady])

sqnr_quant = 10 * log10(p_signal / mean(abs2, quant_error[steady]))
sqnr_total = 10 * log10(p_signal / mean(abs2, total_error[steady]))

println("SQNR по ошибке квантования: ", round(sqnr_quant, digits = 1), " дБ")
println("SQNR по итоговой ошибке:    ", round(sqnr_total, digits = 1), " дБ")
println("Теоретическая оценка:        ", round(sqnr_theory, digits = 1), " дБ")

SQNR по ошибке квантования: 43.8 дБ
SQNR по итоговой ошибке:    42.9 дБ
Теоретическая оценка:        43.9 дБ

Ограничение верхнего уровня добавляет ошибку и снижает SQNR восстановленного сигнала.

Найдём задержку, при которой среднеквадратическая разность сигналов минимальна. Проверим диапазон 0–0.5 с с шагом 1 мс.

τ_grid = 0:0.001:0.5
lag_rms = [rms(d_k[steady] .- hold_at(x_sine, t_k[steady] .- τ)) for τ in τ_grid]
τ_est = τ_grid[argmin(lag_rms)]

println("Задержка по данным: ", τ_est, " с; расчёт для одного фильтра: ", round(τ_filter, digits = 4), " с")
println("СКЗ разности без сдвига: ", round(rms(d_k[steady] .- hold_at(x_sine, t_k[steady])), sigdigits = 4))
println("СКЗ разности со сдвигом: ", round(minimum(lag_rms), sigdigits = 4))

Задержка по данным: 0.171 с; расчёт для одного фильтра: 0.1709 с
СКЗ разности без сдвига: 0.06158
СКЗ разности со сдвигом: 0.002467

Сравним найденную задержку с расчётом по фазе фильтра и наложим сигналы с учётом сдвига.

plot(t_k, hold_at(x_sine, t_k);
    label = "Sine Wave", linewidth = 2,
    xlabel = "Время, с", ylabel = "Амплитуда",
    title = "Выход генератора и выход PCM-декодера",
    legend = :outertopright, size = (1100, 450))
plot!(t_k, d_k; label = "Uniform Decoder", linewidth = 2, linestyle = :dash)
plot!(t_k, hold_at(x_sine, t_k .- τ_est); label = "Sine Wave, сдвиг на τ", linewidth = 1.5, linestyle = :dot)

После компенсации задержки сигналы почти совпадают. Основное различие до сдвига связано с фильтрацией.

Посмотрим, какие коды использует модель.

plot(t_k, c_k;
    seriestype = :steppost, label = "Uniform Encoder", linewidth = 1.5,
    xlabel = "Время, с", ylabel = "Код",
    title = "Коды на выходе PCM-кодера",
    legend = :outertopright, size = (1100, 400))
hline!([2^(B - 1) - 1]; label = "Наибольший код 127", color = :gray, linestyle = :dash)

Для этого входного сигнала используются коды от 0 до 127. Верхняя граница ограничивает восстановление вершины.

Для PWM используем сетку 1 мс и период 0.05 с. Рассчитаем среднее значение за каждый период; шаг заполнения составляет 0.02.

n_per = round(Int, T_pwm / Ts_pwm)
t_all = collect(0:Ts_pwm:last(x_pwm.time) - Ts_pwm)
n_periods = length(t_all) ÷ n_per
t_p = t_all[1:n_per * n_periods]     # только целые периоды

pwm_p  = hold_at(x_pwm,  t_p)
duty_p = hold_at(x_duty, t_p)

pwm_avg  = vec(mean(reshape(pwm_p,  n_per, n_periods); dims = 1))
duty_avg = vec(mean(reshape(duty_p, n_per, n_periods); dims = 1))
t_avg    = T_pwm .* (0:n_periods - 1)
pwm_dev  = pwm_avg .- duty_avg

println("Значения выхода PWM: ", first(sort(unique(pwm_p)), 5))
println("Отсчётов на период: ", n_per, "; периодов: ", n_periods)
println("Среднее входа PWM за расчёт: ", round(mean(duty_p), digits = 5))
println("Среднее выхода PWM за расчёт: ", round(mean(pwm_p), digits = 5))
println("Разность средних за период: среднее ", round(mean(pwm_dev), digits = 4),
        ", от ", round(minimum(pwm_dev), digits = 4), " до ", round(maximum(pwm_dev), digits = 4))
println("Средние выхода за период кратны 1/", n_per, ": ",
        all(isapprox.(pwm_avg .* n_per, round.(pwm_avg .* n_per); atol = 1e-9)))

Значения выхода PWM: [0.0, 1.0]
Отсчётов на период: 50; периодов: 400
Среднее входа PWM за расчёт: 0.49865
Среднее выхода PWM за расчёт: 0.488
Разность средних за период: среднее -0.0106, от -0.1275 до 0.0188
Средние выхода за период кратны 1/50: true

Среднее значение PWM следует за входом. Отклонения связаны с дискретным шагом заполнения и изменением входа в течение периода.

Покажем импульсы на первых 0.5 с моделирования.

win_p = t_p .<= 0.5

plot(t_p[win_p], pwm_p[win_p];
    seriestype = :steppost, label = "Выход PWM", linewidth = 1.2,
    ylims = (-0.05, 1.05),
    xlabel = "Время, с", ylabel = "Уровень",
    title = "Выход PWM и управляющий сигнал на его входе",
    legend = :outertopright, size = (1100, 450))
plot!(t_p[win_p], duty_p[win_p]; label = "Вход PWM", linewidth = 2)

При росте входного сигнала ширина импульсов увеличивается.

Сравним среднее PWM с входом и построим их разность.

p_avg = plot(t_avg, duty_avg;
    label = "Вход PWM, среднее за период", linewidth = 2,
    ylabel = "Уровень", title = "Средние значения за период PWM",
    legend = :outertopright)
plot!(p_avg, t_avg, pwm_avg; seriestype = :steppost, label = "Выход PWM, среднее за период", linewidth = 1.5)

p_dev = plot(t_avg, pwm_dev;
    seriestype = :steppost, label = "Выход − вход", linewidth = 1.5,
    xlabel = "Время, с", ylabel = "Разность",
    legend = :outertopright)
hline!(p_dev, [1 / n_per, -1 / n_per]; label = "±ΔD", color = :gray, linestyle = :dash)

plot(p_avg, p_dev; layout = (2, 1), size = (1100, 700))

Верхний график показывает слежение среднего значения за входом, нижний — отклонение.

Выведем измеренные ошибки, SQNR и задержку.

println("Шаг квантования q:                    ", q)
println("max |quant_error| / q:                ", round(maximum(abs.(quant_error)) / q, digits = 3))
println("СКЗ quant_error / (q/√12):            ", round(rms(quant_error) / (q / sqrt(12)), digits = 3))
println("max |pcm_error| / q:                  ", round(maximum(abs.(pcm_error)) / q, digits = 3))
println("Доля отсчётов с насыщением кодера, %: ", round(100share, digits = 2))
println("Задержка фильтра: расчёт / данные, с: ", round(τ_filter, digits = 3), " / ", τ_est)
println("Среднее (выход − вход) PWM:           ", round(mean(pwm_dev), digits = 4))
@assert maximum(abs.(quant_error)) <= q/2 + 1e-8
@assert all(isfinite, d_k) && all(isfinite, pwm_dev)
string("PCM PWM passed; samples=", length(t_k), "; SQNR=", round(sqnr_quant; digits=1), " dB")

Шаг квантования q:                    0.0078125
max |quant_error| / q:                0.5
СКЗ quant_error / (q/√12):            0.986
max |pcm_error| / q:                  1.0
Доля отсчётов с насыщением кодера, %: 3.8
Задержка фильтра: расчёт / данные, с: 0.171 / 0.171
Среднее (выход − вход) PWM:           -0.0106

Вывод

Модель показывает влияние квантования и ограничения диапазона PCM на точность восстановления сигнала. Ошибка квантования не превышает половины шага, задержка фильтра составляет около 0,17 с, а среднее значение PWM следует за входным сигналом.

PCM и PWM в Engee

Преобразование сигнала с помощью PCM и PWM

Описание модели

Ошибка квантования

Вывод

Type

Short description

Related posts

Tags

Categories

Languages

Formats

Level

Status

Published

Updated

Source