rlevinson

Обратная рекурсия Левинсона — Дурбина.

Библиотека

EngeeDSP

Синтаксис

Вызов функции

r,u,k,eEvolution = rlevinson(a,eFinal) — возвращает автокорреляционную последовательность r для предсказывающего фильтра с коэффициентами, указанными в аргументе a, и мощностью конечной ошибки предсказания eFinal. Также возвращает полиномиальную матрицу u рекурсивного предсказывающего фильтра, коэффициенты отражения k и ошибки предсказания eFinal из каждой итерации обратной рекурсии Левинсона — Дурбина.

Аргументы

Входные аргументы

# a — коэффициенты предсказывающего фильтра
вектор

Details

Коэффициенты предсказывающего фильтра, заданные как вектор.

Фильтр, связанный с a, должен иметь минимальную фазу для генерации допустимой автокорреляционной последовательности r.

Если a не имеет минимальной фазы, функция rlevinson может возвращать значения NaN или Inf, даже если решение существует.

Типы данных	`Float32`, `Float64`
Поддержка комплексных чисел	Да

# eFinal — мощность конечной ошибки предсказания
скаляр

Details

Мощность конечной ошибки предсказания, заданная как скаляр.

Типы данных	`Float32`, `Float64`
Поддержка комплексных чисел	Да

Выходные аргументы

# r — автокорреляционная последовательность
вектор-столбец

Details

Автокорреляционная последовательность, возвращаемая в виде вектора-столбца с тем же количеством элементов, что и в векторе a.

# u — полиномиальная матрица рекурсивного предсказывающего фильтра
матрица

Details

Полиномиальная матрица рекурсивного предсказывающего фильтра, возвращаемая в виде квадратной матрицы, количество строк и столбцов которой равно количеству элементов a.

# k — коэффициенты отражения
вектор-столбец

Details

Список коэффициентов отражения, возвращаемый в виде вектор-столбца длиной , где — количество элементов вектора a.

# eEvolution — рекурсивные ошибки предсказания
вектор-строка

Details

Список рекурсивных ошибок предсказания из каждой итерации обратной рекурсии Левинсона — Дурбина, возвращаемый в виде вектора-строки длиной , где — количество элементов вектора a.

Примеры

Оптимальный порядок авторегрессионной модели

Details

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

Сгенерируем сигнал с частотой 50000 отсчетов. Зададим частоту дискретизации 1 кГц и длительность сигнала 50 секунд. Синусоидальные волны имеют частоты 50 Гц и 55 Гц. Дисперсия шума равна (0.2)².

import EngeeDSP.Functions: rlevinson
using Random

Random.seed!(42)

Ns = 50_000
Fs = 1000
t = range(0, step=1/Fs, length=Ns)


x = sin.(2π * 50 * t) + sin.(2π * 55 * t) + 0.2 * randn(Ns)

Оценим параметры авторегрессионной модели, используя модель с порядком 100.

ar,eFinal = arcov(x,100)

Решим модель, используя обратную рекурсию Левинсона — Дурбина, чтобы оценить автокорреляционную последовательность, рекурсивные предсказывающие фильтры, коэффициенты отражения и ошибки предсказания.

r,u,k,eEvol=rlevinson(ar,eFinal)

Построим график ошибок предсказания. Добавим горизонтальную линию для мощности конечной ошибки. Ошибка уменьшается с увеличением порядка фильтра до достижения мощности конечной ошибки.

rlevinson 1

Оценим спектральную плотность мощности (СПМ) для порядков 1, 5, 25, 50 и 100. Оценим СПМ для исходных авторегрессионных параметров.

N = [1, 5, 25, 50, 100]
nFFT = 8096
psdARpred = zeros(nFFT, 5)


for idx in 1:length(N)
    order = N[idx]
    ARpred = reverse(u[:, order], dims=1)
    psdARpred[:, idx] = 1 ./ abs.(fft(ARpred, nFFT)).^2
end


psdAR = 1 ./ abs.(fft(ar, nFFT)).^2

Построим график оценок СПМ. Оценки СПМ для предсказывающих фильтров постепенно приближаются к исходным параметрам авторегрессии по мере увеличения порядка фильтра.

F = (0:1/nFFT:1/2-1/nFFT) * Fs
pow2db(x) = 10 * log10.(x)

plot(F, pow2db.(psdARpred[1:end÷2, :]),
     xlabel = "Frequency (Hz)",
     ylabel = "Power (dB/Hz)",
     grid = true,
     legend = :best,
     xlims = (35, 70))


plot!(F, pow2db.(psdAR[1:end÷2]),
      linestyle = :dash,
      linecolor = :black,
      label = "Original AR")

order_labels = ["Order = $n" for n in N]

push!(order_labels, "Original AR")


plot!(serieslegend = order_labels)

display(current())

rlevinson 2

Алгоритмы

Обратная рекурсия Левинсона — Дурбина реализует алгоритм пошагового решения симметричной системы линейных уравнений Тёплица.

где , , а символ обозначает комплексное сопряжение.

В приложениях линейного предсказания:

Вектор представляет собой автокорреляционную последовательность входных данных для предсказывающего фильтра, где — элемент с нулевой задержкой. Матрица Тёплица — это автокорреляционная матрица:
Вектор представляет собой коэффициенты полинома предсказывающего фильтра в порядке убывания степеней .
Скалярная величина обозначает порядок предсказывающего фильтра.

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

rlevinson ru

Функция rlevinson решает эту систему уравнений для вычисления автокорреляционной последовательности r, имея вектор a с коэффициентами предсказывающего фильтра. Фильтр должен быть минимально-фазовым, чтобы сгенерировать допустимую автокорреляционную последовательность. Скаляр eFinal представляет собой мощность конечной ошибки предсказания и стремится быть целевой ошибкой оценки рекурсии Левинсона — Дурбина, которую реализует функция rlevinson.

Для получения выходных данных u, k и eEvolution функция rlevinson также выполняет разложение обратной матрицы автокорреляции . Это разложение разрешает матрицы и и позволяет эффективно вычислить обратную матрицу автокорреляции :

Выходной аргумент u возвращает матрицу . Эта матрица содержит рекурсивные полиномы предсказывающего фильтра из каждой итерации обратной рекурсии Левинсона — Дурбина:

Каждый вектор содержит рекурсивные коэффициенты полинома, где — -й коэффициент полинома предсказывающего фильтра -го порядка (т.е., -й шаг рекурсии). Например, коэффициенты полинома предсказывающего фильтра -го порядка можно получить как , где a5 = u[5:−1:1, 5], и этот вектор представляет собой полином . Следовательно, -й столбец содержит вектор коэффициентов входного полинома a, который можно получить, используя (u[p+1:−1:1, p+1])'.

Выходной аргумент k возвращает вектор-столбец с коэффициентами отражения, которые являются сопряженными к недиагональным элементам в первой строке . Другими словами, k возвращается как (u[1, 2:end])':

Выходной аргумент eEvolution возвращает вектор-строку с рекурсивными ошибками предсказания, полученными из диагонали матрицы , за исключением первого элемента:

Каждый элемент перечисляет рекурсивные предсказания полинома предсказывающего фильтра -го порядка (т.е., -й шаг в рекурсии). Например, ошибку предсказания полинома предсказывающего фильтра -го порядка можно получить как e5 = eEvolution(4), и это ошибка предсказания при использовании полинома в качестве предсказывающего фильтра.

Первый элемент на диагонали , , совпадает с первым элементом автокорреляционной последовательности , .
-й элемент на диагонали , , — также последний элемент в eEvolution — совпадает с мощностью конечной ошибки предсказания eFinal.

Литература

Kay, Steven M. Modern Spectral Estimation: Theory and Application. Englewood Cliffs, NJ: Prentice-Hall, 1988.