Text
                    ТИИЭР, т. 69, № 11, ноябрь 1981
5
СТАТЬИ
543.42
Современные методы спектрального анализа: Обзор
С. М. КЕЙ, член ИИЭР; С. Л. МАРПЛ мл., член ИИЭР
Spectrum Analysis—A Modern Perspective
STEVEN M. KAY, member, ieee, and STANLEY LAWRENCE MARPLE, JR., member, ieee
Дан систематический обзор многих новых методов спект-
рального анализа временных рядов, которые были разработаны
в течение двух последних десятилетий. В качестве общей основы,
способствующей пониманию различий между использованными
подходами к спектральному анализу, рассматриваются модели
временных рядов, лежащие в основе каждого метода. Описанные
в статье методы включают классический метод, основанный на
применении периодограмм, метод Блэкмана — Тьюки, метод ав-
торегрессии (максимальной энтропии), метод скользящего сред-
него, метод максимального правдоподобия, метод авторегрессии —
скользящего среднего, метод Прони и метод Писаренко. Приведена
сводная таблица с кратким описанием всех рассмотренных мето-
дов, в которой указаны все основные работы и уравнения,
используемые для вычисления спектральных оценок.
I. ВВЕДЕНИЕ
Оценивание спектральной плотности мощности
(СПМ), или просто спектра, дискретизованных
детерминированных и случайных процессов
обычно выполняется с помощью процедур, использу-
ющих быстрое преобразование Фурье (БПФ). Такой
подход к спектральному анализу эффективен в вы-
числительном отношении и обеспечивает получение
приемлемых результатов для большого класса сиг-
нальных процессов. Однако несмотря на указанные
достоинства подходу, основанному на вычислении
БПФ, присущ ряд принципиальных ограничений.
Наиболее важное из них — это ограничение частот-
ного разрешения, т. е. способности различать спект-
ральные линии двух и более сигналов. Частотное
разрешение в герцах примерно равно величине,
обратной временному интервалу в секундах, на кото-
ром можно получить отсчеты сигнала. Второе огра-
ничение обусловлено неявной весовой обработкой
данных при вычислении БПФ. Взвешивание прояв-
ляется в виде «утечки» в частотной области, т. е.
энергия главного лепестка спектральной линии «уте-
кает» в боковые лепестки, что приводит к наложению
и искажению спектральных линий других присутст-
вующих сигналов. В действительности спектральные
линии слабых сигналов могут маскироваться боко-
выми лепестками спектральных линий более сильных
сигналов. Правильный выбор функции окна, значения
которой спадают на краях, позволяет ослабить утечку
в боковые лепестки, однако лишь за счет снижения
разрешающей способности.
Получена 8 сентября 1980 г., в исправленном виде — 2 июня
1981 г.
Ориг., с. 1380—1419. Статья иаписаиа на основании предвари-
тельной заявки авторов, одобренной редакцией.
S. М. Kay is with the Department of Electrical Engineering, Univer-
sity of Rhode Island, Kingston, RI02881.
S. L. Marple, Jr. is with the Analytic Sciences Corporation, McLean
Operation, 8301 Greensboro Drive, Suite 1200, McLean, VA 22102.
Два указанных ограничения подходов на основе
БПФ особенно сильно проявляются при анализе
коротких записей данных. Однако именно с такими
записями чаще всего и приходится иметь дело на
практике, поскольку многие измеряемые процессы
обладают малой длительностью или же медленно
изменяющимися во времени спектрами, которые можно
считать постоянными только на коротких участках
записей данных. В радиолокации, например, по каж-
дому принятому радиолокационному импульсу можно
получить лишь несколько отсчетов. В гидролокации
движение цели приводит к изменению спектра при-
нимаемого сигнала во времени из-за эффектов Доплера.
В течение двух последних десятилетий предложено
большое число самых различных процедур спектраль-
ного оценивания, которые были разработаны для
того, чтобы так или иначе ослабить ограничения,
присущие подходу на основе БПФ. Сравнение спект-
ральных оценок, показанных на рис. 1, позволяет
проиллюстрировать улучшение, которое можно по-
лучить при нетрадиционном подходе. Представленные
на этом рисунке спектры вычислены по 9 значениям
автокорреляционной функции 11 процесса, состоящего
из двух синусоид равной амплитуды с частотами 3
и 4 Гц и белого шума. На рис. 1а показана спектраль-
ная оценка, полученная на основе 9 известных зна-
чений ДхХ(0), . . ., /?жж(8). Спектр в данном случае
представляет кривую, состоящую из 512 точек, зна-
чения которых были вычислены с помощью 512-то-
чечного БПФ по 9 значениям автокорреляционной
функции, причем отсчеты данных были дополнены
503 нулями. Этот спектр, часто называемый оценкой
СПМ Блэкмана — Тьюки (БТ), характеризуется бо-
ковыми лепестками, некоторые из которых имеют
отрицательную величину, и невозможностью разли-
чить два синусоидальных сигнала.
На рис. 1b показан спектр того же процесса,
полученный методом авторегрессии (АР). Улучшение
разрешения по сравнению со спектром, .доказанным
на рис. 1а, обусловило популярность данного метода
спектральной оценки. И хотя спектральная АР-
оценка первоначально была разработана для обра-
ботки геофизических данных, где она определяется
с помощью метода максимальной энтропии (ММЭ)
Автокорреляционная функция #«(£) случайного стационар-
ного в широком смысле дискретного процесса хп при задержке к
определяется в статье как математическое ожидание произведения
*n+k*n> йли Rxx(k) —E[xn4.gXnl> где предполагается,что процесс
х„ имеет нулевое среднее значение. Здесь звездочкой * обозначена
комплексно-сопряженная величина, так как в общем случае рас-
сматриваются комплекснозначные пропаччл. а Е ( ) — оператор
математическрПЕЗЖйДадйяГ-„.-г,-• »	1

6 ТИИЭР, т. 69, № 11, ноябрь 1981 Действительные частоты («) (Ь) а. ( СПМ шума о в ю Частота, Гц (С) Рис. I. Примеры трех спектральных оценок на основе 9 из- вестных значений автокорреляционной функции про- цесса, состоящего из двух синусоид равной амплитуды и аддитивного белого шума (дисперсия шума составляет 10% мощности синусоиды): (а) оценка СПМ методом Блэкмана — Тьюки; (Ь) авторегрессивная оценка СПМ; (с) оценка СПМ методом гармонического разложения Писаренко. § и я [16, 37—39, 50, 84, 136, 138, 158, 221, 231, 246, 2471, она также стала применяться в радиолокации [75, 92, 99, 116, 125, 126, 216], в гидролокации [122, 198], в области обработки изображений [98], в радиоастро- номии [162, 264, 265], биомедицине [71, 74], океано- графии [96], экологии [88] и пеленгации [70, 128, 233]. АР-подход к спектральному анализу тесно связан с методами кодирования на основе линейного пред- сказания (КЛП), используемыми при обработке ре- чевых сигналов [80, 130, 143, 145]. Оценка СПМ с помощью АР-метода согласует АР-модель и анализи- руемые данные. Свои истоки АР-модели берут в области экономического прогнозирования на основе анализа временных рядов [31, 276] и статистического оценивания [189—191]. В подходе на основе ММ3 делаются различные предположения относительно значений автокорреляционной функции, однако для практических целей спектральные АР- и ММЭ-оценки в случае одномерного анализа стационарных в ши- роком смысле гауссовых процессов можно считать идентичными. Максимальное разрешение двух синусоидальных сигналов, которое обеспечивают две спектральные линии в виде дельта-функций на равномерном спект- ральном фоне, характеризующем уровень белого шума СПМ, достигается в методе гармонического разложе- ния Писаренко (ГРП); см. рис. 1с. Этот метод дает наиболее точную оценку спектра синусоид в присут- ствии шума по крайней мере в тех случаях, когда известны значения автокорреляционной функции. Примеры спектров, приведенные на рис. 1, сви- детельствуют также и о том факте, что разработка различных методов спектрального оценивания в самых различных областях привела к смешению противоре- чивых терминов и к различным точкам зрения на разработку алгоритмов. В связи с этим при напи- сании статьи были поставлены две цели: 1) установить общую терминологию и символику и 2) унифициро- вать подходы и алгоритмы, разработанные в различ- ных областях. В ряде работ рассмотрена степень улучшения спектрального разрешения и обнаружимости сигналов в тех случаях, когда для анализа дискретизованных данных применялись АР-метод и метод Писаренко [36, 206, 250, 2511. Однако, как и можно было ожи- дать, улучшение этих характеристик сильно зависит от величины отношения сигнал/шум. В действитель- ности при достаточно малых отношениях сигнал/шум спектральные оценки, получаемые с помощью со- временных методов, часто не лучше оценок, получа- емых с помощью обычных методов обработки на основе БПФ [122, 150]. Даже в тех случаях, когда более высокая точность спектральной оценки достигается с помощью альтернативной процедуры спектрального оценивания, вычислительные требования такого аль- тернативного метода могут оказаться значительно выше требований к обработке на основе БПФ. Этот факт делает некоторые современные методы спект- ральной оценки малопригодными для работы в ре- альном масштабе времени. В связи с этим третья цель, поставленная перед статьей,— рассмотреть ком- промиссные подходы к выбору различных методов. В частности, будут рассмотрены преимущества и недостатки каждого метода, проанализирована их вычислительная сложность и представлены критерии, с помощью которых можно определить соответствует ли выбранный-метод анализируемому процессу. Для того чтобы читатель составил представление о тех основах, на которых зародились современные методы спектрального анализа, полезно дать короткую историческую справку. Историю знаменитого пре- образования Фурье можно проследить более чем на 200 лет назад [34, 223]. Появление спектрального анализа, основанного на преобразовании Фурье,
СПЕКТРАЛЬНЫЙ АНАЛИЗ 7 можно проследить до Шустера, который первым упот- ребил термин «периодограмма» [218, 219]. Шустер применил ряды Фурье для аппроксимации изменений числа солнечных пятен с целью отыскать «скрытые периодичности» в измеренных данных. Следующий пионерский шаг был сделан в классической статье Норберта Винера по «обобщенному гармоническому анализу» [269], в которой устанавливалась теорети- ческая база для анализа случайных процессов с помощью подхода, основанного на преобразовании Фурье. Главный результат этой статьи состоял во введении автокорреляционной функции случайного процесса и получении выражения, связывающего ее преобразование Фурье со спектральной плотностью мощности. Независимо от Винера аналогичное со- отношение было получено советским математиком А. Я. Хинчиным [127]. Блэкман и Тьюки в своей классической работе [25], опубликованной в 1958 г., обосновали возмож- ность практической реализации автокорреляционного подхода Винера для оценки спектра мощности диск- ретных последовательностей данных. В этом методе сначала вычисляются значения автокорреляционной функции при различных задержках, затем оценки автокорреляционной функции соответствующим об- разом взвешиваются (т. е. обрабатываются с помощью спадающей к краям функции окна), после чего для получения оценки СПМ вычисляются преобразования Фурье этих взвешенных оценок автокорреляционной функции. Метод Блэкмана — Тьюки был наиболее популярным методом спектрального оценивания до тех пор, пока в 1965 г. не появился алгоритм БПФ, разработка которого обычно связывается с именами Кули и Тьюки [53]. Этот вычислительно эффективный алгоритм вновь пробудил интерес к методам спект- рального оценивания на основе периодограммы. При использовании периодограммы спектральная оценка получается в виде квадрата амплитуды выходных величин, вычисляемых посредством прямого приме- нения БПФ к последовательности входных данных (данные могут быть предварительно взвешены). В на- стоящее время именно периодограммы и применяются чаще всего для получения оценок СПМ [17, 24, 32, 105—107, 109]. Стандартные методы спектрального оценивания с использованием БПФ основаны на модели представ- ления данных с помощью рядов Фурье, т. е. анализи- руемый процесс полагается состоящим из некоторого набора гармонически связанных синусоид. В нетех- нических областях в течение многих лет, используются модели на основе других временных рядов. Юл [276] и Уокер 1258] использовали АР-модели для прогно- зирования трендов во временных рядах, характеризу- ющих экономические данные. Прони [202] предложил простую процедуру для согласования экспоненци- альных моделей с данными экспериментов в области химии тазов. Ряд других моделей появился в областях статистического и численного анализа. Именно в этих нетехнических областях моделирования времен- ных рядов и берут свое начало современные методы спектрального оценивания. Широкое применение нетрадиционных методов спектрального оценивания началось в 1960-х годах. Парзен [189] в 1968 г. формально обосновал АР-метод спектрального оценивания. Независимо от него Берг [37] в 1967 г. опубликовал метод максимальной энт- ропии, мотивированный его работой по фильтрации на основе линейного предсказания в геосейсмологи- ческих приложениях. Ван-дер-Бош [255] формально показал, что одномерный метод МЭ эквивалентен АР- методу оценки СПМ. Метод Прони также имеет не- которое математическое сходство с АР-алгоритмами оценивания. Область современных исследований — это модели на основе авторегрессии и скользящего среднего (АРСС). АРСС-модель является обобщением АР-модели. Методы, основанные на АРСС-моделях, могут, по-видимому, обеспечить даже лучшее разре- шение и более высокие другие характеристики, чем АР-методы. Метод гармонического разложения Пи- саренко (ГРП) [194, 195] представляет собой пример метода спектрального оценивания, основанного на частном случае АРСС-модели. Используемый в статье унифицированный подход состоит в кратком рассмотрении каждого метода спектрального оценивания как метода, основанного на согласовании измеренных данных с предполагаемой моделью. Различия характеристик рассматриваемых методов спектрального анализа можно часто соот- нести с тем фактом, насколько хорошо применяемая модель аппроксимирует анализируемый процесс [173]. Различные методы могут давать одинаковые резуль- таты, но какой-либо один из них может требовать меньшего числа параметров модели, а следовательно, такой метод будет более эффективно отображать этот процесс. Для иллюстрации сказанного в разд. III с помощью различных методов вычисляются спектраль- ные оценки процесса, состоящего из нескольких синусоид и окрашенного гауссова шума. Этот процесс имеет и узкополосные, и широкополосные компонен- ты. Рассматриваемый там пример позволяет проил- люстрировать тот факт, почему одни методы спект- ральной оценки позволяют лучше оценивать узко- полосные составляющие спектра, а другие — широ- кополосные. Этот пример весьма убедительно говорит о том, что необходимо правильно понимать свойства использованной модели, прежде чем выносить суж- дение о характеристиках основанного на ней метода спектрального оценивания. Данный обзор разбит на пять разделов. Наиболее крупный из них — разд. II, в котором коротко рас- смотрены все обсуждаемые в статье методы. В разд. III помещены сводная таблица с кратким описанием всех рассмотренных методов и рисунки, которые позволяют наглядно сравнить различные современные методы спектрального оценивания. В разд. IV ко- ротко рассматриваются другие области применения методов спектрального оценивания, обсуждаемых в даннойстатье. Для того чтобы читатель мог быстро отыскать интересующий его вопрос, ниже приведен список названий этих трех разделов статьи и их подразделов. II. Обзор методов спектрального оценивания А. Определения спектральной плотности и ос- новные соотношения В. Традиционные методы (периодограммы, метод Блэкмана — Тьюки) С. Подходы к моделированию и идентификации параметров D. Методы моделирования с использованием ра- циональной передаточной функции
8 ТИИЭР, т. 69, № 11, ноябрь 1981 Е. Авторегрессивное (АР) оценивание СПМ F. Оценивание СПМ на основе скользящего сред- него (СС) G. Оценивание СПМ на основе авторегрессии и скользящего среднего (АРСС) Н. Метод гармонического разложения Писаренко (ГРП) I. Энергетическая спектральная плотность Прони К. Оценивание спектральных линий методом Про- ни L. Спектральное оценивание с помощью метода максимального правдоподобия (ММП) Кейпона III. Краткая сводка основных результатов IV. Другие применения методов спектрального оце- нивания А. Введение В. Экстраполяция и интерполяция временных рядов С. Предобеляющий фильтр D. Сжатие полосы частот Е. Спектральное сглаживание F. Формирование луча G. Решетчатые фильтры Н. Другие применения В статье не обсуждается применение полосно-ограни- ченных методов экстраполяции для спектрального оценивания, поскольку хороший систематический об- зор подобных методов представлен в работе [103]. И, наконец, в кратком заключительном разделе V затрагиваются вопросы, касающиеся направлений исследований современных методов спектрального анализа. II. ОБЗОР МЕТОДОВ СПЕКТРАЛЬНОГО ОЦЕНИВАНИЯ А. Определения спектральной плотности и основные соотношения Традиционные методы спектрального оценивания, которые в настоящее время реализуются на основе БПФ, как правило, характеризуются принятием большого числа различных компромиссных допуще- ний при попытке получить статистически надежные спектральные оценки. Эти компромиссы касаются выбора весовой функции (функции окна) и усреднения дискретизованных данных, полученных из случайного упроцесса, во временнбй и частотной областях, с тем чтобы как-то сбалансировать требования к снижению уровня боковых лепестков, к повышению эффектив- ности усреднения по ансамблю и обеспечить прием- лемое частотное разрешение. Для того чтобы коротко изложить основные положения стандартного спект- рального анализа, рассмотрим сначала случай де- терминированного аналогового сигнала х(0, представ- ляющего собой непрерывную функцию времени. Для большей общности функция x(t) будет здесь и ниже в статье полагаться комплексной. Если x(t) абсо- лютно интегрируема, т. е. если сигнал имеет конеч- ную энергию &: а= f |x(t)|2 dt<«>, (2.1) то существует непрерывное преобразование Фурье (НПФ) X (/) сигнала x(f), которое дается следующим выражением: Х(/) = J x(t) exp (-j'iitft) dt. (2.2) (Заметим, что (2.1) является достаточным, но не необходимым условием существования преобразова- ния Фурье [33].) Квадрат модуля этого преобразо- вания Фурье часто называют спектром S (В сигнала х((): §(/)= W)!2. (2.3) Теорема Парсеваля, записанная в виде |x(t)|2 dt=f |Х(/)|2 df, (2.4) представляет собой одну из формулировок закона сохранения энергии; энергия сигнала во временнбй области равна энергии его преобразования в ча- стотную область J_“ 8 (/)^(- Следовательно, 8 (/) — энергетическая спектральная плотность (ЭСП) в том смысле, что она характеризует распределение энергии сигнала по частоте. Если сигнал х(() дискретизуется с интервалом Ai, в результате чего получается диск- ретная последовательность хп=х(п А/), —оо<п<оо, то эту последовательность отсчетов можно представить как произведение исходной временнбй функции x(t) и бесконечного набора эквидистантных дельта-функ- ций Дирака 6((). Согласно теории распределений [33], преобразование Фурье от этого произведения можно записать в виде х\п = f f Z - лДод, *-00 [ л В -ОО exp (- jlitft) dt = = At £ хп exp (-/2д/лД1). (2.5) Выражение (2.5) соответствует аппроксимации пре- образования (2.2) с помощью прямоугольников; мно- житель А( гарантирует здесь сохранение области интегрирования при переходе от (2.2) к (2.5) при А/ 0. Выражение (2.5) идентично выражению (2.2) для преобразования X (/) на интервале —1/.(2 А()^ ^/^1/(2 А() Гц до тех пор, пока сигнал x(t) имеет ограниченную полосу и все его частотные состав- ляющие расположены в этом интервале. Следова- тельно, непрерывная энергетическая спектральная плотность 8'(/) = 1Х'(/)12 (2.6) для дискретизованных данных, полученных из про- цесса с ограниченной полосой частот, идентична ве- личине (2.3). Если а) последовательность данных обрабатывает- ся лишь с помощью конечного временного окна длительностью от л=0 до n=N— 1 й b) ее преобразо- вание Фурье также дискретизуется на N величин посредством взятия отсчетов на частотах f=m &f,
СПЕКТРАЛЬНЫЙ АНАЛИЗ где т=0, 1, . . М—1 и то из (2.5) ** можно получить более привычное дискретное пре- образование Фурье (ДПФ) [33]: N-1 Хт = At £ х„ ехр (-/2ят Д/лДг) = л-о N-1 = Дг J' х„ exp (-jlitmn/N), л-о m^O,-',N- 1. (2.7) Преобразование (2.7) и соответствующее ему обратное преобразование являются циклическими с периодом N. Следовательно, используя (2.7), мы периодически продолжаем и дискретизованные данные, и значения дискретизованных преобразований, даже если исход- ные непрерывные данные могут и не быть периоди- ческими. Теперь дискретную ЭСП можно определить кзк Sm=l*m|2, (2.8) где, как и раньше, O^rn^N—1. И дискретная ЭСП $т , и непрерывная ЭСП $'т получили название спектральных оценок на основе периодограммы. За- метим, что значения . Sm, и £'(/), вычисленные на частоте f=m!NAf, т=0, 1, . . ., N—1, не совпадают. Sm — это по сути дела дискретизованный вариант спектра, определяемого по свертке спектра X (/) с преобразованием прямоугольного окна, которое содержит отсчеты данных. Следовательно, дискретный спектр Sm, полученный по конечной последователь- ности данных, является искаженной версией непре- рывного спектра $'(/) > соответствующего бесконечной последовательности данных. Если процесс х(() представляет собой не детер- минированный сигнал конечной энергии, а стацио- нарный в широком смысле случайный процесс, то необходимо учитывать ряд других моментов. Энергия такого процесса, как правило, бесконечна, поэтому искомой величиной здесь будет являться частотное распределение мощности (средняя по времени плот- ность энергии). Кроме того, интегралы типа (2.2) для случайных процессов обычно не существуют. В случае стационарного случайного процесса основой спектрального анализа становится не сам случайный процесс, а его автокорреляционная функция Я»(т) = £[х(г + т)х*(г) ], (2.9) которая с помощью теоремы Винера — Хинчина по- средством преобразования Фурье связывается со спектральной плотностью мощности: ?(/)= f Rxx(.r)evp(.-i2vfr)dT. (2.10) Однако на практике статистическая автокорреляцион- ная функция обычно не известна, поэтому весьма часто принимается дополнительное допущение о том, *• Обратное преобразование, соответствующее выражению (2.5), дается выражением хи=Д)2т-аАжехр(+/2ляш/А), а теорема Персеваля записывается в следующей форме: N-1 М-1 Е Е 1хм1*а/. п*0 ЛсО Функция временя х(О — Вре*«вк4я , Уравнение (2.10) «... Спектральная автокорреляция • плотность мощности Рис. 2. Прямой и косвенный методы получения СПМ (в пред- положении свойств стационарности и эргодичности). что случайный процесс эр годичен относительно пер- вого и второго моментов. Это свойство позволяет заменить усреднение по ансамблю усреднением по времени. Тогда для эргодического процесса стати- стическую автокорреляционную функцию можно оп- ределить как лхх(т)“ lim f х(г + т)х*(г) Л. г-»« 2Т J_T (2.11) С помощью (2.11) нетрудно показать [107, 132, 187], что спектр (2.10) можно определить эквивалентным выражением ?(/)’ lim Е т-*- х(Г) ехр (-/2яД) dt (2.12) Оператор математического ожидания здесь необхо- дим, поскольку свойство эргодичности функции Rxx(x) не связано с через преобразование Фурье, т. е. предел в (2.12) без оператора математического ожи- дания не сходится в любом статистическом смысле. На рис. 2 показаны прямые и косвенные подходы к получению СПМ сигнала x(t), основанные на фор- мальных соотношениях (2.10), (2.11) и (2.12). При использовании выражения (2.12) могут воз- никнуть трудности, если оно применяется к конечной совокупности данных без выполнения операций ма- тематического ожидания и предельного перехода. Статистически несостоятельные (неустойчивые) оценки получаются в том случае, когда статистическое ус- реднение не выполняется; иными словами, дисперсия оценки СПМ будет стремиться к нулю, когда интервал усреднения Т неограниченно возрастает. В. Традиционные методы Два метода спектрального оценивания основаны на выполнении операций преобразования Фурье. Популяризации метода оценки СПМ, получаемой посредством косвенного подхода через оценку авто- корреляционной функции, способствовали Блэкман и Тьюки [14]. Второй метод оценки СПМ, основанный на прямом вычислении спектра с помощью алгоритма БПФ, теперь обычно называют периодограммой. При наличии конечной последовательности дан- ных можно оценить лишь конечное число дискретных значений автокорреляционной функции, или лагов. Блэкман и Тьюки предложили спектральную оценку м л ?вт(Л = Дг 22 Л«(т)ехр(-/2я/тД0, (2.13) п—М 2 ТИИЭР № 11
10 ТИИЭР, т. 89, № 11, ноябрь 1981 основанную на использовании имеющихся оценок автокорреляционной функции ,Rxx(/n), где —1/(2 Д/)^ ^/^1/(2Д0, а ~ обозначает оценку. Оценка (2.13) представляет собой дискретный вариант выражения Винера — Хинчина (2.10). Очевидно, что соответст- вующей оценкой автокорреляционной функции, осно- ванной на (2.11), будет несмещенная оценка 1 JV-m-1 = L хп+тхп, (214) П т П-0 где т=0, 1, . . -., М, a M^N—1. Оценки для отрица- тельных задержек определяются из оценок для поло- жительных задержек в соответствии с выражением £„(-m)-&»(»») (2.15) вследствии свойства комплексно-сопряженной сим- метрии автокорреляционной функции стационарного процесса. Вместо оценки (2.14) Дженкинс и Ватте [107] и Парзен [188, 1891 обосновали применение следующей оценки автокорреляционной функции: л, 1 IV-m-l ^xx(m) = — £ xn+mx*, (2.16) п «о определенной при т=0, . . ., М, поскольку она дает меньшую, чем (2.14), среднеквадратичную ошибку для многих конечных наборов данных. Оценка Rxx (m) является смещенной, поскольку E[Rxx(m)]=[(N— — m)/N]Rxx(m). Среднее значение этой оценки равно истинной автокорреляционной функции, взвешенной треугольным окном (иногда его называют окном Бартлетта). Прямой метод спектрального анализа представ- ляет собой современный вариант периодограммы Шу- стера. Дискретный вариант выражения (2.12), для которого данными измерения являются лишь отсчеты х0, . . ., xN-i, имеет следующий вид: ?Пер(-^= N-1------ 2 А* У. х„ ехр (-/2я/пДг) л-о (2.17) и также определяется на частотном интервале —1/(2Д/)</<1/(2 Д(). Заметим, что присутствующая в (2.12) операция вычисления математического ожи- дания здесь впервые опущена. Использование БПФ позволяет вычислять оценку (2.17) на дискретном мно- жестве из N эквидистантных частот Гц, /п=0, 1.....N— 1, Д/=1/УД/: = $пер V»") = I2» где Хт — ДПФ, определяемое выражением (2.7). Оценка 5*т идентична энергетической спектральной плотности Sm, определяемой выражением (2.8), за исключением деления на интервал времени УД/ с, требуемое для того, чтобы &т было оценкой спект- ральной плотности мощности. Полная мощность про- цесса, которая полагается периодической благодаря свойству ДПФ, в соответствии с аппроксимацией интеграла в оценке £пер. с помощью прямоуголь- ников будет равна JV-I . Мощность = 22 (2.19) т =0 Если множитель Д/ ввести в ^т, то получим -9тЫ= (ЛД/)2 |Хт |2 = 1 АГ-1 2 77 Е *л ехр (-/2ятл/ЛГ) " л «о (2.20) Эта величина часто вычисляется как периодограмма, но она надлежащим образом не промасштабирована как СПМ. С помощью (2.20) мы определяем не площадь под кривой графика, а пик на графике СПМ, величина которого равна мощности предполагаемого периоди- ческого сигнала. Благодаря вычислительной эффек- тивности алгоритма БПФ этот метод получил широкое распространение. Когда анализируемый процесс содержит детерми- нированные составляющие, погруженные в случай- ный шум, периодограмма для N отсчетов данных часто вычисляется с помощью выражения (2.18). Как указывалось выше, здесь следует соблюдать осторожность, поскольку при этом можно получить статистически несостоятельные результаты, если ис- пользовать его буквально, т. е. без учета операции математического ожидания. Эта необходимость в некоторого рода усреднении по ансамблю, или сгла- живании выборочного спектра, была продемонстри- рована Оппенгеймом и Шафером [183, с. 5461 и От- весом и Энохсоном [184, с. 328] на примере процессов типа белого шума, при которых дисперсия спект- ральной оценки не уменьшалась даже при исполь- зовании все более длинной и длинной последователь- ности данных. Бартлетт [18] указал статистические трудности, возникающие при использовании выра- жения (2.18), и предложил разбивать данные на сег- менты, вычислять оценки для каждого сегмента и усреднять периодограммы по всем сегментам. Уэлш [262, 263] предложил использовать специальную цифровую процедуру с БПФ, которая включает ус- реднение периодограмм. При других подходах для аппроксимации усред- нения по ансамблю используются окна во временной или частотной области или же одновременно в обеих областях [183, 184]. Для повышения устойчивости оценок и минимизации влияния боковых лепестков Наттолл и Картер [43, 44, 175, 179] предложили ус- реднять перекрывающиеся взвешенные сегменты. Ряд вопросов, касающихся БПФ и его применения для оценивания СПМ, рассмотрен в работах Берг- лэнда [19], Бертрама [22, 23], Бригхэма и Морроу [32, 33], Кокрэна и др. [52], Кули и др. 154, 55], Глис- сона и др. [77], Наттолла и Картера [43,181], Ричардса [207], Райфа и Винсента [208], Уэбба [259] и Юаня [274, 275). В общем случае спектральные оценки ^Вт (/) и ^пер.(Л не идентичны. Однако если используется смещенная оценка автокорреляции (2.16) и число вычисленных дискретных значений автокорреляци- онной функции равно числу отсчетов данных (М= =N— 1), то оценка Блэкмана — Тьюки (БТ) и оценка на основе периодограммы дают одинаковые численные результаты [184]. Следовательно, периодограмму мож- но рассматривать как частный случай процедуры оценивания Блэкмана — Тьюки. Именно по этой
СПЕКТРАЛЬНЫЙ АНАЛИЗ 11 (а) Доли частоты дискретизации (С) Рас. 3. Иллюстрация эффекта маскирования слабого сигнала соседним сильным сигналом при анализе спектра на основе периодограммы (число использованных отсчетов во всех случаях равно 16): (а) СПМ одиночной синусо- иды единичной амплитуды с частотой (в долях частоты дискретизации), равной 0,15, и начальной фазой 45°; (Ь) СПМ (относительно уровня предыдущей СПМ) оди- ночной синусоиды с амплитудой 0,19, частотой 0,24 и начальной фазой 162°; (с) СПМ объединенного сигнала; отметим слабый отклик на частоте, соответствующей слабому сигналу. причине обе эти оценки иногда называют подходом «сужай и преобразуй» (taper and transform,-ТТ) (1161. Многие трудности, связанные с применением ме- тода оценивания СПМ по периодограмме, можно проследить вплоть до тех допущений, которые были сделаны в отношении данных, находящихся вне ин- тервала измерения. Конечную последовательность данных можно рассматривать как полученную из бесконечной последовательности отсчетов с помощью перемещаемой прямоугольной функции окна. Ис- пользование только этих данных неявно подразуме- вает, что неизмеренные данные полагаются равными нулю, но обычно это не так. Умножение реального временного ряда на функцию окна означает, что полное преобразование будет представлять собой свертку требуемого преобразования с преобразова- нием функции окна. Если истинный спектр сигнала сосредоточен в узкой полосе частот, то такая операция свертки приведет к просачиванию мощности в сосед- ние частотные области. Это явление, названное про- сачиванием (или утечкой) мощности, является след- ствием того неявного взвешивания отсчетов последо- вательности, которое «внутренне» присуще вычисле- нию периодограммы. Помимо искажения спектральной оценки утечка ухудшает точность оценивания мощности и обнаружи- мость синусоидальных составляющих анализируемого процесса 1201, 224, 234]. Боковые лепестки из соседних частотных ячеек складываются или вычитаются с главным лепестком отклика в других частотных ячейках спектра, влияя тем самым на оценку мощ- ности в этой ячейке разрешения. В предельном случае боковые лепестки сильных частотных составяющих могут маскировать главные лепестки слабых частот- ных составляющих в соседних ячейках; иллюстрацией этого эффекта может служить рис. 3, на котором боковые лепестки функции (sin л/)/л/ (преобразо- вания прямоугольного окна во временнбй области) ухудшают оценку мощности главного лепестка сла- бого синусоидального сигнала. Взвешивание данных с помощью функции окна является также фундаментальным фактором, опреде- ляющим частотное разрешение периодограммы. Сверт- ка преобразования окна с преобразованием реального сигнала означает, что ширина наиболее узкого спект- рального отклика результирующего преобразования ограничивается шириной главного лепестка преоб- разования окна, которая не зависит от данных. Для прямоугольного окна ширина главного лепестка (а следовательно, разрешение) результирующего пре- образования (sinnf)/iif по уровню 3 дБ приблизи- тельно равна величине, обратной времени наблю- дения N&f с. Можно, конечно, использовать другие окна, но разрешение будет всегда пропорционально величине 1ЛУД/ Гц. Эффекты просачивания, обуслов- ленные взвешиванием, можно ослабить посредством выбора окна с неравномерным взвешиванием. В статье Хэрриса [90] дан хороший обзор свойств различных окон. Наттолл [182] уточнил характеристики боковых лепестков для некоторых из описанных Хэррисом окон. Среди других работ по этому вопросу можно назвать работы 112, 61, 161, 224, 260, 2731. Поэтому никаких результатов, суммирующих относительные качества различных функций окна, мы здесь не при- водим. Плата, которую приходится платить за сни- жение уровня боковых лепестков,— это всегда уши- рение главного лепестка преобразования окна, т. е. ухудшение разрешающей способности спектральных оценок. Широко распространено ошибочное мнение о том, что введение нулей в последовательность данных перед ее преобразованием улучшает разрешение пе- риодограммы. Однако преобразование дополненной
12 ТИИЭР, т. 89, № 11, иоябрь 1981 (с) (d) Рве. 4. Влияние введения нулевых отсчетов в периодограмму для интерполяции огибающей спектра и разрешения неопреде* ленностей. Показаны оценки спектров по 16 отсчетам процесса, состоящего из трех синусоид с частотами (в долях ча- стоты дискретизации), равными 0,1335, 0,1875, 0,3375, и начальными фазами 0°, 90° и 0° соответственно: (а) беа иведения дополнительных нулевых отсчетов;, в спектре присутствуют неопределенности; (Ь) число дополнительных нулевых от- счетов равно числу исходных отсчетов данных; неоднозначности разрешены; (с) число отсчетов данных увеличено в 4 раза за счет введения дополнительных нулевых отсчетов; форма спектра еще более сглажена; (d) число отсчетов данных увеличено в 32 раза за счет введения дополнительных нулевых отсчетов; огибающая спектра является аппроксимацией непрерывного преобразования Фурье. нулями последовательности данных служит только для получения в частотном интервале —ll(2M)^f^ ^1/(2Д() Гц дополнительных интерполированных значений СПМ между значениями СПМ, полученными по преобразованию исходной последовательности (не дополненной нулями). На рис. 4 показаны периодо- граммы М-точечной последовательности данных и этой же последовательности, дополненной N, 7N и 31М нулями. В каждом случае дополнительные зна- чения периодограммы, вычисленной с помощью БПФ по исходной последовательности, постепенно при- ближают ее к непрерывной периодограмме, опреде- ляемой выражением (2.17). Однако ни в одном из случаев введения дополнительных нулей улучшения разрешения основных частотных составляющих не наблюдается (поскольку оно обратно пропорцио- нально интервалу наблюдения). Введение дополни- тельных нулей полезно для а) сглаживания оценки периодограммы за счет получения интерполированных значений, Ь) для устранения потенциальных неопре- деленностей, показанных на рис. 4, и с) для умень- шения влияния ошибки <квантования> на точность оценивания частот спектральных линий. Математи- чески 2М-точечное ДПФ 2М-точечной последова- тельности х„ . . ., xilf-i записывается в виде 2JV-1 Хт=Д/ 22 хп ехр (-/2ятп/2М), (2.21) п *0 где т=0, 1, . . ., 2N—1. Если эта последовательность дополнена N нулями, т. е. хп—0 при n=N, . . ., 2N—1, то выражение (2.21) будет иметь следующий вид: 'О да 1 /2я — n/Nj, (2.22) что совпадает с М-точечным преобразованием (2.7), но вычисляется на интервале —1/(2Д0^/^1/(2Д(), который содержиТ'Вдвое больше частот, чем интервал для (2.7). Исключая вычислительные операции над дополнительно введенными нулями, или вообще не обращая на эти нули внимания, можно получить более эффективный алгоритм вычисления БПФ 11481. Итак, стандартные подходы Блэкмана — Тьюки и на основе периодограммы к спектральному оцени- ванию обеспечивают следующие преимущества: 1) вы- числительную эффективность, если используется толь-
СПЕКТРАЛЬНЫЙ АНАЛИЗ 13 ко несколько дискретных значений автокорреляцион- ной функции (БТ-метод) или если используется БПФ (периодограмма); 2) оценка СПМ прямо пропорцио- нальна мощности синусоидальных процессов; 3) они являются хорошей моделью для некоторых прило- жений (этот вопрос более подробно рассмотрен в следующем подразделе). К недостаткам этих методов следует отнести: 1) подавление главных лепестков слабых сигналов боковыми лепестками сильных сиг- налов; 2) частотное разрешение ограничено длиной имеющейся записи анализируемых данных и не за- висит от характеристик этих данных и их отношения сигнал/шум; 3) искажение спектра из-за просачи- вания энергии боковых лепестков; 4) необходимость какой-либо разновидности псевдоусреднения по ан- самблю для получения статистически состоятельных оценок спектров периодограмм; 5) появление отри- цательных значений СПМ при использовании БТ- метода в случае, когда применяются некоторые ав- токорреляционные оценки последовательностей. С. Подходы к моделированию и идентификации параметров В этом подразделе мы отвлечемся от рассмотрения традиционных ^методов спектрального анализа, ко- торым был посвящен предыдущий подраздел. При стандартных подходах используются операции вы- числения БПФ либо по взвешенным данным, либо по взвешенным значениям автокорреляционной функции. Взвешивание данных или значений автокорреляцион- ной функции с помощью функции окна предполагает неявное допущение о том, что ненаблюдаемые данные или значения автокорреляционной функции вне этого окна равны нулю, что обычно является нереальным допущением. Следствием такого допущения является искажение спектральных оценок. Часто о процессе, из которого берутся отсчеты, мы располагаем доста- точно большими сведениями или по крайней мере можем сделать более разумные допущения, чем просто предположить, что отсчеты равны нулю вне окна. Использование априорной информации (или пред- спектральной оценки с помощью подстановки оценок параметров модели в расчетное выражение для СПМ, соответствующее этой модели. Основная причина наблюдаемого сейчас интереса к подходу, моделирую- щему подход к спектральному оцениванию,— это более высокое частотное разрешение, достижимое с помощью этих современных методов, по сравнению с тем разрешением, которое достижимо с помощью традиционных методов спектрального анализа, о которых говорилось выше. Степень улучшения раз- решения и точности спектральной оценки будет оп- ределяться имеющейся возможностью согласования выбранной модели с малым числом параметров с измеренными данными. Выбор модели для спектраль- ной оценки тесно связан с методами оценивания и идентификации параметров, используемыми в теории линейных систем [271, 272]. Для иллюстрации моделирующего подхода к спект- ральному оцениванию покажем, что оценка СПМ дискретной периодограммы (2.18) эквивалентна сред- неквадратичной аппроксимации данных с помощью гармонической модели, а именно с помощью ряда Фурье. Метод среднеквадратичной аппроксимации с помощью ряда Фурье достаточно хорошо изучен [28], поэтому здесь мы изложим лишь основные идеи этого метода. Если N отсчетов х0, . . ., xN-i непрерыв- ного процесса x(t) моделируются дискретной после- довательностью хп, состоящей из N комплексных синусоид с произвольными частотами Д, . . ., fjy-i, то лг-i х(пДг) = хп= £ ат exp (2.23) т » О где п=0, .... —1. Таким образом, сигнал х(1) на интервале ЛГА/ с представляется периодическими функциями независимо от того, периодичен или нет процесс x(t) сам по себе. Для большей общности веса ат, которые необходимо определить, полагаются комплексными. Выражение (2.23) можно записать в матричной форме Х=ФА, (2.24) где 1 1 ••• 1 .exp (Хо ) exp (Xi) • • • exp (Xjy. j) _ехр (XoUV-11) exp (Xi [АГ- 1]) ••• exp (Xjy^UV- 1])_ положений) может облегчить выбор точной модели процесса, отсчетами которого мы располагаем, или по крайней мере такой модели, которая является хорошей аппроксимацией реального процесса. В этом случае можно, как правило, получить более точную спектральную оценку, определяя параметры выбран- ной модели по результатам измерений. Таким обра- зом, применительно к моделированию спектральный анализ становится трехэтапной процедурой. Первый этап состоит в выборе модели временного ряда. Второй этап состоит в оценивании параметров принятой модели либо с использованием отсчетов имеющихся данных, либо значений автокорреляционной функции (которые либо известны, либо оцениваются по имею- щимся данным). Третий этап состоит в получении и X|=/2n/fAt Если N частот fm заданы, то вектор амплитуд А, определяемый посредством минимизации суммарной среднеквадратичной ошибки оценивания /v-t У. |хи-х„|2, Л “О (2.25) получается из решения уравнения А = (ФнФ)~'Фнх, (2.26) где *о _XN-1
14 ТИИЭР, т. 69, №11, ноябрь 1981 — вектор данных, а верхний индекс Н означает транс- понирование комплексно-сопряженной матрицы. Аппроксимация рядом Фурье требует предвари- тельного выбора частот N синусоид, гармонически связанных между собой; fm=rnbf Гц, где Д/=1/АД/. Как известно [28], такой выбор гармонических частот делает каждый столбец (строку) вектора Ф ортого- нальным всем остальным столбцам (строкам) этого вектора, поэтому <ФНФУ1=^1, (2.27) где I — единичная матрица. Следовательно, урав- нение (2.26) можно записать в следующем виде: A**—#1* (2.28) или 1 ЛГ-1 ат = Tr Z ех₽ (-/2ятл/Я), ™ л«о где /п=0, .... N—1. Мощность синусоидальных со- ставляющих на заданных частотах fm равна \ат I2 = — У хп exp (-jlitmnlN) N (2.29) Выражение (2.29) идентично выражению (2.20). Та- ким образом, спектральная оценка дискретной перио- дограммы может рассматриваться как среднеквад- ратичная аппроксимация данных с помощью гармони- ческого набора комплексных синусоид. Случай, когда частоты заранее не выбираются и необязательно связаны гармоническим соотноше- нием, рассмотрен в подразд. II.J при обсуждении метода Прони. В гармонической модели частоты и число синусоид задается заранее, поэтому необхо- димо оценивать только мощность этих синусоидаль- ных составляющих. Негармоническая модель, ис- пользуемая в методе Прони, требует оценки не только мощности, но и числа синусоид и их частот. В от- ношении гармонической модели необходимо отметить и тот факт, что в этой модели наличие шума никак не учитывается. Иными словами, любой присутст- вующий шум также моделируется с помощью гармо- нически связанных синусоид, и, следовательно, для уменьшения флюктуаций, обусловленных наличием шума, необходимо выполнять усреднение по некоторой совокупности периодограмм, определяемых по ана- лизируемым данным. Главная особенность моделирующего подхода к спектральному оцениванию, которая отличает его от обобщенной задачи идентификации, состоит в том, что для- анализа можно использовать только процесс на выходе модели, а входной процесс пола- гается как бы недоступным для решения общей задачи идентификации. Один из многообещающих аспектов моделирую- щего подхода к спектральному оцениванию обусловлен тем фактом, что можно принять более реалистичные допущения относительно измеряемого процесса вне интервала измерения, а не просто положить, что он равен нулю вне этого интервала или что он перио- дичен. Следовательно, можно отказаться от исполь- зования функции окна и тем самым устранить его искажающее воздействие. В результате точность оценки по сравнению со стандартным подходом на основе БПФ может значительно улучшиться, особенно в случае коротких записей данных. D. Методы моделирования с использованием рациональных передаточных функций Многие детерминированные и случайные процессы с дискретным временем, с которыми приходится иметь дело на практике, хорошо аппроксимируются моделью на основе рациональной передаточной функ- ции. В этой модели входная последовательность {п^} и выходная последовательность {хп}, которые используются для моделирования данных, связаны линейным разностным уравнением вида я р *п = Z blnn-l ~ £ акхп-к- (2.30) 1-0 к-1 Эта наиболее общая линейная модель называется АРСС-моделыо. Интерес к подобным моделям обус- ловлен их связью с линейными фильтрами, обладаю- щими рациональными передаточными функциями. Передаточная функция системы, связывающая вход пп и выход хп для АРСС-процесса (2.30),. опреде- ляется рациональным выражением B(z) Я(г) = ~77, (2.31) Л(г) где г-преобразования АР- и СС-частей процесса даются выражениями р Л(г)|= £ amz~m) т О Q £ bmz~m. т жО Как известно, спектр мощности на выходе линейного фильтра связан со спектром мощности Рх(г) выход- ного случайного процесса следующим соотношением: Рх(г) = Я(г)Я‘(1/г*)Р„(г)= ^^.рп(2). (2.32) Выражение (2.32) обычно вычисляется на единичном круге г=ехр(/2л/Д(), где —1/(2Д()^/^1/(2Д(). Ча- сто входной процесс полагается белым шумом с ну- левым средним значением и дисперсией <т!. СПМ этого шума равна <т’Д(. (Заметим, что мы включили множитель Д( в выражение для спектральной плот- ности мощности шума, а следовательно, после ин- тегрирования Рх(ехр[/2я/Д/]) на интервале —1/ (2Д/)^ <3/^1/(2Д() получим истинное значение мощности аналогового сигнала.) СПМ выходного АРСС-про- цесса будет, следовательно, равна ® АРСС^ = ?«(/) = I® (/)/«(/) I2, (2.33) где в (/)=А (ехр[/2л/Д(]) и ®(/)=В(ехр[/2л/Д(1). Та- ким образом, определение параметров {afe} (называ- емых коэффициентами авторегрессии), параметров {£>*} (называемых коэффициентами скользящего сред- него (СС)) и а1 эквивалентно определению спектра
СПЕКТРАЛЬНЫЙ АНАЛИЗ 15 процесса {хп}. Без потери общности можно положить о,= 1 и 60=1, так как коэффициент усиления любого фильтра можно ввести в о*. Если все коэффициенты {ак}, за исключением Оо=1, стремятся к нулю, то я xn = Z. ьШп-1 (2.34) 1-0 и процесс будет строгим процессом скользящего среднего порядка q, и ?сс(Г) = а’ДИ«(/')1Е. 2. (2.35) Эту модель иногда называют чистонулевой моделью 1266]. Если все коэффициенты {6(}, за исключением 60=1, равны нулю, то р хп ~ ~ 22 + пп (2.36) к-1 и процесс будет строгим процессом авторегрессии порядка р. Этот процесс называется АР-процессом в том смысле, что процесс хп линейно регрессирует сам на себя, а пп представляет собой ошибку. При использовании этой модели текущее значение про- цесса выражается через взвешенную сумму его пре- дыдущих значений и шумовой член. СПМ в этом случае равна a2 At Эту модель иногда называют чистополюсной моделью. Теорема декомпозиции Уолда [266], связывающая АРСС-, СС- и АР-модели, утверждает, что любой стационарный АРСС- или СС-процесс с конечной дисперсией можно представить единственной АР- моделью, возможно, бесконечного порядка; и ана- логично, любой АРСС- или АР-процесс можно пред- ставить как СС-процесс некоторого (возможно, бес- конечного) порядка. Важность этой теоремы обуслов- лена тем фактом, что если мы из этих трех моделей выбираем наименее подходящую модель, то мы все же можем получить приемлемую аппроксимацию, увеличив порядок модели. Следовательно, любую АРСС-модель можно аппроксимировать АР-моделыо более высокого порядка. Поскольку оценивание па- раметров АР-модели приводит к линейным уравне- ниям, то, как будет показано ниже, она имеет вы- числительные преимущества перед методами оце- нивания параметров АРСС- и СС-моделей. Именно по этой причине большая часть усилий по моделиро- ванию на основе рациональных передаточных функ- ций посвящена методам на основе АР-модели. Е. Авторегрессивное оценивание СПМ Введение. Поскольку в этом подразделе мы под- робно остановимся на работах, посвященных авто- регрессивному методу оценивания СПМ, целесообразно коротко представить материал, который будет в нем изложен. Сначала будет показан вывод уравнений Юла — Уокера, которые описывают линейное со- отношение между параметрами АР-процесса и автокор- реляционной функцией. Решение этих уравнений ищется с помощью эффективного в вычислительном отношении алгоритма Левинсона — Дербина; анализ полученного решения позволяет вскрыть некоторые фундаментальные свойства АР-процессов. Затем анализируется связь между АР-моделиро- вайием, теорией линейного предсказания и спект- ральным оцениванием по максимуму энтропии. Далее рассматриваются вопросы выбора порядка АР-модели и ее параметров. Обсуждаются также методы группо- вого оценивания, основанные на теории линейного предсказания, и методы последовательной оценки, основанные на рекурсивных алгоритмах среднеквад- ратичной оценки и адаптивных алгоритмах. Отме- чаются компромиссы между различными подходами. Наконец, рассматриваются некоторые ограниче- ния авторегрессивных методов оценивания, которые сужают область их практического применения. К ним относятся различные вредные эффекты, обусловленные наличием шумов и ложных спектральных линий, а также ряд аномальных эффектов при наличии в данных доминирующих синусоидальных составляю- щих. Описаны некоторые методы ослабления подобных эффектов. Уравнения Юла — Уокера. При использовании авторегрессивной модели оценку спектральной мощ- ности АР-процесса, основанную на выражении (2.37), можно переписать в следующем виде: ?ар(/) = 1Я(ехр [;2w/At ] )|2 ?П(Г) = о2 At (2.38) 1 + £ ак ехр (-/2tr/fcAt) k=l Таким образом, для оценки СПМ необходимо лишь оценить значения {at, а....... ар, о*}. Для этого необходимо определить соотношение между парамет- рами АР-модели и автокорреляционной функцией (известной или оцениваемой). Это соотношение из- вестно под названием уравнений Юла — Уокера [31]. Их вывод проводится следующим образом. Запишем Яхх(£)=£[хп+кх*] = Е р 22 а1&хх(к ~ I) + Elnn+kxni 1-1 Поскольку Н(г) полагается устойчивым каузальным фильтром, то £(пп+кхп) = £ |nn+k У'. flfnn-l L 1-0 = 22 Ь*°28к + 1 = 1—0 = a2h*k = ГО, fc>0; [Лоа2, к = 0. Заметим что 6Ж — это дискретная дельта-функция, т. е. бж=1, если т=0, и 6ж=0, если лпМ). Но Ле=
16 ТИИЭР, т. 69, № 11, ноябрь 1981 = Ит2_»жЯ (г)=1, поэтому р ~ Е -7), Z = 1 RxxW=< fc>0; (2.39) к = 0. р - Е alRxxH)^O2, < ,=1 Выражение (2.39) и представляет собой уравнения Юла — Уокера. Для определения параметров АР- модели необходимо лишь выбрать из (2.39) р урав- нений при fc>0, решить их относительно {ai, аг, . . ., ар} и затем определить а® из (2.39) при k=0. Система уравнений, которая требует наименьшего числа дискретных значений функции автокорреляции, и определяет значения Л=1, 2, . . ., р. Эти уравнения можно записать в матричной форме Необходимо заметить, что коэффициенты {ak, 'Лхх(О) R^-D • *X«(1) *xx(0) e< ’ai' a2 = - 1 T-' OS (2.40) RxxtP - 1) RxxlP - 2) • • «^(0) aP яхх(р) Заметим, что в (2.40) автокорреляционная матрица Rxx является эрмитовой, так как Rxx—Rxx> и тёп- лицевой, так как элементы вдоль любой диагонали одинаковы. Кроме того, эта матрица является поло- жительно определенной (мы полагаем, что процесс хп не является чисто гармоническим), что следует из свойств положительной определенности автокорре- ляционной функции [41, 81, 92]. Следует заметить, что систему (2.40) можно рас- ширить, введя в нее уравнение для о®; в результате получим \«(0) Я„(-1) ••• Яхх(1) RxxW Rxx{p) RXX(.P - 1) • • • что следует из (2.39). Уравнение в такой форме по- надобится нам позднее. Таким образом, для опре- деления АР-параметров и о®, необходимо решить систему уравнений (2.41) относительно р+1 оцени- ваемых значений автокорреляционной функции Rxx (0), • • •> Rxx(p) и использовать свойство Rхх(—пг)= =Rzx(m). Алгоритм Левинсона — Дербина. Алгоритм Ле- винсона— Дербина [56, 60, 142, 270] обеспечивает эффективное решение системы уравнений (2.41). Этот алгоритм требует только р® операций (их число мы будем обозначать через о(р®)) в отличие от алгоритма гауссова исключения, который требует о(р’) опера- ций. И хотя сначала он был предложен только как эффективный алгоритм, в нашем случае он позволяет вскрыть фундаментальные свойства авторегрессив- ных процессов. Этот алгоритм позволяет рекурсивно вычислять набор параметров {au, a?}, {asx, ait, a*}, • • {ар1, api, . . ., app, dp}. Заметим, что дополни- тельный подстрочный индекс у АР-коэффициентов введен здесь для обозначения порядка АР-процесса. Последний набор коэффициентов порядка р представ- ляет собой искомое решение. В частности, этот ре- курсивный алгоритм инициализируется коэффициен- тами «и ^-RxxWIRxxiO), (2.42) О? “(1 - |eu|2)«xx(0), (2.43) а рекурсии для k=2, 3, . . р вычисляются с помощью выражений [к-1 /=1 (2.44) >-1 = (2.45) at = (1 ’ lafcfcl2)Ofc-i- (2.46) ак2, . . ., akh, ok], полученные выше, совпадают с коэффициентами, которые можно получить из-решения системы уравнений (2.41) при р=+. Следовательно, алгоритм Левинсона — Дербина позволяет также по- лучить АР-параметры для всех аппроксимирующих АР-моделей более низкого порядка. Это — весьма полезное свойство данного алгоритма, которое можно использовать в тех случаях, когда требуемый порядок модели заранее не известен, поскольку выражения (2.42) — (2.46) можно при этом использовать для генерации моделей последовательно возрастающего порядка до тех пор, пока ошибка моделирования о| не уменьшится до приемлемой величины. В част- ности, если процесс действительно является АР(р)- процессом (АР-процессом порядка р), то ap+1=apk, где k=l, 2, . . ., р, и, следовательно, ap+i, p+i=0. В общем случае для любого АР (р)-процесса аЛй=0 и о1=ор при k>p. Таким образом, дисперсия вход- ного шума является константной для модели, порядок которой равен или больше порядка правильно вы- бранной модели. Следовательно, момент, начиная с которого величина о| больше не меняется, может служить хорошим индикатором правильного выбора порядка модели. Можно показать, что |алй1 = 1, по- этому [9, 41], а это означает, что о* впервые достигает своего минимума при правильном выборе порядка модели. Ниже мы вернемся к этому вопросу при обсуждении подходов к определению порядка модели. Параметры {an, a22, . . ., арр} часто называют коэффициентами отражения и обозначают через {Кх, К2, • -, Они обладают тем свойством, что для того чтобы {Rxx(0)t Rxx(1), , Rxx(p)} действи- тельно образовывали последовательность значений автокорреляционной функции, т. е., чтобы автокор- реляционная матрица была положительно полуопре- деленной, необходимо и достаточно, чтобы lah)tl = = IKftl=l про й=1, 2, . . ., р [41]. Кроме того, не- обходимым и достаточным условием того, чтобы полюсы А (г) были расположены на единичной ок-
СПЕКТРАЛЬНЫЙ АНАЛИЗ 17 Рис. 5. Решетчатая структура предсказывающего ошибку (обеляющего, или обратного) фильтра. ружности или внутри единичного круга на г-пло- скости, является условие при k=l, 2, . . ., р [41, 140]. Следует заметить, что если при некоторых й|К*1 = 1, то рекурсии (2.42) — (2.45) должны за- канчиваться, поскольку а*=0. Процесс в данном случае является чисто гармоническим (т. е. содержит лишь синусоидальные составляющие). Все эти. свой- ства перечислены в табл. 1. Задача идентификации параметров АР-процесса тесно связана с теорией линейного предсказания. Пусть хп — АР (р)-процесс. Если мы хотим пред- сказать значение хп по предыдущим р отсчетам, т. е. р хп = - X акхп-к. (2.47) fc=l то параметры {«i, а2, . . ., ар) можно выбрать так, чтобы минимизировать мощность ошибки предска- зания Qp, где QP =М1*П - £п|2]. (2.48) Согласно принципу ортогональности [187], - *nWtl =0 при к = п - 1, • • •, л - р или р *xx(fc) = ~ X «!*«(* - О при к = 1, 2, • • • , р. (2.49) 1=1 Минимальная мощность ошибки предсказания равна р GPmin =£Kxn *«)**] =Я*х(0)+ X (2.50) к-1 Это уравнение идентично уравнению (2.39). Следо- вательно, должно быть справедливым, что afc=apfc при k=l, . . ., р и Qpmin—^p, поэтому наилучшим результатом линейного предсказания будет просто xn=—E^+1apfcxn-h. Последовательность ошибок, хо- тя и не коррелирована с этой линейной оценкой, не обязательно будет белым шумовым процессом (она будет таковым, если хп является АР(р)-процессом). В пределе, когда р -* оо, последовательвость ошибок становится белым шумовым процессом. Заметим, что ............} и а} представляют собой параметры оптимального линейного устройства предсказания k-ro порядка и соответствующую мини- мальную мощность ошибки предсказания. Следова- Таблица 1 Свойства АР-процессов (включая чисто гармонические процессы) • Автокорреляционная матрица является положительно опре- деленной: XHRxxX > 0 для всех векторов X » Последовательность коэффициентов отражения удовлетво- ряет условию | К; | < 1, i = 1, 2.р • Нули А (г) лежат внутри единичного круга: | z,-1 < 1, 1=1, 2,.... р • Мощности ошибок предсказания монотонно убывают: а?Э=а|ё=...Э=а2=0 3 ТИИЭРМН тельно, идентификация параметров АР-процесса и линейное предсказание АР (р)-процесса дают иден- тичные результаты, и поэтому теоретические подходы, разработанные в теории идентификации параметров и в теории линейного предсказания, применимы в каждой из этих областей. Теория линейного предсказания приводит к важ- ной интерпретации алгоритма Левинсона — Дербина. Пусть ерп — ошибка предсказания для устройства предсказания р-го порядка. Тогда, используя (2.45), получим р еря=хл+Х аркхп-к = (2.51) к-1 Р-1 ~хп + X (sp-i,fc + ^рар-1,р-кУхп-к + Крхп-р(2.52) к=1 ер-1,п +КрЬр-1,п-1> где р bpn~xn-p + X аркхп-р+к- (2.53) к-1 Член Ьрп представляет собой ошибку предсказания назад, т. е. ошибку предсказания значения хп~р по отсчетам хп-р+1, . . ., хп. Нетрудно видеть, что коэффициенты предсказания, вырабатываемые уст- ройством предсказания назад, являются комплексно- сопряженными коэффициентами предсказания вперед, которые представляют собой последовательность зна- чений стационарной автокорреляционной функции [41]. Аналогично можно показать, что Ьрп ~ bp-i,n-i + ^рер-1,п- (2.54) Выражения (2.52) и (2.54) приводят к так называемой структуре решетчатого фильтра, показанной на рис. 5 [146]. Заметим, что передаточная функция полного фильтра дается выражением р A(z) = 1 + X “рк2~к, к-1 которое является обращением функции Я(г)=1/А (г). Последнее соотношение следует из (2.51). Данный фильтр называется «обратным» фильтром, или фильт- ром, «предсказывак^щим ошибку». Реализация алго- ритма Левинсона — Дербина на основе решетчатого фильтра приводит к важной временной рекурсивной формулировке метода оценивания АР-параметров. Обновляя коэффициенты отражения с помощью вы- ражений (2.52) и (2.54), получим значения АР-пара- метров, выраженные через рекурсию (2.45). Под робности будут даны ниже при рассмотрении после- довательной оценки АР-параметров. При минимизации Е(|ерп1’), определяемой фор- мулой (2.52), относительно величины КР, получим
18 ТИИЭР, т. 69, № 11, ноябрь 1981 Как отмечалось выше, и еРп, и Ьрп имеют одинаковые статистические свойства, поэтому (2.55) можно за- писать в следующем виде: к = _____ Таким образом, К.р является взятым с обратным зна- ком нормированным коэффициентом корреляции меж- ду вр-1 и bp-lt п-1, поэтому МОЖНО ПОЛОЖИТЬ IKpl^l. В литературе по статистике —КР называется частным коэффициентом корреляции, поскольку он представ- ляет собой нормированный коэффициент корреляции между значениями хп и хп-р [31]. Спектральное оценивание на основе метода мак- симальной энтропии [41, 62, 78, 171, .172, 205, 215, 248, 249]. Данный метод основан на экстраполяции сегмента известной автокорреляционной функции на значения временных задержек, где ее вид неизвестен. Подобным образом можно устранить размытие спект- ральной оценки мощности из-за усечения автокорре- ляционной функции. Если значения {/?хх(0), 7?хх(1), . . ., /?хх(р)} известны, то возникает вопрос, как следует выбрать значения {Яхх(р+1), Яхх(р+2), . . .}, для того чтобы гарантировать положительную полуопределенность всей последовательности значе- ний автокорреляционной функции. В общем случае возможно бесконечное множество экстраполяций, при- водящих к желаемой автокорреляционной функции. Берг [37, 41] показал, что возможна такая экстрапо- ляция, при которой временной ряд, характеризуемый экстраполируемыми значениями автокорреляционной функции, будет иметь максимальную энтропию. Такой временной ряд будет максимально случайным рядом, который имеет данный набор из р+1 значений функ- ции автокорреляции. И наоборот, спектральная плот- ность мощности является спектральной плотностью с наиболее плоским (наиболее белым) спектром из всех спектров, для которых известны значения {/?хх(0), /?хх(1), . . ., Rxx(.P)}- Результирующая спектральная оценка называется спектральной оценкой по макси- муму энтропии. Основная причина выбора критерия по максимуму энтропии состоит в том, что такой критерий налагает наименьшее число ограничений на известный временной ряд за счет максимизации его случайности, что приводит к смещенной оценке. В частности, в случае гауссовского случайного процесса энтропия на отсчет пропорциональна вели- чине f 1/2ДГ I ln?x(/)dr, (2.57) где 5*х(/)— СПМ хп. Величина 5*я(/) определяется посредством максимизации интеграла (2.57) при на- ложенных ограничениях, что р+1 известных значе- ний автокорреляционной функции удовлетворяет со- отношению Винера — Хинчина [62]: {1/2ДГ ?х(/) ехр (-/2ff/hAt) df= Rxx(.n), l/2At л = 0, 1, •••,₽• (2.58). Решение находится с помощью метода множителей Лагранжа: а/Дг 1 + £ аРк ех₽ (-/2я/ЬД Г) (2.59) 2 где {ар1.....арр} и о*р — соответственно параметры устройства предсказания р-го порядка и мощность ошибки предсказания. При известных (Рхх(0), Rxx(1), ... ., Rxx(p)} метод спектральной оценки по макси- муму энтропии будет эквивалентен методу спектраль- ной АР-оценки с ^х(/)> даваемой выражением (2.59) и определяемой параметрами, получаемыми из реше- ния (2.41). Выражение для спектральной АР-оценки мощности, полученное на основе максимальной энт- ропии, справедливо лишь для гауссовых случайных процессов и известных значений автокорреляцион- ной функции. Автокорреляционное обобщение АР-оценки. Еще одно выражение для (2.38) дано в работе 159]: = Дг 22 г*х(") ез^Р (_/2я/иД0» (2.60) п 00 где бсх(п) = 4 RxxW> |л|<р; " £ apkrxx(n- к), . к-1 М>р. (2.61) Отсюда нетрудно видеть, что спектральная АР- оценка мощности сохраняет известные значения авто- корреляционной функции й рекурсивно продолжает их за пределы используемого окна. Суммирование в (2.60) идентично суммированию в (2.13) вплоть до значения, соответствующего р-му значению автокор- реляционной функции, но продолжает экстраполяцию автокорреляционной функции в бесконечных преде- лах, а не полагает их равными нулю вне окна. Таким образом, спектр, получаемый посредством АР-оценки, не обладает боковыми лепестками, обусловленными воздействием окна. Кроме того именно подразуме- ваемая выражением (2.61) экстраполяция, опреде- ляет высокие разрешающие свойства спектральной АР-оценки. АР-оценивание параметров. Во многих практиче- ских случаях для процедуры спектрального оцени- вания доступны отсчеты данных, а не значения авто- корреляционной функции. Для получения надежных оценок АР-параметров можно использовать стандарт- ную статистическую теорию оценивания. Для оценки неслучайного множества параметров обычно исполь- зуется оценка по максимуму правдоподобия (ОМП). Однако точную СМП-оценку АР-процесса получить затруднительно [31]. Если №^>р, то приближенную ОМП-оценку можно получить, решая уравнение Юла — Уокера для случая, когда искомая автокор- реляционная функция заменяется соответствующими оценками. Заметим, что при использовании длинных
СПЕКТРАЛЬНЫЙ АНАЛИЗ 19 записей данных АР-оценка дает хорошую спектраль- ную оценку. Однако при использовании коротких записей данных, с которыми приходится иметь дело на практике (и для которых спектральная оценка на основе периодограммы обеспечивает худшее раз- решение), применение подхода Юла — Уокера при- водит к оценкам с худшим спектральным разреше- нием. Для улучшения разрешения приближенного ОМП- подхода в случае коротких записей данных предложен ряд методов группового оценивания, основанных на среднеквадратичном критерии. Для длинных записей данных предложен ряд методов последовательной оценки, с помощью которых производится обновление АР-оценок по мере поступления новых данных. Эти методы особенно полезны для слежения за процес- сами, которые медленно меняются во времени. Сле- дующие два подраздела посвящены рассмотрению группового и последовательного методов АР-оцени- вания. Статистические свойства спектральных АР-оценок исследуются в многих работах [13, 100, 116, 186, 212]. Можно также указать ряд работ, посвященных ана- лизу точности частотного оценивания спектральных АР-оценок [124, 212, 239]. Сакаи [212] эксперимен- тально показал, что частотная дисперсия спектраль- ной АР-оценки обратно пропорциональна и длине записи данных, и квадрату отношения сигнал/шум; Килер [124] также экспериментально показал, что величина этой дисперсии обратно пропорциональна как длине записи данных, так и величине отношения сигнал/шум (а не квадрату этого отношения). Групповая оценка АР-параметров. Несмотря на то, что уравнения Юла — Уокера можно решить, используя оценки значений автокорреляционной функ- ции, предложен ряд процедур на основе критерия среднеквадратичной ошибки, которые позволяют не- посредственно по данным получить более точные оценки АР-параметров. По сравнению с методом Юла — Уокера эти методы часто позволяют получить более точные АР-спектры. Ниже мы рассмотрим два метода среднеквадратичной оценки. В первом из них используется линейное предсказание оценки вперед, во втором — комбинация линейного предска- зания вперед и назад. Предположим, что для отыскания АР-оценок пара- метров используется последовательность данных ха, . . ., Xjy-i. Для предсказания вперед мы будем ис- пользовать обычную форму хп~~ вркхп-к- (2.62) Предсказание вперед осуществляется в том смысле, что значение текущего отсчета данных ищется в виде взвешенной суммы р предыдущих отсчетов. Ошибка линейного предсказания вперед равна л р ерп=хп~ хп = £ -к, (2.63) где Дро = 1. Ошибку ерп можно вычислить для значений п от п=0 до n=N+p—1, если положим, что члены вне изме- рений равны нулю, т. е. хп=0 при п<0 и n>N—1. Мы можем применить взвешивание последовательности данных, для того чтобы распространить индекс ошибки ерп от 0 до N+p—1. Используя матричную запись, имеем Е !>о eN-l _eN+p-l_ ИЛИ (2.64) (2.65) Ошибка предсказания равна просто = 2L 1еР«|2 = £2 £ “ркХп-к п п к*0 (2.66) Диапазон суммирования для &р здесь пока намеренно не определен. Для минимизации величины 6Р произ- водные бр по {орь} приравниваются к нулю, после чего результирующие уравнения решаются относи- тельно АР-параметров: р / \ S °рк [ 2". хп-кхп-I] ~ о, (2.67) к-о \п / В результате получим р / \ ёР = Z арк ( L -кХп ) • (2.68) к=о \п / Выражения (2.67) и (2.68) можно переписать в форме нормальных уравнений: (2Л>9)_ где в соответствии с (2.64) использовано четыре зна- чения индекса 1=1, 2, 3, 4. Заметим, что (2.69) имеет ту же структуру, что и (2.42); однако матричное произведение (XjXi) здесь необязательно тёплицево, как в уравнениях Юла — Уокера. Если выбирается матрица данных Хи то нормаль- ные уравнения (2.69) называются ковариационными уравнениями, которые часто используются для ко- дирования речевых сигналов на основе линейного предсказания (см. работу Макхола [44]). Если вы- брана матрица данных Xt, то результирующие нор- мальные уравнения называются автокорреляцион- ными уравнениями, поскольку матричное произве- дение {XfX^IN сводится точно к уравнениям Юла — Уокера, в которых вместо известной автокорреля- ционной функции используется смещенная оценка автокорреляции (2.16). Заметим, что окно данных уменьшает разрешение АР-спектров, оцениваемых с помощью матрицы данных Xt, что будет проиллю- стрировано в разд. III. Если выбрана матрица Xt, то эти нормальные уравнения называются предвзве-
20 ТИИЭР, т. 69, ^11, ноябрь 1981 Оценить автокорреляцию (уравнение (2.16)) I Инициализировать рекурсию (уравнения (2.42), (2.43)) Увеличить ~ I I порядок Рекурсия Левинсона на I (уравнения (2.44)—(2.46)) Вычислить АР-спектр (уравнение (2.38)) Рис. 6. Краткая запись рекурсивного алгоритма Левинсона для спектрального АР-оценивания. шенными нормальными уравнениями, поскольку всем отсчетам, предшествующим х0, приписываются ну- левые значения. Если выбрана матрица %«, то эти уравнения называются поствзвешенными нормаль- ными уравнениями, поскольку нулевые значения приписываются здесь всем отсчетам, следующим после Хдг-1. Может показаться, что только матрица Ха будет давать нормальные уравнения с тёплицевой струк- турой, которая позволяет получить эффективное ре- курсивное решение (а именно, рекурсию Левинсона); см. рис. 6. Однако даже в том случае, когда произ- ведение матриц (X'lXi) может оказаться не тёплице- вым, каждая из четырех матриц данных Xi имеет тёплицеву структуру. Это свойство позволяет разра- батывать рекурсивные алгоритмы с о(р2) операциями в каждом из названных четырех случаев. Морф и др. [163, 164, 167, 169] подробно рассмотрели рекурсивные алгоритмы для случая ковариации; случай/ пред- взвешивания рассмотрен в работах [66, 168]. К ука- занным работам мы и отсылаем читателя, /которого интересуют подробности, касающиеся указанных ал- горитмов. Здесь же мы их подробно рассматривать не будем, так как методы предсказания вперед и назад, которые обсуждаются ниже, позволяют в большинстве случаев получить более точные спектральные оценки. При использовании спектрального АР-оценивания, основанного только на предсказании вперед, возни- кает множество различных проблем. Решение на основе автокорреляционной матрицы (ХаХг) дает спектры с наименьшим разрешением среди отмечен- ных четырех среднеквадратичных АР-оценок в слу- чае, когда последовательность отсчетов данных имеет небольшую длину. Ухудшение разрешения вызвано здесь неявным взвешиванием в матрице данных Ха. Решение на основе ковариационной матрицы (Х^Хг) дает АР-параметры, результирующие спектры которых имеют, как отмечается некоторыми исследователями [227, 228, 251], большее число ложных максимумов и большее частотное смещение спектральных линий по сравнению с другими методами АР-оценивания. Ковариационные нормальные уравнения, которые также используются в методе Прони (рассмотренном в подразд. II.J), приводят к получению оценок АР- параметров с большей чувствительностью к шуму. Расщепление спектральных линий, т. е. появление двух и более близко расположенных максимумов Дам, где должен находиться только один максимум, при- суще всем четырем рассмотренным здесь подходам на основе предсказания вперед. Причины расщепле- ния спектральных линий анализируются в работе [1201. Если исследуемый процесс стационарен в широком смысле, то коэффициенты оптимального фильтра для предсказания ошибки назад будут идентичны коэф- фициентам оптимального фильтра для предсказания ошибки вперед, но будут при этом комплексно-со- пряженными и инвертированными во времени [41]. Использование ошибок предсказания назад для оце- нивания АР-параметров введено Бергом [37, 41]. Наиболее популярный подход для оценивания АР-параметров по отсчетам данных был предложен Бергом в 1967 г. Алгоритм Берга, который отличается от описанного выше алгоритма на основе максималь- ной энтропии, можно рассматривать как минимизацию методом наименьших квадратов с некоторыми нало- женными ограничениями. В случае процесса, стацио- нарного в широком смысле, ошибка линейного пред- сказания вперед (2.63) определяется для p^n^N—1 (где р — порядок устройства предсказания), а ошибка линейного предсказания назад дается выражением р Ьрп= У" аркхп-р+к (2.70) к=о также для р^п^А—1. Напомним, что аРо полагается равным 1, арк — параметры устройства предсказа- ния порядка р, а хп — отсчеты данных. Для того чтобы получить оценки параметров уст- ройства предсказания (или АР-параметров), Берг минимизирует сумму энергий ошибок предсказания вперед и назад: лг-i xv-i = 1ерд1 У lbpnl2> (2.71) л«р п=р наложив при этом ограничение, что АР-параметры должны удовлетворять рекурсивному уравнению Ле- винсона арк = ®p-l,fc + аррар-1,р -к (2.72) для всех порядков от 1 до р. Это ограничение было обусловлено желанием Берга получить устойчивый АР-фильтр (с полюсами внутри единичного круга). Подставляя рекурсивные выражения (2.52) и (2.54) для решетчатого фильтра в (2.71), получим выражение для &р, зависящее только от коэффициентов отра- жения арр и ошибок предсказания порядка р—1, которые полагаются неизвестными. Следовательно, необходимо лишь оценить а при 1=1, 2, . . ., р, если используется (2.72). Приравнивая производную от &р по ац к нулю и решая результирующее уравнение относительно ац, получим ЛГ-1 “ii-Kt- —---------------------- (2.73) X I2 + 1^-1,kl2) k-i Заметим, что |aif|^l, что нетрудно показать с по- мощью (2.73). Таким образом, совместное исполь- зование (2.72) и (2.73) гарантирует устойчивость чистополюсного фильтра. Рекурсивное выражение
СПЕКТРАЛЬНЫЙ АНАЛИЗ 21 для знаменателя в (2.73) получено Андерсоном [10]: JV-1 |2 Знаменатель (О = V (lfy-i, it-il + 1^-!,*! ) = k-t = Знаменатель(i~ 1) [1 — ~ lbf-1.АГ-/12 " l^-i.il2- (2.74) Процедура Берга для спектрального АР-оценивания схематически представлена на рис. 7. Анализ вычис- лительной сложности этой процедуры показывает, что она требует 3Np—рг—2N—р сложений и 3Np— —p*—N+3p умножений комплексных чисел и р де- лений вещественных чисел. Кроме того, необходима память для хранения 3AZ+p+2 комплексных величин [9, 57, 91, 92]. В работах [166, 177, 178, 210, 225] описаны многоканальные варианты алгоритма Берга. Применение алгоритма Берга связано с решением ряда проблем, включая расщепление спектральных линий и частотное смещение спектральных оценок. Если Sp в (2.71) минимизируется по всем aph, k=\, р, то эти проблемы можно ослабить. Ульрих и Клейтон [25 Н и Наттолл [176] независимо пред- ложили процедуру на основе метода наименьших квадратов для предсказания вперед и назад, в ко- торой ограничение рекурсией Левинсона, наложен- ное Бергом, было устранено. Для того чтобы получить р нормальных урав- нений для алгоритма по методу наименьших квад- ратов (называемого также предсказанием вперед и назад), определим минимум величины бр,приравняв ее производные по всем АР-параметрам от ар1 по арр к нулю: Э&» л р а *2 £ WU/) = 0- (2-75) Oapi /-0 Здесь 1=1, .... р, аРв=1 по определению и tv-p-i . , . »₽(«./)= Е Ufc+p-/*k+p-f + jc* +№+/)> <2.76) к-0 где 0<i, /Хр. Минимальную энергию ошибки пред- сказания можно определить как VfWM (2.77) /-о Выражения (2.75) и (2.77) можно объединить и за- писать в форме матричного уравнения Прямое решение уравнения (2.78) методом гауссова исключения требует вычислительных операций, число которых пропорционально величине о(р*). Бэрродейл и Эриксон [16] проанализировали такое решение и его численную устойчивость. Однако выражение (2.79) Инициализация Знаменатель,, - г £0 Вычислить коэффициенты отражении (уравнение (2.73)) Увеличить порядок Рекурсия Левинсона HaJ (уравнения (2.45), (2.40)) .1 Обновить ошибки предсказания ^уравнения (2.52), (2.53))' I Вычислить АР-спектр (уравнение (2.38)) Рис. 7. Краткая запись алгоритма Берга для спектрального АР-оценивания. имеет структуру, которую можно использовать для построения алгоритма с числом операций, пропорцио- нальным величине о(р’). По существу матрицу Rp можно записать в виде сумм и произведений тёпли- цевой и ганкелевой матриц, что позволяет разрабо- тать рекурсивный алгоритм с числом операций, пропорциональным о(р*); см. работу Марпла [154]. В этой же работе приведена структурная схема дан- ного алгоритма, которую мы здесь не приводим из-за ее сложности. В вычислительном отношении алго- ритм по методу наименьших квадратов почти всегда столь же эффективен, как и алгоритм Берга, и в типичном случае требует примерно на 20% больше вычислений. Однако улучшения, достигаемые при использовании этого подхода, перевешивают незна- чительные затраты на дополнительные вычисления. К этим улучшениям относятся меньшее частотное смещение оценок спектральных составляющих, умень- шение дисперсии частотных оценок и отсутствие рас- щепления наблюдаемых спектральных линий [141, 154, 227]. Последовательное оценивание АР-параметров. Име- ется три основных подхода к временному обновлению оценок АР-параметров по отсчетам данных. Это — рекурсивный метод наименьших квадратов, гради- ентный адаптивный подход и последовательная иден- тификация с использованием решетчатого фильтра. Рекурсивный метод наименьших квадратов [14, 20, 211, 220] внешне напоминает процедуру калманов- ской фильтрации. Исключая член из нормальных уравнений (2.69), получим решение по методу наи- меньших квадратов для вектора оценок АР-парамет- ров Ат для порядка р: (2.80) где всем векторам и матрицам добавлен нижний ин- декс т, для того чтобы указать, что включены вре- менные отсчеты вплоть до отсчета с индексом т. Вектор Ут составлен из отсчетов данных
22 ТИИЭР, т. 69, № 11, ноябрь 1981 а матрица Ят представляет собой модифицированный вариант матрицы данных Xi или Хл: которые соответствуют случаям ковариации или пред- взвешивания, обсуждавшимся при рассмотрении ме- тодов групповой обработки. Добавление нового временного отсчета xm+i можно представить в виде следующего разбиения матриц Ут+1 И где Ят+1=1хт- • Xm-p+J- Вводя матрицу Рт= и подставляя разложенные матрицы Гт+1 И ^т+1, получим Ат +1 = +1 +1 +1 = = ^m + l +^m + lxm + l] = ~^т*1 ^т^т^-т^т +^т + 1хт + 1^ • (2.81) Заметим, ЧТО И = ^m+l^Ge + l> тогда 4m + l =^т + 1[(^т + 1 “+Вт + 1хт-ц1 = = Ат *’^т + 1®т + 1С!ст+1 ~ ^т + 1-^т) • (2.82) Используя лемму об обращении матрицы: (4+ +BCD)-1=A~1—А-'Вх (C-1+DA~1B)-1DA-i, полу- чим для Pm+i другую рекурсивную формулировку: Рщ + 1 ~ (^т* +Вт-цВт + 1) = -Рт - Ртн"+1(1 + Hm + 1PmH%<.iY1Hm.1Pm = = (1-Кт^Нт + 1)Рт, (2.83) где, по определению, хт+1 =л»тя"+1(1 +нт+1Ртн»+1)~1- (2.84) Последовательная рекурсия (2.82) сводится в этом случае к Am + i =Ат +^m+l(*m + i “ ^m + iA») — = Ат +Хт.ц ер'Ш*!* (2.85) Уравнения (2.13) — (2.85) по своей структуре ана- логичны калмановскому фильтру; здесь вектор данных Нт и матрица данных Рт аналогичны вектору ко- вариации и матрице, используемым в формулировке калмановского фильтра. В действительности пред- ставленный выше метод наименьших квадратов можно модифицировать, для того чтобы ввести априорные знания об имеющихся статистиках шумовой состав- ляющей ошибки линейного предсказания. Рекурсивный метод наименьших квадратов можно использовать и для решения других задач оценивания параметров, например для спектрального оценивания в реальном масштабе времени (т. е. до тех пор, пока модель остается линейной по параметрам) [170]. Для того чтобы начать рекурсию, необходимо задать 40 и Ро. Если Ро1 отлична от нулевой матрицы, то к выбору А) следует относиться с должной тщатель- ностью, поскольку оценка Ат будет смещаться в направлении к Ал. Это смещение можно устранить, полагая Рё1 нулевой матрицей, но это требует за- писи уравнений (2.84) и (2.85) относительно матрицы Рё1, а не Ро. Обновление соотношений при использовании ре- курсии по методу наименьших квадратов требует о(р3) операций для каждой новой точки данных, что может оказаться неприемлемо большим. Альтерна- тивным подходом, требующим о(р) операций для каждого последовательного обновления, является под- ход, основанный на использовании адаптивного ли- нейного предсказывающего фильтра [8, 83, 85, 243]. При адаптивном подходе вектор АР-параметров оце- нивается рекурсивно с помощью градиентного метода: 4m + 1 = Ат - MVK(|epm|2), (2.86) где р — размер шага, а у — градиент. Подставляя в (2.86) р ерт =хт + У" аркхт -к к-1 и определяя математическое ожидание, получим ^(l^pml1) = ^хх(О) "* AmRxxAm + 2 Re (Amrxx), где Rxx — автокорреляционная матрица, а гхх — авто- корреляционный вектор: гхх= (/?хх(1) . . . /?Хх(р))т- Градиент от Е(|ертР) равен V£(|epm I2) = 2rxx + 2RXX Ат. (2.87) Если градиентный метод сходится, то vE(|epml)s)=0 и ГХХ = ~&ХхА , что аналогично уравнениям Юла — Уокера (2.40). На практике величины гхх и Rxx не известны, поэтому используются их оценки, что дает Вхх _ tx£ _ ]. Заметим, что ерт=хт+Хт-1Ат; поэтому, используя далее (2.87), уравнение (2.86) можно записать в следующем виде: Аи + i = Ат - 2liepmXm -!. (2.88) Уравнение (2.88) составляет основу алгоритма оценки по методу наименьших средних квадратов (МНСК). Оно имеет структуру, аналогичную уравнению (2.85), но с тем лишь отличием, что р здесь фиксировано, тогда как в (2.85) коэффициент усиления Km+i яв- ляется переменной величиной. Сходимость данного алгоритма гарантируется до тех пор, пока константа р лежит в интервале от 0 до 1/К-»., где к.., — максимальное собственное значение матрицы Rxx. Выбор р требует принятия некоторой) компромисса между скоростью сходи- мости Е [4т1 к 4m и величиной дисперсии в устано-
СПЕКТРАЛЬНЫЙ АНАЛИЗ 23 вившемся режиме (иногда называемой рассогласо- ванием), допускающей достижение сходимости. Сле- довательно, цена, которую приходится платить за значительное уменьшение объема вычислений по сравнению с рекурсивным методом наименьших квад- ратов,— это более медленная сходимость, что тре- бует применения более длинных записей для дости- жения надежных оценок АР-параметров. С целью достичь точности рекурсивного метода наименьших квадратов и вычислительной эффек- тивности адаптивных градиентных методов был пред- ложен третий метод, основанный на использовании решетчатого фильтра. Однако решетчатое рекурсивное соотношение обновляет лишь коэффициенты отраже- ния, для чего требуется о(р) операций. Если оценки АР-параметров должны обновляться при поступлении каждого нового отсчета данных, то необходимо ис- пользовать рекурсию Левинсона (2.45), которая тре- бует о(р’) операций. Следовательно, вычислительная экономия при использовании метода на основе ре- шетчатого фильтра достижима только в том случае, когда оценки АР-параметров обновляются не по каж- дому новому отсчету, а гораздо реже. Алгоритм последовательной оценки, основанный на решетчатых структурах, можно построить для каждого из методов наименьших квадратов, рассмот- ренных при обсуждении группового оценивания [164, 1681. Одной из иллюстраций такого алгоритма, ос- нованного на алгоритме Берга, может служить опи- санный ниже алгоритм. Коэффициент отражения Ктп порядка т и временного индекса (момента вре- мени) т в алгоритме Берга определяется с помощью следующего выражения: Л "2 У? ет Xm„ = —----—-------------------. (2.89) Е (Ibm-U-il2 + l«m-i,d2) i-m Выражение для рекурсивного обновления этого коэф- фициента дается выражением [22] Кт +1, п *1 = Кт +1 ,п ~ — *1',|(l^nnl2 + lfyn,n -112) + 2^ПИЬД|,И -11 Е [l^nil2 + l^m,*-l 121 Z-m Таким образом, объединяя (2.90) с (2.52) и (2.54) при т=1, ..., р и начальном условии e0n=Z>0n=xn, получим алгоритм последовательного временного об- новления коэффициентов отражения. АР-оценивание спектральной мощности. В работе [136] показано, что в отличие от обычных спектраль- ных оценок на основе преобразования Фурье ампли- туды спектральных линий в спектральных АР-оценках нелинейно зависят от мощности в том случае, когда входной процесс представляет собой смесь синусоид и шума. Лакосс [136] показал, что при больших отношениях сигнал/шум амплитуда спектральной ли- нии пропорциональна квадрату мощности, хотя пло- щадь под спектральной линией пропорциональна мощности. Один из методов получения оценки дей- ствительной мощности реальных синусоид в АР- спектре предложен Андерсеном и Джонсеном [108]. Этот метод хорошо работает при больших отношениях сигнал/шум для отдельных компонент процесса. В тер- минах z-преобразования АР-оценка СПМ дается вы- ражением РАр(г)~ MAW)' <291> где Л(г)= 1 + У apkz~k. k-i Если какому-либо максимуму в спектральной АР- оценке соответствует частота ft, т. е. г(=ехр(/2л/Д/), то величина оцениваемой мощности будет приближенно равна: мощность (/,) = 2 X Re I остаток от на z,| I 2 J (2.92) где „ ^Ар(2) . J^Ap(2i)| Остаток от -----= (z - zA —— 2 I 2 Ь-ь/ Заметим, что при этом полагается, что максимумы соответствуют угловым положениям (координате) по- люсов zt. При использовании данного метода возможно появление отрицательных оценок мощности, если спектральные линии расположены очень близко друг к другу. Данный подход тесно связан с процедурой оценивания мощности, используемой в методе Прони, рассмотренном выше в подразд. II.J. Выбор порядка модели [21, 100, ПО, 114, 133, 235, 237, 238]. Поскольку наилучший выбор порядка модели заранее, как правило, не известен, то на практике необходимо использовать несколько поряд- ков модели. В этом случае рассчитывается какой- либо критерий ошибки, по которому и определяется требуемый порядок модели. При слишком низком порядке модели получаются сильно сглаженные спект- ральные оценки, при слишком высоком — в спектре появляются ложные детали. Интуитивно ясно, что следует увеличивать порядок АР-модели до тех пор, пока вычисляемая ошибка предсказания не достиг- нет минимума. Однако во всех обсуждаемых в статье процедурах оценивания по методу наименьших квад- ратов мощности ошибок предсказания монотонно уменьшаются с увеличением порядка модели р. Так, например, в алгоритме Берга и в уравнениях Юла — Уокера используется соотношение S/=S/_ i[ 1 - |в«|2]. (2.93)( До тех пор пока величина |а„I* отлична от нуля (она должна быть равной или меньше 1), мощность ошибки предсказания уменьшается. Следовательно, сама по себе мощность ошибки предсказания не может служить достаточным критерием окончания процедуры оце- нивания. Для выбора порядка АР-модели предложено не- сколько целевых критериев. Акаике [1—7) предложил два критерия. Первым критерием является величина оконечной ошибки предсказания (ООП). Согласно
24 ТИИЭР, т. 69, № 11, ноябрь 1981 этому критерию, порядок АР-процесса выбирается таким образом, чтобы средняя ошибка на одном шаге предсказания была минимальна. Акаике рассматривал ошибку как сумму мощностей в непредсказуемой (или необновляемой) части процесса и как некоторую величину, характеризующую неточность оценивания АР-параметров. Оконечная ошибка предсказания для АР-процесса определяется как с (N + P+ 1\ °°nP=Sp(^J-r;’ (2.94) где N — число отсчетов данных. Заметим, что в (2.94) предполагается использование отсчетов, из которых вычтено выборочное среднее. Член в круглых скобках увеличивает оконечную ошибку предсказа- ния, по мере того как р приближается к N, характе- ризуя тем самым увеличение неопределенности'оценки &р для мощности ошибки предсказания. Выбираемый порядок р равен порядку, при котором величина оконечной ошибки/ предсказания минимальна. Кри- терий на основе оконечной ошибки предсказания исследовался в различных приложениях Гершом и Шарли [73], Джонсом [111], Фрайером и ,,др. [69] и Ульрихом и Бишопом [250]. Для АР-процессов он обеспечивает превосходные результаты. Однако при обработке реальных геофизических данных этот кри- терий, как показали Джонс [111] и Берримэн [21], приводит к выбору слишком малого порядка модели. Второй критерий Акаике основан на методе мак-, симального правдоподобия и получил название ин- формационного критерия Акаике (ИКА). Согласно этому критерию, порядок модели выбирается посред- ством минимизации некоторой теоретико-информаци- онной функции. Если исследуемый процесс имеет гауссовы статистики, то ИКМ определяется выра- жением ИКАр = In (Sp) + 2(р + l)/N. (2.95) Член (р+1) в (2.95) иногда заменяется на р, поскольку величина 2/N — это всего лишь аддитивная кон- станта, которая учитывает вычитание выборочного среднего значения. Второй член в (2.95) характеризует плату за использование дополнительных АР-коэф- фициентов, но это не приводит к значительному уменьшению мощности ошибки предсказания. И здесь выбирается тот порядок модели, при котором ИКА минимален. При N -► оо первый и второй критерии эквивалентны. Кашьяп [1151 утверждает, что ИКА является статистически несостоятельным критерием в том смысле, что вероятность ошибки при выборе правильного порядка модели не стремится к нулю при IV-»- оо. Третий метод выбора критерия предложен Пар- зеном [192] и получил название авторегрессивной передаточной функции критерия (АПФК). Порядок модели р выбирается в этом случае равным порядку, при котором оценка разности среднеквадратичных ошибок между истинным предсказывающим ошибку фильтррм (его длина может быть бесконечной) и оценкой этого фильтра минимальна. Парзен показал, что эту разность можно вычислить, даже если истин- ный предсказывающий ошибку фильтр точно не из- вестен: апфк’-& МН' <2W где &j = (N/(N—И здесь снова р выбирается так, чтобы минимизировать АПФК. Результаты оценивания спектра при использова- нии трех указанных критериев мало отличаются друг от друга, особенно в случае реальных данных, а не моделируемых АР-процессов. Ульрих и Клейтон [251] показали, что в случае коротких записей дан- ных ни один из этих критериев не обеспечивает удов- летворительных результатов. Для гармонических про- цессов в присутствии шума использование первого и второго критериев приводит к заниженной оценке порядка модели, если отношение сигнал/шум велико [92, 138]. При анализе коротких сегментов данных Ульрих и Уи [92] предложили выбирать порядок модели равным от JV/3 до N/2, что во многих случаях приводит к удовлетворительным результатам. За- метим, что выбор порядка модели для данных, полу- ченных из реальных процессов, пока еще требует учета большего числа субъективных факторов, чем для процессов, моделируемых и управляемых с по- мощью ЭВМ. Аномалии и коррекция спектральных АР-оценок. Рядом исследователей отмечается несколько аномалий спектральных АР-оценок. Так, если выбранный по- рядок модели слишком велик по сравнению с числом отсчетов данных, то в спектральных АР-оценках появляются ложные максимумы [212, 250]. В идеаль- ном случае, если значения автокорреляционной функ- ции или, что эквивалентно, коэффициенты отражения были оценены без ошибки, то оценки АР-параметров для АР-модели порядка р будут равны л ррь i=l,2, •••,₽; “pi = I (2.97) 10, i = р + 1, • • •,п, где api — параметры АР-модели порядка р. Однако при наличии ошибки оценивания api=/=0 при i>p. Соответственно появятся п—р «лишних» полюсов. Если оцениваемые значения этих полюсов располо- жены вблизи единичного круга, появятся ложные спектральные максимумы. Именно возможность по- явления ложных максимумов и учитывается в реко- мендации, что максимальный порядок модели не должен превышать величину N/2, где N — длина (число отсчетов) записи данных [251]. В работах [45, 229] отмечается, что при анализе процессов, состоящих из смеси синусоиды с шумом, положение максимума в спектральной АР-оценке сильно зависит от фазы этой синусоиды. Кроме того, в ряде случаев отмечено, что спектральная оценка может иметь два близко расположенных максимума, тем самым ложно указывая на присутствие второй синусоиды. Это явление называется расщеплением спектральных линий (РСЛ) [64]. Зависимость спектральных АР-оценок от фазы уменьшается с увеличением длины записи данных. В количественном отношении эта зависимость раз- лична для различных процедур АР-оценивания. Для алгоритма Берга сдвиг пика максимума достигает 16% [230]. Метод предсказания ошибки вперед и назад имеет наименьшую чувствительность к фазе [154, 251]. Для ослабления этого эффекта предложено два метода. В первом из них зависимость от фазы приписывается взаимодействию между положитель- ными и отрицательными частотными составляющими
СПЕКТРАЛЬНЫЙ АНАЛИЗ 25 вещественной синусоиды, что во многом напоминает зависимость максимумов периодограммы от фазы (119, 231, 234]. Следовательно, вещественнозначный сигнал должен заменяться аналитическим сигналом. Аналитический сигнальный процесс дискретизуется с вдвое меньшей частотой, и, следовательно, в этом случае используется спектральная АР-оценка для комплексных данных. Порядок модели для комплекс- ных данных составляет лишь половину от порядка модели для вещественных данных, поскольку пары сопряженных комплексных полюсов, используемые в случае вещественного сигнала, не обязательны при комплексном подходе. С помощью такого подхода можно ослабить зависимость спектральной оценки Берга от фазы [119]. Еще одна процедура (для спектральной оценки Берга) состоит в использовании оценки А IN-1 V’2 Z ^-^пЬГ-л.п-г / Е + n=i / n»i + (2.98) в которой члены, соответствующие коэффициентам отражения, взвешиваются в соответствии с вещест- венной последовательностью {оп}. Такое взвешивание остаточного временного ряда, предложенное Суинг- лером [228], ослабляет концевые эффекты при анализе коротких записей данных. Результаты моделирования показали, что этот метод ослабляет также и фазовую зависимость. Явление расщепления спектральных линий в АР- оценках спектра при использовании алгоритма Берга впервые подробно рассматривалось в работе [64]. Авторы этой работы отмечали, что расщепление спект- ральных линий наиболее вероятно в тех случаях, когда: 1) отношение сигнал/шум велико, 2) начальная фаза синусоидальных составляющих в некоторое не- четное число раз кратна углу 45°, 3) временная дли- тельность последовательности данных такова, что синусоидальные составляющие имеют нечетное число четвертей периодов, и 4) оцениваемое число АР-пара- метров сравнимо с числом значений данных, исполь- зуемых для оценивания. Спектры, в которых наблю- дается расщепление спектральных линий, могут со- держать большое число ложных линий (максимумов). И снова это явление обусловлено малой длиной записи данных, поскольку число ложных максимумов умень- шается с ростом длины записи данных. Расщепление спектральных линий наблюдается при использовании алгоритмов Берга и Юла — Уокера в случае про- цесса, состоящего из нескольких синусоид и белого шума [120]. Решение этой проблемы предложено Фу- жером [65], который связал расщепление спектраль- ных линий с .тем фактом, что мощность ошибки пред- сказания при использовании оценки Берга для коэф- фициентов отражения на самом деле не минимизи- руется. В его методе модность ошибки предсказания минимизируется посредством одновременного изме- нения всех коэффициентов отражения. Судя по ре- зультатам выполненного им моделирования, этот метод может устранить расщепление спектральных линий по крайней мере в случае одной синусоиды. Еще один метод устранения расщепления спектраль- ных линий для случая одной синусоиды основан на использовании аналитического сигнала [120]. Этот подход дает истинные значения коэффициентов отра- жения для спектральной оценки Берга в случае, когда шум пренебрежимо мал, т. е. при условии, когда расщепление спектральных линий наиболее вероятно. Для спектральной оценки Юла — Уокера при ис- пользовании подхода на основе аналитического сиг- нала и несмещенной оценки автокорреляции полу- чается истинная оценка автокорреляционной функции. Поэтому расщепление спектральных линий в любом случае устраняется. В случае многих синусоид харак- теристики алгоритма Фужера не изучены, а следова- тельно, использование комплексных данных сов- местно с оценкой Берга для коэффициентов отражения все еще может сопровождаться расщеплением спект- ральных линий [94]. При использовании метода ли- нейного предсказания вперед и назад совместно с рекурсивным алгоритмом Марпла расщепление спект- ральных линий не наблюдалось [154]. Помимо фазы АР-оценки спектров могут также искажаться из-за наличия большой постоянной со- ставляющей сигнала или линейного тренда [113], особенно на узкочастотном конце спектральной оценки. Эти составляющие сигнала необходимо устранять до применения методов спектрального АР-анализа. Весьма важная проблема, касающаяся спектраль- ной АР-оценки и ограничивающая ее использование, обусловлена чувствительностью этой оценки к шумам наблюдения, присутствующим во временных рядах [186]. Пример показан на рис. 8. Из этого рисунка нетрудно видеть, что спектральные линии расширя- ются и смещаются со своих правильных позиций Доли частоты дискретизации (Ь) Рис. 8. Спектральные оценки двух синусоид в белом гауссовом шуме при большом (а) и малом (Ь) отношениях сигнал/ /шум. 4 тииэр м 11
26 ТИИЭР, т, 69, № 11, ноябрь 1981 (указаны стрелками). В частности,было показано, что разрешение спектральной АР-оценки двух синусоид равной амплитуды в присутствии белого шума умень- шается с увеличением отношения сигнал/шум [149, 152]. При низких отношениях сигнал/шум разрешение не превышает разрешения периодограммы. Причина этого ухудшения состоит в том, что чистополюсная модель, используемая при спектральном АР-анализе, не обеспечивает правильных результатов при наличии шума наблюдения. Чтобы показать это, предположим, что уп обозначает искаженный шумом АР-процесс хп. Тогда yn=xn + w„, (2.99) где wn — шум наблюдения. Если wn — белый шум с дисперсией oj,, некоррелированный с процессом хп, то Ру(х) = о2 At 2 = [о2+р^,Л(г)Л*(1/г*)] Дг 10()) Л(г)Л*(1/г*) Таким образом, СПМ процесса уп характеризуется нулями и полюсами, т. е. уп есть АРСС (р, р)-процесс. Несостоятельность АР-модели для искаженного шу- мом АР-процесса приводит к искажениям, показанным на рис. 8 [121]. Это явление можно объяснить следую- щим образом. Влияние шума должно уменьшить динамический диапазон СПМ процесса хп. Поскольку предсказывающий ошибку фильтр А (г) «пытается» обелить СПМ, то неудивительно, что при низких отношениях сигнал/шум нули фильтра А (г) распола- гаются вблизи начала координат г-плоскости, т. е. Л(г)«1. Это обусловлено тем, что СПМ процесса уп уже относительно уплощена из-за наличия шума, поэтому последующая операция фильтрации, т. е. использование предсказывающего ошибку фильтра, лишь незначительно обеляет СПМ. Для ослабления степени ухудшения спектральной АР-оценки в присутствии шума предложено четыре общих подхода. Для этой цели можно: 1) использовать спектральную АРСС-оценку, 2) отфильтровать данные для снижения уровня шума, 3) использовать АР-модель высокого порядка, 4) компенсировать влияние шума либо в оценках автокорреляционной функции, либо в оценках коэффициентов отражения. АРСС-подход предполагает, что искаженный шу- мом АР-процесс является обобщенным АРСС(р, р)- процессом, даже если АР- и СС-параметры связаны соотношением (2.100). Наиболее часто применяемый подход основан на использовании модифицированных уравнений Юла — Уокера, как описано в подразд. II.F [58, 72]. Для АРСС(р, р)-процесса это означает решение системы уравнений р Ryy(k) — ~ У ajRyyik- I) (2.101) 1ч при А>=р+1, р+2, . . ., 2р, для того чтобы получить АР-параметры. Хотя этот метод и прост в реализации, проблемы полностью он все же не решает. Прием- лемые результаты получаются при использовании длинных записей данных и (или) больших отношений сигнал/шум. Более подходящее решение проблемы шума состоит в использовании АРСС-оценки по методу максимального правдоподобия. Однако ис- пользование этой процедуры приводит к системе нелинейных уравнений [31]. В случае, когда урав- нения максимального правдоподобия определяются специально для АР-процесса в присутствии белого шума, субоптимальное решение этих уравнений при- водит к итеративной схеме фильтрации [143]. Еще одна схема АРСС-фильтрации рассмотрена в работах [118, 151, 159, 172]. Все эти методы основаны на ите- ративных подходах к расчету фильтра, поскольку спектральная плотность мощности процесса хп, кото- рую мы пытаемся оценить, не известна. Второй метод ослабления влияния шума основан на использовании АР-модели, порядок которой больше порядка истинной АР-модели. Объясняется это тем, что АРСС (р, р)-процесс эквивалентен АР-модели бесконечного порядка, что гарантируется декомпози- цией Уолда. Используя (2.100), положим а2 + a2wA(z)A*(l/z*) = a25(z)fi*(l/z*), где 2?(z)=l+ X bpkz~k, кч поэтому процесс уп можно представить как выход фильтра из общего вида (т. е. с нулями и полюсами) Н(г), на вход которого подан белый шум с диспер- сией о?) (рис. 9). Поделив А (г) на В (г), получим ff(z)=l/(A(z)/B(z))=l/C(z), где A(z) , _ь Таким образом, у„ можно моделировать с помощью АР(оо)-процесса с параметрами {cft}. Ясно, что когда порядок предполагаемой АР-модели возрастает, АР- оценка СПМ будет стремиться к истинному значению СПМ процесса уп- Это свойство вытекает также из формулировки метода максимальной энтропии, по- скольку показано, что ^(fc) = Ry у (fc), I к | С р, где Ryy (й) — автокорреляционная функция, соответ- ствующая АР (р)-модели, а /?„„(£) —истинная ав- токорреляционная функция процесса Таким об- разом, в пределе при р -> «> автокорреляционная функция модели совпадает с истинной автокорреля- ционной функцией, а следовательно, совпадают и спектры. o^(z)B(r? - o'. apktz) А(гЧ Рис. 9. АРСС-модель для АР-процесса в белом шуме.
СПЕКТРАЛЬНЫЙ АНАЛИЗ 27 (а) в тех случаях, когда нули функции В (г) расположены вблизи единичной окружности в г-плоскости. В таких случаях последовательность с* медленно сходится. Поскольку с уменьшением отношения сигнал/шум ну- ли функции В (г) приближаются к границе единичного круга [121], увеличивать порядок модели необходимо при малых отношениях сигнал/шум. Для количест- венной оценки влияния порядка модели на спектраль- ную АР-оценку АР-процесса в присутствии белого шума рассмотрим смесь из двух синусоид равной ам- плитуды и белого шума. Для этого случая в работе [149] показано, что разрешение 6f (в Гц) спектральной АР-оценки в предположении известной автокорреля- ционной функции приближенно равно [р(р + 1)10’31 6,471 X 2ярДг ’ (2.102) -40,0 -50,0 -во,0 -70/) t,At = O, 143 »» At =0,200 N = 100 С/Ш=0 дБ р=32 где р — отношение сигнал/шум. Как и ожидалось, разрешение увеличивается с ростом порядка модели, иллюстрацией чего может служить рис. 10. Заметим, что дополнительные полюсы примерно равномерно распределены внутри единичного круга, что приводит к равноволновой аппроксимации плоского шумового спектра. Предложено большое число схем компенсации влияния шума на значения автокорреляционной функ- ции [123, 151, 159, 214, 240, 268]. Подробности можно найти в указанных работах. Оценка СПМ является частным случаем таких схем. В общем случае схемы компенсации шума могут уменьшить смещение оце- нок, но лишь за счет увеличения дисперсии этих оце- нок. Серьезный недостаток всех этих методов состоит в том, что, вообще говоря, не известно, какую мощ- ность шума необходимо устранить. Следовательно, если дисперсия шума <т£, слишком велика, АР-оценка спектра будет иметь более острые максимумы, чем истинный спектр. Поэтому необходимо быть очень внимательным при использовании этих методов. -во,О--------——-------------------------- ’ 0,0 0,1 0,2 0,3 0,4 0,5 Доли частоты дискретизации Рас. 10. Спектральная оценка Берга для двух синусоид в белом шуме. Показано влияние порядка модели: (а) р=4; (Ь) расположение полюсов при р=4; (с) р— 32; (d) рас- положение полюсов при р=32. Может показаться, что следует выбирать по воз- можности наибольший порядок модели. Однако из-за возможности появления ложных максимумов порядок модели следует ограничивать таким образом, чтобы он не превышал половины числа отсчетов данных, о чем говорилось выше. На практике модель большого порядка необходима F. Оценивание СПМ на основе скользящего среднего Как отмечалось во введении к разд. II, процесс скользящего среднего представляет собой случайный процесс, получаемый на выходе фильтра, передаточ- ная функция которого содержит только одни нули и на выход которого подается белый шумовой процесс: я хп= £ ЬтПп-«1, (2.103) т «о причем £[ля]=0, Е[пп+тп*] = а28т, где 6т=1 при т=0 и 6т=0 в остальных случаях. С помощью (2.103) автокорреляционная функция Со- процесса порядка q определяется следующим обра- зом: С я-к I /«О о, * = 0,1, • • •,<?; (2.104) k>q. 4»
28 ТИИЭР, т. 69, № И, ноябрь 1981 Таким образом, если известны q+1 значений автокор- реляционной функции, то параметры СС-процесса порядка q определяются из решения системы нели- нейных уравнений (2.104)„для чего часто использует- ся так называемый метод моментов [31]. Однако если необходимо определить только спектральную оценку, то решать эти уравнения для определения параметров СС-процесса не обязательно и достаточно определить только автокорреляционную функцию, поскольку ч Е •^хх(т) ехР (-/2а/тДт), т -~q что идентично спектральной оценке Блэкмана — Тью- ки. В этом случае метод моментов для решения задачи спектрального оценивания не применим. Если для параметров СС-процесса мы будем использовать оцен- ку по методу наименьших квадратов, то она будет так- же соответствовать аналогичной оценке автокорреля- ционной функции, поскольку с помощью (2.104) задается однозначное преобразование (предполагает- ся, что нули расположены внутри единичного круга в z-плоскости). В данном случае параметры СС-процес- са удобно рассматривать как промежуточный этап к оцениванию спектра. Однако такой подход не ис- пользуется, так как оценка параметров СС-процесса по методу наименьших квадратов в сильной степени нелинейна [31]. Кроме того, для моделирования СС- процесса необходимо использовать слишком много коэффициентов даже в случае узкополосных спектров, что приводит к получению плохих спектральных оценок в подобных случаях. Когда для анализа до- ступны лишь отсчеты данных, необходимо определять порядок СС-модели. Интуитивный подход, предложен- ный Чжоу [48], основан на использовании несмещен- ных оценок значений автокорреляционной функции (2.14) и проверке того, что эти оценки быстро дости- гают нуля при небольшом числе членов, поскольку из (2.104) мы знаем, что Rxx(m)=0 для т, больших порядка СС-процесса. В противном случае более удоб- ной может оказаться АР- или АРСС-модель. Чжоу предложил также проверять последовательные значе- ния автокорреляционной функции, с тем чтобы опре- делить, достаточно ли близко значение Rxx(q) к нулю с учетом дисперсии значений с индексом, меньшим q. Если такая гипотеза подтверждается, то порядок СС- процесса полагается равным q. Оценки значений авто- корреляционной функции используются в (2.104) для определения параметров СС-процесса. Если порядок q определен, можно осуществить дальнейшее уточнение оценок СС-параметров [236], налагая на оценки ограничение Rxx(m)—Q при m>q. G. Оценивание СПМ на основе авторегрессии и скользящего среднего [И, 31, 35, 72, 79, 86, 87, 89, 117, 165, 242, 266, 267] Уравнения Юла—Уокера. Напомним, что в АРСС- модели предполагается, что временнбй ряд хп можно моделировать как выход фильтра с р полюсами и q нулями, на вход которого подается белый шум: р я хп = ~ Е akxn-k + Ё ^fcnn-k, (2.105) k-t k-o где Rnn(k)=at6k и 50=l- Предполагается также, что полюсы этого фильтра расположены внутри единич- ного круга в z-плоскости. Нули же могут распола- гаться в любом месте этой плоскости. Если параметры АРСС(р, ?)-модели идентифициро- ваны, то спектральная оценка получается с помощью следующего выражения: ?,(/)= |Я(ехр [/2irfAt] |2 ?„(/) = 1 + £ Ьк ехр (->/2а/*Дг) fc-i р 1 + ак ехР (~i2iifk&t) (2.106) Соотношение между АРСС-параметрами и автокор- реляционной функцией нетрудно определить следую- щим образом. Умножая (2.105) на и беря матема- тическое ожидание, получим *«(0 = - Е *)+ Е bkRm{l- к), (2.107) fc-l fc-o где /{„(*) »£(п„*„*_*). Но /?пх(^)=0 при й>0, поскольку будущий вход кау- зального устойчивого фильтра не может влиять на текущий выход, ал, — белый шум. Поэтому С р я Е акРххО к) + Ё bkRnxU ~ к), Л-1 к-О - Е “кЯхЛ-Ъ, к-1 / = <г + 1,<? + 2, •••; При выводе уравнений Юла — Уокера было показа- но, что /Jnx(fc) = a2h!fc, и, следовательно, - Е а^«(/“ *) + <72 £ МГь к-t к-1 *хх(0=< 1 = 0, 1,---,<?; (2.108) - Е “kRxxU- *). fc-i 1 = 9 + 1, q + 2, • • •. Эти нормальные уравнения для АРСС-процесса ана- логичны уравнениям Юла — Уокера для АР-процесСа. Оценивание параметров АРСС-процесса. Предло- жено большое число методов оценивания АРСС-пара- метров, однако все они, как правило, требуют значи- тельных матричных вычислений и (или) итеративных методов оптимизации. Такие методы не пригодны для обработки в реальном масштабе времени. Для того чтобы сделать вычислительную нагрузку более при- емлемой, были разработаны субоптимальные методы. В них обычно используется критерий среднеквадра-
СПЕКТРАЛЬНЫЙ АНАЛИЗ 29 тачной ошибки, й они требуют решения линейных уравнений. Однако эти методы обеспечивают раз- дельную, а не совместную оценку АР- и СС-парамет- ров, которая требуется для оптимального оценивания параметров. Заметим, что АР-параметры можно оце- нить независимо от СС-параметров, если уравнения Юла — Уокера используются в форме (2.108). Основ- ной довод в пользу субоптимальных линейных под- ходов состоит в том, что итеративные методы оптими- зации не гарантируют сходимость или же сходятся к неверным решениям. Примером нелинейных уравне- ний, с которыми приходится иметь дело, могут слу- жить уравнения (2.108). Поскольку импульсный от- клик фильтра является функцией коэффициентов аъ . . ., ар, bi, . . ., bq, уравнения (2.108) нелинейны относительно АРСС-параметров. В качестве примера рассмотрим АРСС(1, 1)-процесс. В данном случае hk = (-Я1 )ки(к) + bi(.-ai)k~1u(k- 1), Кхх(0) = -Я1Кхх(-1) + Р2(1 + Ibi I2 - a*bi), лхх(1) = -oi-R^O) + a2 bi, RxxU) = ~alRxx(l ~ 1) при 1>2, (2.109) где и (k) — единичная ступенчатая функция. Для решения этих уравнений предложено большое число методов, однако известно всего несколько случаев их успешного применения. Более распространенный подход к этой проблеме состоит в использовании уравнений (2.108) при t>q для определения коэффициентов (о», at, . . ., ар), после чего с помощью какого-либо подходящего мето- да определяются коэффициенты (bi, Ьг, . . ., bq) или некоторый эквивалентный набор параметров. Так, например, для определения АР-параметров с помощью (2.108) при l=q+\, <7+2, .... q+p необходимо ре- шить следующее матричное уравнение: Лхх(я) Rxx(Q “О * * ’ “ Р 1) Rxxtq) ••• ^(<7-р + 2) E^fa+p-l) Rxx(q+p~2) ••• Rxx(q) ч___________________У---------------------- Rxx Rxx(.4 + 1) Яхх(« + 2) (2.110) Rxx(.4 + Р) Уравнения (2.110) получили название обобщенных, или модифицированных, уравнений Юла — Уокера. Матрица здесь тёплицева, хотя и несимметричная, а потому не гарантируется, что она будет либо поло- жительно-определенной, либо невырожденной. Для ранения уравнений (2.110) Зохар [278] разработал ал- горитм, требующий о(р*) операций. Для того чтобы выбрать подходящий порядок мо- дели р для авторегрессивной части АРСС-модели, можно использовать свойство [48] |Л»1=0, если размерность матрицы R^ больше, чем порядок АР-модели р. Здесь |/?^| — определитель матрицы R'xx. Указанное свойство означает, что необходимо лишь проверять определитель при i= 1, 2, ... до тех пор, пока он не станет достаточно малым. После на- хождения оценок АР-параметров {ак} параметры СС-процесса можно определить посредством фильтра- ции данных с помощью фильтра А (г), обладающего только одними нулями, где l(z) = 1 + £ fa~k, к-1 в результате чего получается «чистый» СС-процесс. Выполнив эту операцию, можно далее использовать методы, описанные в подразд. II. F. Для определения АР-параметров требуется применить метод спектраль- ной факторизации. Заметим, что для того чтобы устра- нить спектральную факторизацию, для спектрального оценивания можно ограничиться отысканием A(z)x xA*(l/z*) и B(z) B*(l/z*), поскольку спектральная оценка равна [117, 129] УХ(Г) = g24rg(z)J*(l/z*) A(z)A*(l/z*) ж = ехр(/2я/Д»). Если Bk = Z~l [a2415(z)B*(l/z*)], где Z и Z-1 — соответственно прямое и обратное г-пре- образования, то спектральная оценка будет равна У. Якехр(-/2хДДг) х |Л(ехр [/2я/Дг])|2 где В-к=Вк. Для того чтобы получить Вк, заметим, что fifc=Z-‘[A(z)A*(l/z*)Px(z)]. Полагая Afc=Z-1[A(z)A*(l/z*)], получим р Вк = Е AnRxxtt-n), fc = 0, 1,•••,<?• п = -р Для того чтобы определить Вк, необходимо знать {Ях«(0)}, Яхх(1), . . ., 7?хх(р+?)- Для того чтобы га- рантировать неотрицательность спектральных оце- нок, необходимо, чтобы коэффициенты Вк составляли положительно-полуопределенную последовательность. Характеристики модифицированного подхода Юла — Уокера значительно изменяются при его ис- пользовании для моделирования АРСС-процессов. Для некоторых процессов получаемые оценки АРСС-пара- метров будут достаточно точны. Однако для ряда дру- гих процессов их точность будет недостаточной. Рас- смотрим, например, асимптотическую дисперсию (при N -* оо) оценки АР-параметров для реального АРСС(1, 1)-процесса. Оценка, получаемая из решения (2.110), имеет вид £1 =-^(2)1^1). Можно показать [72], что ее дисперсия равна Var(^4 (2.ш) ™ ^ХХ^)
30 ТИИЭР, т. 69. № 11, ноябрь 1981 Но *хх(0) = а2 Г--** - *bi L 1 - “г А«(1) = ~ в1^жх(0) + a2 bi, поэтому (2.Ш) можно переписать в виде 1 Var(ai) = j^ 2b? (1-а?) _ («1 -А)2^]2 0-а,2) /.. Для получения приемлемо точной оценки at можно, например, потребовать, чтобы среднеквадратичная ошибка не превышала 0,1: VVar(aj) < 0,1. Для того чтобы определить условия, удовлетворяющие этому требованию при bt=0,5, была вычислена зави- симость минимального числа отсчетов от величины at, показанная на рис. 11. Нетрудно видеть, что при за- данном значении N статистическая флюктуация будет сильно зависеть от формы спектральной оценки. При ЛГ->- оо АРСС(1, 1)-модель оказывается, по-видимо- му, не пригодной при ai -> 0,5, поскольку это приво- дит к компенсации нулей и полюсов в белом шумовом процессе. Этот пример иллюстрирует тот факт, что следует очень осторожно применять модифицирован- ный подход Юла — Уокера, особенно в тех случаях, когда порядок модели р не известен. Обусловлено это тем, что модифицированная матрица Юла — Уокера сформированная из точно известных значений автокорреляционной функции, будет вырожденной, если ее размерность превышает порядок точной мо- дели. Второй метод оценивания АРСС-параметров осно- ван на использовании равенства Д(г) 1 A(z)~C(z)’ где С(2)=1+ £ Cfc2-*, it-i для приравнивания АРСС-модели АР-модели беско- нечного порядка. Коэффициенты {ск} можно оценить, Ряс. 11. Требуемое число отсчетов в записи данных для точного оценивания АР-параметров АРСС (1,1)-процесса с по- мощью уравнений Юла — Уокера. используя только АР-методы и затем связывая полу- ченные результаты с АРСС-параметрами. В частности, пусть С (z)—1 +2 — оцениваемые АР-пара- метры, где M=p+q. Полагая p>q, получим (791 B(z) 1 Аг) C(z) ’ я л д-о л = 1,2, •••, где ё0=1. Поскольку при п>р значение ап должно полагаться равным нулю, то я л 22 Ьк<т-к = 0> п = р+1,р + 2,- ,р + я. к-0 Это выражение можно переписать в матричной форме: После определения {bu bit .... bq} коэффициенты {а1( а,, . . ., ар} можно найти из решения следующего уравнения: я л £ bkbi-k=4n, for л = 1, 2, • • • ,р, (2.113) к-о где ?о = 1. В матричной форме это уравнение имеет следующий вид: Поскольку C(z)=A (z)/B(z), АР-модель очень боль- шого порядка должна использоваться в том случае, когда нули функции В (г) расположены вблизи единич- ной окружности. В этом случае последовательность с» будет быстро затухать. Данный случай представляет практический интерес в связи с тем фактом, что, если нули функции В (г) расположены вблизи начала коор- динат, они мало влияют на оценки СПМ. В этом случае лк бая АР-модель является приемлемой. Тем не менее с помощью этого метода получен ряд многообещающих результатов [79, 197J. Третий метод, основанный на идентификации по входу и выходу на основе метода наименьших квад- ратов, предложен в работах [131, 155]. Из уравнении (2.107) можно видеть, что нелинейный характер нор- мальных уравнений обусловлен неизвестной взаимной корреляцией между входом и выходом. Если лп не поддаются наблюдению, то значения Rnx(k) оценить нельзя. Если, однако, п* известны, то значения ЯВ*(Л) можно оценить, и поэтому АРСС-параметры
спектральный анализ 31 можно определить как решения системы нелинейных уравнений. На практике значения пп можно оценить по хп с помощью подхода, который будет рассмотрен ниже. Для записи требуемой системы линейных урав- нений перепишем хп из (2.105) в виде Из (2.115) получаем где 2=Нв+р, (2.116) *=[х0 xN.i]T, V=[no ni ~а2 -ар bt b2 — Л-! X-2 • • • X-p n-f Л-2 X0 X-J ••• X_p + I n0 n-1 -1 ~3 * * *N -p -1 П/f «2 -3 ’ ‘ ’ Заметим, что матрица H имеет размерность N X (р+ч)- Уравнение (2.115) имеет стандартную форму для за- дачи линейного оценивания по методу наименьших квадратов, решение которого имеет вид $ = (НнЙ)~1Ннх. (2.117) Данный подход аналогичен методу на основе наимень- ших квадратов для оценивания АР-параметров «чис- того» АР-процесса. Корреляционная матрица ННН АРСС-процесса включает оценку взаимной корреля- ционной функции Япк(А). Для этой цели необходимо задать начальные условия {х-р, х-р+1, . . ., х_ъ и_4, п.,+1, .... п-1} или же положить их равными нулю (так же, как и при оценивании АР-параметров по ме- тоду наименьших квадратов). Для оценки значений пп [155] необходимо моделировать хп с помощью АР- модели большого порядка и положить значения лп равными значениям последовательности ошибок пред- сказания. Оценки АРСС-параметров можно далее улучшить с помощью итеративной процедуры. Данный метод хорошо работает лишь в том случае, когда нули функции В (г) расположены внутри единичного круга. Следует заметить,' что уравнение (2.117) является решением по методу наименьших квадратов, соответ- ствующим^,(2.107), в котором точно известные стати- стические корреляции заменены их оценками. Подход, основанный на методе наименьших квадратов, не ис- пользует дополнительную информацию о том, что Япх(*)=0 при fc>0. Методы оценивания АРСС-параметров продолжают оставаться областью активного исследования, посколь- ку ни один из имеющихся методов пока еще, по-види- мому, не превосходит другие методы как с точки зре- ния достижимых характеристик, так и в смысле более низкой вычислительной сложности. H. Метод гармонического разложения Писаренко Если случайный процесс состоит из смеси сину- соид и аддитивного белого шума, то его можно моде- лировать как частный случай АРСС-процесса. В отли- чие от периодограммы в данной модели предполагает- ся, что в общем случае синусоиды гармонически меж- ду собой не связаны. Математические свойства данно- го АРСС-процесса приводят к анализу собственных значений оценок его параметров. Следовательно, под- ход к анализу этого процесса несколько отличается от подхода к анализу обобщенного АРСС-процесса. Смесь синусоид и аддитивного шума часто исполь- зуется в качестве испытательного процесса для ана- лиза методов оценивания спектра. Для обоснования причины выбора АРСС-процесса в качестве подходя- щей модели для смеси синусоид и белого шума рас- смотрим следующее тригонометрическое тождество: sin (Пп) = 2 cos 12 sin (S2(n _ 1]) - sin (J2[n - 2]), (2.118) где —л<й^л. Если Й=2л/Д(, где —1/2Д«/^1/2Д/, то sin Qn представляет собой синусоиду, дискретизо- ванную с интервалом AZ с. Полагая xn—sin(Qn), вы- = матрица входных/выходных данных. Пдг-д-! ражение (2.118) можно переписать в виде разностного уравнения второго порядка x„ = (2cosl2)x„_1 - х„_2, (2.119) которое позволяет рекурсивно вычислять текущее зна- чение синусоиды по двум предыдущим значениям хп-2 и хп_2. Если берется г-преобразование от (2.119), то X(z)(l - 2 cos 12z-1 + z~2] = D(z), (2.120) где D (z) — полином второго порядка, который ха- рактеризует начальные условия. Он имеет характери- стический полином 1—2 cos Dz-1-bz~s или, что экви- валентно, z’—2 cos Qz+1 с корнями Zi=exp(/2nfA() и z2=zj=exp(—j2nfAt). Эти корни имеют единичный модуль, т. е. |zll=|zal = 1, и частоты синусоид (в Гц) выражаются через эти корни следующим образом: fi= [ arctg (Im{z,}/Re{z,})]/2nAt, i=l,2. (2.121) Заметим, что Д=—Нетрудно видеть, что (2.119) является предельным случаем АР (2)-г.роцесса, у кото- рого дисперсия возбуждающего шума стремится к ну- лю. а полюсы — к единичному кругу. Кроме того, зная лишь два коэффициента и два отсчета, с помощью уравнения (2.119) можно точно предсказать синусои- дальный процесс для всего времени. В общем случае разностное уравнение второго по- рядка с вещественными коэффициентами хп = ~ У, атхп-т (2.122) т =1 может описывать детерминированный процесс, со- стоящий из р вещественных синусоид вида sin(2n/jA0- В этом случае коэффициенты {ат} являются коэффи- циентами полинома z2p +aiz2p~i +• • • + ap-1zp + 1 + apzp +ар^.1г1'~1 + ••• + + а2р z + а2р = У (z - zt)(z - z*) = 0 (2.123) i = i с корнями, имеющими единичный модуль, которые появляются в виде комплексно-сопряженных пар вида Zj=ехр(/2л/(Д(), где ft — произвольные частоты, та-
32 ТИИЭР, т. 69, № 11, ноябрь 1981 кие что —l/2M^fi<.l/2&t, 1=1, . . ., р. Можно пока- зать, что для чисто гармонического процесса а{= =ау>-ь где 1=1, р. Для синусоиды в аддитивном белом шуме шп наблю- даемым процессом является процесс 2₽ л.=*п + *’я=- £ атхп-т+*'п, (2.124) m «1 где £ [шпшп+а]=<46*, £[щп]=0 и Е [хпшп]=0, по- скольку шум полагается не коррелированным с сину- соидами. Подставляя в (2.124) xn_m=yh_m—wn-m, уравнение (2.124) можно переписать в виде 2р гр У. <ЬпУп-т= £ Wn-m, (2.125) т«о т«о где Оо= 1 по определению. Выражение (2.125), впервые полученное Ульрихом и Клейтоном [251], описывает синусоиды в белом шумовом процессе с помощью шума wn и зашумленных наблюдений уп;'оно имеет структуру АРСС(р, р)-процесса. Однако этот АРСС- процесс обладает частным видом симметрии, в резуль- тате чего АР-параметры оказываются идентичными СС-параметрам модели Если автокорреляционная функция процесса уп известна, то АРСС-параметры можно найти из решения собственного уравнения, которое здесь не рассматри- вается. Эквивалентная матричная запись уравнения (2.125) имеет следующий вид [67]: УТА-*ТА, (2.126) = [j’n Уп-i Уп-zplt А^’=[1 01 ••• Ojp^iO2p], — [wn Wn-i Wn_2p]. Умножая обе стороны этого уравнения на вектор Y и беря математическое ожидание, получим Е[УУТ]А =E[YWT]A. (2.127) Полагая — [^n ‘ ‘ -2pi> получим 2?^у(0) * • • Ryy(~ 2р) е[уут]=луу= • •_ ; 1^уу(2р) * * * Луу(О) _ £[УИ'Г] = £Г(Х + JF) И'7’] = £[ WWT] = (2.128) = o»Z. (2.129) Здесь Ryv — тёплицева автокорреляционная матрица Наблюдаемого процесса, а / — единичная матрица. Тот факт, что £[XWr]=0, следует из предположения, что синусоиды не коррелированы с шумом. Выражение (2.127) можно переписать в виде RyyA=al,A, (2.130) которое представляет собой собственное уравнение процесса, в котором дисперсия шума of, является собственным значением автокорреляционной матрицы Вектор АРСС-параметров А является собствен- ным вектором, связанным с собственным значением oj, и промасштабированным так, чтобы первый элемент был равен 1. Уравнение (2.130) позволяет определить значения АРСС-параметров в теш случае, когда из- вестны значения автокорреляционной функции. Знать дисперсию шума oj, не обязательно. Можно показать [149, приложение С], что для процесса, состоящего из р вещественных синусоид и аддитивного белого шума, дисперсия aj, соответствует минимальному собствен- ному значению матрицы Rvv, когда размер этой мат- рицы равен или превышает (2р+1)Х(2р+1) (это ми- нимальное собственное значение повторяется, если размерность этой матрицы превышает 2Р+1). Уравнение (2.130) составляет основу процедуры гармонического разложения, разработанного Писа- ренко [195]. Эта процедура позволяет определить точ- ные частоты и мощности р вещественных синусоид в присутствии белого шума, если точно известны 2р+1 значений автокорреляционной функции, включая ее значение при нулевой задержке. Так как известными полагаются только значения автокорреляционной функции, информация о фазе каждой синусоиды те- ряется. Писаренко отметил применимость тригономет- рической теоремы Каратеодори [82] для разработки метода отыскания не только частот Ог=2л/)Д(, но также мощностей Pt=A}/2 и СПМ шума oJ,A( лишь по 2р+1 известным значениям автокорреляционной функции. Для синусоид в белом шуме автокорреляци- онная функция имеет вид Л^(0) = о» + Л, i=i р RyyW=y cos (2я/,*Д1) при *=#0.(2.131) 1-1 Заметив, что белый шум влияет только на член с ну- левой задержкой, Писаренко пришел к уравнению от- носительно собственных значений (2.130) с помощью подхода, основанного на теореме Каратеодори, а не с помощью рассматриваемого здесь подхода. После определения коэффициентов i корни гп по- линома г2р + a1z2₽ -1 + • • • + вгр-1г + агр = (2.132) формируемые из этих коэффициентов, будут опреде- лять частоты синусоид, поскольку эти корни имеют единичные модули при z„ = exp (/27т/„Дг). (2.133) [См. аргументацию, приведшую к (2.123).] Благодаря структуре уравнения (2.128) оказывается, что поли- ном (2.132) должен быть симметричным, т. е. i—aip-t. Ни одного рекурсивного метода решения уравне- ния (2.130) порядка р на основе решения этого урав- нения для порядка р—1 не известно. Если число сину- соид не известно, но значения автокорреляционной функции известны точно, то независимо от решения уравнения (2.130) для последовательно возрастающих порядков должны вычисляться до тех пор, пока мини- мальное собственное значение перестанет изменяться при переходе к следующему более высокому порядку. Это будет означать достижение правильного порядка. В этой точке минимальное собственное значение равно дисперсии шума. На практике, как правило, известны оценки значений автокорреляционной функции, п<> этому число синусоид р необходимо выбирать в каче- стве значения того порядка, при котором минимальное собственное значение уравнения (2.130) лишь незна- чительно отличается от минимального собственного
СПЕКТРАЛЬНЫЙ АНАЛИЗ 33 значения для порядка р—1. Вычислительные требо- вания для решения уравнения (2.130) можно несколь- ко ослабить, используя теплицеву структуру матрицы Rvw Минимальное собственное значение и связанный с ним собственный вектор можно определить с помощью классического метода, в котором последовательность векторов А(к + 1) = Я~^А(к), к = 0, (2.134) сходится в пределе к собственному вектору минималь- ного собственного значения при некотором заданном начальном значении 4(0). Уравнение (2.134) можно переписать в виде ЯууА(к + 1) = А(к), (2.135) которое можно решить относительно неизвестного вектора А (£4-1), задав вектор Л (А). Методы типа гаус- сова исключения требуют о(р’) операций. Однако известны алгоритмы [63, 277], в которых за счет использования тёплицевой структуры матрицы Ryy решение уравнения (2.135) получается за о(р’) опе- раций. Хорошим начальным вектором является еди- ничный вектор Дг(0)=[1, .... 1]; в этом случае собст- венный вектор Д(оо) получается уже после несколь- ких итераций. После того как вектор А найден, мини- мальное собственное значение (и, следовательно, оцен- ка дисперсии шума) определяется с помощью следую- щего выражения: л Tg л = (2-136) А А Заметим, что для использования в (2.134) вектор A(k) каждый раз повторно масштабируется с помощью рэлеевского коэффициента в (2.136) до тех пор, пока не будет достигнута сходимость к ^min- После определения частот по корням полинома для А можно определить мощности синусоид. С помощью (2.131) значения автокорреляционной функции от jRyy(l) до Ryy(ty можно записать в матричной форме я> = г, (2.137) где cos • • cos (2п/рДг) ' _СО8(2Я/1РДГ) ГЛ 1 • • • cos (2я/ррД^)_ “Л^(1)4 UpJ L*w(P)_ Матрица F состоит из членов, зависящих от частот синусоид, которые определяются посредством нахож- дения корней полинома. Мощности этих синусоид находятся посредством решения системы уравнений (2.137) относительно вектора мощности Р. Мощность шума можно также определить из уравнения р о»=Rj)jl(o) - ^2 л- <•1 Краткая запись метода гармонического разложения Писаренко (ГРП) представлена на рис. 12. Поскольку порядок обычно не известен, определе- ние порядка должно проверяться минимальным собст- венным значением, используемым в нескольких реше- ниях собственного уравнения (2.130). Однако это не только требует больших вычислительных затрат, но р-1 -------------- Вычислить 2р-|-1 смещенных оценок автокорреляция (уравнение ^(2.16)) Определить минимальное собственное значение А.т^п р * р*1 н связанный с ним собственный вектор Л (уравнения (2.135), (2.136)) Отличается ли найденное минимальное собственное значение Л_1п от его н предыдущей оценки? пет I Да Определить корни полинома по собственному вектору А (уравнение (2.132)) Вычислить частоты синусоид по найденным корням (уравнение (2.133)) Определить мощности синусоид (уравнение (2.137)) Рис. 12. Краткая запись процедуры спектрального разложения Писаренко. часто также не ясно, когда достигается минимальное собственное значение, поскольку обычно используют- ся оценки значений автокорреляционной функции, а не их точные значения. Если выбранный порядок слишком велик, то по собственному уравнению поряд- ка 2р будет вычислено р синусоид, даже если на самом деле их число меньше р. Если порядок слишком мал, то найденные спектральные составляющие будут по- являться на неправильных частотах. Использование смещенных оценок значений автокорреляционной функции гарантирует положительную полуопреде- ленность тёплицевой автокорреляционной матрицы. Однако неявное взвешивание смещенной оценки авто- корреляционной функции треугольным окном, о чем говорилось в подразд. II.В, существенно ухудшает точность оценок частоты и занижает мощность реаль- ных присутствующих сигналов. Это может также при- вести к появлению ложных составляющих при спек- тральном разложении. Можно было бы, конечно, ис- пользовать несмещенные оценки значений автокорре- ляционной функции типа (2.14), но в этом случае не гарантируется положительная неопределенность ав- токорреляционной матрицы, как того требует разло- жение Писаренко. Это может привести к отрицатель- ным собственным значениям и к неимеющим физиче- ского смысла частотным оценкам. В уравнении (2.130) вместо р пытались использовать нетёплицевы положи- тельно-определенные автокорреляционные матричные формы, такие как XfXi, рассмотренные в подразд. II.Е, однако при этом получались корни полиномов с модуляциями отличными от 1, что, по-видимому, не позволяло улучшить результаты. Для окрашенного шума, который приводит к конечному числу значений автокорреляционной функции с задержками, боль- шими нуля, можно применить модифицированный ме- тод Писаренко [214]. Подход к анализу собственных значений, лежащий в основе метода Писаренко, можно связать с идеей выделения максимума информации о сигнале посред- ством обработки с целью отыскания наибольших соб- ственных значений и соответствующих им собствен- ных векторов матрицы оценок корреляций [135, 185, 252]. 5 ТИИЭР№ 11
34 ТИИЭР, т. 69, № 11, ноябрь 1981 I. Энергетическая спектральная плотность Прони [36, 49, 101, 102, 156, 160, 199, 200, 217, 253, 261] Метод Прони, а именно метод моделирования по- следовательности эквидистантных отсчетов данных с помощью линейных комбинаций экспоненциальных функций, не является методом спектрального оцени- вания в обычном смысле, однако в этом подразделе мы дадим его спектральную интерпретацию. Гаспар Фриш (барон Пронй) [202] пришел к выводу, что законы, описывающие расширение различных газов, можно представить с помощью суммы экспоненциальных функций. Он предложил метод для интерполяции данных своих измерений, основанный на согласовании параметров экспоненциальной модели с измеренными точками и вычислении интерполированных значений посредством оценки экспоненциальной модели в этих точках. Современный вариант метода Прони сильно отличается от своего оригинала вследствие изменений, которым он подвергся за прошедшие с того времени годы. Исходная процедура точно согласует экспонен- циальную кривую, содержащую р экспоненциальных членов (каждый член характеризуется двумя пара- метрами — амплитудой At и показателем степени at, где Л<ехр(а{0), с 2р результатами измерения данных. Этот подход обсуждался Хильдебрандом в работе [95]. В случае, когда требуется лишь аппрокси- мировать последовательность из N отсчетов данных с помощью р экспоненциальных функций, где N>2p, как правило, используется процедура оценивания по методу наименьших квадратов. Эта процедура и полу- чила название обобщенного метода Прони. Модель, предлагаемая в обобщенном методе Прони, представляет собой набор из р экспоненциальных функций с произвольными амплитудами, фазами, час- тотами и коэффициентами затухания. Функция диск- ретного времени р X bmz%., (2.138) m = l является моделью, которая используется для аппрок- симации измеренных значений данных х0, .... xN~i. Для общности величины Ьт и гт полагаются комплекс- ными: ьт =Ат ехр(/вт), zm = ехр [(ат +/2irfm)At], (2.139) где Ат — амплитуда, 0т — фаза Ь радианах, ат — коэффициент затухания, fm — частота осцилляций в герцах и Ai — интервал дискретизации в секундах. Отыскание параметров {Ат, 0m, ат, fm} и р, миними- зирующих квадратичную ошибку & = *£ |хп-£„|2, (2.140) П=0 представляет собой трудную нелинейную задачу ап- проксимации по методу наименьших квадратов. Ре- шение включает итеративный процесс, посредством которого первоначально принятые оценки значений этих параметров постепенно улучшаются. Макдунаф и Хаггинс [157] и Хольц [97] разработали такие ите- ративные схемы для решения уравнения (2.140). Аль- тернативное субоптимальное решение, которое не ми- нимизирует ошибку (2.140), но все же обеспечивает удовлетворительные результаты, основано на метода Прони. Метод Прони решает две системы линейных уравнений для последовательных значений оценок с использованием промежуточного этапа нахождения корней полиномов, который переносит проблему нели- нейности в процедуру нахождения корней полинома. Ключевым моментом метода Прони является тог факт, что функция (2.138) является однородным реше- нием линейного разностного уравнения с постоянными коэффициентами, форму которого можно определить следующим образом. Определим полином ¥(z) как р р Ф(г) = П (z - zk) = £ а/2₽-', а0 = 1. (2.141) i=0 Таким образом, полином V (г) состоит из комплексных экспонент zft, определяемых выражением (2.139), ко- торые используются в качестве его корней, и комп- лексных коэффициентов а^. С помощью (2.138) хпт можно записать в виде р хп-т = £ (2.142) Z = 1 где О^п—m^.N—1. Умножая (2.142) на ат и сумми- руя р+1 предшествующих произведений, получим р р р £ атх„-т = Х ь1 £ amz"’mr (2.143) m=o Z = 1 m«O где р^п^У—1. Если в (2.143) подставить z?-m= =z?~₽zf~m, то р р р £ W»-»-£ brf-p £ ^””=0. (2.144) m-o 1=1 т-0 Равенство нулю в (2.144) следует из того факта, что вторая сумма в этом выражении дает полином ’F(zi), определенный в (2.141) и вычисленный в точке, соот- ветствующей одному из его корней. Выражение (2.144) дает рекурсивное разностное уравнение ₽ *n=-£ “m*n-m> (2.145) т = 1 определенное при р^п^У—1. [Сравните это уравне- ние с уравнением (2.122) для процедуры гармониче- ского разложения Писаренко.] Таким образом, экс- поненциальные параметры находятся посредством оп- ределения корней полинома (2.141) с помощью коэф- фициентов От- Для формулировки обобщенного метода Прони положим сначала, что разность между реальными измеренными данными хп и аппроксимацией хп равна еп, так что хп=хп+еп, (2.146) где —1. Подставляя (2.145), получим урав- нение р хп=~ £ атхп-т+еп = т=1 Р Р = ~ £ атхп-т + £ атеп-т> (2.147) т=1 т=0 определенное при р^п^У—1, где использовано со- отношениехп_т—хп_т—еп_т. С помощью (2.147) мож- но показать, что альтернативной моделью суммы экс-
СПЕКТРАЛЬНЫЙ АНАЛИЗ 35 понеит и аддитивного шума является АРСС-модель с идентичными АР- и СС-параметрами, возбуждаемая шумовым процессом еп. В отличие от метода Писарен- ко коэффициенты а{ не ограничиваются здесь получе- нием корней полинома с единичным модулем (затуха- ние отсутствует). И хотя истинная оценка этих пара- метров по методу наименьших квадратов получается посредством минимизации величины квадратов дает известное решение ЭМФ^ФГЧ^Х. (2.151) Для уменьшения объема вычислений, гребуемых для (2.151), можно использовать следующее полезное со- отношение: это приводит к системе нелинейных уравнений, реше- ние которой получить трудно. Альтернативная про- цедура, получившая название обобщенного подхода Прони [266], определяет ошибку еп=£втеи_т, „=₽,- Л-1,(2.148) т»0 поэтому Р х„=-Е “т*п-т+еп. (2.149) т-1 Далее минимизируется величина 2«-plenl2, а не Следовательно, обобщенная процедура оце- нивания параметров Прони сводится к процедуре оце- нивания АР-параметров, в которой можно использо- вать алгоритм оценки ковариации (2.69) по методу наименьших квадратов с матрицей XiXi. Заметим, что небелый случайный входной процесс еп получает- ся из СС-процесса, возбуждаемого ошибкой аппрокси- мации еп, как это следует из (2.148). Кроме того, еп — это разность между хп и его линейным предсказанием, основанным на р предыдущих отсчетах, тогда как — это разность между хп и его экспоненциальной аппрок- симацией. Число экспонент р определяется с помощью методов выбора порядка АР-модели, рассмотренных в подразд. II.Е. Другая схема определения р основана на анализе собственных значений матрицы X^Xi [252] и тесно связана с анализом собственных значений не- тёплицевой матрицы с помощью метода Писаренко, рассмотренным в подразд. П.Н. После того как величины z{ определены из корней полинома, выражение (2.138) сводится к системе ли- нейных уравнений относительно неизвестных пара- метров Ьт, матричная форма которого имеет следую- щий вид: ФЭ-Х, (2.150) где 1 1 ••• 1 ф = zi *2 2р В=1Ь1 ••• ЬР]Т, Х=[х0 хлг_1]г. Заметим, что матрица Ф является матрицей Вандер- монда, аналогичной матрице (2.24), с тем лишь отли- чием, что члены zt теперь имеют затухание и произ- вольные, а не гармонически связанные частоты. Ми- нимизация величины 2 (*—*)* по методу наименьших где (2.152) Определение параметров аг посредством оценивания по методу наименьших квадратов, определение корней полинома и последующее решение уравнений относи- тельно параметров bf (или остатков) и составляет об- общенный метод Прони. Для определения амплитуды Ai, фазы 0Ь коэффициента затухания а( и частоты ft по оценкам величин zt и bt используются простые вы- ражения ^ = |Z>il, et = arctg [Im (Z»/)/Re (*,)], a( = ln |Z(|/Ar, ft = arctg [Im(z,)/Re(г()]/2тгДг. (2.153) Обычно метод Прони завершается вычислением экспоненциальных параметров, определенных выра- жениями (2.153). Как таковой, метод Прони наиболее часто применяется при анализе переходных процес- сов, например, для определения резонансных мод при воздействии электромагнитных импульсов [196]. Од- нако с его помощью можно выполнить «спектральный анализ» следующим образом. И хотя можно опреде- лить много различных спектров, один «спектр» может оказаться полезным, если сделать предположение, что модель процесса имеет симметрию, что проиллюст- рировано для одной затухающей вещественной сину- соиды на рис. 13. Предполагаемая аппроксимирующая функция принимает форму х(Г)= 2 Ат exp(am|t|)exp(;[2ff/mr + em]), (2.154) m-1 определенную при —оо<«со. Для вещественного процесса x(t) в (2.154) потребуются комплексно-сопря- женные пары типа exp/(2nfmf-|-0m) и exp — / (2nfmt+ +0m). Далее предполагается, что коэффициенты зату- хания отрицательны, поэтому получаются затухаю- Рнс. 13. Симметричная экспоненциальная модель огибающей.
36 ТИИЭР, т. 69, № 11, ноябрь 1981 Энергетическая спектральная плотность Частота Рнс. 14. Узкополосная и широкополосная спектральные оценки- Прони. щие экспоненциальные функции. Одна из причин вы- бора симметричной огибающей состоит в том, что при а=0 аппроксимирующая функция x(t) будет содер- жать незатухающие синусоидальные составляющие, определенные на интервале —оо</<оо. В результате невзвешанные синусоиды точно моделируются при таком подходе. Поскольку функция (2.154) имеет конечную энер- гию, то детерминированное выражение для ее энерге- тической спектральной плотности, основанное на пре- образовании Фурье для (2.154), имеет вид $ Прони СП = 1*(Л12, (2.155) где л р 2а™ Х</) = £1 Ат (/0Й,) [«^+(2^/-/ml)2] ’ (2156) Последнее выражение и представляет собой один из возможных «спектров» Прони. Заметим, что спек- тральная оценка (2.155) сохраняет максимумы, вели- чины которых линейно пропорциональны этой энер- гии, в отличие от максимумов АР-спектра, которые нелинейно связаны с мощностью [136]. Спектр Прони обладает способностью получать узкополосные или широкополосные спектральные линии, форма которых зависит от величины коэффициента затухания (см. рис. 14). Колоколообразная кривая имеет ширину (на уровне —3 дБ), равную a/я Гц, поэтому разрешение зависит от коэффициента затухания. Заметим, что для выбора порядка модели р метод Прони требует 2р параметров для описания спектра, что вдвое больше числа параметров, требуемых для спектральной АР-оценки. Кроме того, метод Прони обеспечивает получение информации о фазе, которая теряется при спектральном АР-оценивании. Процедура спектраль- ного оценивания по методу Прони кратко представле- на на рис. 15. Метод Прони является адаптивной процедурой в том смысле, что она подстраивает параметры модели, С помощью любого алгоритма АР-оценивания оценить коэффициенты полинома и порядок модели I Найти корни полинома для определения частот и коэффициентов затухания (уравнения (2.141), (2.153)) Определить амплитуды и фазы (Уравнения (2.151), (2.153) Вычислить преобразование экспоненциальной модели и взять модуль для получения спектра (уравнения (2.155), (2.156)) Ряс. 15. Процедура спектрального оценивания по методу Прони. а именно частоту, фазу, амплитуду и коэффициент за- тухания экспоненциальных функций модели для их согласования с реальными данными. В отличие от нее периодограмма использует фиксированное число не- затухающих синусоид с фиксированными частотами. Применение метода Прони связано с рядом труд- ностей, одну из которых необходимо четко осознавать. Это — проблема, связанная с определением числа экспоненциальных членов, которая аналогична про- блеме выбора порядка модели при АР-оценивании, поэтому здесь приложимы те же самые соображения. Однако поскольку здесь вычисляется 2р параметров, максимальный порядок ограничен величиной p^.NI2, тогда как при спектральном АР-оценивании возможен порядок p>NI2 (хотя это и не рекомендуется). В ряде случаев точность оценивания полюсов по методу Прони будет сильно зависеть от шума [252—254]. Наличие шума может также привести к слишком боль- шим коэффициентам затухания. К. Оценивание спектральных линий по методу Прони Для процессов, состоящих из р вещественных не- затухающих (а=0) синусоид и шума, разработан частный вариант метода Прони. Основной подход описан Хильдебрандом [95]. В рассматриваемом слу- чае выражение (2.138) можно записать в виде р р хп ~ ^mzmn] = £ 008 (2я’ЛплА1 + 0т), т = 1 т = 1 (2.157) где 6т=Атехр(/0т)/2 и гт=ехр(/2л/тД(). Заметим, что гт являются корнями единичного модуля с произ- вольными частотами, которые появляются в виде комплексно-сопряженных пар до тех пор, пока fm=^) или 1/2Д(. Таким образом, необходимо решить урав- нение (2.141) для отыскания корней полинома р гр Ф(г) = П ’ г<)(* " 2*) = £ в*г2р'* = 0, (2.158) /»1 fc=O где Оо=1, a at— вещественные коэффициенты. По- скольку эти корни имеют единичный модуль и появ- ляются в виде комплексно-сопряженных пар, уравне- ние (2.158) должно быть инвариантным относительно подстановки z~J вместо z: /1 \ гр гр z*P* -) = z2p £ akz ip = £ акгк = 0. (2.159) Сравнивая (2.158) и (2.159), можно видеть, что Я/= ~а2р-/ ПРИ /—0. Р> где Оо=аар=1. Следовательно, требование на комплексно-сопряженные пары корЛей единичного модуля реализуется посредством наложе- ния на коэффициенты полинома свойства симметрии относительно центрального элемента. Для процесса порядка 2р ошибку линейного предсказания, анало- гичную (2.149), можно переписать в виде 𠈄= £ Om(xn-tm+хп-т)< (2.160) т »0 что позволяет вдвое уменьшить требуемое число коэф- фициентов. Здесь применимы все подходы к миними- зации по методу наименьших квадратов. Отличие те- перь состоит в том, что матрица данных X, которая
СПЕКТРАЛЬНЫЙ АНАЛИЗ 37 характеризует уравнение (2.160), имеет тёплицево- ганкелевскую структуру, а не просто тёплицеву, как в АР-случае. Например, используя ошибку е„, изме- няющуюся от п=р до N, получим ep+i eN-p (2.161) где Тёплицева матрица данных Ганкелева матрица данных поэтому минимум квадратичной ошибки определяет вещественные коэффициенты аи ..., ар_и ар/2, аналогичные решениям при групповом АР-оце- нивании. Заметим, что последний коэффициент равен ар/2, а не ар. Множитель 1/2 обусловлен симметрией матрицы X, что приводит к тому, что последний коэф- фициент учитывается дважды. Эти коэффициенты ис- пользуются для записи полинома порядка 2р+ 1 с сим- метричными коэффициентами (2.158). И хотя корни единичного модуля приводят к симметричному поли- ному, справедливость обратного утверждения вовсе не обязательна. Симметричные коэффициенты лишь га- рантируют, что если встретился корень zit то обяза- тельно должен быть и корень zf1; требовать для этого |z(|=1 не обязательно. На практике [226] корни с от- личными от 1 модулями встречаются очень редко, но когда они появляются, то они, как правило, соответст- вуют частотам ft—Q или ft—1/2 At. Данный алгоритм завершается определением амплитуды и фазы из уравнения (2.151), размерность которого можно умень- шить вдвое, объединяя связанные комплексные пары. Спектр будет состоять из дельта-функций, представ- ляющих синусоиды, и затухающих экспоненциальных функций для тех редких случаев пар корней с отлич- ными от 1 модулями. Метод гармонического разложения Прони, описан- ный выше, имеет несколько преимуществ перед мето- дом ГРП. Во-первых, он не требует оценки значений автокорреляционной функции. Метод Прони дает меньшее число ложных спектральных линий, чем ме- тод Писаренко, поскольку порядок модели можно точ- нее определить посредством контроля величины оста- точной квадратичной ошибки в частном методе Прони. Во-вторых, оценки частоты и мощности имеют меньшее смещение по сравнению с аналогичными оценками, получаемыми с помощью метода Писаренко [153, 226]; см. рис. 16 для сравнения спектров, получаемых с помощью этих методов. Метод Прони требует реше- ния только двух систем линейных уравнений и опре- деления корней полинома, тогда как метод Писаренко требует более сложного в вычислительном отношении решения собственного уравнения. L. Спектральное оценивание е помощью метода максимального правдоподобия Кейлона [42, 137, 203] При спектральном оценивании с помощью метода максимального правдоподобия (ММП), который перво- начально был предложен для пространственно-времен- ного анализа сигналов сейсмических решеток [42], оценивание СПМ осуществляется посредством изме- рения мощности сигналов на выходах набора узкопо- лосных фильтров [136]. Следует заметить, что назва- ние этого метода спектрального оценивания — ММП — на самом деле неправильно в том смысле, что получае- мая с его помощью спектральная оценка не Является истинной оценкой СПМ, соответствующей максималь- ному правдоподобию. Спектральную оценку с помощью ММП иногда называют спектральной оценкой Кейпона [92]. Название ММП мы сохраняем здесь только по историческим причинам. Различие между спектраль- ным оцениванием с помощью ММП и стандартным методом спектрального оценивания Блэкмана — Тью- ки (на основе периодограммы) состоит в том, что фор- ма частотной характеристики узкополосных фильт- ров, используемых в ММП, вообще говоря, различна для каждой частоты, тогда как в методе Блэкмана — Тьюки она фиксированна. Эти фильтры адаптивно подстраиваются под процесс, для которого ищется СПМ. В частности, эти фильтры относятся к типу филь- тров с конечной импульсной характеристикой (КИХ) и имеют р весовых коэффициентов (отводов у линии задержки) Л = [*о«1 •••вР-1]г. (2.162) Коэффициенты фильтра выбираются таким образом, чтобы на рассматриваемой частоте ft его отклик был равен единице (т. е. входная синусоида с этой часто- той остается неизменной на выходе фильтра), а дис- персия выходного процесса была минимальна. Следо- вательно, этот фильтр должен подстраиваться для режекции составляющих спектра, удаленных от час- тоты /о, поэтому выходная мощность будет обуслов- лена частотными составляющими, близкими к частоту ft. Для получения такого фильтра необходимо мини- мизировать дисперсию о1 выходного процесса, опре- деляемую выражением <?~АНЯХХЛ, (2.163) при единичном отклике фильтра на частоте ft (т. е. при наложенном ограничении, что синусоида с часто- той ft проходит через фильтр без искажений) ЕНА = 1, (2164) где Rxx — ковариационная матрица хя, Е — вектор, определяемый выражением Е= [1 ехр(/2я/оДО '’ ехр(/2тг[р - 1 ]/оДО]Т,
38 ТИИЭР, т. 69, М 11, ноябрь 1981 Доли частоты дискретизации (а) (Ь) ‘°0 0,1 0,2 0,3 0.4 0,5 Доли частоты дискретизации (С) (d) (®) (f) Доли частоты дисретизации 0) S-” и О 0,1 0,2 0,3 0,4 0,5 Доли частоты дисретизации (j) (к) (1) Рис. 16. Спектральные оценки для одной и той же 64-точечной последовательности, полученные с помощью различных методов: (а) истинная СПМ; (Ь) СПМ на основе периодограммы; БПФ, число отсчетов удвоено за счет введения дополнительных нуле- вых отсчетов; (с) СПМ Блэкмана — Тьюки; (d) АР-оценка СПМ на основе подхода Юла — Уокера; (е) АР-оценка СПМ с помощью алгоритма Берга; (f) АР-оценка СПМ по методу наименьших квадратов или алгоритма предсказания вперед и назад; (g) ОС оценка СПМ; (h) АРСС-оценка СПМ на основе обобщенных уравнений Юла — Уокера; (!) СПМ на осно- ве метода спектрального разложения Писаренко; (j) энергетическая спектральная плотность Прони; (к) частный вари- ант метода Прони на основе подхода Хильдебранда; (1) СПМ Кейпона (метод максимального правдоподобия). а надстрочный индекс Н означает транспонированную комплексно-сопряженную величину. Решение для ве- совых коэффициентов фильтра имеет, как нетрудно по- казать, следующий вид [203]: Минимальная дисперсия выходного процесса в этом случае равна а °min t? ’ Д Дуу Д (2.166) . R~'XE opt EhR-^E ' (2.165) Нетрудно видеть, что отклик оптимального фильтра на частоте f=ft равен единице и что характеристики
Таблица 2 Сводка современных методов спектрального оценивания Метод Основные работы Джскретный (д) или неп- рерывный (Н) (спектр) Алгоритм Вычислительная сложность (прибли- зительное число операций сложения и умножения) Модель Преимущества и недостатки Примечания Периодо- - грамма (ва- риант на основе БПФ) **> (24, 46, 207, 218, 262] д Уравне- ние (2.18) N log,N Сумма гармони- чески связан- ных синусо- ид Выход прямо пропорцио- нален мощности Наиболее эффективен в вычислительном отно- шении Разрешение приближенно равно величине, обрат- ной интервалу наблю- дения Плохие характеристики при малой длине запи- си данных Просачивание искажает спектр и маскирует слабые сигналы Гармоническое согласова- ние по методу наимень ших квадратов Для стабилизации спект- ра необходимо приме- нять какую-либо разно- видность статистичес- кого усреднения (на- пример, метод Уэлша) Взвешивание может уменьшить боковые ле- пестки, но лишь за счет ухудшения разрешения Блэкмана— Тьюки *♦> (БТ) [25, 33,' 204] н Уравне- ния (2.13), (2.16) Оценка значений автокорреляци- онной функции: NM Оценка СПМ*>: MS Идентична СС- модели со взвешиванием значений ав- токорреляци- онной функ- ции В вычислительном отно- шении наиболее эффек- тивен при М Разрешение приближенно равно величине, обрат- ной интервалу наблю- дения Просачивание искажает' спектр и маскирует сла- бые сигналы При некоторых видах взвешивания (типах окон) и оценок авто- корреляции (например, несмещенных) в спектре могут появиться отри- цательные значения СПМ Авторегрес- сивный (АР) на ос- нове урав- нений Юла— Уокера [59, 108, 250, 276] н Рис. 6**» Оценка значений автокорреляци- онной функции: NM АР-коэф^ициенты: Оценка СПМ •>: MS Авторегрессив- ный процесс (только с од- ними полю- сами) Требуется определять по- рядок модели Более высокое разреше- ние, чем у БПФ или БТ, но не такое хоро- шее, как у других АР-методов Возможно расщепление спектральных линий Неявное взвешивание ис- кажает спектр Боковые лепестки отсут- ствуют Модель применима для описания сейсмических и речевых сигналов, а также радиолокацион- ных сигналов, отражен- ных от облака пассив- ных отражателей Гарантируется минималь- но-фазовый (устойчи- вый) линейный пред- сказывающий фильтр, если вычисляются сме- щенные оценки значе- ний автокорреляцион- ной функции АР-процесс, связанный с анализом на основе ли- нейного предсказания и адаптивной фильтра- цией Максимумы (пики) в спек- тре моделируются луч- ше, чем минимумы (впадины) Авторегрес- сивный (АР) на ос- нове алго- ритма Бер- га [9,15,37, 41,57,64, 65) н Рис. 7 АР-коэффициенты: ЛГЛН-М* Оценка СПМ *>: MS Авторегрессив- ный процесс (только с од- ними полю- сами) Высокое разрешение при низком уровне шума Хорошая точность спек- тральных оценок в слу- чае коротких записей данных Смещение спектральных оценок спектральных линий по частоте Неявное взвешивание от- сутствует Боковые лепестки отсут- ствуют Требуется определять по- рядок модели Гарантируется устойчи- вый линейный предска- зывающий фильтр Применима адаптивная фильтрация Используется рекурсив- ный подход на основе метода наименьших квадратов с наложен- ными ограничениями
Продолжение Метод । Основные работы дискретный (Д) или неп- рерывный (Н) спектр Алгоритм Вычислительная сложность (прибли- зительное число операций сложения и умножения) Модель Преимущества и недостатки Примечания Авторегрес- сивный (АР) на ос- нове мето- да наи- меньших квадратов или линей- ного пред- сказания вперед и назад 154, 227] н Струк- турная схема алгорит- ма при- ведена в работе [154] АР-коэффициенты: NM-\-M* Оценка СПМ *>: MS Авторегрессив- ный процесс (только с од- ними полю- сами) Более острые спектраль- ные линии для узкопо- лосных процессов, чем у других АР-оценок Расщепление спектраль- ных линий отсутствует Уменьшенное смещение частотных оценок Требуется определять по- рядок модели Боковые лепестки отсут- ствуют Не гарантируется устой- чивый линейный пред- сказывающий фильтр, хотя в большинстве случаев получается ус- тойчивый фильтр На основе точного рекур- сивного решения по методу наименьших квадратов без наложен- ных ограничений Скользящего среднего (СС) [47, 236] н Уравне- ния (2.104), (2.35) СС-коэффициенты системы нелиней- ных уравнений: NM Оценка СПМ •>: MS Процесс сколь- зящего сред- него (только с одними ну- лями) Широкие спектральные линии (низкое разре- шение) Требуется определять по- рядок модели Присутствуют боковые лепестки Обобщенная форма метода Блэкмана—Тьюки АРСС (вари- ант Юла— Уокера [31, 97] н Уравне- ния (2.108) (уравне- ния Юла — Уокера высокого порядка) Оценка значений автокорреляци- онной функции: Коэффициенты: М* Оценка СПМ MS АРСС-процесс (рациональ- ная переда- точная функ- ция; порядок АР-части про- цесса не равен порядку СС- части) Требуется определять по- рядок АР- и СС-частей процесса Моделируются все про- цессы с рациональными передаточными функци- ями Для получения хороших результатов требуются точные оценки значений автокорреляционной функции Писаренко (метод гар- мониче- ского раз- ложения, ГРП) [195, 251] д Рис. 12 Оценка значений автокорреляци- онной функции: Собственные ура- внения: от Afa до Л13 Корни полинома: число операций зависит от при- меняемого алго- ритма Мощности: М* Частный слу- чай АРСС- процесса с равными АР- и СС-коэффи- циентами Сумма негармо- нических не- затухающих синусоид в белом шуме Требуется определять по- рядок модели Плохо работает при вы- соком уровне шума Собственное уравнение и определение корней в вычислительном отно- шении не эффективны ,Для получения хороших результатов требуются точные оценки значений автокорреляционной функции Если порядок модели вы- бран слишком боль- шим, могут появиться ложные спектральные линии Прони (обоб- щенный) [95, 97, 252] н Рис. 15 АР-коэффициенты: Корни полинома: число операций зависит от при- меняемого алго- ритма Амплитудные коэф- фициенты: ЛР Оценка СПМ *>: MS Сумма негармо- нически свя- занных сину- соид АРСС с равны- ми коэффи- циентами и порядками [p = q) АР- и СС-частей мо- дели Требуется определять по- рядок модели Выход линейно зависит от мощности Требуется определять корни полинома Разрешение такое же хо- рошее, как у АР-мето- дов, иногда лучше Боковые лепестки отсут- ствуют Для получения экспонен- циальных параметров используются оценки по методу наименьших квадратов Первый этап тот же, что и при АР-оценивании по методу наименьших квадратов Прони (ме- тод гармо- нического разложе- ния) [95, 226] д Уравне- ния (2.161), К2Л58), (2.151), (2.153) Коэффициенты: ЛР Корни: число операций зависит от при- меняемого алго- ритма Амплитудные ко- эффициенты: АР Сумма негармо- нически свя- занных сину- соид Требуется определять по- рядок модели Выход линейно зависит от мощности Требуется определять корни полинома Разрешение такое же хо- рошее, как у АР-мето- дов, иногда лучше Боковые лепестки отсут- ствуют Используетсяоценивание по методу наименьших квадратов
СПЕКТРАЛЬНЫЙ АНАЛИЗ 41 П родомсение Метод Основные работы Дискретный <Д) или не- 1 прерывный | (Н) спектр Алгоритм Вычислительная сложность (прибли- зительное число операций сложения и умножения) Модель Преимущества и недостатки Примечания Кейлона (ме- тод мак- симально- го правдо- подобия, ММП) [40, 136, 193, 203) н Уравне- ния (2.167), (2.16) Оценки значений автокорреляци- онной функции: NM Обращение матриц: М» Оценка СПМ *>: MS Для каждой спектральной компоненты формируется отдельный по- лосовой фильтр Разрешение лучше, чем у метода Блэкмана— Тьюки В статистическом смысле ММП-спектры менее вариабельны, чем АР- спектры Спектральное оценивание по методу максималь- ного правдоподобия связано с АР-спектр амн (см. уравнение (2.168)) ♦ Для вычисления S=2V значений СПМ можно использовать БПФ. ** Машинные программы приведены в сборнике Programs for Digital Signal Processing, IEEE Press, 1979, выпущенном Комитетом по обра- ботке цифровых сигналов Общества акустики, речи и обработки сигналов ИИЭР. /V —число отсчетов данных; 3 —число рассчитанных точек спектра (обычно о»Л4); М — порядок модели (или число дискретных значений автокорреляционной функции). фильтра зависят от автокорреляционной функции вы- ходного процесса. Поскольку минимальная дисперсия выходного процесса зависит от мощности составляю- щих спектра вблизи частоты /0. то величину CTmin М можно интерпретировать как оценку СПМ. Таким об- разом, оценка СПМ с помощью ММП определяется как ?ММП(/о) = рЯр-1 • (2.167) " ^хх^ Для вычисления этой спектральной оценки необхо- димо знать только оценку автокорреляционной мат- рицы. На практике оценивание с помощью ММП обеспе- чивает более высокое разрешение, чем методы спек- трального оценивания на основе периодограммы и метода Блэкмана — Тьюки, но меньшее, чем АР-ме- тоды спектрального оценивания [136]. В случае, когда необходимо оценивать автокорреляционную функцию, спектральное оценивание с помощью ММП обеспечи- вает (как показали эксперименты и расчеты для длинных записей данных) меньшее значение диспер- сии, чем у спектральной АР-оценки [13]. Следует за- метить, что для узкополосных процессов, автокорре- ляционная функция которых известна, максимум (пик) АР-спектра пропорционален квадрату мощности ана- лизируемого процесса, тогда как для ММП-спектра он пропорционален мощности [136, 137]. Кроме того, мощность спектральной АР-оценки можно найти, определяя площадь пика (спектральной линии), пло- щадь же под пиком ММП-спектра пропорциональна корню квадратному из этой мощности. Спектральные АР- и ММП-оценки связаны следую- щим соотношением [40]: 1 1 р 1 ---Ё ’ (2.168) где ^др’(/) — АР-оценка СПМ для модели m-го по- рядка, а 5*мп(/) — ММП-оценка СПМ; обе эти оценки получены на основе известной автокорреляционной матрицы порядка р [40, 203]. Таким образом, более низкое разрешение ММП-оценки можно объяснить эффектом «усреднения посредством параллельной ре- зистивной цепи», т. е. эффектом, обусловленным объе- динением спектральной АР-оценкой малого порядка 6 ТИИЭР № 11 с наименьшим разрешением и спектральной АР- оценкой высокого порядка с наиболее высоким раз- решением. Отметим также тот интересный факт, что обратное преобразование Фурье для оценки которое дает оценку автокорреляционной функции, не идентично автокорреляционной функции, исполь- зуемой для получения СПМ. С другой стороны, обрат- ное преобразование Фурье для АР-оценки СПМ дает идентичную автокорреляционную функцию во всем диапазоне известных задержек, как нетрудно видеть из выражения (2.61). XII. КРАТКАЯ СВОДКА ОСНОВНЫХ РЕЗУЛЬТАТОВ В табл. 2 сведены и коротко описаны все 11 наибо- лее часто используемых методов спектрального оце- нивания, которые были рассмотрены в данной статье. В этой таблице коротко описаны основные свойства, указаны номера уравнений, используемых для вы- числения каждой спектральной оценки, а также при- ведены номера основных работ по каждому методу, что должно помочь читателю быстро реализовать лю- бой избранный им метод. На рис. 16 показаны типичные оценки спектров, получаемых с помощью всех 11 методов, представлен- ных в табл. 2. Каждая спектральная оценка основана на 64-точечной последовательности вещественных от- счетов, полученной из процесса, состоящего из трех синусоид и окрашенного шумового процесса, прлу- ченного посредством фильтрации белого гауссовою шума. В табл. 3 представлены значения использован- ных отсчетов данных. Истинная СПМ показана на рис. 16а. Частота изменяется от 0,0 до 0,5 и выражена в долях частоты дискретизации входного процесса. Три синусоиды имели частоты 0,10, 0,20 и 0,21 и отно- шение сигнал/шум +10, +20 и +30 дБ соответствен- но, где отношение сигнал/шум определяется как от- ношение мощности каждой синусоиды к полной мощ- ности шума в полосе пропускания. Полоса шумового процесса центрирована относительно частоты 0,35. Описанный входной сигнал был выбран для демонст- рации поведения каждого метода спектрального оце- нивания при наличии узкополосного и широкополос- ного процессов. Рисунок 16 предназначен лишь для иллюстрации свойств каждого метода, в частности в случае коротких записей данных, а потому его не еле-
42 ТИИЭР, т. 69, № 11, ноябрь 1981 Таблица 3 Список значений отсчетов данных Х( 1)- 1.291061 Х( 2)- -2.006368 Х( 3>- -1.691316 Х( 4)- 1.243138 Х( 5)- 1.641072 Х( 6)- -0.008688 Х( 71- -1.659390 Х< 8)- -1.111467 Х< 9)- 0.985908 Х(10)- 1.991979 Х(11)- -0.046613 Х(12>- -1.649269 Х<13>* -1.040810 Х(14)- 1.054665 Х(15)- 1.855816 ХС16)- -0.951182 Х<17>- -1.476495 Х(18)- -0.212242 Х<19)- 0.780202 ХС201- 1.416003 Х<21>- 0.199202 Х(22)- -2.027026 Х(2Э>- -0.483577 ХС24)- 1.664913 Х(25)- 0.614114 Х(2б)> -0.791469 Х(27>- -1.195311 Х<28>- 0.119801 Х(29)- 0.807635 Х(30>- 0.895236 Х<31)- -0.012734 Х(32)- -1.763842 Х(33>- 0.309840 Х<34>- 1.212892 Х<35>- -0.119905 Х(Э6>- -0.441686 Х<37)- -0.879733 Х<38)> 0.306181 Х(Э9)- 0.795431 ХС40)- 0.189598 ХС417- -0.342332 Х<42)- -0.328700 Х(43)- 0.197881 Х(44>- 0.071179 Х(45)- 0.185931 Х(46>- -0.324595 Х<47)- -0.366092 Х(48)- 0.368467 Х(49>- -0.191935 Х(50)- 0.519116 Х(51>" 0.008320 Х<52>- -0.425946 Х(53)- 0.651470 Х(54)- -0.639978 Х(55>- -0.344389 Х<56>- 0.814130 Х<57)- -0.385168 Х(58>- 0.064218 Х(5У>- -0.380008 Х(60>- -0.163008 Х(61>- 1.180961 Х(62)- 0.114206 Х<63)- -0.667626 Х<64>- -0.814997 дует рассматривать как основу для сравнения относи- тельных характеристик этих методов. На рис. 16b показана периодограмма 64-точечной последовательности, значения которой приведены в табл. 3. Периодограмма получена с помощью БПФ, причем число отсчетов входной последовательности было удвоено за счет введения 64 дополнительных нулевых отсчетов. Номинальное разрешение (в Гц) 64-точечной последовательности равно 0,015625 час- тоты дискретизации, поэтому синусоиды с относитель- ными частотами 0,20 и 0,21 разделены расстоянием, меньшим интервала разрешения. Ясно, что показан- ная на рисунке периодограмма не позволяет разре- шить эти две синусоидальные составляющие. Более слабую синусоиду можно вицеть среди боковых лепест- ков (взвешивание данных не применялось). Присутст- вие окрашенного шума также указывается дискретны- ми спектральными линиями в верхней части частот- ного диапазона. На рис. 16с показана БТ-оценка спек- тра, основанная на 16 значениях оценок автокорреля- ционной функции, что составляет примерно 20% числа отсчетов данных, в соответствии с рекомендацией Блэкмана — Тьюки. Несколько АР-оценок СПМ показано на рис. 16d— 16f. И хотя все это АР-оценки различаются только спо- собом оценивания АР-коэффициентов, результирую- щие спектры весьма различны. По 64 отсчетам данных были вычислены 16 коэффициентов для АР-оценки и всех оставшихся методов, рассмотренных в статье. АР-подход на основе уравнений Юла — Уокера, ко- торый требует оценивания значений автокорреляци- онной функции, не позволяет разрешить две близко расположенные синусоиды (он имеет наихудшее раз- решение из всех АР-методов) и обеспечивает незначи- тельную информацию относительно состава спектра по обе стороны от главного максимума. АР-оценка СПМ, основанная на алгоритме Берга, показанная на рис. 16е, имеет острые максимумы на частотах всех трех синусоид, хотя одна из них с частотой 0,10 едва видна при избранном масштабе. Она также показы- вает мощность на высокочастотном конце спектра, хотя она и не столь гладка и широкополосна, как должна быть на самом деле. Эта оценка иллюстрирует «остроконечную» природу АР-спектра. Более точную оценку для частот трех синусоид дает метод спектраль- ного АР-оценивания на основе линейного предсказа- ния вперед и назад (или на основе наименьших квад- ратов), как показано на рис. 16f. В то же время сле- дует отметить сравнимость результатов АР-рценок, показанных на этих двух рисунках. На рис. 16g показана СС-оценка СПМ. Она идентич- на оценке СПМ методом Блэкмана — Тьюки, посколь- ку в ней также используются только оценки значений автокорреляционной функции. Широкие максимумы СС-оценки контрастируют с узкими максимумами АР-оценок. СС-оценка не позволяет разрешить две близкие по частоте синусоиды. Спектральная линия на частоте 0,10 так же широка, как и спектральные линии на высокочастотном конце спектра, что затрудняет обнаружение узкополосных составляющих. На рис. 16h показана АРСС (8, 8)-оценка СПМ, основанная на модифицированных уравнениях Юла — Уокера с использованием смещенных оценок значений автокор- реляционной функции, вычисленных по отсчетам дан- ных. Это — не очень хорошая спектральная оценка, хотя АРСС(16, 16)-оценка спектра, которая здесь не показана, позволяет разрешить три синусоидальные составляющие входного процесса. Метод спектрального разложения Писаренко при- водит, подобно БПФ-Ьериодограмме, к получению ди- скретного спектра (рис. 161). Две близкие по частоте синусоиды разрешаются, но частоты и мощности весь- ма не точны. Более слабая синусоида с частотой 0,10 имеет спектральную линию вблизи этой частоты, од- нако появилось много других спектральных линий, что затрудняет выбор истинных сигналов из ложных спектральных составляющих. Широкополосная часть спектра моделировалась посредством помещения не- скольких спектральных линий в область спектра ши- рокополосного процесса. Таким образом, метод Пи- саренко не позволяет достаточно точно моделировать широкополосные процессы, хотя он и показывает, что в данном частотном диапазоне присутствует некоторая мощность. Оценка энергетической спектральной плотности, основанная на обобщенном методе Прони, показана Таблица 4 Значения оценок параметров, полученные с помощью метода Прони Номер сигнала Оценки параметров 1 ЯГ- 0.0924044 PHASE- 127.8651 ОРРР- 0.0001263 ERE»- 0.108033 2 ATF- 1.2756225 PHASE- 150.5105 DHf>- 0.0006901 FREO- 0.201151 3 ЯГ- 0.0447961 PHASE- 190.7096 DW- 0.0092935 FREO- 0.208949 4 ЯГ- 0.2363050 PHASE- 208.9569 W*P- -0.0592305 FREO- 0.275976 5 (HP* 0.5034477 PHASE- 51.5229 DWF- -0.0383019 FREO- 0.318000 6 ЯИР- 0.2562665 PHASE- 101.3553 DAMP- -0.0320243 FREO- 0.357704 ? ppp- 0.2442874 PHASE- 103.9500 DPTF- -0.8119484 FREO- 0.480836 8 яг- 0.1313714 PHASE- 51.7046 D»P- -0.8738351 FREO- 0.451616 Прилипание. АМР — амплитуда, PHASE— фаза, DAMP—затухание, FREQ—частота.
СПЕКТРАЛЬНЫЙ АНАЛИЗ 43 Таблица 5 Список оценок параметров спектральных линий, полученных с помощью метода Прони Номер сигнала Оценки параметров 1 ЯТ* 8.1833*4 PHASE- 125.27?? ОЯР* 2 ЯР* 1.8259654 PHASE- 156.4714 DAMP- Э ЯР* 1.8396414 PHASE* 178.5962 РЯР* 4 ЯР* 8.6564418 WASE* 36.3465 РЯР* 5 ЯР* 6.2168961 PHASE* 65.7856 РЯР- 6 ЯР* 8,1166532 РНЙ6Е- 174.4637 РЯР* 7 ЯР* 6.1971165 PHASE* 95.5934 РЯР* в ЯР* 6.8276794 PHASE* 32.9251 РЯР* Примечание.*АМР — амплитуда, PHASE—фаза. DAMP — затухание, FREQ —частота. на рис. 16j. Три синусоидальные составляющие имеют очень острые максимумы на частотах, соответствую- щих этим синусоидам, а в высокочастотной части спек- тра присутствует широкий отклик. В табл. 4 представ- лены оценки реальных параметров, полученные с помощью обобщенного метода Прони. Действительные амплитуды синусоид на частотах 0,1, 0,2 и 0,21 равны соответственно 0,1, 1,0 и 1,0. Наиболее точные оценки мощностей и частот этих трех синусоид получены с помощью метода спектрального разложения, который является частным случаем метода Прони (рис. 16k). Этот результат неудивителен, поскольку данный ме- тод основан на оценке по методу наименьших квадра- тов, который предполагает использование синусои- дальной модели. Этот метод дает дискретный спектр, поэтому широкополосные процессы моделируются плохо, хотя в получаемом спектре и присутствуют спектральные линии, указывающие наличие энергии в данной области. Действительные оценки парамет- ров, полученные с помощью этого метода, приведены в табл. 5. Спектр, полученный с помощью метода максималь- ного правдоподобия, показан на рис. 161 и имеет весь- ма сглаженную форму. Он не позволяет разрешить две близкие по частоте синусоиды. Сглаженный характер ММП-оценки спектра, эквивалентной некоторой сред- ней оценке для всех АР-методов спектрального оцени- вания с порядком модели от 1 до 16, типичен для дан- ного метода. Если в результате спектрального оценивания не- обходимо получить более точные оценки частот «зашум- ленных» синусоид, а не просто оценить форму спектра, то можно использовать некоторые полученные не- давно Тафтсам и Кумаресаном результаты [244, 245, 134, 135], которые рассматривали возможности улуч- шения методов линейного предсказания, анализа соб- ственных значений и максимального правдоподобия для уменьшения дисперсии частотных оценок и уве- личения разрешения. Однако следует заметить, что они рассматривали только процессы, состоящие из синусоид и шума. IV. ДРУГИЕ ПРИМЕНЕНИЯ МЕТОДОВ СПЕКТРАЛЬНОГО ОЦЕНИВАНИЯ А. Введение В предыдущих разделах статьи были рассмотрены теория и приложения современных методов спектраль- ного анализа. Однако многие из представленных там 6* результатов с успехом .приложимы и для решения задач, отличных от задач спектрального оценивания. Поскольку эти приложения достаточно интересны специалистам во многих областях, то в данном разделе мы коротко остановимся на некоторых из них. Заметим лишь, что обсуждаемые здесь вопросы отнюдь не со- ставляют исчерпывающий перечень всех возможных приложений методов спектрального анализа, а только характеризуют наиболее типичные приложения в бо- лее или менее близких областях. В. Экстраполяция и интерполяция временных рядов Теоретические основы современных методов спек- трального анализа позволяют применить их и для других приложений. Одним из очевидных приложений является экстраполяция временных рядов с неиз- вестной СПМ. Если, например, временной ряд пред- ставляет собой АР(р)-процесс, то параметры опти- мального устройства линейного предсказания будут АР-параметрами. Оценки последних получаются по имеющимся данным так, как показано в подразд. II.Е. Если же временной ряд не является АР-процессом, то число параметров оптимального линейного устрой- ства предсказания будет в общем случае бесконечным. Теоретически с ростом числа параметров устройства предсказания величина ошибки экстраполяции долж- на уменьшаться. При наличии конечного набора дан- ных, по которым должны оцениваться параметры устройства предсказания, ошибку предсказания мож- но минимизировать посредством выбора наибольшего возможного порядка модели при наложенном ограни- чении, что параметры устройства предсказания долж- ны оцениваться точно [51]. _ Несмотря на то что методы, описанные в под- разд. П.Е, позволяют получить только одношаговое устройство предсказания, предсказанные значения от- счетов можно использовать, как если бы они были частью исходного набора данных, и продолжить экст- раполяцию к следующему отсчету (27, 291. Было даже -.Предсказание вперед Яп \ по отсчетам набора 1 \ р Предсказание назад \ \ по отсчетам набора 2 Набор 1 «Плохие» Набор 2 данные Рис. 17. Интерполяция «плохих» отсчетов данных. предложено использовать увеличенное множество ис- ходных и экстраполированных данных совместно с периодограммой или спектральной БТ-оценкой для улучшения разрешения [68]. Помимо экстраполяции с помощью устройства предсказания вперед и назад можно выполнить интерполяцию, как показано на рис. 17. Эта возможность имеет большую ценность для замены плохих отсчетов данных [176]. С. Предобеляющий фильтр Предобеляющий фильтр является естественным об- общением использования параметров, получаемых из спектральных оценок. Так, например, при спектраль-
44 ТИИЭР, т. 69, №11, ноябрь 1981 ном АР-анализе величина А (г), определяемая вы* ражением (2.91), является обеляющим фильтром. Вы- ход этого фйльтра (ошибка предсказания) представ- ляет собой белый шум, если наблюдаемый процесс является АР (р)-процессом и оцениваемые коэффи- циенты устройства предсказания действительно яв- ляются АР-параметрами. В случае, когда временной ряд не является АР-процессом, форма оценки спект- ральной плотности мощности выходного временного ряда будет все еще более плоской, чем СПМ входного процесса, и «приближенно» белой. Это свойство важ- но, в частности, при проектировании детекторов (об- наружителей) сигналов в присутствии окрашенного шума с неизвестным спектром. Типичным примером может служить обнаружение сигналов, отраженных от цели, в присутствии фонового отражения от облака пас- сивных отражателей. Оптимальный детектор представ- ляет собой предобеляющий фильтр с последующим фильтром, согласованным с сигналом на выходе пред- обеляющего фильтра [256]. Поскольку спектр сигна- лов, отраженных от облака пассивных отражателей, обычно изменяется во времени, параметры предобе- ляющего и согласованного фильтров должны обнов- ляться. Успех применения предобеляющей схемы будет зависеть от временного изменения спектра сиг- налов, отраженных от облака пассивных отражателей, и возможности оценивания параметров этого спектра, прежде чем они изменятся [301. Близко связанный подход состоит в предваритель- ном обелении временного ряда для уменьшения сме- щения стандартных спектральных оценок. Можно показать [107], что z-l/2Af £[&(/)] = I 9x(v)W(f-v)dv, (4.1) •'-l/jAf где Px(f) — спектральная оценка БТ-типа, 9*,$) — истинная СПМ, W (f) — спектральное окно, требуемое для умень- шения дисперсии оценки [ \tW(f)df= =11. Если значение 5*х(/) слабо зависит от частоты, то /-1/2АГ E[?x(f>] * 9Х(П W(f-v)dv=3>x(f). (4.2) •'-1/2ДГ Таким образом, для уменьшения смещения следует попытаться обелить входные данные [213, 232]. После обеления находится БТ-оценка 5\(f) профильтрован- ного временного ряда еп. Поскольку еп является выхо- дом предобеляющего фильтра A (z), ее СПМ равна 5*,(f)=|A (ехр [/2лД^])|>5*ж(/). Спектральная оценка хп будет в этом случае определяться выражением $ (/) =------, (4.3) х |А(ехр[/2я/Д1])|2 где предполагается использование предобеляющего фильтра только с одними нулями. Интересно заметить, что данный подход дает спектральную оценку, которая является стандартной спектральной АР-оценкой, в которой СПМ белого шума ojAf заменена оценкой СПМ остаточного временного ряда 5*e(f). D. Сжатие полосы частот Одной из важных проблем при анализе речевых сигналов является проблема сжатия полосы частот. Если избыточность речи можно уменьшить, то многие речевые сигналы можно передавать по каналу с фикси- рованной шириной полосы частот или хранить в за- поминающих устройствах большой емкости. Один из наиболее часто применяемых методов основан на при- менении дифференциальной импульсно-кодовой мо- дуляции (ДИКМ) [104]. Основу ДИКМ составляет передача только той информации, которую нельзя предсказать; такую информацию часто называют об- новлением процесса [112]. В действительности, если речевой сигнал точно предсказуем по набору пред- шествующих отсчетов то приемник, получив эти от- счеты, может точно восстановить полный сигнал (пред- полагается, что шум в канале полностью отсутствует). Иными словами, передача может прерываться! На практике при использовании ДИКМ отсчеты речевого сигнала анализируются в передатчике с целью опре- деления параметров устройства предсказания. В этом случае передаются только временной ряд ошибки пред- сказания и параметры устройства предсказания, что позволяет восстановить речевой сигнал в приемнике. Уменьшение ширины полосы частот оказывается воз- можным, поскольку дисперсия временного ряда ошиб- ки предсказания меньше дисперсии речевого сигнала [251], т. е. a* = Rxxm П (1 - К?) < йхх(0). (4.4) i-1 Следовательно, для кодирования остаточного времен- ного ряда достаточно нескольких уровней квантова- ния. Заметим, что для максимального сжатия полосы частот параметры устройства предсказания должны непрерывно обновляться вследствие статистического характера изменений речевого сигнала; иными сло- вами, должны отслеживаться вокализованные и не- вокализованные участки речи. Методом, обеспечивающим наибольшее уменьшение полосы частот, является.метод кодирования на основе линейного предсказания (КЛП), в котором для пред- ставления речевого сигнала используется АР-модели- рование [145]. Если речевой сигнал допускает точное моделирование в виде выхода чистополюсного фильтра, возбуждаемого чисто белым шумом в случае невокали- зованной речи или некоторой импульсной последова- тельностью в случае вокализованной речи, то речевой сигнал можно представить с помощью небольшого числа параметров. Таким образом, необходимо пере- давать или запоминать параметры модели или период повторения импульсной последовательности. Следова- тельно, синтез речевого сигнала можно выполнить на основе соответствующей модели для каждого звука речи. Е, Спектральное сглаживание Стандартная периодограмма и БТ-анализ приводят к спектральным оценкам, которые характеризуются большим числом «холмов и впадин», поскольку пре- образование Фурье для случайного процесса с нуле- вым средним значением также является случайным про- цессом с нулевым средним значением. Взвешивание
СПЕКТРАЛЬНЫЙ АНАЛИЗ 45 значений автокорреляционной функции или сглажива- ние спектральным окном будет существенно умень- шать флюктуации, но все же их не устраняет. Для сглаживания этих флюктуаций можно использовать устройство спектральной АР-оценки, поскольку на спектральную АР-оценку р-го порядка налагается ограничение, состоящее в том, что она должна иметь р или менее максимумов (или минимумов). При малых р получается сглаженная оценка. Покажем теперь, что спектральная АР-модель может точно характеризо- вать максимумы периодограммы, но не ее минимумы [26, 76, 145]. Рассмотрим оценку АР-параметров, най- денную посредством минимизации критерия ошибки (2.64) при использовании матрицы Х$Хг (автокорре- ляционных нормальных уравнений): 6=^7 f 1е„12, (4.5) 4* П«-«О где полагается, что Хп, n=0, 1, ..., W—1, доступно наблюдению, а хп=0 находится вне этого интервала. Тогда, согласно теореме Парсеваля. 1 гi,2At &=—— |£(exp[/2ff/At])|2d/, где £(z) = 4r £ enz~n . п * Поскольку Х(г)=Дг £ x„z-n, п = то Е становится равным 1 Г1/2Af Е~ТГ~ I |Х(ехр [;'2я/Дг])|2 х х|Л(ехр [/21Г/ДГ] )|2 df = -фЫ i / = о2Д^ / —|Х(ехр[/2я/Дг])|2/^(/)4Г, Л/2ДГ / (4.6) где IX (exp [J2nfAfl)|’M(Д( — периодограмма, а £ж(/)= =<т‘Д(/|Л (ехр [/2л/Д/])|*. Следовательно, когда т4"|Х(ехр [/2я/Дг] )|2 Частота. кГц Pic. 18. СПМ на основе периодограммы я сглаженная АР(28)- оценка СПМ речевого сигнала [144]. г x(t.O) Волновой фронт плоской волны х(1Дх) x(t,2Ax) J Датчики x(t,[M-llAx) Рис. 19. Геометрия линейной решетки. IX (ехр [/2л/Д/])|’/ЛГД( мало, то существует лишь не- значительный вклад в эту ошибку, поэтому согласова- ние не обязательно. Этот результат говорит о том, что спектральная АР-оценка согласует максимумы, но не минимумы периодограммы. Если необходимо охарак- теризовать максимумы спектра, то необходимо лишь взять обратное преобразование Фурье для отыскания первых р+1 значений автокорреляционной функции, которые затем используются для определения АР(р)- модели. Пример такого подхода показан на рис. 18 для периодограммы речевого сигнала. F. Формирование луча При формировании луча необходимо получить оценку пространственной структуры случайного про- странственного поля. Если для получения отсчетов этого поля используется линейная решетка датчиков, то получается векторный временной ряд {x(t, 0), x(t, Дх), . . ., x(t, [М—1]Дх)}, где x(t, (Дх) — непре- рывный сигнал на выходе i-ro датчика, —1, а М — число датчиков решетки (рис. 19). Поле можно представить как [42] x(t, /Дх) = f ГФ(/, кх) ехр - kxibx)] df dkx, (4.7) т. e. как сумму бесконечного числа монохроматиче- ских плоских волн со случайными амплитудами ¥ (/, kx). Частота изменения поля во времени обозначена здесь через /, а пространственная частота в направле- нии оси х — через kx. Составляющая пространствен- ной частоты kx, характеризующая число волн, являет- ся величиной, обратной длине волны монохроматиче- ской плоской волны, распространяющейся в направ- лении оси х. Поскольку jfeI=(f/c)sin 0, где 0 — угол, показанный на рис. 19, то Е [|¥ (7, fex)|a] представляет собой мощность на частоте f, принимаемую с направ- ления, характеризуемого углом 0. В соответствии с (4.7) обратное преобразование Фурье равно М-1! г \ '^(f,kx)= £ ( |х(Г, /Дх) ехр (-/2я/г) dr) ехр (/2mtxiAx) = г «о w / М-1 “ Z XzGAx) ехр (/21гЛжгДх). (4.8) 1 = 0 Выражение (4.8) представляет собой преобразование Фурье между пространственным «временным» рядом Xz(iAx), где Дх — расстояние между отсчетами, и его спектром ¥ (f, kx). Если поле в пространстве полагает- ся однородным, т. е. Е[x(t, iДх) х*(Г,/Дх)] = /(Г, [/ - }] Дх), (4.9) то пространственный «временной ряд» является ста- ционарным в широком смысле, и оценка £(|¥(f, Ml’)
46 ТИИЭР, т. 69, № 11, ноябрь 1981 для всех kx на заданной временной частоте f будет ана- логична одномерной временной спектральной оценке мощности. В данном случае применим любой из опи- санных в статье методов, если временная запись дан- ных заменяется пространственной. записью данных {х(<0, 0), х(1в, Дх), . . ., x(f0, [Af—1]Дх)} в некоторый момент времени ta. Заметим, что в пространственной области возможно некоторое дополнительное усред- нение, которое невозможно во временной области. На- пример, оценку пространственной «автокорреляции» можно выбрать равной 1 N-1 М-к-1 Ах(*:Дх)= — — Z Е х*(ЛД1,1Дх)х N(M -к) „ _0 хх(пД1, (i+Л)Дх), fc>0, (4.10) где х(пД(, (Дх) полагается стационарным во времени на интервале 0^nAt^(M—1)А/. Эта оценка включает операцию дополнительного временного усреднения. G. Решетчатые фильтры Минимально-фазовый решетчатый фильтр, описан- ный в подразд. П.Е, обладает тем свойством, что его коэффициенты ограничены по величине единицей. Это — весьма ценное свойство в том случае, когда эти коэффициенты необходимо проквантовать для передачи или хранения в запоминающем устройстве [257]. Решетчатую структуру можно использовать для синтеза устойчивого минимально-фазового КИХ- фильтра, устойчивого чистополюсного БИХ-фильтра или устойчивого чистонулевого БИХ-фильтра, нули которого расположены внутри единичного круга в г-плоскости. Н. Другие применения К другим возможным применениям описанных в статье методов относятся: 1) коррекция в цифровых системах связи [147], 2) анализ переходных процессов [252, 254], 3) проектирование цифровых фильтров [183], 4) обратная свертка на основе предсказания [209], 5) кепстральный анализ [139]. За подробностями интересующийся читатель отсы- лается к указанным работам. V. ЗАКЛЮЧЕНИЕ Современные методы спектрального анализа осно- ваны на моделировании данных с помощью неболь- шого числа (набора) параметров. В случае, когда модель точно описывает данные, можно получить спектральные оценки, характеристики которых пре- восходят характеристики классических оценок на основе периодограммы или метода Блэкмана — Тью- ки. Улучшение характеристик проявляется в улучше- нии разрешения и отсутствии боковых лепестков. Сле- дует особо подчеркнуть тот факт, что помимо точной модели данных спектральные оценки должны основы- ваться на хороших оценках параметров модели. Обыч- но это приводит к оценкам параметров на основе мак- симального правдоподобия. Если модель не совсем точна, как, например, в случае использования АР- модели для АР-процесса с аддитивным шумом наблю- дения, то получаются плохие (смещенные) спектраль- ные оценки. Если модель точна, но используются плохие статистические оценки параметров, как, на- пример, в случае спектральных АРСС-оценок на осно- ве модифицированных уравнений Юла — Уокера, то также получаются плохие (с большой дисперсией) спектральные оценки. Для получения спектральных АР-оценок по методу максимального правдоподобия разработаны эффектив- ные в вычислительном отношении процедуры. Эти ме- тоды обычно требуют ненамного больше вычислений, чем стандартные методы спектрального оценивания на основе преобразования Фурье. Однако спектральное АРСС-оценивание по методу максимального правдо- подобия включает решение нелинейных уравнений, поэтому для этого метода спектрального оценивания пока еще не существует эффективных вычислительных процедур. Поскольку благодаря своей большей устой- чивости спектральные АРСС-оценки более желатель- ны, чем спектральные АР-оценки в тех случаях, когда характеристики данных не известны, будущие исследо- вания должны, быть направлены на повышение вычис- лительной эффективности спектрального АРСС-оце- нивания по методу максимального правдоподобия. В настоящее время предложено большое число ал- горитмов спектрального оценивания, лишь небольшая часть которых рассмотрена в данной обзорной статье. К сожалению, статистические свойства в случае ко- нечных записей данных проанализированы лишь у нескольких алгоритмов. Сравнение различных кон- курирующих алгоритмов было основано на ограни- ченном числе результатов машинного моделирования, поэтому результаты такого сравнения могут быть оши- бочными. По этой причине будущие исследования должны быть также направлены на получение более полных статистических описаний современных мето- дов спектрального анализа. В заключение следует заметить, что современные методы спектрального оценивания при правильном их использовании являются чрезвычайно ценным сред- ством для анализа данных. Цель авторов настоящей статьи состояла в представлении различных методов с единой точки зрения и единой терминологии. Мы надеемся, что такой подход поможет пользователям при выборе метода спектрального оценивания, наи- более полно удовлетворяющего их нуждам. ЛИТЕРАТУРА Работы, которые, в приведенном ниже списке лите- ратуры помечены звездочкой, помещены в сборнике Modern. Spectrum Analysis («Современный спектраль- ный анализ»), New York, 1978, выпущенный под об- щей редакцией Д. Дж. Чайлдерса издательством IEEE Press в серии сборников избранных статей. [ 1 ] Н. Akaike, “Fitting autoregressive models for prediction/* Ann. Inst. Statist. Math.,vo\. 21, pp. 243-247, 1969. [2] , “Power spectrum estimation through autoregression model fitting/* Ann. Inst. Statist. Math., vol. 21, pp. 407-419, 1969. [3] —, “On a semi-automatic power spectrum estimation proce- dure,” in Proc. 3rd Hawaii Int. Conf. System Sciences, Part 2, pp. 974-977, 1970. [4] —, “Statistical predictor identification/* Ann. Inst. Statist. Math., vol. 22, pp. 203-217, 1970. 15] , “Autoregressive model fitting for control/* Ann. Inst. Statist. Math., voL 23, pp. 163-180, 1971. [6] —, “Use of an information theoretic quantity for statistical model identification.” hi Proc^ 5th Hawaii Int. Conf. System Sciences, pp. 249-250, Jan. 11^13/1971.
СПЕКТРАЛЬНЫЙ АНАЛИЗ 47 * [71 -, “A new look at the statistical model identification,” IEEE Trans. Au tom. Contr., vol. AC-19, pp. 716-723, Dec. 1974. • [81 M. A. Alam, “Adaptive spectral estimation,” in Proc. Joint Automat. Contr. Conf., vol. 1, San Francisco, CA, June 22-24, 1977, pp. 105-112. * [9] N. O. Andersen, “On the calculation of filter coefficients for maximum entropy spectral analysis,” Geophys., vol. 39, pp. 69-72, Feb. 1974. [10] Андерсон H. Замечания о характеристиках алгоритмов максимальной энтропии. ТИИЭР, 1978, т. 66, № 11, с. 342— 343. (11] Т. W. Anderson, “Estimation by maximum likelihood in auto- regressive moving average models in the time and frequency domains,” Dep. Statistics, Tech. Rep. 20, Contract NR-042-034, Stanford University, Stanford, CA, June 1975. 112] H. Babic and G. C. Temes, “Optimum low-order windows for discrete Fourier transform systems,” IEEE Trans. Acoustics, Speech, Signal Process., vol. ASSP-24, pp. 512-517, Dec. 1976. (13) A. B. Baggeroer, “Confidence intervals for regression (MEM) spectral estimates,” IEEE Trans. Inform. Theory, vol. IT-22, pp. 534-545, Sept. 1976. [14] Бэнес. Замечания к сообщению «Гармонический анализ с помощью метода калмановской фильтрацию. ТИИЭР, 1973, т. 61. №12. с. 111. [15] Т. Е. Barnard, “The maximum entropy spectrum «nd the Burg technique,” Tech. Rep. 1, Advanced Signal Processing, Texas Instruments, prepared for Office of Naval Research, ALEX (03)- TR-75-01, June 25, 1975. (16] I. Barrodale and R. E. Erickson, “Algorithms for least-squares linear prediction and maximum entropy spectral analysis—Part I: Theory and Part II: FORTRAN Program,” Geophys., vol. 45, pp. 420-446, Mar. 1980. [17] M. S. Bartlett, “Periodogram analysis and continuous spectra,” Biometrika, vol. 37, pp. 1-16, June 1950. (18] M. S. Bartlett and J. Medhi, “On the efficiency of procedures for smoothing periodograms from time series with continuous spectra,” Biometrika, vol. 42, pp. 143-150, 1955. [19] G. D. Bergland, “A guided tour of the fast Fourier transform,” IEEE Spectrum, vol. 6, pp. 41-52, July 1969. (20] A. J. Berkhout and P. R. Zaanen, “A comparison between Weiner filtering, Kalman filtering, and deterministic least squares estimation,” Geophysical Prospecting, vol. 24, pp. 141-197, Mar. 1976. [21] J. G. Berryman, "Choice of operator length for maximum en- tropy spectral analysis,” Geophys., vol. 43, pp. 1384-1391, Dec. 1978. [22] S. Bertram, “On the derivation of the fast Fourier transform,” IEEE Trans. Audio Electroacoust., vol. AU-18, pp. 55-58, Mar. 1970. [23] --, “Frequency analysis using the discrete Fourier transform,” IEEE Trans. Audio Electroacoust., vol. AU-18, pp. 495-500, Dec. 1970. *[24] C. Bingham, M. D. Godfrey, and I. W. Tukey, “Modem tech- niques of power spectrum estimation,” IEEE Trans. Audio Electroacoust., vol. AU-15, pp. 56-66, June 1967. [25] R. B. Blackman and J. W. Tukey, The Measurement of Power Spectra From the Point of View of Communications Engi- neering. New York: Dover, 1959. [26] S. F. Boll, “A priori digital speech analysis,” Ph.D. dissertation, Dep. Elec. Eng., Utah Univ., Mar. 1973. [27] S. B. Bowling, “Linear prediction and maximum entropy spec- tral analysis for radar applications,” M.I.T. Lincoln Lab., Project Rep. RMP-122 (ESD-TR-77-113), May 24, 1977. [28] P. Bloomfield, Fourier Analysis oif Time Series: An Introduction. NewYork: Wiley, 1976. [29] S. B. Bowling and S. Lai, “The use of linear prediction for the interpolation and extrapolation of missing data prior to spectral analysis,” in Proc. RADC Workshop on Spectrum Estimation, pp. 39-50, Oct. 3-5, 1979. [30] D. E. Bowyer, P. K. Rajasekaran, and W. W. Gebhart, “Adaptive clutter filtering using autoregressive spectral estimation,” IEEE Trans. Aerospace Elec. Syst., vol. AES-15, pp. 538-546, July 1979. [31] Бокс Дж., Дженкинс Г. Анализ временных рядов: Прогноз и управление. Выл. 1 и 2. М.: Мир, 1974. [32] Бригхэм, Морроу. Быстрое преобразование Фурье. ТИИЭР, 1967, т. 55, № 10. с. 21—29. [33] Е. Oran Brigham, The Fast Fourier Transform. Englewood Cliffs, NJ: Prentice-Hall, 1974. [34] Брнллинджер. Фурье-анализ стационарных процессов. ТИИЭР, 1974, т. 62, № 12, с. 15—33. [35] S. Bruzzone and М. Kaveh, “On some suboptimum ARMA spectral estimates,” IEEE Trans. Acoustics, Speech, Signal Process., vol. ASSP-28, pp. 753-755, Dec. 1980. [36] H. P. Bucker, “Comparison of FFT and Prony algorithms for bearing estimation of narrow-band signals in a realistic ocean environment,” J. Acoust. Soc. Amer., vol. 61, pp. 756-762, Mar. 1977. •[37] J. P. Burg, “Maximum entropy spectral analysis,” in Proc. 37th Meeting Society of Exploration Geophysicists (Oklahoma City, OK), Oct. 31, 1967. •[38] -----, “A new analysis technique for time series data,” NATO Advanced Study Institute on Signal Processing with Emphasis on Underwater Acoustics, Enschede, The Netherlands, Aug. 12-23, 1968. [39] ----, “New concepts in power spectral estimation,” in Proc. 40th Annu. Int. Society of Exploration Geophysicists (SEG) Meeting (New Orleans, LA) Nov. 11, 1970. •[40] -----, “The relationship between maximum entropy and maxi- ' mum likelihood spectra,” Geophys., vol. 37, pp. 375-376, Apr. 1972. [41] , “Maximum entropy spectral analysis,” Ph.D. dissertation, Dep. Geophysics, Stanford Univ., Stanford, CA, May 1975. [42] Кейпон. Пространственно-временндй спектральный анализ с высоким разрешением. ТИИЭР, 1967, т. 57, № 8, с. 69— 79. V [43] Картер Дж. К., Наттолл А. X., Юэнь Ч. К- О методе взве- шенного усреднения перекрывающихся сегментов для оце- нивания спектра. ТИИЭР, 1980, т. 68, № 10, с. 216—218. [44] G. С. Carter and А. Н. Nuttall, «А brief summary of a ge- neralized framework for power spectral estimations. Sig- nal Processing, vol. 2, pp. 387—390, Oct. 1980. [45] W. Y. Chen and G. R. Stegen, “Experiments with maximun entropy power spectra of sinusoids,” J. Geophysical Res., vol. 79, pp. 3019-3022, July 10, 1974. [46] D. Childers and A. Durling, Digital Filtering and Signal Pro- cessing. St. Paul, MN: West Publishing, 1975. [47] J. C. Chow, “On the estimation of the order of a moving-average process,” IEEE Trans. Automat. Contr., vol. AC-17, pp. 386- 387, June 1972. [48] , “On estimating the orders of an autoregressive moving- average process with uncertain observations,” IEEE Trans. Automat. Contr., vol. AC-17, pp. 707-709, Oct. 1972. [49] C. W. Chuang and D. L. Moffatt, “Natural resonances of radar targets via Peony’s method and target discrimination,” IEEE Trans. Aerospace Electron. Syst., vol. AES-12, pp. 583-589, Sept. 1976. [50] J. F. Claerbout, Fundamentals of Geophysical Data Processing. NewYork: McGraw-Hill: 1976. [51] W. S. Cleveland, “Htting time series models for prediction,” Technometrics, pp. 713-723, Nov. 1971. [52] W. T. Cochran et al., “What is the fast Fourier transform?,” IEEE Trans. Audio Electroacoust., vol. AU-15, pp. 45-55 June 1967. [53] J. W. Cooley and J. W. Tukey, “An algorithm for machine cal- culation of complex Fourier series,” Math. Comput., vol. 19, pp. 297-301, Apr. 1965. [54] J. W. Cooley, P.A.W. Lewis, and P. D. Welch, “The fast Fourier transform algorithm: Programming consider-tions in the cal- culation of sine, cosine, and Laplace transforms,” J. Sound Vibration, vol. 12, pp. 315-337, July 1970. [55] , “The application of the fast Fourier transform algorithm to the estimation of spectra and cross-spectra,” J. Sound Vibra- tion, vol. 12, pp. 339-352, July 1970. [56] G. Cybenko, “Round-off error propagation in Durbin's, Levin- son'S, and Trenchb Algorithms,” Rec. 1979 IEEE Int. Conf. Acoustics, Speech, and Signal Processing, pp. 498-501. [57] Датта A. K-, ХайкннС. Замечания к сообщению «Комп- лексная форма метода максимальной энтропии для оценки спектральной плотности». ТИИЭР, 1977, т. 65, № 8, с. 134— 135. [58] W. J. Done, “Estimation of the parameters of an autoregressive process in the presence of additive white noise,” Computer Science Dep., Univ. Utah, Rep. UTEC-CSC-79-021, Dec. 1978 (order NTIS no. ADA 068749). V [59] Дуброфф. Эффективная автокорреляционная функция спек* тров максимальной энтропии. ТИИЭР, 1975, т. 63, № 11, с. 96—9Z- [60] J. Durbin, “The fitting of time series models,” Rev. Inst. Int. de Stat., vol. 28, pp. 233-244, 1960. [61] A. Eberhard, “An optimal discrete window for the calculation of power spectra,” IEEE Trans. Audio Electroacoust., vol. AU-21, pp. 37-43, Feb. 1973. •[62] J. A. Edward and M. M. Fltelson, “Notes on maximum entropy processing,” IEEE Trans. Inform. Theory, vol. IT-19, pp. 232- 234, Mar. 1973. [63] D. C. Farden, “Solution of a Toeplitz set of linear equations,” IEEE Trans. Antennas Propagat., vol. AP-24, pp. 906-908, Nov. 1976.
48 ТИИЭР, т. 69, № 11, ноябрь 1981 [64] Р. F. Fougere, Е. J. Zawalick, and Н. R. Radoski, “Spontaneous line splitting in maximum entropy power spectrum analysis/* Physics Earth Planetary Interiors, vol. 12, pp. 201-207, Aug. 1976. [65] P. F. Fougere, “A solution to the problem of spontaneous line splitting in maximum entropy power spectrum analysis,” J, Geophysical Res., vol. 82, pp. 1051-1054, Mar. 1, 1977. [66] B. Friedlander, M. Mort, T. Kailath, and L. Ljung, “New inver- sion formulas for matrices classified in terms of their distance from Toeplitz matrices/* Linear Algebra and Its Applications, vol. 27, pp. 31-60, 1979. [67] O. L. Frost, “Power spectrum estimation,** in Proc. 1976 NATO Advanced Study Institute on Signal Processing with Emphasis on Underwater Acoustics, Portovener (LaSpezia), Italy, Aug. 30-Sept 11, 1976. [68] O. L. Frost and T. M. Sullivan, “High resolution two dimensional spectral analysis/* in Conf. Rec. 1979 ICASSP, pp. 673-676, 1979. [69] J. Fryer, M. E. Odegard, and G. H. Sutton, “Deconvolution and spectral estimation using final prediction error/* Geophys., vol. 40, pp. 411-425, June 1975. V [70] ГеЙбриел У. Ф. Спектральный анализ и методы сверхраз- решения с использованием адаптивных решеток. ТИИЭР, 1980, т. 68, №6/с. 19—32. [71] W. Gersch, “Spectral analysis of EEG’s by autoregressive de- composition of time series,** Math. Biosci., vol. 7, pp. 205-222, 1970. [72] , “Estimation of the autoregressive parameters of a mixed autoregressive moving-average time series/* IEEE Trans. Auto- mat. Contr., vol. AC-15, pp. 583-588, Oct. 1970. [73] W. Gersch and D. R. Sharpe, “Estimation of power spectra with finite-order autoregressive models/* IEEE Trans. Automat. Contr., vol. AC-18, pp. 367-369, Aug. 1973. [74] W. Gersch and J. Yonemoto, “Automatic classification of EEG’s: A parametric model new features for classification approach/* in Proc. 1977 Joint Automatic Control Conf. (San Francisco, CA), June 22-24, 1977, pp. 762-769. [75] J. Gibson, S. Haykin, and S. B. Kesler, “Maximum entropy (adaptive) filtering applied to radar clutter/* in Rec. 1979 IEEE Int. Conf. Acoustics, Speech, and Signal Processing, pp. 166-169. [76] J. D. Gibson et al., “Digital speech analysis using sequential estimation techniques/* IEEE Trans. Acoustics, Speech, Signal Process., pp. 362-369, Aug. 1975. [77] T. H. Glisson, С. I. Black, and A. P. Sage, “The digital computa- tion of discrete spectra using the fast Fourier transform, ” IEEE Trans. Audio Electroacoust., vol. AU-18, pp. 271-287, Sept. 1970. [78] J. Grandell, M. Hamrud, and P. Toll, “A remark on the corre- spondence between the maximum entropy method and the autoregressive models/* IEEE Trans. Inform. Theory, vol. IT-26, pp. 750-751, Nov. 1980. [79] D. Graupe, D. J. Krause, and J. B. Moore, “Identification of autoregressive moving-average parameters of time series,” IEEE Trans. Automat. Contr., vol. AC-20, pp. 104-107, Feb. 1975. [80] A. H. Gray, Jr. and D. Y. Wong, “The Burg algorithm for LPC speech analysis/synthesis,” IEEE Trans. Acoustics, Speech, Signal Process., vol. ASSP-28, pp. 609-615, Dec. 1980. [81] R. M. Gray, “Toeplitz and circulant matrices: A review/* Infor- mation Systems Laboratory, Center for Systems Research, Stan- ford University, Tech. Rep. No. 6502-1, June 1971. [82] O. Grenander and G. Szego, Toeplitz Formsand Their Applica- tions. Berkeley, CA: Univ. California Press, 1958. *[83] L. J. Griffiths, “Rapid measurement of digital instantaneous frequency/* IEEE Trans. Acoustics, Speech, Signal Process., vol. ASSP-23, pp. 207-222, Apr. 1975. [84] L. J. Griffiths and R. Prieto-Diaz, “Spectral analysis of natural seismic events using autoregressive techniques/’ IEEE Trans. Geosci. Electron., vol. GE-15, pp. 13-25, Jan. 19-77. [85] L. J. Griffiths, “Adaptive structures for multiple-input noise canceling applications/* in Conf. Rec. ICASSP 79, pp. 925-928. [86] P. R. Gutowski, E. A. Robinson, and S. Treitel, “Novel aspects of spectral estimation,” in Proc. 1977 Joint Automatic Con- trol Conf., vol. 1 (San Francisco, CA), June 22-24, 1977, pp. 99-104. *[87] ----, spectral estimation: Fact or fiction?/* IEEE Trans. Geosci. Electron., vol. GE-16, pp. 80-84, Apr. 1978. [88] C. S. Hacker, “Autoregressive and transfer function models of mosquito populations,’’ in Time Series and Ecological Processes, H. M. Shugat, Ed. Boulder, CO: SIAM, 1978, pp. 294-303. [89] Хеннан Э. Многомерные временные ряды. М.: Мир, 1974. 190] Хэррис Ф. Дж. Использование окон при гармоническом анализе методом дискретного преобразования Фурье. ТИИЭР, 1978, т. 66, № 1, с. 60—96. [91] Хайкин, Кеслер. Комплексная форма метода максимальной энтропии для оценки спектральной плотности. ТИИЭР, 1976, т. 64, №5, с. 313—314. [92] S. S. Haykin, Ed., Nonlinear Methods of Spectral Analysis. New York: Springer-Verlag, 1979. [93] Хайкин С., Рейлли Дж. Моделирование отклика линейной антенной решетки на падающие плоские волны на основе смешанного процесса авторегрессии — скользящего сред- него. ТИИЭР, 1980, т. 68, № 5, с. 92-93. [94] R. W. Herring, “The cause of line splitting in Burg maximum- entropy spectral analysis,” IEEE Trans. Acoustics, Speech, Signal Process., vol. ASSP-28, pp. 692-701, Dec. 1980. [95] F. B. Hildebrand, Introduction to Numerical Analysis. New York: McGraw-Hill, 1956, ch. 9. [96] S. Holm and J. M. Hovem, “Estimation of scalar ocean wave spectra by the maximum entropy method,” IEEE J. Ocean Eng., vol. OE-4, pp. 76-83, July 1979. [97] H. Holtz, “Prony’s method and related approaches to exponen- tial approximation,” Aerospace Corp., Rep. ATR-73(99901-5. June IS, 1973. [98] M. Hsu, “Maximum entropy principle and its application to spectral analysis and image reconstruction,” Ph.D. dissertation, Ohio State Univ., 1975. [99] F. M. Hsu and A. A. Giordano. “Line tracking using autore- gressive spectral estimates,” IEEE Trans. Acoustics, Speech, Signal Process., vol. ASSP-25,pp. 510-519, Dec. 1977. [100] M. Huzii, “On a spectral estimate obtained by an autoregressive model fitting,” Ann. Inst. Statist. Math., vol. 29, pp. 415-431, 1977. [101] L. B. Jackson and F. K. Soong, “Observations on linear estima- tion,” in Rec. 1978 IEEE Int. Conf. Acoustics, Speech, and Signal Processing, pp. 203-207. [102] L. B. Jackson, D. W. Tufts, F. K. Soong, and R. M. Rao, “Fre- quency estimation by linear prediction,” in Rec. 1978 IEEE Int. Conf. Acoustics, Speech, and Signal Processing, pp. 352-356. [103] A. K. Jain and S. Ranganath, “Extrapolation algorithms for discrete signals with a application in spectral estimation,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, pp. 830-845, Aug. 1981. / [Khj Джайант. Цифровое кодирование речевых сигналов. Кван- тизаторы для ИКМ, ДИКМ и ДМ. ТИИЭР, 1974, т. 62, № 5, с. 83—107. [105] G. М. Jenkins, “General considerations in the analysis of spec- tra,” Technometrics, vol. 3, pp. 133-166, May 1961. [106] , “A survey of spectral analysis,” Appl. Stat., vol. 14, pp. 2- 32, 1965. [107] G. M. Jenkins and D. G. Watts, Spectral Analysis and Its Appli- cations. San Francisco, CA: Holden-Day, 1968. [108] S. J. Johnsen and N. Andersen, “On power estimation in maxi- mum entropy spectral analysis,” Geophys., vol. 43, pp. 681- 690, June 1978. [ 109 ] R. H. Jones, “A reappraisal of the periodogram in spectral analy- sis,” Technometrics, vol. 7, pp. 531-542, November 1965. •[110 ] --, “Identification and autoregressive spectrum estimation,” IEEE Trans. Automat. Contr., vol. AC-19, pp. 894-898, Dec. 1974. [Ill] , “Autoregression order selection,” Geophys., vol. 41, pp. 771-773, Aug. 1976. [112] T. Kailath, “A view of three decades of linear filtering theory,” IEEE Trans. Informa. Theory, vol. IT-20, pp. 146-181, Mar. 1974. [113] R. P. Kane and N. B. Trivedi, “Effects of linear trend and mean value on maximum entropy spectral analysis,” Institute de Pes- quisas Espaciais, Rep. INPE-1568-RPE/069, Sao Jose dos Com- pos, Brasil (Order NTIS no. N79-33949). [114] R. P. Kane, “Maximum entropy spectral analysis of some artifi- cial. samples,” J. Geophys. Res., vol. 84, pp. 965-966, Mar. 1, 1979. [115] R. L. Kashyap, “Inconsistency of the AIC rule for estimating the order of autoregressive models,” IEEE Trans. Automa. Contr., vol. AC-25, pp. 996-998, Oct. 1980. [116] M. Kaveh and G. R. Cooper, “An empirical investigation of the properties of the autoregressive spectral estimator,” IEEE Trans. Informa. Theory, vol. IT-22, pp. 313-323, May 1976. [ 117 ] M. Kaveh, “High resolution spectral estimation for noisy signals,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-27, pp. 286-287, June 1979. [118] S. M. Kay, “Improvement of autoregressive spectral estimates in the presence of noise,” in Rec. 1978 Int. Conf. Acoustics, Speech, and Signal Processing, pp. 357-360. [119] , “Maximum entropy spectral estimation using the analytical signal,” IEEE Trans. Acoust., Speech Signal Process., vol. ASSP- 26, pp. 467-469, Oct. 1978. [120] S. M. Kay and S. L. Marple, Jr., “Sources of and remedies for spectral line splitting in autoregressive spectrum analysis,” in Rec. 1979 Int. Conf. Acoustics, Speech, and Signal Processing, pp.151-154.
СПЕКТРАЛЬНЫЙ АНАЛИЗ 49 [121] S. М. Kay, “The effects of noise on the autoregressive spectral estimator,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-27, pp. 478-485, Oct. 1979. [122] ----, “Autoregressive spectral analysis of narrowband processes in white noise with application to sonar,” Ph.D. dissertation, Georgia Institute of Technology, Mar. 1980. [ 123] , “Noise compensation for autoregressive spectral estimates,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-28, pp. 292-303, June 1980. [124] R. J. Keeler, "Uncertainties in adaptive maximum entropy fre- quency estimators,” NOAA Tech. Rep. ERL-105-WPL53, Feb. 1979 (Order NTIS no. N79-32489). [125] S. B. Kesler and S. S. Haykin, “The maximum entropy method applied to the spectral analysis of radar clutter,” IEEE Trans. Inform. Theory, vol. IT-24, pp. 269-272, Mar. 1978. [126] ----, “Maximum entropy estimation of radar clutter spectra,” in Natl. Telecommunications Conf. Rec. (Birmingham, AL), pp. 18.5.1-18.5.5, Dec. 3-6, 1978. [127] A. Ya. Khinchin, “Korrelationstheorie der Stationaren Sto- chastischen Prozesse,” Math. Annalen, vol. 109, pp. 604-615, 1934. [128] W. R. King, “Maximum entropy spectral analysis in the spatial domain,” Rome Air Development Center: Tech. Rep. RADC- TR-78-160, July 1978 (Order NTIS no. ADA 068558). [129] J. F. Kinkel et al., “A note on covariance-invariant digital filter design and autoregressive moving average spectrum analysis,” IEEE Trans. Acoust., Speech, Signal Process., pp. 200-202, Apr. 1979. [130] H. Kobatake et al., “Linear predictive coding of speech signals in a high ambient noise environment,” in Conf. Rec. 1978ICASSP, pp. 472-475. [131] I. S. Konvalenka er al., “Simultaneous estimation of poles and zeros in speech analysis and ITIF-Iterative inverse fitting algo- rithm,” IEEE Trans. Acoust., Speech, Signal Process., pp. 485- 491, Oct. 1979. [132] L. H. Koopmans, The Spectral Analysis of Time Series. New York: Academic Press, 1974. [133] F. Kozin and F. Nakajima, “The order determination problem for linear time-varying AR models,” IEEE Trans. Automat. Contr., vol. AC-25, pp. 250-257, Apr. 1980. [134] Кумаресан P., Тафтс Д. У. Улучшенный метод спектраль- ного разрешения, III : Эффективная реализация. ТИИЭР, 1980, т. 68, № 10, с. 218—220. [135] R. Kumaresan and D. W. Tufts, <Data-adaptive principal component signal processing», in Proc. 19th IEEE Conf. Decision and Control (Albuquerque, NM), Dec. 10—12, 1980, pp. 949—954. *[136] R. T. Lacoss, “Data adaptive spectral analysis methods,” Geo- physics, vol. 36. pp. 661-675, Aug. 1971. [137] -----, “Autoregressive and maximum likelihood spectral analysis methods.” in Proc. 1976 NATO Conf. Signal Processing. *[138] T. E. Landers and R. T. Lacoss, "Some geophysical applications of autoregressive spectral estimates,” IEEE Trans. Geosci. Elec- Iron., vol. GE-15, pp. 26-32, Jan. 1977. [139] T. Landers, “Maximum entropy cepstral analysis,” in Proc. 1978 RADC Spectral Estimation Workshop, pp. 245-258. [140] Лэнг С. У., Макклеллан Дж.'X. Простое доказательство устойчивости моделей линейного прогнозирования, со- держащих только плюсы. ТИИЭР, 1979, т. 67, № 5, с. 185— 186. [141] S. W. Lang and J. Н. McClellan, «Frequency estimation with maximum entropy spectral estimators», IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-28, pp. 716— 724, Dec. 1980. [ 142] N. Levinson, “The Wiener (root mean square) error criterion in filter design and prediction,” J. Math. Phys., vol. 25, pp. 261- 278, 1947. [143] J. S. Lim, “All pole modeling of degraded speech,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-26, pp. 197-209, June 1978. [144] Макхол. Линейное предсказание. Обзор. ТИИЭР, 1975, т. 63, №4, с. 20—44. [145] J. Machoul, «Spectral linear predication: Properties and ap- plications», IEEE Trans, Acoust., Speech, Signal Proces- sing, vol. ASSP-25, pp. 423—428, Oct. 1977. * [ 146 ] -, “Stable and efficient lattice methods for linear prediction,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-25, pp. 423-428, Oct. 1977. [147] J. Makhoul and R. Viswanathan, “Adaptive lattice methods for linear prediction,” in Conf. Rec. 1978 ICASSP, pp. 83-86. [ 148] J. D. Market, “FFT pruning,” IEEE Trans. Audio Electroacoust., vol. AU-19, pp. 305-311, Dec. 1971. [149] S. L. Marple, Jr., “Conventional Fourier, autoregressive, and special ARMA methods of spectrum analysis,” Engineer’s Dis- sertation, Stanford Univ., Stanford, CA, Dep. Elec. Eng Dec 1976. [150] , “Resolution of conventional Fourier, autoregressive and special ARMA methods of spectral analysis,” in Rec. 1977 Int. Conf. Acoustics, Speech and Signal Processing (Hartford, CT), May 9-11, 1977, pp. 74-77. [151] --, “High resolution autoregressive spectrum analysis using noise power cancellation,” in Rec. 1978 IEEEInt. Conf. Acous- tics, Speech and Signal Processing., pp. 345-348. [152] —“Frequency resolution of high-resolution spectrum analysis techniques,” in Proc. 1978 RADC Spectrum Estimation Work- shop, pp. 19-35. [153] -—, “Spectral line analysis by Pisarenko and Prony methods,” in Rec. 1979 IEEE Int. Conf. Acoustics, Speech, and Signal Pro- cessing, pp. 159-161. [154] --, “A new autoregressive spectrum analysis algorithm,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-28, pp. 441-454, Aug. 1980. [155] D. Q. Mayne et al., “Linear estimation of ARMA systems,” Dep. Computing and Control, Imperial College of Science and Technology, London, England, Research Rep. 77113. [156] R. N. McDonough, ""Representations and analysis of signals. Part XV. Matched exponents for the representation of signals,” Ph.D. dissertation, Dep. Elec. Eng., Johns Hopkins Univ., Balti- more, MD, Apr. 30, 1963. [157] R. N. McDonough and W. H. Huggins, “Best least-squares repre- sentation of signals by exponentials,” IEEE Trans. Automat Contr., vol. AC-13, pp. 408-412. Aug. 1968. *[158] R. N. McDonough, “Maximum entropy spatial processing of array data,” Geophysics, vol. 39, pp. 843-851 Dec 1974. ’ [159] J. M. Melsa and J. D. Tomcik, “Linear predictive coding with additive noise for application to speech digitization,” in Proc. 14th Allerton Conf. Circuits and Systems Theory, pp. 500-508. Sept. 27,1976. [160] R. Mittra and L. W. Pearson, “A variational method for efficient determination of SEM poles,” IEEE Trans. Antennas Propagat., vol. AP-26. pp. 354-358, Mar. 1978. [161] A. Mohammed and R. G. Smith, "Data windowing in spectral analysis,” Defense Research Establishment Atlantic, Darts- mouth, Nova Scotia, Rep. 7512, June 1975. [162] D. R. Moorcroft, “Maximum entropy spectral analysis of radio- auroral signals,” Department of Physics and Centre for Radio Science, Univ. Western Ontario, London, Ontario, Canada. [163] M. Morf, “Fast algorithms for multivariable systems,” Ph.D. dis- sertation, Stanford University, Stanford, CA, Dep. Elec. Eng., .Aug. 1974. [164] M. Morf, T. Kailath, and L. Ljung, “Fast algorithms for recursive ; identification,” in Proc. 1976 IEEE Conf. Decision and Control (Clearwater, FL), Dec. 1-3, 1976, pp. 916-921. •[165] M. Morf, D. T. Lee, J. R. Nickoils, and A. Vieira, “A classifica- tion of algorithms for ARMA models and ladder realizations,” in Rec. 1977 IEEE Int. Conf. Acoustics, Speech, and Signal Pro- cessing. (Hartford, CT), May 9-11, 1977, pp. 13-19. [166] M. Morf, A. Vieira, D. T. Lee, and T. Kailath, “Recursive multi- channel maximum entropy method,” in Proc. 1977 Joint Auto- matic Control Conf. (San Francisco, CA), pp. 113-117, June 22-24, 1977. •[167] M. Morf, B. Dickinson, T. Kailath, and A. Vieira, “Efficient so- lution of covariance equations for linear prediction,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-25, pp. 429- 433, Oct. 1977. [168] M. Morf, A. Vieira, and D. T. Lee, “Ladder forms for identifica- tion and speech processing,” in Proc. 1977 Conf. Decision and Control (New Orleans, LA), pp. 1074-1078, Dec. 1977. *[169] M. Morf, A. Vieira, D. T. Lee, and T. Kailath, “Recursive multi- channel maximum entropy spectral estimation,” IEEE Trans. Geosci. Electron., vol. GE-16,j>p. 85-94, Apr. 1978. v [170] Нейдорфер. Варианта методов гармонического анализа. ТИИЭР, 1973, т. 61, № 11, с. 184—185. •[171] W. I. Newman, "Extension to the maximum entropy method,” IEEE Trans. Inform. Theory, vol. IT-23, pp. 89-93, Jan. 1977. [172] ----, “Extension to the maximum entropy method II,” IEEE Tmns. Inform. Theory, vol. IT-25, pp. 705-708, Nov. 1979. V [173] Ннтцберг P. О возможности определения спектральной плотности. ТИИЭР, 1979, т. 67, №3, с. 120—121. [174] В. Noble and J. W. Daniel, Applied Linear Algebra, 2nd ed. Englewood Cliffs, NJ: Prentice-Hall, 1977. [175] A. H. Nuttall, “Spectral estimation by means of overlapped fast Fourier transform processing of windowed data,” NUSC Tech. Rep. 4169, New London, CT, Oct. 13,1971. [176] , “Spectral analysis of a univariate process with bad data points, via maximum entropy, and linear predictive techniques,” Naval Underwater Systems Center, Tech. Rep. 5303, New Lon- don, CT, Mar. 26,1976. [177] , “FORTRAN program for multivariate linear predictive spectral analysis, employing forward and backward averaging,” 7 ТИИЭР J* 11
50 ТИИЭР, т. 69, № 11, ноябрь 1981 Naval Underwater Systems Center, Tech. Document 5419, New London, CT, May 19, 1976. [178] , “Multivariate linear predictive spectral analysis employing weighted forward and backward averaging: A generalization of Burg’s algorithm,” Naval Underwater Systems Center, Tech. Rep. 5501, New London, CT,'Oct. 13,1976. [179] , “Probability distribution of spectral estimates obtained via g overlapped FFT processing of windowed data,” Naval Under- I water Systems Center, Tech. Rep. 5529, New London, CT, * Dec. 3,1976. [180] —, “Positive definite spectral estimate and stable correlation recunion for multivariate linear predictive spectral analysis,” NUSC Tech. Rep. 5729, Naval Underwater Systems Center, New London, CT, Nov. 14,1977. [181] A. H. Nuttall and G. C. Carter, “A generalized framework for power spectral estimation,” IEEE Trans. Acoust., Speech, Sig- nal Process., vol. ASSP-28, pp. 334-335, June 1980. [182] A. H. Nuttall, “Some windows with very good sidelobe behav* ior,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP- 29, pp. 84-89, Feb. 1981. Ц83] Оппенгейм А. В., Шафер P. В. Цифровая обработка сиг- налов. M.: Связь, 1979. [184] R. К. Otnes and L. Enochson, Digital Time Series Analysis. New York: Wiley, 1972. [185] N. L. Owsley, “Adaptive data orthogonalization,”in Proc. 1978 IEEE ICASSP (Tulsa, OK), pp. 109-112, April 1978. [186] M. Pagano, “Estimation of models of autoregressive signal plus white noise,” Ann. Statistics, vol. 2, pp. 99-108,1974. [187] A. Papoulis, Probability, Random Variables, and Stochastic Pro- cesses. New York: McGraw-Hill, 1965. [188] E. Parzen, “Mathematical considerations in the estimation of spectra,” Technometrics, vol. 3, pp. 167-190, May 1961. [189] —, “Statistical spectral analysis (single channel case) in 1968,” Dep. Statistics, Stanford Univ., Stanford, CA, Tech. Rep. 11, June 10, 1968. [190] ---, “Multiple time series modelling,” Dep. Statistics, Stanford, ' Univ., Stanford, CA, Tech. Rep. 12, July 8, 1968 (Contract Nonr-225-80). [191 ] —, “Multiple time series modelling,” Multivariate Analysis II, P. R. Krishnaiah, Ed. New York: Academic Press, pp. 389- 409,1969. [192] , “Some recent advances in time series modeling,” IEEE Trans. Automat. Contr., vol. AC-19, pp. 723-730, Dec. 1974. [193] J. V. Pendrel and D. E. Smylie, “The relationship between maxi- mum entropy and maximum likelihood spectra,” Geophysics, vol. 44, pp. 1738-1739, Oct. 1979. [194] V. F. Pisarenko, “On the estimation of spectra by means of non- linear functions of the covariance matrix,” Geophysical J. Royal Astronomical Soc., vol. 28, pp. 511-531, 1972. [195] —, “The retrieval of harmonics from a covariance function,” Geophysical J. Royal Astronomical Soc., vol. 33, pp. 347-366, 1973. [196] A. J. Poggio, M. L. Van Blaricum, E. K. Miller, and R. Mittra, “Evaluation for a processing technique for transient data,” IEEE Trans. Antennas Propagat., vol. AP-26, pp. 165-173, Jan. 1978. [197] J. Ponnusamy et al., “Identification of complex autoregressive processes,” in Rec. 1979JICASSP, pp. 384-387. [198] G. Prado and P. Moroney, “Linear predictive spectral analysis for Doppler sonar applications,” Charles Stark Draper Labora- tory, Tech. Rep. R-1109, Sept. 1977. [199] , “The application of linear predictive methods of spectral analysis to frequency estimation,” in Rec. 1978 NAECON. [200] , “The accuracy of center frequency estimators using linear predictive methods,” in Rec. 1978 IEEE Int. Conf. Acoustics, Speech and Signal Processing, pp. 361-364. [201 ] J. F. Prewitt/Amplitude bias in the Fourier transforms of noisy signals,” IEEE Trans. Antennas Propagat., vol. AP-26, pp. 730- 731, Sept. 1978. [202] G.R.B. Prony, “Essai experimental et analytique, etc.,” Paris, J. de L’Ecole Poly technique, vol. 1, cahier 2, pp. 24-76,1795. [203] L. C. Pusey, “High resolution spectral estimates,” Lincoln Lab- oratory, M.I.T., Tech. Note 1975-7, Jan. 21,1975. [204] С. M. Rader, “An improved algorithm for high speed autocorre- lation with applications to spectral estimation,” IEEE Trans. Audio Electroacoust., vol. AU-18, pp. 439-442, Dec. 1970. [205] H. R. Radoski, P. F.Fougere, and E. J. Zawaiick, “A compari- son of power spectral estimates and applications of the maxi- mum entropy method,” J. Geophysical Res., vol. 80, pp. 619- 625, Feb. 1,1975. '[206 ] H. R. Radoski, E. J. Zawaiick, and P. F. Fougere, “The superior- ity of maximum entropy power spectrum techniques applied to geomagnetic-micropulsations,” Phys. Earth Planetary Interiors, vol. 12, pp. 208-216, Aug. 1976. [207] P. I. Richards, “Computing reliable power spectra,” IEEE Spec- trum, vol. 4, pp. 83-90. Jan. 1967. [208] D. C. Rife and G. A. Vincent, “Use of the discrete Fourier trans- form in the measurement of frequencies and levels of tones,” Bell Syst. Tech. J., vol. 49, pp. 197-228, Feb. 1970. [209] E. A. Robinson, “Predictive decomposition of time series with application to seismic exploration,” Geophysics, Vol. 32, pp. 418-484, June 1967. [210] —, Multichannel Time Series Analysis with Digital Computer Programs. San Francisco, CA: Holden-Day, 1967. [211] Сейдж Э. П., Меле Дж. Теория оценивания и ее применение в связи и управлении. М.: Связь, 1976, гл. 6. (212] Н. Sakai, “Statistical properties of AR spectral analysis,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-27, pp. 402-409, Aug. 1979. [213] M. R. Sambur, “A preprocessing filter for enhancing LPC analy- sis/synthesis of noisy speech,” in Rec. 1979 ICASSP, pp. 971- 974. [214] E. H. Satorius and J. T. Alexander, “High resolution spectral analysis of sinusoids in correlated noise,” in Rec. 1978 ICASSP (Tulsa, OK), pp. 349-351, Apr. 10-12, 1978. [215] E. H. Satorius and J. R. Zeidler, “Maximum entropy spectral analysis of multiple sinusoids in noise,” Geophysics, vol. 43, pp. 1111-1118, Oct. 1978. [216] J. H. Sawyers, “Applying the maximum entropy method to adaptive digital filtering,” in Rec. 12th Asilomar Conf. Circuits, Systemsand Computers, pp. 198-202. [217] D. H. Scahuhert, “Application ofProny’s method to time-domain reflectometer data and equivalent circuit synthesis,” IEEE Trans. Antennas Propagat., vol. AP-27, pp. 180-184, Mar. 1979. [218] A. Schuster, “On the investigation of hidden periodicities with application to a supposed 26 day period of meteorological phe- nomena,” Terrestrial Magnetism, vol. 3, pp. 13-41, Mar. 1898. [219] , “The periodogram of magnetic declination as obtained from the records of the Greenwich Observatory during the years 1871-1895,” Trans. Cambridge Philosophical Soc., vol. 18, pp. 1 107-135, 1899- , [220] Шарма, Махаланабис. Гармонический анализ с помощью метода калмановской фильтрации. ТИИЭР, 1973, т. 61, № 3, с. 162—163. [221 ] D. Е. Smylie, G.K.C. Clarke, and Т. J. Ulrych, “Analysis of irreg- ularities in the earth’s rotation,” in Methods in Computational Physics, vol. 13. B. A. Bolt, Ed. New York: Academic Press, 1973,pp. 391-430. [222] M. D. Srinath and M. M. Viswanathan/“Sequential algorithm for identification of parameters of an autoregressive process,” IEEE Trans. Automat. Contr., vol. AC-20, pp. 542-546, Aug. 1975. [223] G. G. Stokes, “Note on the search for periodic inequalities,” Proc. Royal Soc. London, vol. 29, pp. 122-123, May 29,1879. (2241 Стрейдер II H. P. Влияние субгармоник на коэффициента дискретного преобразования Фурье. ТИИЭР, 1980, т. 68, №2, с. 104—105. [225] О. N. Strand, “Multichannel maximum entropy spectral analy- sis,” IEEE lions. Automat. Contr., vol. AC-22, pp. 634-640, Aug. 1977. [226] T. M. Sullivan, O. L. Frost, and J. R. Treichler, “High resolution signal estimation,” ARGO Systems, Internal Rep. June 16,1978. [227] D. N. Swingler, “A comparison between Burg’s maximum en- tropy method and a nonrecursive technique for the spectral analysis of deterministic signals,” J. Geophysical Res., vol. 84, pp.679-685, Feb. 10,1979. V [228] Суинглер Д. H. Модифицированный алгоритм Берга для спектрального анализа максимальной энтропии. ТИИЭР, 1979, т. 67, №9, с. 223—225. [229] D. N. Swingler, «Frequency errors in MEM processing», IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-28, pp. 257—259, Apr. 1908. [230] ---, “Comments on maximum entropy spectral estimation us- ing the analytical signal,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-28, pp. 259-260, Apr. 1980. [231] D. J. Thomson, M. F. Robbins, C. G. Maclennan, and L. J. Lan- zerotti, “Spectral and windowing techniques in power spectral analyses of geomagnetic data,” Phys. Earth Planetary Interiors, vol. 12, pp. 217-231, Aug. 1976. [232] D. J. Thomson, “Spectrum estimation techniques for charac- terization and development of WT4 waveguide-I,” Bell Syst. Tech. J., vol. 56, pp. 1769.-1815, Nov. 1977. [233] T. Thorvaldsen, A. T. Waterman Jr., and R. W. Lee, “Maximum entropy angular response patterns of microwave transhorizon signals,” IEEE Trans. Antennas Propagat., vol. AP-28, pp. 722- 724, Sept. 1980. r [234] K. Toman, “The spectral shift of truncated sinusoids, J. Geo- physical Res., ml. 70, pp. 1749-1750, Apr. 1, 1965. *[235] H. Tong, “Autoregressive model fitting with noisy data by Akaike’s information criterion,” IEEE Trans. Inform. Theory, vol. IT-21, pp. 476-480, July 1975.
СПЕКТРАЛЬНЫЙ АНАЛИЗ 51 [236] —, “Fitting a smooth moving average to noisy data,” IEEE Trans. Inform. Theory, vol. IT-22, pp. 493-496, July 1976. •[237] ----, “More on autoregressive model fitting with noisy data by Akaike’s information criterion,” IEEE Trans. Inform. Theory, vol. IT-23, pp. 409-410, May 1977. [238] H. Tong, “Final prediction error and final interpolation error: A paradox?” 1ЕЕЁ Trans. Inform. Theory, vol. IT-25, pp. 758- 759, Nov. 1979. (239) J. P. Toomey, “High-resolution frequency measurement by lin- ear prediction,” IEEE Trans. Aerospace Electron. Sy st., vol. AES-16, pp. 517-525, July 1980. [240] J. Treichler, “7-LMS and its use in a noise-compensating adap- tive spectral analysis technique,” in Rec. 1979 ICASSP, pp. 933-936. [241] W. Trench, “An algorithm for the inversion of finite Hankel matrices,” /. Soc. Indust. Appl. Math., vol. 13, no. 4, pp. 1102- 1107, Dec. 1965. [242] S. A. Tretter and K. Steiglitz, “Power spectrum identification in terms of rational models,” IEEE Trans. Automat. Contr., vol. AC-12, pp. 185-188, Apr. 1967. 1243] Тафтс Д. УГриффитс Л. Дж., Уидроу Б., Гловер Дж., Маккул Дж., Трейчлер Дж. Адаптивная селекция и спек- ?>альный анализ. ТИИЭР, 1977, т. 65, № 1, с. 199—203. афте Д. У., Кумаресан Р. Улучшенные методы спектраль- ного разрешения. ТИИЭР, 1980, т. 68, № 3, с. 137—138. (245] D. W. Tufts and R. Kumaresan, «Improved spectral re- solution II», in Proc. 1980 ICASSP (Denver, CO), vol. 2, pp. 592—597, Apr. 9—11, 1980. [246] T. J. Ulrych, “Maximum entropy power spectrum of long period geomagnetic reversals,” Nature, vol. 235, pp. 218-219, Jan. 28, 1972. (247] ----, “Maximum entropy power spectrum of truncated sinu- soids,” J. Geophysical Research, vol.'77, pp. 1396-1400, March 10, 1972. [248] T. J. Ulrych, D. E. Smylie, O. G. Jensen, and G.K.C. Clarke, “Predictive filtering and smoothing of short records by using maximum entropy,” J. Geophysical Res., vol. 78, pp. 4959- 4964, Aug. 10, 1973. [249] T. J. Ulrych and O. G. Jensen, “Cross-spectral analysis using maximum entropy,” Geophysics, vol. 39, pp. 353-354, June ... 1224. ’[250] T. J. Ulrych and T. N. Bishop, “Maximum entropy spectral anal- ysis and autoregressive decomposition,” Rev. Geophysics Space Phys., vol. 13, pp. 183-200, Feb. 1975. [251] T. J. Ulrych and R. W. Clayton, “Time series modelling and maximum entropy,” Phys. Earth Planetary Interiors, vol. 12, pp. 188-200, Aug. 1976. [252] M. L. Van Blaricum and R. Mittra, “Problems and solutions as- sociated with Prony’s method for processing transient data,” IEEE Trans. Antennas Propagat., vol. AP-26, pp. 174-182, Jan. 1978. _ _ , [253] M. L. Van Blaricum, “A review of Prony’s method techniques for parameter estimation,” in Proc. Rome Air Development Center Spectrum Estimation Workshop, Griffis Air Force Base, May 24-26, 1978. pp. 125-139. [254] M. L. Van Blaricum and R. Mittra, Correction to “Problems and : solutions associated with Prony’s method for processing tran- sient data,” IEEE Trans. Antennas Propagat., vol. AP-28, p. 949, Nov. 1980. • [255 ] A. VanDenBos, "Alternative interpretation of maximum entropy spectral analysis,” IEEE Trans. Inform. Theory, vol. IT-17, pp. 493-494, July 1971. /[256] Ван-Трис Г. Теория обнаружения, оценок и модуляций. Т. 1. М.: Сов. радио, 1972. [257] R. Viswanathan and J. Makhoul, “Quantization properties of transmission parameters in linear predictive systems,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-23,- pp. 309-321, June 1975. [258] G. Walker, “On periodicity in series of related terms,” Proc. Roy. Soc. London, Series A, vol. 131, pp. 518-532, 1931. [259] C. Webb, "Practical use of the fast Fourier transform (FFT) al- gorithm in time series analysis,” ARL, Univ. Texas, Austin, TX, ARL-TR-70-22, June 22, 1970. [260] Уэбстер P. Дж. Регулирование просачивания спектральных составляющих ДПФ. ТИИЭР, 1980, т. 68, № 10, с. 203—205. [261 ] L. Weiss and R. N. McDonough, “Prony’s method, Z-transforms, and Pade approximation,” SIAM Rev., vol. 5, pp. 145-149, Apr. 1963. *[262] P. D. Welch, “The use of fast Fourier transform for the estima- tion of power spectra: A method based on time averaging over short, modified periodograms,” IEEE Trans. Audio Electro- acoust., vol. AU-15, pp. 70-73, June 1967. [263] P. D. Welch, “On the variance of time and frequency averages over modified periodograms,” in Rec. 1977 IEEE Int. Conf. Acoustics, Speech, and Signal Processing (Hartford, CT), May 9-11, 1977, pp. 58-62. *[264] S. J. Wemecke and L. R. D’Addario, “Maximum entropy image reconstruction,” IEEE Trans. Comput., vol. C-26, pp. 351-364, Apr. 1977. [265] S. J. Wemecke, “Two-dimensional maximum entropy recon- struction of radio brightness,” Radio Sci., vol. 12, pp. 831-844, Sept.-Oct. 1977. [266] E. C. Whitman, “The spectral analysis of discrete time series in terms of linear regressive models,” Naval Ordnance Labs Rep. NOLTR-70-109, White Oak, MD, June 23.1974. [267] P. Whittle, “On the fitting of multivariate autoregressions, and the approximate canonical factorization of a spectral density matrix,” Biometrika, vol. 50. on. 129-134,1963. [268] Уидроу, Гловер мл., Маккул, Кауниц и др. Адаптивные компенсаторы помех. Принципы построения и примене- ния. ТИИЭР, 1975, т. 63, № 12, с. 69—98. [269 ] N. Wiener, “Generalized harmonic analysis,” Acta Mathematics, vol. 55, pp. 117-258, 1930. [270] R. A. Wiggins and E. A. Robinson, “Recursive solution to the multi-channel filtering problem,” J. Geophysical Res., vol. 70 . pp. 1885-1891, Apr. 1965. [271] Уилски А. С. Взаимосвязь между теорией цифровой обработ- ’ ки сигналов и теорией управления и оценивания. ТИИЭР, 1978, т. 66, №9, с. 5—33. [272] A. S. Willsky, Digital Signal Processing and Control and Estimation theory: Points of Tangency, Areas of Intersection, and Parallel Directions. Cambridge, MA: M.I.T. Press, 1980. [273] T. H. Wonnacott, “Spectral analysis combining a Bartlett win- dow with an associated inner window,” Technometrics, vol. 3, j [274] Юань С. К- Сравнение пяти методов вычисления спектра мощности случайного процесса с использованием сегмен- тации данных. ТИИЭР, 1977, т. 65, № 6, с. 200—202. [275] , “Comments on modem methods for spectral estimation,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-27, pp. 298-299, June 1979. [276] G. U. Yule, “On a method of investigating periodicities in dis- turbed series, with special reference to Wolfer’s sunspot num- bers,” Philosophical Trans. Royal Soc. London, Series A, vol. 226, pp. 267-298, July 1927. [277] S. Zohar, “The solution of a Toeplitz set of linear equations,”/. Assoc. Comput. Mach., vol. 21, pp. 272-276, Apr. 1974. [278] , “FORTRAN subroutines for the solution of Toeplitz sets of linear equations,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-27, pp. 656-658, Dec. 1979. Стивен M. Kefi родился в 1951 г. Степень бакалавра в области электротехники получил в Технологическом институте Стивенса, Хобокен, шт. Нью-Джерси, в 1972 г., степени магистра и доктора философии в области электротехники — в Колумбийском уни- верситете, шт. Нью-Йорк, и Технологическом институте шт. Джор- джия, Атланта, соответственно в 1972 и 1973 гг. С 1972 по 1975 г. работал в фирме Bell Telephone Laboratories, Хоумдел, шт. Нью- Джерси, где принимал участие в планировании передач для рече- вых систем связи и в моделировании и непосредственном испыта- нии алгоритмов обработки .речевых сигналов. С 1975 по 1977 г. посещал Технологический институт шт. Джорджия, где изучал теорию связи и методы цифровой обработки сигналов. С 1977 по 1980 г. работал в Отделении сигналов ,Ц|Пя m™?"* тДЬ фирмы ~ Raytheon Company, Портсмут, шт. Род-Айленд, где принимал участие в исследовании авторегрессивных методов спектральной оценки и проектировании гидролокационных систем. В настоя- щее время занимает должность доцента кафедры электротехники Университета шт. Род-Айленд, Кингстон. Сфера научных инте- ресов включает спектральный анализ, теорию обнаружения и оценивания и адаптивную фильтрацию. Стэнли Лоренс Марпл мл. родился в 1947 г. Степень бакалав- ра и магистра в области электротехники получил в Университете Райса, Хьюстон, шт. Техас, в 1970 г., диплом инженера-электри- ка — в Стэнфордском университете, Станфорд, шт. Калифорния, в 1976 г. С 1971 по 1977 г. работал в фирме ArgoSystems, Inc., Саннивейл, шт. Калифорния, где принимал участие в разработке аппаратурных средств для цифровой обработки радиолокацион- ных сигналов и анализа связных сигналов. С 1978 по. 1979 г. ра- ботал в должности старшего научного сотрудника в фирмах Signal Science, Inc., и Advent Systems, Inc., Маунтин-Вью, шт.Ка- лифорния. В настоящее время сотрудник фирмы The Analytic Sciences Corporation (TASC), Мак-Лин, шт. Виргиния. Сфера на- учных интересов включает приложения методов цифровой обра- ботки сигналов, связь, радиолокацию и обработку изображений.