Инвертируйте рекурсию Левинсона-Дербина
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.
