Перекрестный спектр и квадрат когерентности

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

Создайте двумерные временные ряды. Отдельные ряды состоят из двух синусоид с частотами 100 и 200 Гц. Ряды встраиваются в аддитивный белый Гауссов шум и производятся на уровне 1 кГц. Синусоиды в x-ряду оба имеют амплитуды, равные 1. Синусоида на 100 Гц в y-ряду имеет амплитуду 0.5, и синусоида на 200 Гц в y-ряду имеет амплитуду 0.35. Синусоиды на 200 Гц и на 100 Гц в y-ряду изолированы фазой $$\pi /4$$ радианы и $$\pi /2$$ радианы, соответственно. Можно думать о y-ряде как о поврежденном шумом выходе линейной системы с входом x. Установите генератор случайных чисел на настройки по умолчанию для восстанавливаемых результатов.

rng default

Fs = 1000;
t = 0:1/Fs:1-1/Fs;

x = cos(2*pi*100*t) + sin(2*pi*200*t) + 0.5*randn(size(t));
y = 0.5*cos(2*pi*100*t - pi/4) + 0.35*sin(2*pi*200*t - pi/2) + 0.5*randn(size(t));

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

Чтобы предотвратить получение оценки квадрата когерентности, которая является тождественно 1 для всех частот, необходимо использовать усредненное средство оценки когерентности. И Перекрытое усреднение сегмента валлийцев (WOSA) и методы мультизаострения являются соответствующими. mscohere реализует средство оценки WOSA.

Установите длину окна на 100 выборок. Эта длина окна содержит 10 периодов синусоиды на 100 Гц и 20 периодов синусоиды на 200 Гц. Используйте перекрытие 80 выборок с Окном Хэмминга по умолчанию. Введите частоту дискретизации явным образом, чтобы получить выходные частоты в Гц. Постройте квадрат когерентности. Квадрат когерентности больше 0.8 на уровне 100 и 200 Гц.

[Cxy,F] = mscohere(x,y,hamming(100),80,100,Fs);

plot(F,Cxy)
title('Magnitude-Squared Coherence')
xlabel('Frequency (Hz)')
grid

Получите перекрестный спектр X и Y с помощью cpsd. Используйте те же параметры, чтобы получить перекрестный спектр, который вы использовали в оценке когерентности. Пропустите перекрестный спектр, когда когерентность будет мала. Постройте фазу перекрестного спектра и укажите на частоты со значительной когерентностью между этими двумя разами. Отметьте известные задержки фазы между синусоидальными компонентами. На уровне 100 Гц и 200 Гц, задержки фазы, оцененные от перекрестного спектра, близко к истинным значениям.

[Pxy,F] = cpsd(x,y,hamming(100),80,100,Fs);

Pxy(Cxy < 0.2) = 0;

plot(F,angle(Pxy)/pi)
title('Cross Spectrum Phase')
xlabel('Frequency (Hz)')
ylabel('Lag (\times\pi rad)')
grid

Смотрите также

cpsd | mscohere | pwelch