Исходный массив данных, представляющих собой реализацию из 2048 отсчетов многомерного случайного сигнала записан в файле «EEG1(2,3,4).txt» текстового формата. Количество колонок массива соответствует числу компонент случайного процесса (16), количество строк - числу отсчетов (2048). Интервал дискретизации, использованный при регистрации сигнала, равен 10 мс.

Задача 1. Оценка статистических характеристик реализации

случайного процесса.

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

Задача 2. Оценка плотности распределения реализации

случайного процесса.

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

Задача 3. Оценка корреляционных характеристик реализации

случайного процесса.

Оценить и построить корреляционные функции заданных в соответствии с номером варианта компонент СП.

Оценить и построить взаимные корреляционные функции компонент СП.

Задача 4. Оценка спектральных характеристик реализации

случайного процесса.

Исследовать спектральные характеристики (СПМ) случайного процесса методом Уэлча.

Длительности интервалов (сегментов сигнала) выбрать из таблицы.

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

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

Вариант №

Номера

каналов

Спектр.

окно

Длител.

интервала

3,4,5

Чебышева

128

1. Основные теоретические сведения

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

1.1 Оценка моментов

Оценка математического ожидания:

где w(x) -плотность распределения случайной величины.

Для стационарного эргодического процесса:

для дискретного сигнала:

В качестве оценки математического ожидания используют:

Качество оценки определяется степенью её разброса вокруг точного значения. Количественно оценка описывается доверительным интервалом значений вокруг точной величины, в который оценка попадает с заданной вероятностью. Под доверительной вероятностью понимается площадь под кривой внутри границ доверительного интервала. Чем уже интервал, тем лучше оценка. Количественной мерой ширины интервала служит дисперсия оценки:

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

Часто используется оценка максимального правдоподобия. Такая оценка основана на рассмотрении совместной плотности вероятности, как функции оцениваемого параметра.

Оценкой максимального правдоподобия называют:

- уравнение правдоподобия. Его решение: .

Пусть - независимые случайные величины, тогда

Функция правдоподобия:

Для нормального закона распределения: .

Оценкой дисперсии для нормального случайного процесса является:

1.2 Оценивание закона распределения

Существует два класса законов распределения:

· Интегральный

· Дифференциальный

Для оценки сигнала используется дифференциальный закон (плотность распределения).

Известно два способа оценки:

· Непараметрический (когда тип исходного распределения неизвестен)

а Гистограммный

а Парзена

а Разложения на базисные функции

а Полигонов Смирнова

а К-ближайших соседей

· Параметрический (известен закон, надо определить параметры)

Мы имеем дело с априорно непрерывной плотностью распределения вероятности, отличной от нуля на всем рассматриваемом интервале х. Объем выборки должен быть велик.

Наиболее простой и часто используемый метод - метод гистограмм.

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

Размеры интервалов делаются одинаковыми. Для выбора количества интервалов и их ширины нет общих правил.

Достоинства такой гистограммы:

+ Простота оценивания

+ Ясный физический смысл

Недостатки:

- При увеличении объема выборок, но неизменном количестве интервалов, оценка не сходится к точному значению закона распределения.

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

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

Второе условие: закон распределения на всем интервале отличен от нуля.

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

Существует два способа выполнения этих условий:

1. Сжатие подобласти таким образом, чтобы был обратно пропорционален корню квадратному из N (метод Парзена).

2. Подобласти так сжимаются, чтобы (метод к-ближайших соседей)

В случае обработки цифрового сигнала общий объем является целой степенью 2. Максимальное значение сигнала , где r-число разрядов используемого кода. Поэтому границы интервалов выбираются так, чтобы они совпадали с уровнями квантования, их количество 8-20.

1.3 Корреляционный анализ

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

1.3.1 Корреляционная функция

Корреляционная функция (КФ; английский термин -- correlation function, CF) детерминированного сигнала с конечной энергией представляет собой интеграл (в бесконечных пределах) от произведения двух копий сигнала, сдвинутых друг относительно друга на время ф:

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

1. Значение КФ при равно энергии сигнала, то есть интегралу от его квад-рата:

2. КФ является четной функцией своего аргумента ф: .

3. Значение КФ при ф = 0 является максимально возможным значением:

C ростом абсолютного значения ф КФ сигнала с конечной энергией затухает:

5. Если сигнал s(t) не содержит особенностей в виде дельта-функций, его КФ не может иметь разрывов (то есть обязана быть непрерывной функцией).

6. Если сигнал -- напряжение, то размерность его КФ равна В2 * с.

1.3.2 Взаимная корреляционная функция

Если КФ показывает степень сходства между сдвинутыми копиями одного и того же сигнала, то взаимная корреляционная функция (ВКФ; английский тер-мин -- cross-correlation function, CCF) позволяет измерить аналогичную величи-ну для сдвинутых экземпляров двух разных сигналов.

Общий вид формулы КФ сохраняется, но под интегралом стоит произведение двух разных сигналов, один из которых задержан на время ф.

Свойства ВКФ несколько отличаются от свойств КФ:

1. , где и -- энергии сигналов s1(t) и s2(t).

2. , то есть изменение знака ф равносильно взаимной перестанов-ке сигналов.

3. Значение ВКФ приф = 0 ничем не выделяется; максимум может быть распо-ложен в любом месте оси .

4. С ростом абсолютного значения т ВКФ сигналов с конечной энергией зату-хает:

5. Если сигналы s1(t) и s2(t) не содержат особенностей в виде дельта-функций, их ВКФ не может иметь разрывов (то есть обязана быть непрерывной функцией).

6. Если сигналы -- напряжение, то размерность их ВКФ равна В2 * с.

Для периодических сигналов понятие ВКФ обычно не применяется, хотя оно мо-жет быть введено в случае, если сигналы s1(t) и s2(t) имеют одинаковый период.

1.3.3 Связь между корреляционными функциями и спектрами

сигналов

Поскольку как корреляционные функции, так и спектры являются интегральны-ми преобразованиями анализируемых сигналов, логично предположить, что эти характеристики как-то связаны друг с другом. Для выявления этой связи подвергнем взаимную корреляционную функцию преобразованию Фурье, считая, что сигналы s1(t) и s2(t) имеют спектральные функции и :

Полученный результат очень прост: ВКФ связана преобразованием Фурье с так называемым взаимным спектром сигналов. Взаимный спектр для сигналов s1(t) и s2(t) представляет собой произведение их спектральных функций, одна из которых подвергнута комплексному сопряжению:

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

Приняв s1(t) = s2(t) = s(t), получаем аналогичный результат для КФ:

Итак, КФ сигнала связана преобразованием Фурье с квадратом модуля спектраль-ной функции, или с энергетическим спектром сигнала.

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

1.4 Спектр дискретного случайного процесса

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

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

Выполнять вычисления непосредственно по данной формуле неудобно, поэтому попробуем привести ее к более приемлемому виду. Для этого раскроем выраже-ние для квадрата модуля:

Суммируемые слагаемые зависят от разности индексов k и m, поэтому можно преобразовать двойную сумму в одиночную:

поскольку при любом l:

окончательно получаем

Это выражение представляет собой дискретный аналог теоремы Винера-Хинчина: «Спектр дискретного случайного процесса является преобразованием Фурье от его корреляционной функции».

1.4.1 Непараметрические методы расчета

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

Периодограмма

Этим термином -- периодограмма (periodogram) -- называется оценка спектральной плотности мощности, полученная по N отсче-там одной реализации случайного процесса согласно определению (1) (есте-ственно, не путем взятия предела, а усреднением конечного числа слагаемых). Таким образом, периодограмма рассчитывается по следующей формуле:

Деление на частоту дискретизации необходимо для получения оценки спек-тральной плотности мощности аналогового случайного процесса, восстановлен-ного по отсчетам x(k).

Если при расчете спектра используется весовая функция (окно) с коэффициен-тами w(k), формула (2) слегка модифицируется -- вместо числа отсчетов N b знаменателе должна стоять сумма квадратов модулей коэффициентов окна. По-лученная оценка спектра мощности называется модифицированной периодограм-мой (modified periodogram):

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

Метод Уэлча.

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

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

Итак, вычисления при использовании метода Уэлча (он называется еще методом усреднения модифицированных периодограмм -- averaged modified periodogram method) организуются следующим образом:

1. Вектор отсчетов сигнала делится на перекрывающиеся сегменты. Как прави-ло, на практике используется перекрытие на 50 %. Строго говоря, оптималь-ная степень перекрытия зависит от используемой весовой функции.

2. Каждый сегмент умножается на используемую весовую функцию.

3. Для взвешенных сегментов вычисляются модифицированные периодограммы.

4. Периодограммы всех сегментов усредняются.

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

1.4.2 Весовые функции (окна)

Для уменьшения растекания спектра при ДПФ применяются весовые функции (weighting functions), которые также называют окнами (windows). В этом случае перед расчетом ДПФ сигнал умножается на весовую функцию w(k), которая должна спадать к краям сегмента. Формула прямого ДПФ при использова-нии весовых функций принимает следующий вид:

Роль весовой функции в этой формуле можно рассматривать с различных точек зрения. Сначала проанализируем ситуацию во временной области. Если мы ис-пользуем весовую функцию, которая имеет максимум в середине (при k = N/2) и плавно спадает к краям (k = 0 и k = N-1), то это приведет к ослаблению эффектов, связанных с возникновением скачков сигнала при периодическом повторе-нии анализируемой конечной последовательности, и, таким образом, к уменьше-нию растекания спектра.

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

Если трактовать ДПФ как фильтрацию, при использовании весовой функции w(k) получаются частотные характеристики фильтров следующего вида:

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

Текст программы.

-----------------------------------

kurs.m

-----------------------------------

load EEG1.txt %загрузка исходного файла

k=input('Введите номер канала k='); %выбор каналов

l=input('Введите номер канала l=');

m=input('Введите номер канала m=');

x(:,1)=EEG1(:,k); %запись выбранных каналов

x(:,2)=EEG1(:,l); %в отдельную матрицу

x(:,3)=EEG1(:,m);

clc;

st %переход в меню

-----------------------------------

st.m

-----------------------------------

%меню выбора задания

function st

disp('Выберите задание')

disp('<a href="matlab:z1(x)">Задание 1. Оценка статистических характеристик реализации случайного процесса</a>')

disp('<a href="matlab:z2(x)">Задание 2. Оценка плотности распределения реализации случайного процесса</a>')

disp('<a href="matlab:z3(x)">Задание 3. Оценка корреляционных характеристик реализации случайного процесса</a>')

disp('<a href="matlab:z4(x)">Задание 4. Оценка спектральных характеристик реализации случайного процесса</a>')

-----------------------------------

z1.m

-----------------------------------

%Функция расчета математического ожидания и дисперсии

%для выбранных компонент реализации СП

function z1(x)

format short g

%Следующий цикл формирует 4 массива:

%mx - c мат.ожиданием

%varx - c дисперсией

%sostmx - вариативность мат. ожидания

%sostvarx - вариативность дисперсии

%1 столбец массива - временные интервалы

%2,3,4 столбцы соответствуют разным каналам

for ii=1:3

dlit=64;

for jj=1:4

for kk=1:4

a=(kk-1)*dlit+1;

b=kk*dlit;

tmpmx(kk)=mean(x(a:b,ii))

tmpvarx(kk)=var(x(a:b,ii))

end;

mx(jj,1)=dlit;

mx(jj,ii+1)=mean(tmpmx);

varx(jj,1)=dlit;

varx(jj,ii+1)=mean(tmpvarx);

sostmx(jj,1)=dlit;

sostmx(jj,ii+1)=var(tmpmx);

sostvarx(jj,1)=dlit;

sostvarx(jj,ii+1)=var(tmpvarx);

dlit=dlit*2;

end;

end

%Вывод данных

disp('-----------------------------------------------------');

disp(' | N канала ');

disp(' Длит. ------------------------------------')

disp(' интервала | 1 | 2 | 3 ');

disp('-----------------------------------------------------');

disp('Математическое ожидание');

disp(mx);

disp('Дисперсия');

disp(varx);

disp('Дисперсия оценки математического ожидания');

disp(sostmx);

disp('Дисперсия оценки дисперсии');

disp(sostvarx);

grafik=input('Хотите построить графики дисперсии оценок? (1 - Да, 2 - Нет) ');

if grafik == 1

figure;

subplot(2,1,1);

plot(sostmx(:,2:4));

title('Графики дисперсии оценки математического ожидания для разных временных интервалов','FontName','Times New Roman');

k=legend('Канал 1','Канал 2','Канал 3');

set(k,'FontName','Times New Roman');

subplot(2,1,2);

plot(sostvarx(:,2:4));

title('Графики дисперсии оценки дисперсии для разных временных интервалов','FontName','Times New Roman');

k=legend('Канал 1','Канал 2','Канал 3');

set(k,'FontName','Times New Roman');

end;

-----------------------------------

z2.m

-----------------------------------

%Функция построения гистограмм для заданных компонент

function z2(x)

disp('Введите количество интервалов (не менее 5)');

disp('0 - автоматический выбор');

nbins=input('Количество интервалов = ');

figure;

for ii=1:3

subplot(3,1,ii);

if nbins<5

hist(x(:,ii));

else

hist(x(:,ii),nbins);

end;

title(['Гистограмма',num2str(ii),'-го канала'],'FontName','Times New Roman');

end

-----------------------------------

z3.m

-----------------------------------

%Функция оценки корреляционных характеристик

function z3(x)

N=50; %размер корреляционной матрицы

[temp1, R1]=corrmtx(x(:,1) ,N); %Расчет АКФ (корреляционных

[temp2, R2]=corrmtx(x(:,2) ,N); %матриц)

[temp3, R3]=corrmtx(x(:,3) ,N);

k=0:N;

subplot(3,1,1);

stem(k, R1(1,:));

title('АКФ 1-го канала','FontName','Times New Roman');

subplot(3,1,2);

stem(k, R2(1,:));

title('АКФ 2-го канала','FontName','Times New Roman');

subplot(3,1,3);

stem(k, R3(1,:));

title('АКФ 3-го канала','FontName','Times New Roman');

R4 = xcorr(x(:,1),x(:,2)); %Расчет ВКФ

R5 = xcorr(x(:,1),x(:,3));

R6 = xcorr(x(:,2),x(:,3));

figure;

subplot(3,1,1);

stem(R4(2048:2098));

title('ВКФ 1-го и 2-го каналов','FontName','Times New Roman');

subplot(3,1,2);

stem(R5(2048:2098));

title('ВКФ 1-го и 3-го каналов','FontName','Times New Roman');

subplot(3,1,3);

stem(R6(2048:2098));

title('ВКФ 2-го и 3-го каналов','FontName','Times New Roman');

-----------------------------------

z4.m

-----------------------------------

%Функция расчета СПМ методом Уэлча

function z4(x)

nxt=1;

ur=input('Введите уровень боковых лепестков =');

dlit=input('Введите длительность сигнала =');

while nxt>0

long=input('Введите длину окна n=');

w=chebwin(long,ur); %окно Чебышева

figure;

plot(w);

figure;

for ii=1:3

subplot(3,2,(2*ii-1))

pwelch(x(1:dlit,ii),boxcar(dlit),[],[],100); %СПМ без выделяющей функции

title(['СПМ ',num2str(ii),' канала без выделяющей функции'],'FontName','Times New Roman');

subplot(3,2,2*ii)

pwelch(x(1:dlit,ii),w,[],[],100); %Сглажено с помощью окна Чебышева

title(['СПМ ',num2str(ii),' канала с использованием окна Чебышева'],'FontName','Times New Roman');

end;

nxt=input('Изменить длину окна и продолжить расчет? (1-Да, 0-Нет) ');

end;

Результаты выполнения программы

Оценка статистических характеристик реализации случайного процесса

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

Оценка плотности распределения реализации случайного процесса.

Оценка производится методом гистограмм. Рассмотрим гистограммы при количестве интервалов n =15.

Сигнал имеет распределение, близкое к Гауссовому в диапазоне примерно -500..500. Математическое ожидание стремится к нулю. Такое распределение характерно для помех в каналах связи.

Оценка корреляционных характеристик реализации случайного процесса

Отсчеты сигнала коррелированны слабо.

Оценка спектральных характеристик реализации случайного процесса

Длина окна n=32

Длина окна n=64

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

Список литературы

1. Дженкинс Г., Ваттс Д. Спектральный анализ и его приложения. Т. 1. М.: Мир, 1971. 316 с.

2. Марпл-мл. С.Л. Цифровой спектральный анализ и его приложения. М.: Мир, 1990. 584 с.

3. Сергиенко А.Б. Цифровая обработка сигналов. СПб.: Питер, 2002.

4. Шелухин О.И., Лукьянцев Н.Ф. Цифровая обработка и передача речи. М.: Радио и связь, 2000.