Методы спектрального оценивания стационарных случайных процессов, основанные на быстром преобразовании Фурье (БПФ), хорошо известны и широко применяются в инженерной практике. К их недостаткам следует отнести, в частности, высокую дисперсию (низкую точность) оценки при недостаточно длительном интервале наблюдения за процессом, что визуально обычно проявляется в сильной «изрезанности» графика спектральной плотности мощности(СПМ). Одним из альтернативных методов спектрального оценивания является авторегрессионный метод, рассмотренный на примере ниже, который в инженерной практике известен гораздо меньше. Метод во многих случаях позволяет сравнительно просто получить гораздо более качественную оценку СПМ (рис.1), а иногда и более глубокие сведения об исследуемом случайном процессе.
Рис.1 Классическая и авторегрессионная оценка СПМ «короткого» процесса
Для демонстрационных целей был синтезирован дискретно-временной сигнал (последовательность) x[i]. Сигнал смоделирован при помощи ARMA-модели (цифрового фильтра), имитирующей свойства механической системы (1) — перемещение материальной точки x(t) в «одномассовом» осцилляторе с параметрами m=1 кг, c= 100 Н/м, k=2,5 кг/с, и силовым возмущением — гауссовым «белым»(с учетом дискретизации) шумом f(t) с дисперсией 1 Н2, интервал дискретизации по времени Δt=0,12 с.
Построена модель (2). Способ построения модели уже рассматривался ранее здесь .
x[i] — 0.6388· x[i-1] + 0.7408· x[i-2] = 0.009667·f[i-1] (2)
С помощью (2) синтезирована последовательность в 50 тыс. отсчетов, для чего использован генератор нормально-распределенной случайной величины randn( ) общеизвестной программной среды.
После завершения моделирования процесса x[i], количественные параметры модели (2) предполагаются неизвестными — для исследования доступен только сам процесс и, в какой-то мере, сведения о свойствах модели в самых общих чертах.
Было проведено спектральное оценивание 50000-точечной последовательности методом Уэлча, размер сегмента принят равным 256 отсчетам, применено окно Хэмминга и 60% перекрытие сегментов. Среднеквадратичное отклонение такой оценки, исходя из того, что последовательность имеет длину около 200 неперекрывающихся сегментов, может быть примерно оценено, как ~7%.
Далее, предполагая, что в реальных условиях в эксперименте для исследования доступна гораздо менее длинная последовательность, проведены исследования только по первым 500 отсчетам этого сигнала.
Получена оценка методом Уэлча с теми же параметрами. СКО такой оценки ~70%, заметна очень сильная «изрезанность» графика (рис.2).
Рис.2 Оценивание СПМ «длинного» и «короткого» процессов классическими методом
Исходя из того, что примерный вид функции (графика) СПМ процесса нам известен (например, исходя из известной физической природы процесса — одномассовый осциллятор под белым шумом, либо из оценивания аналогичных процессов, для которых доступны более длинные реализации), принято решение об оценивании с помощью модели авторегрессии второго порядка (AR(2), или =ARMA(2,0)).
Определение порядка модели — весьма важный момент, ошибка в порядке может повлечь очень грубые ошибки в результатах оценивания. Существуют методы, здесь пока не рассматриваемые, помогающие в определении порядка модели на основании только самого анализируемого процесса.
Оценивание параметров модели поизведено с помощью известных уравнений Юла-Уолкера для авторегрессионного процесса (несущественно модифицированных с целью некоторого упрощения структуры скрипта ):
Как видно из уравнений, для определения параметров понадобятся только три первых члена авторегрессионной последовательности Rxx[0], Rxx[1], Rxx[2], которые и были оценены по исходной 500-точечной последовательности x[i] корелограммным методом, СКО такой оценки ~4,5%.
(Кстати, видно, что «минусы» перед a1, a2 2и т.д, крайне неудобны. Они и появились-то из-за преимущественно «предсказательного» использования ARMA-моделей в экономике, в более ранних «инженерных» источниках их нет. Уже сомневаюсь, что надо было здесь использовать такое понимание AR-коэффициентов.)
Корреляционная матрица в (3) на практике всегда имеет строгое диагональное преобладание | Rxx[0] | >| Rxx[i] |, в том числе по причине присутствия шумов наблюдения, вследствие чего трудностей с ее обращением (нахождением решения(3)) не возникает.
(Для пояснения вопроса о величине статистической ошибки моделирования интересно упомянуть, например, оценку Rxx[0] =2.2606e-04 м2, полученную по 500 отсчетам, в сравнении с полученными корелограммной оценкой дисперсии по 50000 отсчетам, = 2.4238e-04 м2 и оценкой по подынтегральной площади СПМ, полученной методом Уэлча по 50000 отсчетам (рис.2), = 2.4232e-04 м2)
После подстановки найденных оценок Rxx[i] имеем:
Определены следующие параметры модели a0=11325.9; a1=7090.1; a2=-8411.5; Как видно из (3), дисперсией гипотетического вхоящего белого шума здесь задались =1, определив вместо нее коэффициент усиления a0. Авторегрессионная оценка СПМ построена путем преобразования Фурье над последовательностью коэффициентов a0, a1, a2:
Рис.3 Классическая и авторегрессионная оценка СПМ «короткого» процесса
Таким же образом, по выражению, аналогичному (5), был ранее построен и «теоретический» график СПМ, только коэффициенты модели там, естественно, были взяты иные (из (2)).
Из графика видно, что AR-оценка СПМ получилась весьма близка к теоретически ожидаемой. Помимо графика, есть возможность попытаться оценить некоторые аналитические характеристики процесса и связанной с ним механической системы. В данном случае это «полюса» модели, численно характеризующие частоты «резонансных» пиков модели и связанные с ними «добротности».
Из (5) находим соотношение для поиска разрывов передаточной функции нашей модели, используя преобразование Лапласа (заменяя jω на λ=-ε+ jω):
Для полученной AR- модели таким способом вычислены λ1,2= -1.5427 ± j· 10.1514, что весьма близко к исходной модели, использованной для генерации процесса
λ1,2теор=-1.2500 ± j · 9.9216 (т.е положения резонансного пика соответственно, 1,615 Гц (в теории) и 1,579 Гц (определено)).
Рис.4 О понятии «полюсов»
Несколько замечаний и рекомендаций в заключение.
Рис.1 Классическая и авторегрессионная оценка СПМ «короткого» процесса
Для демонстрационных целей был синтезирован дискретно-временной сигнал (последовательность) x[i]. Сигнал смоделирован при помощи ARMA-модели (цифрового фильтра), имитирующей свойства механической системы (1) — перемещение материальной точки x(t) в «одномассовом» осцилляторе с параметрами m=1 кг, c= 100 Н/м, k=2,5 кг/с, и силовым возмущением — гауссовым «белым»(с учетом дискретизации) шумом f(t) с дисперсией 1 Н2, интервал дискретизации по времени Δt=0,12 с.
Построена модель (2). Способ построения модели уже рассматривался ранее здесь .
x[i] — 0.6388· x[i-1] + 0.7408· x[i-2] = 0.009667·f[i-1] (2)
С помощью (2) синтезирована последовательность в 50 тыс. отсчетов, для чего использован генератор нормально-распределенной случайной величины randn( ) общеизвестной программной среды.
После завершения моделирования процесса x[i], количественные параметры модели (2) предполагаются неизвестными — для исследования доступен только сам процесс и, в какой-то мере, сведения о свойствах модели в самых общих чертах.
Было проведено спектральное оценивание 50000-точечной последовательности методом Уэлча, размер сегмента принят равным 256 отсчетам, применено окно Хэмминга и 60% перекрытие сегментов. Среднеквадратичное отклонение такой оценки, исходя из того, что последовательность имеет длину около 200 неперекрывающихся сегментов, может быть примерно оценено, как ~7%.
Далее, предполагая, что в реальных условиях в эксперименте для исследования доступна гораздо менее длинная последовательность, проведены исследования только по первым 500 отсчетам этого сигнала.
Получена оценка методом Уэлча с теми же параметрами. СКО такой оценки ~70%, заметна очень сильная «изрезанность» графика (рис.2).
Рис.2 Оценивание СПМ «длинного» и «короткого» процессов классическими методом
Исходя из того, что примерный вид функции (графика) СПМ процесса нам известен (например, исходя из известной физической природы процесса — одномассовый осциллятор под белым шумом, либо из оценивания аналогичных процессов, для которых доступны более длинные реализации), принято решение об оценивании с помощью модели авторегрессии второго порядка (AR(2), или =ARMA(2,0)).
Определение порядка модели — весьма важный момент, ошибка в порядке может повлечь очень грубые ошибки в результатах оценивания. Существуют методы, здесь пока не рассматриваемые, помогающие в определении порядка модели на основании только самого анализируемого процесса.
Оценивание параметров модели поизведено с помощью известных уравнений Юла-Уолкера для авторегрессионного процесса (несущественно модифицированных с целью некоторого упрощения структуры скрипта ):
Как видно из уравнений, для определения параметров понадобятся только три первых члена авторегрессионной последовательности Rxx[0], Rxx[1], Rxx[2], которые и были оценены по исходной 500-точечной последовательности x[i] корелограммным методом, СКО такой оценки ~4,5%.
(Кстати, видно, что «минусы» перед a1, a2 2и т.д, крайне неудобны. Они и появились-то из-за преимущественно «предсказательного» использования ARMA-моделей в экономике, в более ранних «инженерных» источниках их нет. Уже сомневаюсь, что надо было здесь использовать такое понимание AR-коэффициентов.)
Корреляционная матрица в (3) на практике всегда имеет строгое диагональное преобладание | Rxx[0] | >| Rxx[i] |, в том числе по причине присутствия шумов наблюдения, вследствие чего трудностей с ее обращением (нахождением решения(3)) не возникает.
(Для пояснения вопроса о величине статистической ошибки моделирования интересно упомянуть, например, оценку Rxx[0] =2.2606e-04 м2, полученную по 500 отсчетам, в сравнении с полученными корелограммной оценкой дисперсии по 50000 отсчетам, = 2.4238e-04 м2 и оценкой по подынтегральной площади СПМ, полученной методом Уэлча по 50000 отсчетам (рис.2), = 2.4232e-04 м2)
После подстановки найденных оценок Rxx[i] имеем:
Определены следующие параметры модели a0=11325.9; a1=7090.1; a2=-8411.5; Как видно из (3), дисперсией гипотетического вхоящего белого шума здесь задались =1, определив вместо нее коэффициент усиления a0. Авторегрессионная оценка СПМ построена путем преобразования Фурье над последовательностью коэффициентов a0, a1, a2:
Рис.3 Классическая и авторегрессионная оценка СПМ «короткого» процесса
Таким же образом, по выражению, аналогичному (5), был ранее построен и «теоретический» график СПМ, только коэффициенты модели там, естественно, были взяты иные (из (2)).
Из графика видно, что AR-оценка СПМ получилась весьма близка к теоретически ожидаемой. Помимо графика, есть возможность попытаться оценить некоторые аналитические характеристики процесса и связанной с ним механической системы. В данном случае это «полюса» модели, численно характеризующие частоты «резонансных» пиков модели и связанные с ними «добротности».
Из (5) находим соотношение для поиска разрывов передаточной функции нашей модели, используя преобразование Лапласа (заменяя jω на λ=-ε+ jω):
Для полученной AR- модели таким способом вычислены λ1,2= -1.5427 ± j· 10.1514, что весьма близко к исходной модели, использованной для генерации процесса
λ1,2теор=-1.2500 ± j · 9.9216 (т.е положения резонансного пика соответственно, 1,615 Гц (в теории) и 1,579 Гц (определено)).
Рис.4 О понятии «полюсов»
Несколько замечаний и рекомендаций в заключение.
- «Избыточный» (слишком большой) порядок AR-модели обычно гораздо менее опасен, чем недостаточный, с точки зрения риска получения оценки СПМ с грубыми ошибками.
- Как правило, AR-моделирование позволяет довольно точно определить резонансные частоты jωk и гораздо мене точно — ширины соответствующих им «пиков» -εk
- ARMA — модель может получиться гораздо меньшего порядка (размера), чем AR-модель, к чему вроде бы следует стремиться для повышения точности модели, по мнению многих источников. Однако оценивание MA-части модели гораздо более затруднительно и может вообще включать в себя первым этапом получение AR-модели большого порядка с целью ее дальнейшего преобразования в MA-часть. В связи с этим источниками высказывается также альтернативное мнение о целесообразности применения для целей спектрального оценивания именно AR-моделей, пусть и большего порядка.
- Для очень коротких, а также для нестационарных процессов вместо матрицы оценок автокорреляционной функции в (3) обычно используют матрицу ковариаций.
- Для подробного изучения вопроса авторегрессионного спктрального оценивания можно рекомендовать С.Л. Марпл-мл. «Цифровой спектральный анализ и его приложения», М., Мир, 1990
