Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

信号の類似度の測定

この例では、信号の類似度の測定方法を示します。これは次のような問題の解決に役立ちます。長さやサンプル レートが異なる信号を比較するにはどうするのか。測定データにあるのが信号かノイズにすぎないかを確かめるにはどうするのか。2 つの信号に関連はあるのか。2 つの信号間の遅れを測定するにはどうするのか (また、信号はどう配置するのか)。2 つの信号の周波数成分はどうやって比較するのか。類似度は、1 つの信号の異なる部分で検出されることで、その信号が周期的かどうかの判断も可能にします。

サンプル レートが異なる信号の比較

オーディオ信号データベースやパターン マッチング アプリケーションで、再生中の曲名を特定する場合を考えてみます。メモリ占有量を少なくするため、一般的にデータは低サンプル レートで保存されます。

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.

1 番目および 2 番目のサブプロットは、データベースのテンプレート信号を示します。3 番目のサブプロットはデータベースでの探索対象となる信号を示しています。時系列を見るだけでは、この信号は 2 つのテンプレートのいずれとも一致しないように見えます。さらに詳しく調べると、実際には信号の長さやサンプル レートが異なっていることがわかります。

[Fs1 Fs2 Fs]
ans = 1×3

        4096        4096        8192

信号の長さが異なると 2 つの信号の違いの計算が妨げられますが、これは信号の共通部を抽出することによって簡単に対処できます。さらには必ずしも長さを同じにする必要はありません。長さが異なる信号間の相互相関は可能ですが、必ず同じサンプル レートをもっていることが必要です。これを行う最も確実な方法は信号を低いサンプル レートでリサンプリングすることです。関数 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 とテンプレート T1T2 を相互相関させ、一致しているかどうかを判断することができます。

[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.

1 番目のサブプロットは信号 S とテンプレート T1 の相関が低いことを示しますが、2 番目のサブプロットにおける高いピークは 2 番目のテンプレートの中に信号があることを示しています。

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

相互相関のピークは、テンプレート T2 の 61 ms 後からその信号が存在することを示唆しています。つまり、SampleDiff によって示されるように、テンプレート T2 は信号 S を 499 サンプル分先行することになります。この情報は、信号の整列に使用できます。

信号間の遅延の測定および信号の配置

橋の両側で車によって生じる振動を記録している異なるセンサーからデータを収集している状況を考えてみましょう。信号を解析する場合、それらを整列させる必要があります。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 を使用して、2 つの信号間の遅延を検出することもできます。

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 に向かうコヒーレンスの値は対応する周波数成分が相関していることを示します。2 つの信号とそれぞれのパワー スペクトルを見てみます。

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 は 2 つの信号間のスペクトル コヒーレンスを計算します。sig1sig2 の 35 Hz および 165 Hz 付近に 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.

35 Hz 成分の間の位相遅れは -90°に近く、165 Hz 成分の間の位相遅れは -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

自己共分散における主要な変動と小さな変動に注目してください。主要な変動と小さな変動が等間隔に現れています。等間隔であるかどうかを確認するため、後続の各ピークの位置の差を計算してプロットします。

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

小さなピークは週 7 回、主要なピークは週 1 回のサイクルで現れています。7 日単位の日程表により温度制御されているビルのデータであるということを考えれば、これは理にかなっています。最初の 7 日周期は、ビルの温度に週単位の周期性があり、温度は週末に低くなり、平日は元に戻るということを示しています。1 日周期の動作は、日単位の周期性を、夜間は温度が低く日中に高くなる動作もまたあることを示しています。

参考

| | | | | |