gccphat

一般化相互相関

構文

tau = gccphat(sig,refsig)

tau = gccphat(sig,refsig,fs)

[tau,R,lag] = gccphat(___)

[___] = gccphat(sig)

[___] = gccphat(sig,fs)

説明

tau = gccphat(sig,refsig) は、信号 sig と基準信号 refsig との間の時間遅延 tau を計算します。sig と refsig は、どちらも複数のチャネルを含めることができます。関数は、信号と基準信号が単一のソースから到達すると仮定します。遅延を推定するために、gccphat は sig と refsig との相互相関におけるピーク位置を検出します。相互相関は、一般化相互相関位相変換 (GCC-PHAT) アルゴリズムを使用して計算されます。時間遅延は、既定のサンプリング周波数 1 ヘルツに対応するサンプル間隔の倍数です。

tau = gccphat(sig,refsig,fs) は、信号のサンプリング周波数を指定します。時間遅延は、サンプリング周波数に対応するサンプル間隔の倍数です。すべての入力信号のサンプル レートは同じでなければなりません。

[tau,R,lag] = gccphat(___) は、これまでの構文のいずれかの引数を使用して、さらに相互相関値と相関タイム ラグを返します。ラグは、サンプリング間隔の倍数です。相互相関チャネル数は、sig のチャネル数に等しくなります。

[___] = gccphat(sig) または [___] = gccphat(sig,fs) は、sig 内のすべてのチャネル ペア間の推定遅延と推定相互相関を返します。sig の列数が M である場合、結果の tau と R の列数は M² になります。これらの構文では、基準信号入力が使用されません。tau と R の最初の M 列には、第 1 チャネルを基準として使用する遅延と相互相関が含まれます。2 番目の M 列には、第 2 チャネルを基準として使用する遅延と相互相関が含まれます。これ以降も同様です。

例

2 つの信号と基準信号との間の相互相関

ゴングの音声信号を読み込みます。まず、ゴングの信号を基準信号として使用します。次に、信号を 2 回複製し、5 秒および 25 秒の時間遅延を導入します。サンプル レートは、既定値の 1 ヘルツのままにします。gccphat を使用して、遅延信号と基準信号との間の時間遅延を推定します。

load gong;
refsig = y;
delay1 = 5;
delay2 = 25;
sig1 = delayseq(refsig,delay1);
sig2 = delayseq(refsig,delay2);
tau_est = gccphat([sig1,sig2],refsig)

tau_est = 1×2

     5    25

信号と基準信号との間の相互相関

ゴングの音声信号を読み込みます。ゴングの信号を基準信号として使用します。次に、信号を複製し、5 ミリ秒の時間遅延を導入します。8192 Hz のサンプル レートを使用します。gccphat を使用して、遅延信号と基準信号との間の時間遅延を推定します。

load gong;
delay = 0.005;
refsig = y;
sig = delayseq(refsig,delay,Fs);
tau_est = gccphat(sig,refsig,Fs)

tau_est = 
0.0050

3 つの信号と基準信号との相互相関のプロット

サンプル レートが 8192 ヘルツの楽音信号を読み込みます。次に、信号を 3 回複製し、信号間に時間遅延を導入します。遅延信号と基準信号との間の時間遅延を推定します。相関値をプロットします。

load handel;
dt = 1/Fs;
refsig = y;

信号の遅延バージョンを 3 つ作成します。

delay1 = -5.2*dt;
delay2 = 10.3*dt;
delay3 = 7*dt;
sig1 = delayseq(refsig,delay1,Fs);
sig2 = delayseq(refsig,delay2,Fs);
sig3 = delayseq(refsig,delay3,Fs);

遅延信号と基準信号との相互相関を計算します。

[tau_est,R,lags] = gccphat([sig1,sig2,sig3],refsig,Fs);

gccphat 関数は、最も近いサンプル間隔までの遅延を推定します。

disp(tau_est*Fs)

    -5    10     7

相関関数をプロットします。

plot(1000*lags,real(R(:,1)))
xlabel('Lag Times (ms)')
ylabel('Cross-correlation')
axis([-5,5,-.4,1.1])
hold on
plot(1000*lags,real(R(:,2)))
plot(1000*lags,real(R(:,3)))
hold off

複数の信号の相互相関のプロット

サンプル レートが 8192 ヘルツの楽音信号を読み込みます。次に、信号を 2 回複製し、2 つの信号と基準信号との間に時間遅延を導入します。すべての信号ペア間の時間遅延を推定し、相互相関関数をプロットします。

load handel;
dt = 1/Fs;
refsig = y;

信号の遅延バージョンを 3 つ作成します。

delay1 = -5.7*dt;
delay2 = 10.2*dt;
sig1 = delayseq(refsig,delay1,Fs);
sig2 = delayseq(refsig,delay2,Fs);

すべての信号と他の信号との相互相関を計算します。

[tau_est,R,lags] = gccphat([refsig,sig1,sig2],Fs);

時間遅延をサンプル間隔単位で表示します。このアルゴリズムは、最も近いサンプル間隔に量子化された時間遅延を推定します。3 つの信号の相互相関では、考えられるすべての信号ペアに対して 1 つずつ、合計 9 個の時間遅延が生成されます。

disp(tau_est*Fs)

     0    -6    10     6     0    16   -10   -16     0

信号とその信号自体との相互相関では、ラグがゼロになります。

相関関数をプロットします。

for n=1:9
    plot(1000*lags,real(R(:,n)))
    if n==1
        hold on
        xlabel('Lag Times (ms)')
        ylabel('Correlation')
        axis([-5,5,-.4,1.1])
    end
end
hold off

入力引数

sig — センサー信号

N 行 1 列の複素数値の列ベクトル | N 行 M 列の複素数値の行列

センサー信号。N 行 1 列の列ベクトルまたは N 行 M 列の行列として指定します。N は時間サンプル数、M はチャネル数です。sig が行列である場合、各列は異なるチャネルです。

例: [0,1,2,3,2,1,0]

データ型: single | double
複素数のサポート: あり

refsig — 基準センサー信号

N 行 1 列の複素数値の列ベクトル | N 行 M 列の複素数値の行列

基準信号。N 行 1 列の複素数値の列ベクトル、または N 行 M 列の複素数値の行列として指定します。refsig が列ベクトルである場合、sig 内のすべてのチャネルは、相互相関を計算する際に refsig を基準信号として使用します。

refsig が行列である場合、refsig のサイズは sig のサイズと一致しなければなりません。gccphat 関数は、sig 内と refsig 内の対応するチャネルの間の相互相関を計算します。これらの信号は、異なるソースから取得できます。

例: [1,2,3,2,1,0,0]

データ型: single | double
複素数のサポート: あり

fs — 信号のサンプル レート

1 (既定値) | 正の実数値のスカラー

信号のサンプル レート。正の実数値のスカラーとして指定します。すべての信号のサンプル レートは同じでなければなりません。サンプル レートの単位はヘルツです。

例: 8000

データ型: single | double
複素数のサポート: あり

出力引数

tau — 時間遅延

1 行 K 列の実数値の行ベクトル

時間遅延。1 行 K 列の実数値の行ベクトルとして返されます。K の値は、入力引数の構文によって異なります。

• 基準信号 refsig を使用する場合、K の値は sig の列次元 M と等しくなります。tau の各エントリは、sig 内と refsig 内の対応する信号のペアに対する推定遅延を指定します。

• 基準信号を使用しない場合、K の値は sig の列次元の 2 乗 M² と等しくなります。tau の各エントリは、sig 内の対応する信号のペアに対する推定遅延を指定します。

単位は秒です。

R — 信号間の相互相関

(2N-1) 行 K 列の複素数値の行列

異なるセンサーの信号間の相互相関。(2N-1) 行 K 列の複素数値の行列として返されます。

• 基準信号 refsig を使用する場合、K の値は sig の列次元 M と等しくなります。各列は、sig 内と refsig 内の対応する信号のペア間の相互相関です。

• 基準信号を使用しない場合、K の値は sig の列次元の 2 乗 M² と等しくなります。各列は、sig 内の対応する信号のペア間の相互相関です。

lag — 相互相関ラグ時間

(2N-1) の実数値の列ベクトル

相関ラグ時間。(2N-1) の実数値の列ベクトルとして返されます。lag の各行には、R の対応する行のラグ時間が含まれます。ラグ値は、サンプリング間隔の倍数に制限されます。ラグの単位は秒です。

詳細

一般化相互相関

一般化相互相関を使用して、2 つの異なるセンサーにおける信号の到着時間差 (TDOA) を推定できます。

TDOA を推定するには、あるセンサーで受信した信号と 2 番目のセンサーで受信した信号との間の相互相関関数におけるピーク位置を検出します。相互相関関数のピークは、時間遅延に対応します。phased.GCCEstimator System object と phased.TDOAEstimator System object は、どちらも一般化相互相関位相変換 (GCC-PHAT) アルゴリズムを使用して相関ピークを計算します。ピーク位置は、サンプル間隔の倍数に制限されます。サンプル間隔は、信号のサンプル レートの逆数です。信号のサンプル レートは、両方の System object に共通する SampleRate プロパティで指定できます。

2 つの信号 x₁(t) と x₂(t) との間の時間遅延を計算するには、まずノイズが存在する状態での時間領域表現から始めます。

$$\begin{array}{l}{x}_1(t)=s(t)+n_1(t)\\ x_2(t)=\alpha s(t-D)+n_2(t)\end{array}$$

ここで、n₁(t) と n₂(t) は各信号に加えられるノイズです。D は時間遅延です。

周波数領域では、信号が次のように表されます。

$$\begin{array}{l}{X}_1(\omega )=S(\omega )+N_1(\omega )\\ X_2(\omega )=\alpha S(\omega )e^{-j\omega D}+N_2(\omega )\end{array}$$

時間領域と周波数領域での 2 つの信号間の相互相関関数を定義します。

$$\begin{array}{c}{R}_{x_1,x_2}(\tau )=\text{E}[x_1(t)x_2(t-\tau )]\\ R_{x_1,x_2}(\omega )=X_1(\omega )X_2^{*}(\omega )\\ =\alpha {\left|S(\omega )\right|}^2{e}^{j\omega \tau }\end{array}$$

時間遅延を特定するには、相互相関関数のピーク τ_max の位置を特定します。S/N 比 (SNR) が大きい場合、相関ピーク τ_max は実際の時間遅延 D に対応します。相関関数のピークが鋭くなるほど、パフォーマンスは向上します。入力信号を白色化する重み付け関数を使用することで、相互相関ピークをシャープ化できます。この手法は、"一般化相互相関" (GCC) と呼ばれます。ある特定の重み付け関数では、信号のスペクトル密度をスペクトル振幅で正規化します。これにより、"一般化相互相関位相変換" ("GCC-PHAT") 法が導かれます。ラグ τ における一般化相互相関位相変換 (GCC-PHAT) 関数は、次のように定義されます。

$$\begin{array}{c}{\tilde{R}}_{x_1{x}_2}(\tau )=\frac{1}{2\pi }{\displaystyle {\int }_{-\pi }^{\pi }\frac{R_{x_1{x}_2}(\omega )}{\left|R_{x_1{x}_2}(\omega )\right|}}{e}^{i\omega \tau }d\omega \\ =\frac{1}{2\pi }{\displaystyle {\int }_{-\pi }^{\pi }{e}^{i\omega (D+\tau )}}d\omega \\ =\delta (D-\tau )\end{array}$$

ここで、δ はディラックデルタ関数です。この関数は、到着時間差 (TDOA) である D=τ のときに非ゼロになります。通常、ラグはサンプル単位で計算されるため、秒単位に変換するにはサンプル レートで除算します。時間遅延は、2 つの信号間の相互相関を最大化するタイム ラグを求めることで推定できます。

TDOA から、2 つのセンサーを結ぶ直線に対する平面波のブロードサイド到来角を推定できます。距離 L だけ離れた 2 つのセンサーでは、ブロードサイド到来角ブロードサイド角度と時間ラグとの間に次の関係があります。

$$\sin \beta =\frac{c\tau }{L}$$

ここで、c は媒質中での伝搬速度です。

センサー ペアを 1 組だけ使用する場合、推定できるのはブロードサイド到来角のみです。しかし、URA 内など、同一線上にない複数のセンサー ペアを使用すると、最小二乗推定を使用して、平面波の到来方位角および到来仰角を推定できます。N 個のセンサーがある場合、k^th センサーに到達する信号の、j^th センサーに対する遅延時間 τ_kj は、次のように記述できます。

$$\begin{array}{l}c{\tau }_{kj}=-\left({\overrightarrow{x}}_k-{\overrightarrow{x}}_j\right)\cdot \overrightarrow{u}\\ \overrightarrow{u}=\cos \alpha \sin \theta \widehat{i}+\sin \alpha \sin \theta \widehat{j}+\cos \theta \widehat{k}\end{array}$$

ここで、u は平面波の単位伝搬ベクトルです。角度 α と θ は、伝搬ベクトルの方位角と仰角です。すべての角度とベクトルは、ローカル座標軸に対して定義されます。最小二乗法を使用して最初の方程式を解くことで、単位伝搬ベクトルの 3 つの成分を求めることができます。次に、2 番目の方程式を解くことで、方位角と仰角を求めることができます。

参照

[1] Knapp, C. H. and G.C. Carter, “The Generalized Correlation Method for Estimation of Time Delay.” IEEE Transactions on Acoustics, Speech and Signal Processing. Vol. ASSP-24, No. 4, Aug 1976.

[2] G. C. Carter, “Coherence and Time Delay Estimation.” Proceedings of the IEEE. Vol. 75, No. 2, Feb 1987.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

可変サイズの入力はサポートされていません。

バージョン履歴

R2015b で導入

参考

phased.GCCEstimator