mscohere
크기 제곱 코히어런스
구문
설명
cxy = mscohere(x,y)x와 y의 크기 제곱 코히어런스 추정값 cxy를 구합니다.
x와y가 모두 벡터인 경우 길이가 같아야 합니다.두 신호 중 하나는 행렬이고 다른 하나는 벡터인 경우, 벡터의 길이는 행렬에 포함된 행 개수와 동일해야 합니다. 이 함수는 벡터를 확장하고 열마다 크기 제곱 코히어런스 추정값을 가진 행렬을 반환합니다.
x와y가 행 개수는 같지만 열 개수가 다른 행렬인 경우,mscohere는 다중 코히어런스 행렬을 반환합니다.cxy의 m번째 열은 모든 입력 신호와 m번째 출력 신호 간 상관도의 추정값을 포함합니다. 자세한 내용은 크기 제곱 코히어런스 항목을 참조하십시오.x와y가 같은 크기의 행렬인 경우,mscohere는 각 열에 대해 연산을 수행합니다. 즉,cxy(:,n) = mscohere(x(:,n),y(:,n))입니다. 다중 코히어런스 행렬을 구하려면'mimo'를 인수 목록에 추가하십시오.
mscohere(___)에 출력 인수를 지정하지 않으면 현재 Figure 창에 크기 제곱 코히어런스 추정값을 플로팅합니다.
예제
두 개의 유색 잡음 시퀀스 사이의 코히어런스 추정값을 계산하고 플로팅합니다.
백색 가우스 잡음으로 구성된 신호를 생성합니다.
r = randn(16384,1);
첫 번째 시퀀스를 생성하기 위해 신호에 대역통과 필터를 적용합니다. 0.2π rad/sample에서 0.4π rad/sample 사이의 정규화 주파수를 통과하는 16차 필터를 설계합니다. 저지대역 감쇠량을 60dB로 지정합니다. 원래 신호에 필터를 적용합니다.
dx = designfilt('bandpassiir','FilterOrder',16, ... 'StopbandFrequency1',0.2,'StopbandFrequency2',0.4, ... 'StopbandAttenuation',60); x = filter(dx,r);
두 번째 시퀀스를 생성하기 위해 0.6π rad/sample에서 0.8π rad/sample 사이의 정규화 주파수를 저지하는 16차 필터를 설계하십시오. 통과대역 리플을 0.1dB로 지정합니다. 원래 신호에 필터를 적용합니다.
dy = designfilt('bandstopiir','FilterOrder',16, ... 'PassbandFrequency1',0.6,'PassbandFrequency2',0.8, ... 'PassbandRipple',0.1); y = filter(dy,r);
x와 y의 크기 제곱 코히어런스를 추정합니다. 512개 샘플을 갖는 해밍 윈도우를 사용합니다. 인접 세그먼트 간에 500개 샘플이 중첩되도록 지정하고 2048개 DFT 점을 지정합니다.
[cxy,fc] = mscohere(x,y,hamming(512),500,2048);
코히어런스 함수를 플로팅하고 필터의 주파수 응답을 겹쳐 놓습니다.
[qx,f] = freqz(dx); qy = freqz(dy); plot(fc/pi,cxy) hold on plot(f/pi,abs(qx),f/pi,abs(qy)) hold off
2채널 랜덤 신호 x를 생성합니다. 두 채널에 저역통과 필터를 적용하고 이를 모두 더하여 또 다른 신호 y를 생성합니다. 사각 윈도우를 사용하여 설계된 0.3π의 차단 주파수를 갖는 30차 FIR 필터를 지정합니다.
h = fir1(30,0.3,rectwin(31)); x = randn(16384,2); y = sum(filter(h,1,x),2);
x와 y의 다중 코히어런스 추정값을 계산합니다. 1024개 샘플을 갖는 핸 윈도우를 사용하여 신호에 윈도우를 적용합니다. 인접 세그먼트 간에 512개 샘플이 중첩되도록 지정하고 1024개 DFT 점을 지정합니다. 추정값을 플로팅합니다.
noverlap = 512;
nfft = 1024;
mscohere(x,y,hann(nfft),noverlap,nfft,'mimo')
코히어런스 추정값을 필터에 대한 주파수 응답과 비교합니다. 코히어런스가 뚝 떨어지는 지점은 주파수 응답의 영점에 해당합니다.
[H,f] = freqz(h); hold on yyaxis right plot(f/pi,20*log10(abs(H))) hold off
x와 y의 일반 크기 제곱 코히어런스 추정값을 계산하고 플로팅합니다. 이 추정값은 어떠한 채널에 대해서도 1에 도달하지 않습니다.
figure mscohere(x,y,hann(nfft),noverlap,nfft)
2초 동안 각각 1kHz로 샘플링된 다중채널 신호 두 개를 생성합니다. 첫 번째 신호인 입력값은 각각 120Hz, 360Hz, 480Hz의 주파수를 갖는 세 개의 정현파로 구성됩니다. 두 번째 신호인 출력값은 각각 120Hz와 360Hz의 주파수를 갖는 두 개의 정현파로 구성됩니다. 이 중 한 정현파는 첫 번째 신호보다 π/2만큼 지연된 상태입니다. 다른 정현파는 지연이 π/4입니다. 두 신호 모두 백색 가우스 잡음에 묻혀 있습니다.
fs = 1000; f = 120; t = (0:1/fs:2-1/fs)'; inpt = sin(2*pi*f*[1 3 4].*t); inpt = inpt+randn(size(inpt)); oupt = sin(2*pi*f*[1 3].*t-[pi/2 pi/4]); oupt = oupt+randn(size(oupt));
모든 입력 신호와 각 출력 채널 간 상관도를 추정합니다. 길이가 100인 해밍 윈도우를 사용하여 데이터에 윈도우를 적용합니다. mscohere는 출력 채널당 하나씩 코히어런스 함수를 반환합니다. 코히어런스 함수는 입력값과 출력값이 공유하는 주파수에서 최댓값에 도달합니다.
[Cxy,f] = mscohere(inpt,oupt,hamming(100),[],[],fs,'mimo'); for k = 1:size(oupt,2) subplot(size(oupt,2),1,k) plot(f,Cxy(:,k)) title(['Output ' int2str(k) ', All Inputs']) end
입력 신호와 출력 신호를 서로 바꾸고 다중 코히어런스 함수를 계산합니다. 동일한 해밍 윈도우를 사용합니다. 480Hz에서는 입력값과 출력값 간에 상관관계가 없습니다. 따라서 세 번째 상관 함수에는 피크가 없습니다.
[Cxy,f] = mscohere(oupt,inpt,hamming(100),[],[],fs,'mimo'); for k = 1:size(inpt,2) subplot(size(inpt,2),1,k) plot(f,Cxy(:,k)) title(['Input ' int2str(k) ', All Outputs']) end
mscohere의 플로팅 기능을 사용하여 계산을 반복합니다.
clf
mscohere(oupt,inpt,hamming(100),[],[],fs,'mimo')
두 번째 신호의 일반 코히어런스 함수와 첫 번째 신호의 처음 두 채널을 계산합니다. 피크가 아닌 값은 다중 코히어런스 함수와 다릅니다.
[Cxy,f] = mscohere(oupt,inpt(:,[1 2]),hamming(100),[],[],fs); plot(f,Cxy)
최대 코히어런스를 갖는 점에서 상호 스펙트럼 각도를 계산하여 위상차를 구합니다.
Pxy = cpsd(oupt,inpt(:,[1 2]),hamming(100),[],[],fs); [~,mxx] = max(Cxy); for k = 1:2 fprintf('Phase lag %d = %5.2f*pi\n',k,angle(Pxy(mxx(k),k))/pi) end
Phase lag 1 = -0.51*pi Phase lag 2 = -0.22*pi
1kHz에서 각각 1초 동안 샘플링된 두 개의 정현파 신호를 생성합니다. 각 정현파는 250Hz의 주파수를 갖습니다. 이 중 한 신호가 다른 신호의 위상을 π/3라디안만큼 지연시킵니다. 두 신호 모두에 단위 분산의 백색 가우스 잡음을 포함시킵니다.
fs = 1000; f = 250; t = 0:1/fs:1-1/fs; um = sin(2*pi*f*t)+rand(size(t)); un = sin(2*pi*f*t-pi/3)+rand(size(t));
mscohere를 사용하여 신호의 크기 제곱 코히어런스를 계산하고 플로팅합니다.
mscohere(um,un,[],[],[],fs)
플롯 제목, x축의 레이블, y축의 제한을 수정합니다.
title("Magnitude-Squared Coherence") xlabel("f (Hz)") ylim([0 1.1])
gca를 사용하여 현재 좌표축에 대한 핸들을 구합니다. 눈금 표시의 위치를 변경합니다. y축의 레이블을 제거합니다.
ax = gca; ax.XTick = 0:250:500; ax.YTick = 0:0.25:1; ax.YLabel.String = [];
핸들에 대한 Children 속성을 호출하여 플로팅된 선의 색과 너비를 변경합니다.
ln = ax.Children; ln.Color = [0.8 0 0]; ln.LineWidth = 1.5;
또는, set과 get을 사용하여 line 속성을 수정할 수 있습니다.
set(get(gca,"Children"),Color=[0 0.4 0],LineStyle="--",LineWidth=1)
입력 인수
출력 인수
세부 정보
참고 문헌
[1] Gómez González, A., J. Rodríguez, X. Sagartzazu, A. Schumacher, and I. Isasa. “Multiple Coherence Method in Time Domain for the Analysis of the Transmission Paths of Noise and Vibrations with Non-Stationary Signals.” Proceedings of the 2010 International Conference of Noise and Vibration Engineering, ISMA2010-USD2010. pp. 3927–3941.
[2] Kay, Steven M. Modern Spectral Estimation. Englewood Cliffs, NJ: Prentice-Hall, 1988.
[3] Rabiner, Lawrence R., and Bernard Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975.
[4] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.
[5] Welch, Peter D. “The Use of 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, 1967, pp. 70–73.
