Инвертируйте рекурсию Левинсона-Дербина
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 (e).
[r,u] = rlevinson(a,efinal) возвращает верхнюю треугольную матрицу U в разложение UDU*
где
и E является диагональной матрицей с элементами, возвращенными в выходе e (см. ниже). Это разложение разрешает эффективную оценку инверсии матрицы автокорреляции, R −1.
Выходная матрица u содержит полином предсказательного фильтра, a, от каждой итерации реверса рекурсия Левинсона-Дербина
где ai (j) является j th коэффициент i th полином предсказательного фильтра порядка (т.е. шаг 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, NJ: Prentice Hall, 1988.
