Обратная рекурсия Левинсона-Дурбина
r = rlevinson(a,efinal)
[r,u] = rlevinson(a,efinal)
[r,u,k] = rlevinson(a,efinal)
[r,u,k,e] = rlevinson(a,efinal)
Обратная рекурсия Левинсона-Дурбина реализует понижающий алгоритм для решения следующей симметричной системы Теплица линейных уравнений для r, где r = [r (1 )... r ( p + 1)] и r (i)* обозначает комплексный сопряженный с r (i).
r = rlevinson(a,efinal) решает вышеописанную систему уравнений для r вектора a, где a = [1 a (2 )... a (p + 1)]. В приложениях линейного предсказания r представляет автокорреляционную последовательность входа фильтра ошибок предсказания, где r (1) является элементом нулевой задержки. Рисунок ниже показывает типовой фильтр этого типа, где H (z) является оптимальным линейным предиктором, x (n) является входным сигналом, является предсказанным сигналом, и e (n) является ошибкой предсказания.
Входной вектор a представляет полиномиальные коэффициенты этого фильтра ошибки предсказания в нисходящих степенях z.
Фильтр должен быть минимально-фазовым, чтобы сгенерировать допустимую автокорреляционную последовательность. efinal - скалярная степень ошибки предсказания, которая равна отклонению сигнала ошибки предсказания, σ2(<reservedrangesplaceholder0>).
[r,u] = rlevinson(a,efinal) возвращает верхнюю треугольную матрицу U из UDU* разложение
где
и E является диагональной матрицей с элементами, возвращаемыми в выходе e (см. ниже). Это разложение позволяет эффективно оценивать обратную матрицу автокорреляции R−1.
Выходные матричные u содержит полином предсказательного фильтра, a, от каждой итерации обратной рекурсии Левинсона-Дурбина
где ai (j) - j-й коэффициент полинома предсказательного фильтра i-го порядка (т.е. шаг, i в рекурсии). Для примера полином 5 предсказательных фильтров го порядка
a5 = u(5:-1:1,5)'
Обратите внимание, что u(p+1:-1:1,p+1)' - вектор коэффициента входа полинома a.
[r,u,k] = rlevinson(a,efinal) возвращает вектор k длины p + 1, содержащей коэффициенты отражения. Коэффициенты отражения являются сопряженными значениями в первой строке u.
k = conj(u(1,2:end))
[r,u,k,e] = rlevinson(a,efinal) возвращает вектор длины p + 1, содержащий ошибки предсказания из каждой итерации обратной рекурсии Левинсона-Дурбина: e(1) - ошибка предсказания из модели первого порядка, e(2) - ошибка предсказания из модели второго порядка и так далее.
Эти значения ошибки предсказания образуют диагональ матрицы, E в UDU* разложение R−1.
Оцените спектр двух синусоид в шуме с помощью авторегрессивной модели. Выберите лучший порядок модели из группы моделей, возвращенных обратной рекурсией Левинсона-Дурбина.
Сгенерируйте сигнал. Задайте частоту дискретизации 1 кГц и длительность сигнала 50 секунд. Синусоиды имеют частоты 50 Гц и 55 Гц. Шум имеет отклонение 0,2 ².
Fs = 1000; t = (0:50e3-1)'/Fs; x = sin(2*pi*50*t) + sin(2*pi*55*t) + 0.2*randn(50e3,1);
Оцените параметры авторегрессивной модели.
[a,e] = arcov(x,100); [r,u,k] = rlevinson(a,e);
Оцените спектральную плотность степени для порядков 1 , 5 , 25 , 50 и 100.
N = [1 5 25 50 100]; nFFT = 8096; P = zeros(nFFT,5); for idx = 1:numel(N) order = N(idx); ARtest = flipud(u(:,order)); P(:,idx) = 1./abs(fft(ARtest,nFFT)).^2; end
Постройте график оценок PSD.
F = (0:1/nFFT:1/2-1/nFFT)*Fs; plot(F, 10*log10(P(1:length(P)/2,:))) grid legend([repmat('Order = ',[5 1]) num2str(N')]) xlabel('Frequency (Hz)') ylabel('dB') xlim([35 70])
[1] Кей, Стивен М. Современная спектральная оценка: теория и применение. Englewood Cliffs, Нью-Джерси: Prentice Hall, 1988.
