ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

信号の類似度の測定

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

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

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

% Load data
load relatedsig.mat

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 (secs)')
linkaxes(ax(1:3),'x')
axis([0 1.61 -4 4])

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 sampling rate by rational factor
T2 = resample(T2,P2,Q2);        % Change sampling 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(secs)') 
axis(ax(1:2),[-1.5 1.5 -700 700 ])

1 番目のサブプロットは信号とテンプレート 1 との相関性が低いことを示しますが、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')

関数 finddelay を使用して、2 つの信号間の遅延を検出することもできます。

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

t21 は s2 が s1 よりも 350 サンプル遅れており、t31 は s3 が s1 よりも 150 サンプル進んでいることを示しています。これで、この情報を使用し、信号の時間をシフトして 3 つの信号を整列させることができます。最初の信号を遅延させて 2 つの信号を整列させる関数 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')

信号の周波数成分の比較

パワー スペクトルは各周波数に含まれるパワーを表示します。スペクトル コヒーレンスは信号間の周波数領域の相関を示します。0 に向かうコヒーレンスの値は対応する周波数成分が非相関であることを、また、1 に向かうコヒーレンスの値は対応する周波数成分が相関であることを示します。2 つの信号とそれぞれのパワー スペクトルを見てみます。

Fs = FsSig;         % Sampling 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 (secs)')
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)')

関数 mscohere は 2 つの信号間のスペクトル コヒーレンスを計算します。sig1 と sig2 の 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])

35 Hz 成分の間の位相遅れは -90°に近く、165 Hz 成分の間の位相遅れは -60°に近くなっています。

信号内の周期性の検出

冬期のオフィス ビルにおける温度測定の結果を見てみましょう。測定は約 16.5 週間、30 分間隔で行われました。

% Load Temperature Data
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

華氏 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')

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

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')

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

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