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

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

Создайте двухмерные временные ряды. Отдельные ряды состоят из двух синусоид с частотами 100 и 200 Гц. Серия встраивается в аддитивный белый Гауссов шум и дискретизируется с частотой дискретизации 1 кГц. Синусоиды в x-ряду обе имеют амплитуды, равные 1. Синусоида 100 Гц в y-серии имеет амплитуду 0,5, а синусоиду 200 Гц в y-серии - амплитуду 0,35. Частоты 100 Гц и 200 Гц в синусоиды 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 для всех частот, необходимо использовать среднюю оценку когерентности. Подходящими являются как перекрывающиеся средние уровни Welch (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