Main Content

신호 유사성 측정하기

이 예제에서는 신호 유사성을 측정하는 방법을 보여줍니다. 이는 다음과 같은 문제의 해결에 도움이 됩니다. 길이가 서로 다르거나 샘플 레이트가 서로 다른 신호를 비교하려면 어떻게 해야 합니까? 측정값에 신호가 있는지 또는 잡음만 있는지 확인하려면 어떻게 해야 합니까? 두 신호가 관련이 있습니까? 두 신호 간의 지연을 측정하려면 어떻게 해야 합니까? 또한 이를 정렬하려면 어떻게 해야 합니까? 두 신호의 주파수 성분을 비교하려면 어떻게 해야 합니까? 유사성은 하나의 신호의 여러 부분에서 검출될 수도 있는 것으로, 그 신호가 주기성을 갖는지 여부를 결정하기도 합니다.

샘플 레이트가 서로 다른 신호 비교하기

오디오 신호의 데이터베이스와 패턴 일치 응용 프로그램을 가지고 재생 중인 곡명을 알아내는 상황이라고 가정하겠습니다. 일반적으로 데이터는 메모리를 덜 차지하도록 낮은 샘플 레이트로 저장됩니다.

load relatedsig

figure
ax(1) = subplot(3,1,1);
plot((0:numel(T1)-1)/Fs1,T1,"k")
ylabel("Template 1")
grid on
ax(2) = subplot(3,1,2); 
plot((0:numel(T2)-1)/Fs2,T2,"r")
ylabel("Template 2")
grid on
ax(3) = subplot(3,1,3); 
plot((0:numel(S)-1)/Fs,S)
ylabel("Signal")
grid on
xlabel("Time (s)")
linkaxes(ax(1:3),"x")
axis([0 1.61 -4 4])

Figure contains 3 axes objects. Axes object 1 with ylabel Template 1 contains an object of type line. Axes object 2 with ylabel Template 2 contains an object of type line. Axes object 3 with xlabel Time (s), ylabel Signal contains an object of type line.

첫 번째 서브플롯과 두 번째 서브플롯에서는 데이터베이스의 템플릿 신호를 보여줍니다. 세 번째 서브플롯은 데이터베이스에서 검색하려는 신호를 보여줍니다. 시계열만 봐도 신호가 두 템플릿 중 어떤 템플릿과도 일치하지 않습니다. 더 정밀하게 검사하면 신호가 실제로 길이와 샘플 레이트에서 각각 다른 것을 알 수 있습니다.

[Fs1 Fs2 Fs]
ans = 1×3

        4096        4096        8192

길이가 다르므로 두 신호 간의 거리를 계산할 수 없지만, 이 문제는 신호의 공통 부분을 추출함으로써 손쉽게 해결할 수 있습니다. 사실, 길이를 항상 같게 만들 필요는 없습니다. 상호상관은 길이가 다른 신호 간에 수행할 수 있습니다. 단, 이 경우 샘플 레이트는 반드시 동일해야 합니다. 이를 수행하기 위한 가장 안전한 방법은 더 낮은 샘플 레이트로 신호를 리샘플링하는 것입니다. resample 함수는 리샘플링 과정 동안 안티에일리어싱(저역통과) FIR 필터를 신호에 적용합니다.

[P1,Q1] = rat(Fs/Fs1);          % Rational fraction approximation
[P2,Q2] = rat(Fs/Fs2);          % Rational fraction approximation
T1 = resample(T1,P1,Q1);        % Change sample rate by rational factor
T2 = resample(T2,P2,Q2);        % Change sample rate by rational factor

측정값에서 신호 찾기

이제 xcorr 함수를 사용하여 신호 S와 템플릿 T1, T2의 상호상관을 분석하여 일치하는 부분이 있는지 확인할 수 있습니다.

[C1,lag1] = xcorr(T1,S);        
[C2,lag2] = xcorr(T2,S);        

figure
ax(1) = subplot(2,1,1); 
plot(lag1/Fs,C1,"k")
ylabel("Amplitude")
grid on
title("Cross-Correlation Between Template 1 and Signal")
ax(2) = subplot(2,1,2); 
plot(lag2/Fs,C2,"r")
ylabel("Amplitude") 
grid on
title("Cross-Correlation Between Template 2 and Signal")
xlabel("Time(s)") 
axis(ax(1:2),[-1.5 1.5 -700 700])

Figure contains 2 axes objects. Axes object 1 with title Cross-Correlation Between Template 1 and Signal, ylabel Amplitude contains an object of type line. Axes object 2 with title Cross-Correlation Between Template 2 and Signal, xlabel Time(s), ylabel Amplitude contains an object of type line.

첫 번째 서브플롯은 신호 S와 템플릿 T1의 상관성이 낮음을 나타내는 반면, 두 번째 서브플롯의 높은 피크는 신호가 두 번째 템플릿에 존재함을 나타냅니다.

[~,I] = max(abs(C2));
SampleDiff = lag2(I)
SampleDiff = 499
timeDiff = SampleDiff/Fs
timeDiff = 0.0609

상호상관의 피크는 신호가 템플릿 T2에 존재하고 61ms 후 시작됨을 의미합니다. 다시 말해서, 템플릿 T2SampleDiff에 표시된 대로 499개 샘플만큼 신호 S보다 앞섭니다. 이 정보를 활용하여 신호를 정렬할 수 있습니다.

신호 간의 지연을 측정하고 신호 정렬하기

다리 양쪽에서 자동차로 인해 발생하는 진동을 여러 다른 센서에서 녹음한 데이터를 수집하는 상황을 가정해 보겠습니다. 신호를 분석할 때 신호를 정렬해야 할 수 있습니다. 동일한 샘플 레이트로 작동하고 동일한 이벤트로 인해 발생하는 신호를 측정하는 3개의 센서가 있다고 가정합니다.

figure
ax(1) = subplot(3,1,1);
plot(s1)
ylabel("s1")
grid on
ax(2) = subplot(3,1,2); 
plot(s2,"k")
ylabel("s2")
grid on
ax(3) = subplot(3,1,3); 
plot(s3,"r")
ylabel("s3")
grid on
xlabel("Samples")
linkaxes(ax,"xy")

Figure contains 3 axes objects. Axes object 1 with ylabel s1 contains an object of type line. Axes object 2 with ylabel s2 contains an object of type line. Axes object 3 with xlabel Samples, ylabel s3 contains an object of type line.

finddelay 함수를 사용하여 두 신호 간의 지연을 검출할 수도 있습니다.

t21 = finddelay(s1,s2)
t21 = -350
t31 = finddelay(s1,s3)
t31 = 150

t21s2s1보다 350개 샘플만큼 지연됨을 나타내고, t31s3s1보다 150개 샘플만큼 앞서 있음을 나타냅니다. 이제, 이 정보를 활용하여 시간을 기준으로 신호를 이동하여 3개 신호를 정렬할 수 있습니다. 또한 alignsignals 함수를 사용하여 더 빠른 신호를 지연시켜 신호를 정렬할 수 있습니다.

s1 = alignsignals(s1,s3);
s2 = alignsignals(s2,s3);

figure
ax(1) = subplot(3,1,1);
plot(s1)
grid on 
title("s1")
axis tight
ax(2) = subplot(3,1,2);
plot(s2)
grid on 
title("s2")
axis tight
ax(3) = subplot(3,1,3); 
plot(s3)
grid on 
title("s3")
axis tight
linkaxes(ax,"xy")

Figure contains 3 axes objects. Axes object 1 with title s1 contains an object of type line. Axes object 2 with title s2 contains an object of type line. Axes object 3 with title s3 contains an object of type line.

신호의 주파수 성분 비교하기

파워 스펙트럼은 각 주파수에 존재하는 전력을 표시합니다. 스펙트럼 코히어런스는 신호 간의 주파수 영역 상관관계를 나타냅니다. 0 방향의 코히어런스 값은 대응하는 주파수 성분이 상관 관계가 없음을 나타내는 반면, 1 방향의 값은 대응하는 주파수 성분이 상관 관계가 있음을 나타냅니다. 두 신호와 각각 해당하는 파워 스펙트럼이 있다고 가정하겠습니다.

Fs = FsSig;         % Sample Rate

[P1,f1] = periodogram(sig1,[],[],Fs,"power");
[P2,f2] = periodogram(sig2,[],[],Fs,"power");

figure
t = (0:numel(sig1)-1)/Fs;
subplot(2,2,1)
plot(t,sig1,"k")
ylabel("s1")
grid on
title("Time Series")
subplot(2,2,3)
plot(t,sig2)
ylabel("s2")
grid on
xlabel("Time (s)")
subplot(2,2,2)
plot(f1,P1,"k")
ylabel("P1")
grid on
axis tight
title("Power Spectrum")
subplot(2,2,4)
plot(f2,P2)
ylabel("P2")
grid on
axis tight
xlabel("Frequency (Hz)")

Figure contains 4 axes objects. Axes object 1 with title Time Series, ylabel s1 contains an object of type line. Axes object 2 with xlabel Time (s), ylabel s2 contains an object of type line. Axes object 3 with title Power Spectrum, ylabel P1 contains an object of type line. Axes object 4 with xlabel Frequency (Hz), ylabel P2 contains an object of type line.

mscohere 함수는 두 신호 간의 스펙트럼 코히어런스를 계산합니다. 이 함수의 결과를 통해 sig1sig2가 약 35Hz와 165Hz의, 2개의 상관된 성분을 가진다는 것을 확인할 수 있습니다. 스펙트럼 코히어런스가 높은 주파수에서, 상관된 성분 간의 상대적인 위상은 상호 스펙트럼 위상으로 추정될 수 있습니다.

[Cxy,f] = mscohere(sig1,sig2,[],[],[],Fs);
Pxy = cpsd(sig1,sig2,[],[],[],Fs);
phase = -angle(Pxy)/pi*180;
[pks,locs] = findpeaks(Cxy,MinPeakHeight=0.75);

figure
subplot(2,1,1)
plot(f,Cxy)
title("Coherence Estimate")
grid on
hgca = gca;
hgca.XTick = f(locs);
hgca.YTick = 0.75;
axis([0 200 0 1])
subplot(2,1,2)
plot(f,phase)
title("Cross-Spectrum Phase (deg)")
grid on
hgca = gca;
hgca.XTick = f(locs); 
hgca.YTick = round(phase(locs));
xlabel("Frequency (Hz)")
axis([0 200 -180 180])

Figure contains 2 axes objects. Axes object 1 with title Coherence Estimate contains an object of type line. Axes object 2 with title Cross-Spectrum Phase (deg), xlabel Frequency (Hz) contains an object of type line.

35Hz 성분 간의 위상 지연은 -90도에 가까우며, 165Hz 성분 간의 위상 지연은 -60도에 가깝습니다.

신호의 주기성 찾기

겨울 동안 사무실 건물에서 측정한 일련의 온도 값이 있다고 가정하겠습니다. 온도 값은 약 16.5주 동안 30분 간격으로 측정되었습니다.

load officetemp.mat  

Fs = 1/(60*30);                          % Sample rate is 1 sample every 30 minutes
days = (0:length(temp)-1)/(Fs*60*60*24); 

figure
plot(days,temp)
title("Temperature Data")
xlabel("Time (days)")
ylabel("Temperature (Fahrenheit)")
grid on

Figure contains an axes object. The axes object with title Temperature Data, xlabel Time (days), ylabel Temperature (Fahrenheit) contains an object of type line.

화씨 70도대의 저온에서는 신호의 작은 변동을 분석하기 위해 평균을 제거해야 합니다. xcov 함수는 상호상관을 계산하기 전에 신호의 평균을 제거하고 교차공분산을 반환합니다. 교차공분산을 정확히 추정하기 위해 최대 지연을 신호의 50%로 제한합니다.

maxlags = numel(temp)*0.5;
[xc,lag] = xcov(temp,maxlags);         

[~,df] = findpeaks(xc,MinPeakDistance=5*2*24);
[~,mf] = findpeaks(xc);

figure
plot(lag/(2*24),xc,"k",...
     lag(df)/(2*24),xc(df),"kv",MarkerFaceColor="r")
grid on
xlim([-15 15])
xlabel("Time (days)")
title("Auto-Covariance")

Figure contains an axes object. The axes object with title Auto-Covariance, xlabel Time (days) contains 2 objects of type line. One or more of the lines displays its values using only markers

자기공분산(Auto-covariance)에서 주요 변동과 작은 변동을 확인해 보십시오. 주요 피크와 작은 피크는 등거리로 나타납니다. 이를 확인하기 위해 피크 사이의 거리를 계산하고 플로팅합니다.

cycle1 = diff(df)/(2*24);
cycle2 = diff(mf)/(2*24);

subplot(2,1,1)
plot(cycle1)
ylabel("Days")
grid on
title("Dominant Peak Distance")
subplot(2,1,2)
plot(cycle2,"r")
ylabel("Days")
grid on
title("Minor Peak Distance")

Figure contains 2 axes objects. Axes object 1 with title Dominant Peak Distance, ylabel Days contains an object of type line. Axes object 2 with title Minor Peak Distance, ylabel Days contains an object of type line.

mean(cycle1)
ans = 7
mean(cycle2)
ans = 1

작은 피크(Minor Peak)는 주 7회, 주요 피크(Dominant Peak)는 주 1회 나타나고 있습니다. 이는 온도가 조절되는 건물에서 한 주를 기준으로 하여 측정된 데이터를 가져온 것이란 것을 생각해보면 이해할 수 있습니다. 첫 번째 7일 주기는 건물의 온도에 주 단위의 주기성이 있어 온도가 주말 동안 낮았다가 주중에 다시 정상으로 돌아오는 것을 나타냅니다. 1일 주기 동작은 일 단위의 주기성이 있어 야간에는 온도가 낮았다가 주간에는 온도가 높아지는 것을 나타냅니다.

참고 항목

| | | | | |