1.7. Численное представление случайных сигналов

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

Основу получения случайных сигналов составляют генераторы случайных числовых последовательностей. Такие генераторы входят в математическое обеспечение компьютерных математических программ, предназначенных для физических исследований. Например, в пакете программ Mathcad есть ряд генераторов случайных последовательностей, описанных в [1] (стр. 299-300). Каждое число случайной последовательности является статистически независимым от других. Такая последовательность представляет собой "белый шум". Закон распределения и его параметры определяются выбором конкретной программы. В программу, описанную в [1], входит набор более 10 различных законов распределения, включающих нормальное, Пуассона, равномерное и ряд других.

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

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

(кликните для просмотра скана)

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

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

Результат показан на графике (см. рис. 1.15). Сплошная линия - гистограмма, относящаяся к сигналу, а пунктир - гистограмма сигнала на выходе согласованного фильтра. Графики даны в логарифмическом масштабе в децибелах. Формула, позволяющая получить такой масштаб, приведена на рисунке. Она позволяет логарифмировать данные, превышающие установленный порог (десять в степени Для всех остальных данных, в числе которых могут быть и нули, принято значение логарифма установленного порога.

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

Случайный сигнал с заданной формой спектра. Остановимся на методах получения случайных последовательностей с заданной формой спектра [2].

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

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

(кликните для просмотра скана)

Далее полученная форма спектра преобразуется к той, в которой спектр отображается в пакете программ Mathcad (высокие частоты - в центре, а низкие - по краям). Ниже приведена программа, обеспечивающая вставление в нужное место дискретной составляющей.

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

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

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

Таким образом получается спектр, включающий только целые частоты в силу того, что задаются значения частот лишь в дискретных точках отсчета. Каждая спектральная компонента обладает целочисленным значением частоты. В этом можно убедиться с помощью программы, показанной на рис. 1.17, которая вначале повторяет программу, приведенную на рис. 1.16, кроме того, что в этой программе формируется только одна реализация ансамбля случайных функций с заданной формой спектра. Сформированный спектр показан на рисунке в логарифмическом масштабе в виде функции (левый график). На правом графике - спектр, получающийся при сдвиге спектра левого графика на 0,43 отсчетной точки. Такой сдвш спектра выполнен в функции Крохотное смещение частоты вызывает резкое изменение формы спектра, что характерно для целых частот. Такой фокус мы уже демонстрировали на рис. 1.1 и 1.2. Теперь то же самое - на спектре произволь ной формы.

Зададимся вопросом: каким образом "размыть" спектр так, чтобы получались дробные частоты. Это может быть поучительно, позволит подчеркнуть особенностг численного спектрального анализа, и по-настоящему "размытые" спектры понадобятся нам в дальнейшем.

(кликните для просмотра скана)

Программа, позволяющая "размывать" спектры в пределах одной отсчетной точки, т. е. формировать дробные частоты, показана на рис. 1.18. Эта программа полностью повторяет предыдущую (см. рис. 1.17), для увеличенного в 8 раз массива данных. В результате мы вновь получаем "целочисленный" спектр, который без изменения формы превращается в дробный посредством чрезвычайно простого приема, показанного в программе под названием "искомая функция". Прием заключается в том, что, вводя новую переменную обладающую в 8 раз меньшей областью значений, мы сокращаем область существования нашей функции в 8 раз, причем число значений спектра функции тоже станет в 8 раз меньше, а сама форма спектра при этом не изменится, останется прежней, какой она была у первоначальной функции, занимавшей в 8 раз большую область. Это свойство преобразования Фурье, не имеющее отношения к тому, делается это преобразование численно или иным способом. А теперь учтем специфику численного преобразования. Она состоит в том, что одну и ту же форму спектра надо отобразить сначала большим числом отсчетов, а затем числом отсчетов, уменьшенным в 8 раз. Если в первом случае все точки спектра точно попали на свои отсчеты, то во втором случае это уже невозможно, значения спектральных компонент попадают с неизбежностью в промежуточные точки, порождая дробные частоты. В этом позволяет убедиться приведенная программа. Теперь небольшое смещение положения спектра как целого уже никак не сказывается на виде спектра.

На следующих двух программах показано как целый и "размытый" узкополосные процессы проходят через фильтр, сформированный в виде "спектрального окна". На рис. 1.19 показана программа прохождения через фильтр спектра, использующего только целые частоты. Вначале эта программа повторяет предыдущую, показанную на рис. 1.18, за исключением расширения области в 8 раз. Эта область значений не расширяется, в результате чего спектр получается только целым. Сформированный узкополосный спектр смещается по частоте при изменении значения переменной Фильтр "спектральное окно" формируется точно так, как это показано и описано в программе на рис. 1.8. Полученный выход фильтра в виде функции времени и смещения частоты усредняется по интенсивности по времени Получившийся результат изменения интенсивности входного и выходного сигналов при смещении частоты узкополосного сигнала показан на рисунках. На других двух рисунках показаны спектр входного сигнала при одном фиксированном значении его сдвига и форма частотного окна фильтра.

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

На рис. 1.20 показано то же самое, но при "размытом" спектре узкополосного сигнала. Изменение спектра сигнала достигнуто путем расширения области значений в 8 раз. Теперь сигнал частично проходит через "спектральное окно", даже находясь далеко от него. Это обстоятельство следует учитывать при фильтрации сигналов.

(кликните для просмотра скана)

(кликните для просмотра скана)

(кликните для просмотра скана)

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

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

Рассмотрим решение такой задачи на конкретном примере восстановления формы морского ветрового волнения по его спектру. Такая задача решена в программе, изображенной на рис. 1.21. Вначале идут диапазонные переменные по две для двух измерений спектра и функции. Затем приведена формула морского волнения с включаемыми в нее многочисленными параметрами, которые можно изменять в зависимости от условий. Этот спектр приводится к виду, в котором принято представлять спектры в программе Mathcad, когда высокие частоты помещаются внутрь спектральной области. Все частоты являются целыми положительными числами. Далее формируется матрица случайных чисел, имеющих некоторые характеристики в виде среднего значения и дисперсии. В данном случае использовано нормальное распределение вероятностей. Спектры случайных реализаций чисел умножаются на двумерный спектр морского волнения, чтобы придать последнему случайный характер. Для получения реализаций самого случайного морского волнения теперь достаточно использовать обратное преобразование Фурье. И вот одна из реализаций морского волнения перед вами! Для получения следующей реализации, независимой от существующей, достаточно еще раз вычислить программу по команде: "calculate worksheet".

Сопоставление статистик. Одна из основных целей этого раздела (он основан на результатах работы [2]) - показать, насколько приведенные выше элементарные расчеты способны существенно облегчить численное моделирование сложных физических задач.

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

(кликните для просмотра скана)

На нем приведена картина поля, полученного в результате точного расчета. Чтобы иметь возможность рассчитывать точно, надо существенно упрощать задачу, отбрасывая ряд параметров, условий, сокращая дистанции и т.п. На этом рисунке поле сосчитано не в волноводе, а в свободном пространстве в непосредственной близости (на расстоянии нескольких длин волн) от рассеивающей поверхности ограниченных размеров (квадрат со стороной 128 отсчетных точек). На поверхности поле принимает всего два значения нуль или единица по случайному закону (использован генератор случайного процесса в пакете Mathcad 6.0 plus). С помощью этих упрощений расчет удалось провести точно. На рис. 1.22 приведен угловой, или пространственный, спектр рассеянного волнового поля. Показано только рассеянное поле. Прямая волна, вызывающая рассеяние, исключена.

Рис. 1.22. Картина волнового поля, полученная точным расчетом

Как видно из рисунка, спектр поля сильно изрезан максимумами и минимумами. Видно также, что угловой спектр ограничен. Ограничение связано с тем, что поле распространяется лишь в пределах ±90° от нормали к поверхности. Поле, рассеянное под более крутыми углами, не является распространяющимся и быстро затухает.

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

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

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

Во-первых, в силу физических условий интерференционная изрезанность отдельных реализаций получается максимальной, ограничиваемой только общей шириной спектра сигнала. Используемый нами метод построения случайных сигналов обеспечивает получение именно таких реализаций. Во-вторых, реальный процесс является суммой большого числа слагаемых, что при условии их статистической независимости и в силу центральной предельной теоремы теории вероятностей обеспечивает вполне определенные статистические свойства совокупности интерференционных картин, рассматриваемых как случайный процесс [5]. Такой процесс должен иметь нормальное или гауссово распределение, которое (в силу равенства нулю среднего значения) определяется всего одним параметром - дисперсией. Моделируемые нами реализации получаются в результате преобразования Фурье, которое фактически представляет собою суммирование. Следовательно, в силу центральной предельной теоремы можно ожидать получение нормального распределения. Однако применение центральной предельной теоремы требует отдельного строгого математического обоснования, которым мы пока не располагаем, но надеемся, что физические условия при максимальной изрезанности интерференционной картины обеспечивают нормализацию рассматриваемого нами процесса.

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

Для вычисления распределения поля с учетом интерференции (точное решение задачи) используется метод примененный и изложенный в работе [6]. Вычисляется сигнал, рассеянный ограниченным участком неровной поверхности, в массиве данных 128x128 числовых значений. Плоскость, в которой вычисляется поле, отстоит от исходной плоскости на несколько длин волн. Поверхность освещается плоской волной, падающей на поверхность под малым углом к ней. Этот расчет включает три операции: 1) задание комплексной амплитуды волнового монохроматического поля на выбранном участке поверхности; 2) получение двумерного спектра этого поля; 3) умножение полученного спектра на комплексную двумерную функцию, называемую частотной характеристикой свободного пространства. Вид этой функции известен из теории. Этот расчет назовем сокращенно расчетом А.

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

Картина спектра при расчете А представлена выше на рис. 1.22. На рис. 1.23 показан спектр, получающийся в результате расчета Б.

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

Рис. 1.23. Картина волнового поля, полученная упрощенно

Рис. 1.24. Эмпирические распределения вероятностей значений полей, показанных на рис. 1.22, 1.23, и нормальное распределение вероятностей

Из полученных точного (А) и модельного (Б) спектров полей выберем по 1024 значения действительных частей комплексных спектров (они могут иметь нормальное распределение с нулевым средним) с одинаковыми номерами. Эти последовательности уравняем по интенсивностям и исследуем на соответствие нормальному закону распределения с помощью критерия согласия Колмогорова [7]. Для этого эмпирические гистограммы чисел в последовательностях А и Б проинтегрируем по интервалам значений гистограммы до изменяющегося верхнего предела. В результате получим эмпирические функции распределения последовательностей А и Б, приведенные на рис. 1.24. На том же рисунке показана функция распределения нормально распределенной случайной величины, полученная путем интегрирования выражения

где - среднеквадратическое отклонение, найденное по эмпирическим данным.

Определим значения максимальных отклонений эмпирических функций распределения последовательностей А и Б от теоретической. Согласно теореме Колмогорова [7, стр. 284] и таблице [7, стр. 384], максимальные отклонения эмпирической функции распределения от теоретической могут превысить полученные нами значения для А и Б с вероятностями 0,97 и 0,999 соответственно. Эти вероятности так близки к единице, что такие события осуществляются почти всегда. Значит, практически нельзя получить лучшего совпадения эмпирических законов распределения с теоретическим, чем то, которое представлено на рис. 1.24. На этом основании критерий согласия Колмогорова позволяет считать оба эмпирических закона соответствующими гауссову закону распределения. Заметим, что исходная поверхность, которая задавалась численно, имела не гауссово, а равномерное распределение.

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