cpsd

상호 파워 스펙트럼 밀도

구문

pxy = cpsd(x,y)

pxy = cpsd(x,y,window)

pxy = cpsd(x,y,window,noverlap)

pxy = cpsd(x,y,window,noverlap,nfft)

pxy = cpsd(___,'mimo')

[pxy,w] = cpsd(___)

[pxy,f] = cpsd(___,fs)

[pxy,w] = cpsd(x,y,window,noverlap,w)

[pxy,f] = cpsd(x,y,window,noverlap,f,fs)

[___] = cpsd(x,y,___,freqrange)

cpsd(___)

설명

pxy = cpsd(x,y)는 수정된 주기도를 평균화하는 Welch 스펙트럼 추정 방법을 사용하여 두 개의 이산시간 신호 x와 y에 대한 상호 파워 스펙트럼 밀도(CPSD)를 추정합니다.

• x와 y가 모두 벡터인 경우 길이가 같아야 합니다.

• 두 신호 중 하나는 행렬이고 다른 하나는 벡터인 경우, 벡터의 길이는 행렬에 포함된 행 개수와 동일해야 합니다. 이 함수는 벡터를 확장하고 열마다 상호 파워 스펙트럼 밀도 추정값을 가진 행렬을 반환합니다.

• x와 y가 행의 개수는 동일하지만 열의 개수는 다른 행렬인 경우, cpsd는 입력 열의 모든 조합에 대한 상호 파워 스펙트럼 밀도 추정값을 포함하는 3차원 배열 pxy를 반환합니다. pxy의 각 열은 x의 열에 해당하고, 각 페이지는 y의 열에 해당합니다. 즉, pxy(:,m,n) = cpsd(x(:,m),y(:,n))입니다.

• x와 y가 같은 크기의 행렬인 경우, cpsd는 각 열에 대해 연산을 수행합니다. 즉, pxy(:,n) = cpsd(x(:,n),y(:,n))입니다. 다중 입력/다중 출력 배열을 구하려면 인수 목록에 'mimo'를 추가하십시오.

x와 y가 실수인 경우, cpsd는 단측 CPSD를 반환합니다. x 또는 y가 복소수인 경우, cpsd는 양측 CPSD를 반환합니다.

pxy = cpsd(x,y,window)는 window를 사용하여 x와 y를 세그먼트로 나누고 윈도우를 적용합니다.

pxy = cpsd(x,y,window,noverlap)은 인접 세그먼트 간에 noverlap개의 샘플이 중첩되도록 적용합니다.

pxy = cpsd(x,y,window,noverlap,nfft)는 nfft개 샘플링 점을 사용하여 이산 푸리에 변환을 계산합니다.

pxy = cpsd(___,'mimo')는 상호 파워 스펙트럼 밀도 추정값으로 구성된 다중 입력/다중 출력 배열을 계산합니다. 이 구문은 위에 열거된 구문에 사용할 수 있습니다.

[pxy,w] = cpsd(___)는 상호 파워 스펙트럼 밀도가 추정된 정규화 주파수로 구성된 벡터 w를 반환합니다.

[pxy,f] = cpsd(___,fs)는 상호 파워 스펙트럼 밀도가 추정된 주파수(샘플 레이트 fs의 값으로 표현됨)로 구성된 벡터 f를 반환합니다. fs는 cpsd에 대한 6번째 숫자형 입력값이어야 합니다. 샘플 레이트를 입력하고 위에 열거된 옵션 인수의 디폴트 값을 그대로 사용하려면 이 인수를 빈 값 []로 지정하십시오.

[pxy,w] = cpsd(x,y,window,noverlap,w)는 w에 지정된 정규화 주파수에서의 상호 파워 스펙트럼 밀도 추정값을 반환합니다.

[pxy,f] = cpsd(x,y,window,noverlap,f,fs)는 f에 지정된 주파수에서의 상호 파워 스펙트럼 밀도 추정값을 반환합니다.

[___] = cpsd(x,y,___,freqrange)는 freqrange로 지정된 주파수 범위에 대한 상호 파워 스펙트럼 밀도 추정값을 반환합니다. freqrange에 유효한 옵션은 'onesided', 'twosided', 'centered'입니다.

cpsd(___)에 출력 인수를 지정하지 않으면 현재 Figure 창에 상호 파워 스펙트럼 밀도 추정값을 플로팅합니다.

예제

유색 잡음 신호의 상호 파워 스펙트럼 밀도

두 개의 유색 잡음 신호를 생성하고 이에 대한 상호 파워 스펙트럼 밀도를 플로팅합니다. 길이가 1024인 FFT와 세그먼트 간에 50% 중첩이 있는 500개 점을 갖는 삼각 윈도우를 지정합니다.

rng default
r = randn(2e4,1);

hx = fir1(30,0.2,rectwin(31));
x = filter(hx,1,r);

hy = ones(1,10);
y = filter(hy,1,r);

win = triang(500);
nov = 250;

cpsd(x,y,win,nov,1024)

SISO 및 MIMO 상호 파워 스펙트럼 밀도

1초 동안 1kHz로 샘플링된 2채널 정현파 2개를 생성합니다. 첫 번째 신호의 채널은 200Hz와 300Hz의 주파수를 가집니다. 두 번째 신호의 채널은 300Hz와 400Hz의 주파수를 가집니다. 두 신호 모두 단위 분산 백색 가우스 잡음에 묻혀 있습니다.

fs = 1e3;
t = (0:1/fs:1-1/fs)';

q = 2*sin(2*pi*[200 300].*t);
q = q+randn(size(q));

r = 2*sin(2*pi*[300 400].*t);
r = r+randn(size(r));

두 신호의 상호 파워 스펙트럼 밀도를 계산합니다. 256개 샘플로 구성된 바틀렛 윈도우를 사용하여 신호를 세그먼트로 나누고 세그먼트에 윈도우를 적용합니다. 인접 세그먼트 간에 128개 샘플이 중첩되도록 지정하고 2048개 DFT 점을 지정합니다. cpsd의 내장된 기능을 사용하여 결과를 플로팅합니다.

cpsd(q,r,bartlett(256),128,2048,fs)

기본적으로 cpsd는 동일한 크기의 행렬 입력값에 대해 열 단위로 연산을 수행합니다. 각 채널은 원래 정현파의 주파수에서 피크에 도달합니다.

계산을 반복하되, 이번에는 'mimo'를 인수 목록에 추가합니다.

cpsd(q,r,bartlett(256),128,2048,fs,'mimo')

'mimo' 옵션을 사용하여 호출되는 경우, cpsd는 모든 입력 열 조합에 대한 상호 파워 스펙트럼 밀도 추정값을 포함하는 3차원 배열을 반환합니다. q의 두 번째 채널과 r의 첫 번째 채널에 대한 추정값은 300Hz의 공통 주파수에서 더 높은 피크에 도달하는 것을 보여줍니다.

지연된 정현파의 상호 스펙트럼 위상

296ms 동안 1kHz로 샘플링된 두 개의 100Hz 정현파 신호를 생성합니다. 정현파 중 하나는 다른 정현파보다 2.5ms만큼 지연되고, 이는 π/2의 위상 지연에 상응합니다. 두 신호 모두 분산이 1/4²인 백색 가우스 잡음에 묻혀 있습니다.

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

x = cos(2*pi*t*100)+0.25*randn(size(t));
tau = 1/400;
y = cos(2*pi*100*(t-tau))+0.25*randn(size(t));

상호 파워 스펙트럼 밀도의 크기를 계산하고 플로팅합니다. cpsd에 대한 디폴트 설정을 사용합니다. 크기는 신호 간에 유의미한 코히어런스가 있는 주파수에서 피크에 도달합니다.

cpsd(x,y,[],[],[],Fs)

크기 제곱 코히어런스 함수와 상호 스펙트럼의 위상을 플로팅합니다. 코히어런스가 높은 주파수 지점에서의 세로 좌표는 정현파 간 위상 지연에 해당합니다.

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

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

subplot(2,1,1)
plot(F,Cxy)
title('Magnitude-Squared Coherence')

subplot(2,1,2)
plot(F,angle(Pxy))

hold on
plot(F,2*pi*100*tau*ones(size(F)),'--')
hold off

xlabel('Hz')
ylabel('\Theta(f)')
title('Cross Spectrum Phase')

지수 시퀀스에 대한 상호 파워 스펙트럼 밀도

$$N$$개 샘플로 구성된 두 개의 지수 시퀀스 $$x_a=a^n$$과 $$x_b=b^n$$($$n\ge 0$$임)을 생성합니다. $$a=0.8$$, $$b=0.9$$ 및 작은 $$N$$값을 지정하여 유한한 크기의 효과를 확인합니다.

N = 10;
n = 0:N-1;

a = 0.8;
b = 0.9;

xa = a.^n;
xb = b.^n;

정규화 주파수의 전체 구간 $$[-\pi ,\pi ]$$에 대해 시퀀스의 상호 파워 스펙트럼 밀도를 계산하고 플로팅합니다. 길이가 $$N$$이고 세그먼트 간 중첩이 없는 사각 윈도우를 지정합니다.

w = -pi:1/1000:pi;
wind = rectwin(N);
nove = 0;

[pxx,f] = cpsd(xa,xb,wind,nove,w);

두 시퀀스의 상호 파워 스펙트럼은 큰 $$N$$에 대한 해석적 표현을 갖습니다.

$$R(\omega )=\frac{1}{1-a{e}^{-j\omega }}\frac{1}{1-b{e}^{j\omega }}.$$

이 표현을 $$2\pi N$$으로 나누어 상호 파워 스펙트럼 밀도로 변환합니다. 결과를 비교합니다. cpsd 결과의 리플은 윈도우 적용으로 인한 것입니다.

nfac = 2*pi*N;

X = 1./(1-a*exp(-1j*w));
Y = 1./(1-b*exp( 1j*w));
R = X.*Y/nfac;

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

$$N=25$$를 사용하여 계산을 반복합니다. $$N$$이 100 정도의 작은 숫자인 경우 이 곡선은 6개 유효 자릿수까지 일치합니다.

N = 25;
n = 0:N-1;
xa = a.^n;
xb = b.^n;

wind = rectwin(N);

[pxx,f] = cpsd(xa,xb,wind,nove,w);
R = X.*Y/(2*pi*N);

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

발신음 인식

상호 파워 스펙트럼 밀도를 사용하여 많이 손상된 음색을 식별할 수 있습니다.

디지털 전화에서 숫자나 기호를 누를 때 생성되는 소리 신호는 두 개의 다른 그룹에서 가져온 주파수를 갖는 정현파의 합계입니다. 각각의 음색 쌍은 저주파수 무리의 주파수 하나(697Hz, 770Hz, 852Hz 또는 941Hz)와 고주파수 무리의 주파수 하나(1209Hz, 1336Hz 또는 1477Hz)를 포함합니다.

모든 기호에 대응되는 신호를 생성합니다. 1/2초 동안 4kHz로 각 음색을 샘플링합니다. 참조 테이블을 준비합니다.

fs = 4e3;
t = 0:1/fs:0.5-1/fs;

nms = ["1";"2";"3";"4";"5";"6";"7";"8";"9";"*";"0";"#"];

ver = [697 770 852 941];
hor = [1209 1336 1477];

v = length(ver);
h = length(hor);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones(:,idx) = tone;
    end
end

각 신호의 Welch 주기도를 플로팅하고 성분 주파수에 대한 주석을 표시합니다. 200개 샘플로 구성된 해밍 윈도우를 사용하여 신호를 중첩되지 않은 세그먼트로 나누고 세그먼트에 윈도우를 적용합니다.

[pxx,f] = pwelch(tones,hamming(200),0,[],fs);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        ax = subplot(v,h,idx);
        plot(f,10*log10(pxx(:,idx)))
        ylim([-80 0])
        title(nms(idx))
        tx = [ver(k);hor(l)];
        ax.XTick = tx;
        ax.XTickLabel = int2str(tx);
    end
end

숫자 8을 누를 때 생성된 신호가 잡음이 있는 채널을 통해 전송됩니다. 수신 신호가 너무 손상되었기 때문에 그림만 보아서는 숫자를 식별할 수 없습니다.

mys = sum(sin(2*pi*[ver(3);hor(2)].*t))'+5*randn(size(t'));

% To hear, type soundsc(mys,fs)

손상된 신호와 기준 신호의 상호 파워 스펙트럼 밀도를 계산합니다. 형태 인자 β = 5에다 512개 샘플로 구성된 카이저 윈도우를 사용하여 신호에 윈도우를 적용합니다. 각 스펙트럼의 크기를 플로팅합니다.

[pxy,f] = cpsd(mys,tones,kaiser(512,5),100,[],fs);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        ax = subplot(v,h,idx);
        plot(f,10*log10(abs(pxy(:,idx))))
        ylim([-80 0])
        title(nms(idx))
        tx = [ver(k);hor(l)];
        ax.XTick = tx;
        ax.XTickLabel = int2str(tx);
    end
end

손상된 신호에 포함된 숫자는 가장 높은 피크와 RMS 값이 가장 높은 스펙트럼을 갖습니다.

[~,loc] = max(rms(abs(pxy)));

digit = nms(loc)

digit = 
"8"

입력 인수

x, y — 입력 신호

벡터 | 행렬

입력 신호로, 벡터나 행렬로 지정됩니다.

예: cos(pi/4*(0:159))+randn(1,160)은 백색 가우스 잡음을 포함하는 정현파를 지정합니다.

데이터형: single | double
복소수 지원 여부: 예

window — 윈도우

정수 | 벡터 | []

윈도우로, 정수로 지정되거나 행 벡터 또는 열 벡터로 지정됩니다. window를 사용하여 신호를 여러 세그먼트로 나눕니다.

• window가 정수이면 cpsd는 x와 y를 길이가 window인 세그먼트로 나누고 각 세그먼트에 해당 길이의 해밍 윈도우를 적용합니다.

• window가 벡터이면 cpsd는 x와 y를 벡터와 길이가 같은 세그먼트로 나누고 window를 사용하여 각 세그먼트에 윈도우를 적용합니다.

x와 y의 길이를 noverlap개의 샘플이 중첩된 정수 개수의 세그먼트로 정확히 나눌 수 없는 경우 신호가 이에 맞게 적절하게 잘립니다.

window를 빈 값으로 지정하면, cpsd는 해밍 윈도우를 사용하여 x와 y를 noverlap개의 중첩 샘플을 가지는 8개의 세그먼트로 나눕니다.

사용 가능한 윈도우 목록은 윈도우 항목을 참조하십시오.

예: hann(N+1) 및 (1-cos(2*pi*(0:N)'/N))/2는 모두 길이가 N + 1인 핸 윈도우를 지정합니다.

데이터형: single | double

noverlap — 중첩 샘플의 개수

양의 정수 | []

중첩(Overlap) 샘플의 개수로, 양의 정수로 지정됩니다.

• window가 스칼라이면 noverlap은 window보다 작아야 합니다.

• window가 벡터이면 noverlap은 window의 길이보다 작아야 합니다.

noverlap을 빈 값으로 지정하면 cpsd는 세그먼트 간에 50% 중첩을 만드는 값을 사용합니다. 세그먼트 길이가 지정되지 않은 경우 이 함수는 noverlap을 ⌊N/4.5⌋로 설정합니다. 여기서 N은 입력 신호와 출력 신호의 길이입니다.

데이터형: double | single

nfft — DFT를 적용할 지점의 개수

양의 정수 | []

DFT를 적용할 지점의 개수로, 양의 정수로 지정됩니다. nfft를 빈 값으로 지정하면 cpsd는 파라미터를 max(256,2^p)로 설정합니다. 여기서, 길이가 N인 입력 신호의 경우 p = ⌈log₂ N⌉입니다.

데이터형: single | double

fs — 샘플 레이트

양의 스칼라

샘플 레이트로, 양의 스칼라로 지정됩니다. 샘플 레이트는 단위 시간당 샘플 개수입니다. 시간 단위가 초이면 샘플 레이트의 단위는 Hz입니다.

w — 정규화 주파수

벡터

정규화 주파수로, 적어도 2개의 요소를 가진 행 벡터나 열 벡터로 지정됩니다. 정규화 주파수는 rad/sample을 단위로 사용합니다.

예: w = [pi/4 pi/2]

데이터형: double | single

f — 주파수

벡터

주파수로, 적어도 2개의 요소를 가진 행 벡터나 열 벡터로 지정됩니다. 주파수는 단위 시간당 주기를 단위로 사용합니다. 단위 시간은 샘플 레이트 fs로 지정됩니다. fs의 단위가 샘플/초이면 f는 Hz 단위입니다.

예: fs = 1000; f = [100 200]

데이터형: double | single

freqrange — 상호 파워 스펙트럼 밀도 추정값을 구할 주파수 범위

'onesided' | 'twosided' | 'centered'

상호 파워 스펙트럼 밀도 추정값을 구할 주파수 범위로, 'onesided', 'twosided' 또는 'centered'로 지정됩니다. 디폴트 값은 실수 값을 갖는 신호의 경우 'onesided'이고, 복소수 값을 갖는 신호의 경우 'twosided'입니다.

• 'onesided' — 실수 값을 갖는 두 개의 입력 신호 x와 y의 상호 파워 스펙트럼 밀도에 대한 단측 추정값을 반환합니다. nfft가 짝수이면 pxy가 nfft/2 + 1개의 행을 가지고 구간 [0,π] rad/sample에 대해 계산됩니다. nfft가 홀수이면 pxy가 (nfft + 1)/2개의 행을 가지고 구간은 [0,π) rad/sample이 됩니다. fs를 지정하면 이에 대응되는 구간은 nfft가 짝수인 경우 [0,fs/2] cycles/unit time이고 nfft가 홀수인 경우 [0,fs/2) cycles/unit time입니다.

• 'twosided' — 실수 값 또는 복소수 값을 갖는 두 개의 입력 신호 x와 y의 상호 파워 스펙트럼 밀도에 대한 양측 추정값을 반환합니다. 이 경우, pxy는 nfft개의 행을 가지고 구간 [0,2π) rad/sample에 대해 계산됩니다. fs를 지정하면 구간은 [0,fs) cycles/unit time이 됩니다.

• 'centered' — 실수 값 또는 복소수 값을 갖는 두 개의 입력 신호 x와 y의 상호 파워 스펙트럼 밀도에 대해 중심이 맞춰진 양측 추정값을 반환합니다. 이 경우, pxy는 nfft개의 행을 가지고, nfft가 짝수이면 구간 (–π,π] rad/sample에 대해 계산되며, nfft가 홀수이면 구간 (–π,π) rad/sample에 대해 계산됩니다. fs를 지정하면 이에 대응되는 구간은 nfft가 짝수인 경우 (–fs/2, fs/2] cycles/unit time이고 nfft가 홀수인 경우 (–fs/2, fs/2) cycles/unit time입니다.

주파수로 구성된 벡터 w 또는 f와 freqrange를 입력값으로 지정하면 cpsd 함수는 freqrange를 무시합니다.

출력 인수

pxy — 상호 파워 스펙트럼 밀도

벡터 | 행렬 | 3차원 배열

상호 파워 스펙트럼 밀도로, 벡터, 행렬 또는 3차원 배열로 반환됩니다.

w — 정규화 주파수

벡터

정규화 주파수로, 실수 값 열 벡터로 반환됩니다.

• pxy가 단측값인 경우 w의 구간은 nfft가 짝수이면 [0,π]이고, nfft가 홀수이면 [0,π)입니다.

• pxy가 양측값인 경우 w의 구간은 [0,2π)입니다.

• pxy가 DC 기준인 경우 w의 구간은 nfft가 짝수이면 (–π,π]이고, nfft가 홀수이면 (–π,π)입니다.

데이터형: double | single

f — 주파수

벡터

주파수로, 실수 값 열 벡터로 반환됩니다.

데이터형: double | single

세부 정보

상호 파워 스펙트럼 밀도

상호 파워 스펙트럼 밀도는 단위 주파수당 전력 분포이고 다음과 같이 정의됩니다.

$$P_{xy}(\omega )={\displaystyle \sum _{m=-\infty }^{\infty }{R}_{xy}(m)e^{-j\omega m}.}$$

상호상관 시퀀스는 다음과 같이 정의됩니다.

$$R_{xy}(m)=E\left\{x_{n+m}{y}_n^{\ast }\right\}=E\left\{x_n{y}_{n-m}^{\ast }\right\},$$

여기서 x_n과 y_n은 평균 0의 결합 정상 확률 과정이고 각각 –∞ < n < ∞, $$-\infty <n<\infty$$이며, E {· }는 기대값 연산자입니다.

알고리즘

cpsd는 수정된 주기도를 평균화하는 Welch 방법을 사용한 스펙트럼 추정입니다.

참고 문헌

[1] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

[2] Rabiner, Lawrence R., and B. Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975, pp. 414–419.

[3] Welch, Peter D. “The Use of the Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms.” IEEE^® Transactions on Audio and Electroacoustics, Vol. AU-15, June 1967, pp. 70–73.

확장 기능

C/C++ 코드 생성
MATLAB® Coder™를 사용하여 C 코드나 C++ 코드를 생성할 수 있습니다.

GPU 배열

Parallel Computing Toolbox™를 사용해 GPU(그래픽스 처리 장치)에서 실행하여 코드 실행 속도를 높일 수 있습니다.

이 함수는 GPU 배열을 완전히 지원합니다. 자세한 내용은 GPU에서 MATLAB 함수 실행하기 (Parallel Computing Toolbox) 항목을 참조하십시오.

버전 내역

R2006a 이전에 개발됨

R2024a: 단정밀도 가변 크기 창 입력값에 대한 코드 생성 지원

cpsd 함수는 코드 생성 시 단정밀도 가변 크기 창 입력값을 지원합니다.

R2023b: 단정밀도 데이터 사용

cpsd 함수는 단정밀도 입력값을 지원합니다.

참고 항목

mscohere | pburg | pcov | peig | periodogram | pmcov | pmtm | pmusic | pwelch | pyulear | tfestimate

도움말 항목

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