このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
信号の類似度の測定
この例では、信号の類似度の測定方法を示します。これは次のような問題の解決に役立ちます。長さやサンプル レートが異なる信号を比較するにはどうするのか。測定データにあるのが信号かノイズにすぎないかを確かめるにはどうするのか。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])
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
とテンプレート 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])
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")
関数 finddelay
を使用して、2 つの信号間の遅延を検出することもできます。
t21 = finddelay(s1,s2)
t21 = -350
t31 = finddelay(s1,s3)
t31 = 150
t21
は s2
が s1
から 350 サンプル遅れていることを示し、t31
は s3
が s1
より 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")
信号の周波数成分の比較
パワー スペクトルは各周波数に含まれるパワーを表示します。スペクトル コヒーレンスは信号間の周波数領域の相関を示します。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)")
関数 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 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
小さなピークは週 7 回、主要なピークは週 1 回のサイクルで現れています。7 日単位の日程表により温度制御されているビルのデータであるということを考えれば、これは理にかなっています。最初の 7 日周期は、ビルの温度に週単位の周期性があり、温度は週末に低くなり、平日は元に戻るということを示しています。1 日周期の動作は、日単位の周期性を、夜間は温度が低く日中に高くなる動作もまたあることを示しています。
参考
alignsignals
| cpsd
| finddelay
| findpeaks
| mscohere
| xcov
| xcorr