Main Content

상호 스펙트럼과 크기 제곱 코히어런스

이 예제에서는 상호 스펙트럼을 사용하여 이변량 시계열에서 정현파 성분 간의 위상 지연을 구하는 방법을 보여줍니다. 이 예제에서는 또한 크기 제곱 코히어런스(Magnitude-Squared Coherence)를 사용하여 사인파 주파수에서 유의미한 주파수 영역 상관관계를 식별합니다.

이변량 시계열을 생성합니다. 각각의 시계열은 100Hz 주파수의 사인파와 200Hz 주파수의 사인파로 구성됩니다. 시계열은 가산성 백색 가우스 잡음(AWGN)에 묻혀 있고 1kHz로 샘플링되었습니다. x 시계열의 사인파는 모두 진폭이 1입니다. y 시계열의 100Hz 사인파는 진폭이 0.5이고, y 시계열의 200Hz 사인파는 진폭이 0.35입니다. y 시계열의 100Hz 사인파와 200Hz 사인파는 각각 π/4라디안과 π/2라디안만큼 위상 지연됩니다. 입력값 x의 선형 시스템이 잡음으로 인해 손상되어 y 시계열로 출력된다고 간주해 볼 수 있습니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 디폴트 설정으로 구성합니다.

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개 샘플로 설정합니다. 이 윈도우 길이는 100Hz 사인파의 10개 주기와 200Hz 사인파의 20개 주기를 포함합니다. 기본 해밍 윈도우와 함께 80개 샘플의 중첩을 사용합니다. Hz 단위의 출력 주파수를 얻도록 샘플 레이트를 명시적으로 입력합니다. 크기 제곱 코히어런스를 플로팅합니다. 크기 제곱 코히어런스가 100Hz와 200Hz에서 0.8보다 큽니다.

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

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

Figure contains an axes object. The axes object with title Magnitude-Squared Coherence, xlabel Frequency (Hz) contains an object of type line.

cpsd를 사용하여 xy의 상호 스펙트럼을 구합니다. 동일한 파라미터를 사용하여 코히어런스 추정값에 사용된 상호 스펙트럼을 구합니다. 코히어런스가 적은 경우 이 상호 스펙트럼을 무시합니다. 상호 스펙트럼의 위상을 플로팅하고 두 개의 시간값 사이에 유의미한 코히어런스를 가지는 주파수를 나타냅니다. 정현파 성분 간에 알려진 위상 지연의 위치를 나타냅니다. 100Hz 및 200Hz에서 상호 스펙트럼을 통해 추정된 위상 지연이 실제 값에 근접합니다.

[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

Figure contains an axes object. The axes object with title Cross Spectrum Phase, xlabel Frequency (Hz), ylabel Lag ( times pi blank rad) contains an object of type line.

참고 항목

| |