Преобразование сигнала с помощью PCM и PWM
В этом примере сравним импульсно-кодовую модуляцию (PCM) и широтно-импульсную модуляцию (PWM) в модели Engee. Измерим ошибки квантования и кодирования, оценим задержку фильтров и проверим зависимость средней величины PWM от входного сигнала.
Дискретизация и квантование
При дискретизации берём отсчёты . Для восстановления ограниченного по полосе сигнала частота дискретизации должна превышать удвоенную верхнюю частоту спектра.
Квантователь округляет отсчёты с шагом :
Для равномерно распределённой ошибки её среднеквадратическое значение оценивается как .
Кодирование и декодирование
Для битов и диапазона шаг равен . Знаковый код принимает значения от до , а декодер умножает код на . Верхний восстановимый уровень — .
Отношение мощности сигнала к мощности ошибки квантования:
Широтно-импульсная модуляция
PWM-сигнал переключается между 0 и 1. Его среднее значение за период равно коэффициенту заполнения:
При дискретном шаге заполнение изменяется с шагом .
Описание модели
Модель pcm_pwm_analysis показана на рисунке ниже.
.png)
Модель содержит ветви 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)
Максимальный знаковый код — 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), " с")
Один фильтр даёт задержку около 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
Сигналы записаны на разных временных сетках. Приведём их к общей сетке для сравнения.
Функция 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)))
Рассчитанная ошибка 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))
Ошибка не выходит за ±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)
.png)
Гистограмма не обязана быть равномерной: за 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))
Код 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)
.png)
До уровня 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)
.png)
Ошибка квантования ограничена ±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 восстановленного сигнала.
Найдём задержку, при которой среднеквадратическая разность сигналов минимальна. Проверим диапазон 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))
Сравним найденную задержку с расчётом по фазе фильтра и наложим сигналы с учётом сдвига.
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)
.png)
После компенсации задержки сигналы почти совпадают. Основное различие до сдвига связано с фильтрацией.
Посмотрим, какие коды использует модель.
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)
.png)
Для этого входного сигнала используются коды от 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.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)
.png)
При росте входного сигнала ширина импульсов увеличивается.
Сравним среднее 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))
.png)
Верхний график показывает слежение среднего значения за входом, нижний — отклонение.
Выведем измеренные ошибки, 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")
Вывод
Модель показывает влияние квантования и ограничения диапазона PCM на точность восстановления сигнала. Ошибка квантования не превышает половины шага, задержка фильтра составляет около 0,17 с, а среднее значение PWM следует за входным сигналом.