6.18. Энергетический спектр случайных сигналов

При спектральных измерениях случайных последовательностей встречается ряд трудностей. Во-первых, обычно интерес представляет энергетический спектр (определение энергетического спектра еще не было дано). Во-вторых, измерения должны быть статистически устойчивыми в том смысле, что, какая бы оценка ни рассматривалась, при увеличении интервала измерения она должна сходиться к вполне конкретной величине.

Простое и удобное условие сходимости сводится к тому, чтобы дисперсия оценки стремилась к нулю при увеличении интервала измерения. Известно, что при работе с энергетическим спектром весьма удобно пользоваться корреляционной функцией сигнала. В связи с этим в данном разделе будут кратко рассмотрены понятия корреляционной функции и энергетического спектра и связывающие их соотношения. Затем будут проанализированы с точки зрения эффективности вычислений два хорошо известных метода расчета энергетического спектра случайных сигналов.

До сих пор использовалось довольно общее определение спектра. Он определялся как значение z-преобразования в одной или ряде точек в z-плоскости. В связи с этим возникает следующий вопрос: можно ли энергетический спектр сигнала определить в таком же общем виде? По-видимому, разумным подходом является следующий: если задана спектральная плотность мощности на некоторой кривой в z-плоскости (не обязательно на единичной окружности), то, интегрируя эту плотность вдоль кривой, можно получить общую мощность сигнала. Используя теорему Парсеваля, можно убедиться в том, что такой подход возможен. Так, для ограниченной случайной последовательности , состоящей из N отсчетов, теорема Парсеваля записывается следующим образом:

Соотношение (6.75) является упрощенным вариантом теоремы о комплексной свертке. Если последовательность имеет конечную длину, то интегрирование может проводиться по любому контуру, охватывающему начало координат, поэтому спектральную плотность мощности можно определить как величину, пропорциональную произведению на этом произвольном контуре, т. е.

Из этого соотношения следует, что для определения спектральной плотности мощности в точке плоскости z необходимо вычислять значения z-преобразования заданного сигнала при и при Например, если требуется измерить спектральную плотность на окружности радиуса , необходимо вычислить значения z-преобразования последовательности и на этой окружности, и на окружности радиуса . При оба эти значения -преобразования становятся комплексно сопряженными, так что по существу в качестве спектральной плотности мощности измеряется величина .

Чтобы определить коэффициент пропорциональности в соотношении (6.76), воспользуемся статистическим методом. Запишем

Найдем среднее значение от этого выражения

Рассмотрим простейший пример анализа «белого» шума, для которого среднее значение произведения равно нулю всегда, исключением случая когда оно равно, предположим, К. Таким образом,

Поскольку желательно, чтобы среднее значение сходилось к постоянной величине при увеличении N, спектральную плотность следует определить следующим образом:

Например, при N = 5, согласно формуле (6.80), имеем

Объединяя члены с одинаковыми степенями z, получаем

Из выражений (6.81) и (6.82) легко получить формулу для произвольного N:

где

. Выражение (6.84) определяет автокорреляционную функцию последовательности конечной длины .

Из формулы (6.83) следует, что спектральная плотность мощности равна z-преобразованию автокорреляционной функции сигнала . Хотя выше предполагалось, что случайный сигнал имеет конечную длительность, как правило, последовательность не ограничена. В этих случаях основная формула для вычисления автокорреляционной функции имеет вид

На практике спектральная плотность мощности чаще всего вычисляется на единичной окружности. В этом случае формула (6.83) принимает вид

т.е. равно преобразованию Фурье от . Соответственно может быть получена обратным преобразованием Фурье от , т. е.

Соотношения (6.86) и (6.87) можно обобщить на случай, когда желательно измерить взаимный энергетический спектр или функцию взаимной корреляции двух случайных колебаний и :

(6.88)

Фиг. 6.30. Два метода выполнения спектрального анализа.

Функцию взаимной корреляции измеряют в случаях, когда последовательности связаны операциями фильтрации и задержки. Например, в радиолокационной системе может представлять собой отфильтрованный, зашумленный и задержанный зондирующий сигнал . В таких случаях положение максимума взаимной корреляционной функции дает величину задержки входного сигнала.

Итак, выше было дано определение спектральной плотности мощности, установлена взаимосвязь между спектральной плотностью мощности и автокорреляционной функцией и показано, что средние значения этих оценок сходятся при больших к некоторому определенному постоянному числу, по крайней мере в случае белого шума. Однако хорошо известно, что дисперсия оценки будет сходиться к нулю лишь при выполнении дальнейшего усреднения. Существуют два хорошо известных метода измерения спектральной плотности шума с использованием БПФ, для которых удовлетворяются условия сходимости дисперсии к нулю. Первый метод основан на вычислении корреляционной функции с помощью алгоритма БПФ, а второй состоит в усреднении последовательности непосредственных измерений спектральной плотности. На фиг. 6.30 иллюстрируются все этапы выполнения спектрального анализа этими методами.

В методе 1 БПФ используется непосредственно для вычисления оценок взаимной корреляционной функции при L задержках, где 2L — размер БПФ.

При вычислении по значениям спектральной плотности мощности для конечного числа частот необходимо применять сглаживающее окно , чтобы уменьшить нежелательные эффекты, связанные с конечной длиной выборки, поскольку вместо бесконечной корреляционной последовательности используются только L ее значений. С помощью L-точечного БПФ можно вычислить значения спектральной плотности мощности на L равноотстоящих частотах на единичной окружности.

В методе 2 алгоритм БПФ используется для непосредственного вычисления спектральной плотности мощности в отличие от первого метода, в котором сначала вычисляется корреляционная функция. Поэтому каждая из последовательностей до выполнения преобразования сглаживается с помощью окна . Для последующего вычисления на основе спектральной плотности мощности корреляционной функции можно использовать алгоритм обратного БПФ. Так как спектральная плотность мощности была рассчитана только на L дискретных частотах (вместо бесконечного числа частот, что требуется при теоретическом рассмотрении), после обратного БПФ будет получена корреляционная функция, искаженная эффектом наложения.

Ниже будут рассмотрены вопросы практического применения методов 1 и 2, однако сначала покажем, как можно использовать ДПФ для выполнения корреляционного анализа. Рассмотрим периодические последовательности (с периодом в L отсчетов) ДПФ, которых соответственно равны

(6.91)

и

Выше было показано, что умножение на соответствует свертке периодических последовательностей , т. е.

Покажем теперь, что круговую свертку можно элементарно найти следующим образом:

Для проверки соотношения (6.94) покажем, что обратное ДПФ произведения равно правой части (6.94).

Используя формулы (6.91) и (6.92), получим

что и требовалось доказать. [Заметим, что во всех приведенных выкладках последовательности считаются действительными.]

1. Метод 1 - корреляционный.

На фиг. 6.31 иллюстрируется один из способов вычисления первых пяти отсчетов автокорреляционной функции заданной -точечной последовательности (N = 20). Сама последовательность показана на фиг. 6.31 вверху. Первый блок чисел размером (9x5) соответствует круговой корреляции двух девятиточечных последовательностей, а именно . Пять из полученных отсчетов являются правильными; они показаны справа. Второй блок начинается с аналогичная процедура вычисления дает отсчеты от до . Чтобы найти результирующие значения первых пяти отсчетов автокорреляционной функции, необходимо просуммировать по i следующим образом:

Следующий шаг состоит в том, чтобы показать, как выполнять вычисления с использованием алгоритма БПФ.

Фиг. 6.31. Иллюстрация вычисления круговой корреляции.

Каждую из четырех круговых корреляционных функций девятиточечных последовательностей можно вычислить с помощью ДПФ по следующей схеме:

1. Вычислить

и

2. Вычислить

3. Вычислить

4. Вычислить

5. Первые пять отсчетов являются искомым результатом и также удовлетворяют соотношению

Описанная процедура эффективнее, чем предварительный расчет каждой из частичных корреляционных, функций с помощью обратного ДПФ от и последующее суммирование этих частичных корреляционных функций. В первом случае для каждого из показанных на фиг. 6.31 блоков чисел размером рассчитываются два ДПФ, что дает , а одно «общее» ДПФ искомой корреляционной функции получается путем накопления частичных ДПФ. «Общая» корреляционная функция вычисляется посредством единственного обратного ДПФ от полученного путем накопления «полного» ДПФ. Следовательно, общее число выполняемых преобразований равно , где К — количество обрабатываемых блоков данных. Во втором случае, когда накапливаются частичные корреляционные функции, преобразование приходится выполнять раз.

Из описанной выше методики вычислений видно, что размер ДПФ должен приблизительно вдвое превышать требуемое число отсчетов корреляционной функции, а требуемое количество преобразований равно (здесь , причем N — общее число отсчетов сигнала, а М — требуемое число отсчетов корреляционной функции). Затем, дополнив нулями и выполнив единственное ДПФ, можно получить оценку энергетического спектра с произвольным частотным разрешением. Если является хорошей оценкой корреляционной функции с малой дисперсией, то и преобразование от также будет достаточно хорошей оценкой спектра.

Для дальнейшего повышения эффективности вычислений можно использовать прием, предложенный Рэйдером, который отметил, что вычисляемые на первом зтапе значения можно выразить через . Для этого следует прежде всего несколько изменить характер показанного на фиг. 6.31 секционирования последовательности, предусмотрев, чтобы секции содержали четное число отсчетов, желательно кратное 2, и перекрытие секций в точности было равно 2:1. Эти условия вполне выполнимы, потому что для большинства алгоритмов БПФ требуется или желательно иметь секции именно такой длины, а обеспечение перекрытия, в точности равного 2:1, означает, что в данном случае размер преобразования возрастает на один отсчет, что приводит конечно, к незначительному ухудшению эффективности вычислений.

(см. скан)

Фиг. 6.32. Секционирование при вычислении оценок автокорреляционной функции.

На фиг. 6.32 изображена 20-точечная последовательность с отсчетами от до , разделенная на пять секций, и показано, для каких секций рассчитывается круговая автокорреляция. Справа показаны эквивалентные БПФ-операции. Идея излагаемого метода, как уже было отмечено, состоит в том, что любой коэффициент ДПФ можно представить в функции от . Если, например,

то

так как круговой сдвиг на при ДПФ эквивалентен умножению самого ДПФ на Следовательно,

В более общем виде

Теперь, используя пару преобразований (6.93), описанную выше процедуру спектральных измерений можно изменить следующим образом:

1. Вычислить все от до , где — ДПФ последовательностей, показанных на фиг. 6.32.

2. Вычислить

причем

3. Вычислить

4. Вычислить

После определения можно вычислить спектральные оценки, рассчитав ДПФ от Важно отметить, что количество преобразований, которые в рассматриваемом методе составляют основу всех вычислений, удалось уменьшить до К, где К — количество секций.

Для того чтобы дисперсия оценки сходилась к нулю, нужно, чтобы аргумент функции корреляции был значительно меньше интервала измерения. В результате оценка спектра вычисляется через преобразование усеченной последовательности , а, как известно, усечение во временной области может вызвать чрезмерные пульсации в частотной области; в связи с этим приходится использовать окна.

Проиллюстрируем сказанное на примере. Пусть интервал измерения содержит 1024 отсчета и известно, что отсчеты малы при , причем требуется найти 64 спектральные оценки на 64 частотах. Процедура вычислений состоит в следующем:

1. Разделить секцию из 1024 отсчетов на подсекции длиной по 32 отсчета.

2. Рассчитать 32-точечные БПФ этих подсекций с перекрытием 2:1, т. е. сдвигаясь каждый раз на 16 отсчетов.

3. Используя описанный выше метод, найти для в пределах от 0 до 16.

4. Чтобы получить оценку спектра, умножить на соответствующее окно . Затем дополнить нулями на интервале от до и вычислить БПФ сформированного массива.

В общем случае, когда требуется вычислить L отсчетов корреляционной функции по N заданным отсчетам сигнала (N L), рассматриваемый корреляционный метод спектральных измерений можно сформулировать следующим образом. Искомая функция взаимной корреляции имеет вид

Простая методика вычисления по формуле (6.97) с использованием Р-точечного БПФ (Р = 2L) представлена на фиг. 6.33. Из последовательности образуются подпоследовательности согласно правилу

а из последовательности формируются подпоследовательности по несколько иному правилу:

Для каждой пары Р-точечных подпоследовательностей вычисляется последовательность представляющая собой периодическую корреляцию . Она может быть эаписана в виде

Фиг. 6.33. Секционирование при вычислении функции взаимной корреляции.

Следует отметить, что только отсчетов вычисленной корреляционной последовательности удовлетворяют формуле (6.100); остальные отсчетов соответствуют ненужной периодической корреляции. Суммируя отдельные корреляционные последовательности, находим

(6.101)

(6.102)

что идентично искомому результату (6.97) для . Таким образом, корреляционные функции можно эффективно вычислять с использованием БПФ при любом числе задержек L.

Как указывалось выше, при вычислении спектральной плотности мощности по корреляционной последовательности с использованием окна могут возникнуть трудности. Если ДПФ окна не является положительным, то существует вероятность того, что вычисленная спектральная плотность мощности будет отрицательной, что весьма нежелательно. Причина этого нежелательного результата состоит в том, что вычисленная спектральная плотность мощности равна круговой свертке преобразования окна и преобразования найденной корреляционной последовательности, поэтому, если преобразование окна не является положительным на всех частотах, существует вероятность того, что из-за флуктуационных ошибок оценок корреляционной функции результирующая свертка на некоторых частотах может дать отрицательные значения спектра мощности. Существуют, однако, окна, преобразования которых положительны на всех частотах, например треугольное окно. Именно такие окна и следует применять в тех случаях, когда другие окна приводят к нежелательному результату.

Еще одна трудность применения корреляционного метода встречается при оценке взаимной спектральной плотности по функции взаимной корреляции с использованием окна. Как видно из формулы (6.100), отсчеты функции корреляции здесь измеряются только при положительных значениях . При использовании этого метода для оценки автокорреляционной функции ее отсчеты при отрицательных получаются автоматически благодаря симметрии автокорреляционной функции. Однако для функции взаимной корреляции эта симметрия не существует. Таким образом, неясно, как использовать окна для сглаживания взаимно корреляционной последовательности при оценке взаимной спектральной плотности. Существует несколько возможных способов сглаживания; они иллюстрируются на фиг. 6.34. Трудность заключается в том, что весьма непросто придать физический смысл значению взаимно корреляционной функции при поскольку последняя зависит от задержки между и . Так, если , то будет иметь пик при . И в этом случае можно было бы рассматривать точку как некоторое начало координат На фиг. 6.34 иллюстрируются три возможных способа взвешивания с использованием окна до выполнения ДПФ, позволяющих получить оценки , имеющие физический смысл.

Фиг. 6.34. Варианты окон для вычисления взаимного спектра по взаимной корреляции.

На фиг. 6.34, а окно (в данном примере треугольное) симметрично относительно точки причем предполагается, что размер БПФ при вычислении равен На фиг. 6.34, б размер БПФ также равен но окно симметрично относительно точки , которой придается физический смысл начала координат. Вариант взвешивания, иллюстрируемый на фиг. 6.34, в, по-видимому, наиболее разумный. Здесь оцениваются одновременно и причем вычисляется точно так же, как и но меняются местами. В зтом случае используется Р-точечное БЙФ и снова предполагается, что окно симметрично относительно точки . Конечно, время вычисления в третьем случае в два раза больше, чем в двух Предыдущих, но часто эти затраты времени не имеют смысла. Более того, поскольку используется Р-точечное БПФ, в третьем случае достигается более высокое разрешение оценки взаимной спектральной плотности по частоте.

На фиг. 6.35-6.37 иллюстрируется применение корреляционного метода для измерения автокорреляционной функции и спектральной плотности мощности окрашенного шума, а именно белого шума, пропущенного через КИХ-фильтр нижних частот с равновеликими пульсациями частотной характеристики и импульсной характеристикой, состоящей из 39 отсчетов.

Автокорреляционная функция шума на выходе такого фильтра теоретически отлична от нуля только на интервале . В действительности измеренные значения функции автокорреляции вне зтого интервала в точности не равны нулю, хотя и близки к нему. Аналогично идеальная функция спектральной плотности мощности шума на выходе фильтра нижних частот равна

где - частотная характеристика фильтра нижних частот. На фиг. 6.35 изображено несколько графиков энергетического спектра (в логарифмическом масштабе) при различных значениях числа отсчетов N, по которым вычисляется корреляционная функция. Графики соответствуют значениям N = 128, 256, 1024, 4096 и 16384 соответственно, причем число задержек автокорреляционной функции равно 59. При получении оценки спектральной плотности мощности для взвешивания автокорреляционной функции использовалось окно Хэмминга. Чтобы получить достаточно хорошее частотное разрешение, автокорреляционная функция дополнялась нулями и выполнялось 1024-точечное БПФ. Графики для разных N показывают, что функция быстро сходится к .

Фиг. 6.35. Спектральная плотность мощности шума на выходе фильтра нижних частот, полученная корреляционным методом (в логарифмическом масштабе).

Для сравнения на фиг. 6.35, е показана идеальная характеристика На фиг. 6.35, ж изображен график спектральной плотности мощности, полученный с помощью БПФ по 512 точкам автокорреляционной функции, взвешенной треугольным окном. Поскольку автокорреляционная функция при представляет собой чисто случайную шумовую последовательность, полученная функция спектральной плотности будет совпадать с идеальной характеристикой фильтра с наложенным на нее шумом.

Фиг. 6.36. Спектральная плотность мощности шума на выходе полосового фильтра, полученная корреляционным методом.

На фиг. 6.36 и 6.37 приведены аналогичные результаты прохождения шума через полосовой КИХ-фильтр с 15-точечной импульсной характеристикой и через КИХ-дифференциатор с 27-точечной импульсной характеристикой соответственно. Очевидно, что корреляционный метод обеспечивает достаточно хорошую сходимость оценок, так что в рассмотренных простых примерах получаемые спектры мощости оказываются достаточно хорошей аппроксимацией теоретических спектров.

2. Метод 2 - метод модифицированных периодограмм

Интуиция подсказывает, что достаточно хорошую оценку спектра можно получить, используя гребенку фильтров, как показано на фиг. 6.38.

В соответствии с изложенным в разд. 6.13 гребенку фильтров можно реализовать с помощью скользящего БПФ. Если, кроме того, прореживать выходную последовательность (фиг. 6.38), то при условии, что фильтры гребенки имеют достаточно узкие полосы, будет меняться сравнительно медленно. Однако использование прореживания позволяет вместо скользящего БПФ воспользоваться скачущим БПФ; такая процедура спектральных измерений носит название метода модифицированных периодограмм.

Для иллюстрации применения этого метода рассмотрим сначала случай , т. е. нужно получить спектр мощности стационарного случайного процесса, представленного типичной последовательностью . На фиг 6.39 показана такая последовательность и иллюстрируется способ ее разбиения на подпоследовательности длиной по L отсчетов. Подпоследовательности сдвинуты относительно друг друга на D отсчетов, причем в приведенном на фиг. 6.39 примере , так что связаны с соотношением

Фиг. 6.37. Спектральная плотность мощности шума на выходе дифференциатора, полученная корреляционным методом.

Фиг. 6.38. Один канал спектроанализатора для измерения шумов.

(6.103)

где К — число подпоследовательностей, используемых для вычисления спектральной плотности мощности. Для каждой из взвешенных подпоследовательностей рассчитываются коэффициенты БПФ, , по формуле

(6.104)

где — соответстьующее окно. Величина , называемая периодограммой, вычисляется по формуле

где

а

Оценка спектральной плотности мощности находится по формуле

(6.108)

Фиг. 6.39. Секционирование входной последовательности при вычислении спектральной плотности.

Таким образом, оценка спектральной плотности мощности равна среднему значению периодограмм отдельных подпоследовательностей. Уэлч показал, что математическое ожидание этой оценки описывается формулой

(6.109)

где

(6.110)

а — истинная спектральная плотность мощности анализируемого процесса с отчетами . Из формулы (6.109) следует, что математическое ожидание искомой оценки равно свертке истинной спектральной плотности мощности с квадратом модуля фурье-преобразования от последовательности окна. Кроме того, Уэлч показал, что в случае, когда последовательность формируется из реализации гауссовского случайного процесса, а функция является достаточно гладкой в диапазоне частот, где значения ДПФ окна достаточно велики, дисперсия оценки описывается формулой

где

(6.113)

В качестве иллюстрации применения формулы (6.112) рассмотрим случай , для которого при . Для этого случая

(6.114)

где . Теперь рассмотрим случай , причем считаем, что используется треугольное окно, изображенное на фиг. 6.40. Для этого случая [см. формулу (6.113)] и при . Следовательно,

(6.115)

Фиг. 6.40. Треугольное окно.

Но , поэтому соотношение (6.115) можно переписать в виде

(6.116)

Итак, используя всего лишь сдвиг подпоследовательности на половину ширины окна, а не на полное окно, можно уменьшить дисперсию оценки почти в два раза (но за счет удвоения времени вычисления).

При вычислении взаимного спектра обе последовательности разбиваются на подпоследовательности , как показано на фиг. 6.39. Каждая из подпоследовательностей перед вычислением преобразований взвешивается окном. Периодограмма в этом случае имеет вид

(6.117)

где U определяется согласно формуле (6.107). Оценка ваимного спектра находится по формуле

Фиг. 6.41. Спектральная плотность мощности шума на выходе фильтра нижних частот, полученная с помощью периодограмм.

Фиг. 6.42. Спектральная плотность мощности шума на выходе колосового фильтра, полученная с помощью периодограмм.

Снова можно показать, что математическое ожидание зтой оценки равно свертке истинной взаимной спектральной плотности мощности с квадратом ДПФ от окна.

На фиг. 6.41-6.43 иллюстрируется практическое использование метода модифицированных периодограмм для оценки спектральной плотности мощности окрашенного шума применительно к трем примерам, приведенным ранее на фиг. 6.35 — 6.37. В каждом из этих примеров спектральная плотность оценивалась на 513 равноотстоящих частотах от до , т. е. с использованием L-точечного БПФ (L = 1024).

Фиг. 6.43. Спектральная плотность мощности шума на выходе дифференциатора, полученная с помощью периодограмм.

Для взвешивания входной последовательности применялось окно Хэмминга. С целью уменьшения дисперсии оценки при фиксированном числе входных отсчетов был использован сдвиг подпоследовательностей на D = 512 отсчетов. Общее число подпоследовательностей равно 32, т. е. всего обрабатывалось N = 16384 отсчетов. На фиг. 6.41 изображен энергетический спектр шума на выходе фильтра нижних частот (в логарифмическом масштабе), а на фиг. 6.42 и 6.43 — энергетические спектры шумовых последовательностей на выходах полосового фильтра и дифференциатора соответственно.

3. Выводы

Итак, для оценки энергетического спектра или корреляционных функций случайных процессов можно использовать два эффективных метода. В зависимости от конкретной задачи каждый из них может иметь преимущества относительно другого. Общим для обоих методов является то, что среднее значение оценки спектра сходится к определенной величине, а дисперсия оценки стремится к нулю, когда число обрабатываемых отсчетов стремится к бесконечности.