Signal Processing Toolbox

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

信号の類似度の測定

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

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

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

% Load data
load relatedsig.mat;

figure
ax(1) = subplot(311);
plot((0:numel(T1)-1)/Fs1,T1,'k'); ylabel('Template 1'); grid on
ax(2) = subplot(312);
plot((0:numel(T2)-1)/Fs2,T2,'r'); ylabel('Template 2'); grid on
ax(3) = subplot(313);
plot((0:numel(S)-1)/Fs,S,'b'); 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 =

        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(211);
plot(lag1/Fs,C1,'k'); ylabel('Amplitude'); grid on
title('Cross-correlation between Template 1 and Signal')
ax(2) = subplot(212);
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));
timeDiff = lag2(I)/Fs
timeDiff =

    0.0609

相互相関のピークは、テンプレート T2 の 61 ms 後からその信号が存在することを示唆しています。

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

橋の両側で車によって生じる振動を記録している異なるセンサーからデータを収集している状況を考えてみましょう。信号を解析する場合、それらを整列させる必要があります。3 つのセンサーは同じサンプルレートで、同じ事象によって生じる信号を測定しているとします。

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

s1、s2 間、および s1、s3 間の相互相関の最大値は時間の進みと遅れを示しています。

[C21,lag1] = xcorr(s2,s1);
[C31,lag2] = xcorr(s3,s1);

figure
subplot(211); plot(lag1,C21/max(C21)); ylabel('C21');grid on
title('Cross-Correlations')
subplot(212); plot(lag2,C31/max(C31)); ylabel('C31');grid on
xlabel('Samples')

[~,I1] = max(abs(C21));     % Find the index of the highest peak
[~,I2] = max(abs(C31));     % Find the index of the highest peak
t21 = lag1(I1)              % Time difference between the signals s2,s1
t31 = lag2(I2)              % Time difference between the signals s3,s1
t21 =

  -350


t31 =

   150

t21 は s2 が s1 よりも 350 サンプル遅れており、t31 は s3 が s1 よりも 150 サンプル進んでいることを示しています。これで、この情報を使用して 3 つの信号を配置することができます。

s2 = [zeros(abs(t21),1);s2];
s3 = s3(t31:end);

figure
ax(1) = subplot(311); plot(s1); grid on; title('s1'); axis tight
ax(2) = subplot(312); plot(s2); grid on; title('s2'); axis tight
ax(3) = subplot(313); 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(221); plot(t,sig1,'k'); ylabel('s1');grid on
title('Time Series')
subplot(223); plot(t,sig2); ylabel('s2');grid on
xlabel('Time (secs)')
subplot(222); plot(f1,P1,'k'); ylabel('P1'); grid on; axis tight
title('Power Spectrum')
subplot(224); 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(211);
plot(f,Cxy); title('Coherence Estimate');grid on;
set(gca,'xtick',f(locs),'ytick',.75);
axis([0 200 0 1])
subplot(212);
plot(f,phase); title('Cross Spectrum Phase (deg)');grid on;
set(gca,'xtick',f(locs),'ytick',round(phase(locs)));
xlabel('Frequency (Hz)');
axis([0 200 -180 180])

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

信号内の周期性の検出

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

load officetemp.mat  % Load Temperature Data
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
set(gca,'Xlim',[-15 15])
xlabel('Time (days)')
title('Auto-covariance')

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

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

subplot(211); plot(cycle1); ylabel('Days'); grid on
title('Dominant peak distance')
subplot(212); plot(cycle2,'r'); ylabel('Days'); grid on
title('Minor peak distance')

mean(cycle1)
mean(cycle2)
ans =

     7


ans =

    1.0000

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

製品評価版の入手または製品の購入