ドキュメンテーション

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

xcorr

構文

r = xcorr(x,y)
r = xcorr(x)
r = xcorr(___,maxlag)
r = xcorr(___,scaleopt)
[r,lags] = xcorr(___)

説明

r = xcorr(x,y) では、2 つの離散時間シーケンス (x および y) の相互相関が返されます。相互相関は、xy のシフトされた (ラグのある) コピーとの間の相似性をそのラグの関数として測定します。x および y の長さが異なる場合、両方が同じ長さ N になるように、関数により短いほうのベクトルの最後にゼロが追加されます。

r = xcorr(x) では x の自己相関列が返されます。x が行列の場合、r は、その列に x の列のすべての組み合わせに対する自己相関列と相互相関列を含む行列です。

r = xcorr(___,maxlag) のラグ範囲は –maxlagmaxlag に制限されています。この構文では、1 つまたは 2 つの入力シーケンスを受け入れます。maxlag の既定値は N – 1 です。

r = xcorr(___,scaleopt) は、さらに相互相関または自己相関の正規化オプションを指定します。'none' (既定) 以外のオプションでは xy の長さを同じにする必要があります。

[r,lags] = xcorr(___) は、相関の計算対象になるラグをもつベクトルを返します。

すべて折りたたむ

別の場所にある 2 つのセンサーは、車が橋を渡るときに車の振動を測定します。

信号およびサンプルレート を読み込みます。時間ベクトルを作成し、信号をプロットします。センサー 2 の信号のほうが、センサー 1 の信号よりも早く到着します。

load sensorData

t1 = (0:length(s1)-1)/Fs;
t2 = (0:length(s2)-1)/Fs;

subplot(2,1,1)
plot(t1,s1)
title('s_1')

subplot(2,1,2)
plot(t2,s2)
title('s_2')
xlabel('Time (s)')

2 つの測定値の相互相関は、遅延と等しいラグで最大になります。

相互相関をプロットします。この遅延をサンプル数として秒単位で表します。

[acor,lag] = xcorr(s2,s1);

[~,I] = max(abs(acor));
lagDiff = lag(I)
lagDiff = -350
timeDiff = lagDiff/Fs
timeDiff = -0.0317
figure
plot(lag,acor)
a3 = gca;
a3.XTick = sort([-3000:1000:3000 lagDiff]);

2 つの信号を並べてもう一度プロットします。MATLAB® で使用される 1 ベースのインデックスを考慮して、ラグ差には 1 を追加します。

s1al = s1(-lagDiff+1:end);
t1al = (0:length(s1al)-1)/Fs;

subplot(2,1,1)
plot(t1al,s1al)
title('s_1, aligned')

subplot(2,1,2)
plot(t2,s2)
title('s_2')
xlabel('Time (s)')

音声の録音には壁面からの反射によるエコーが含まれることがあります。自己相関を使用したフィルター処理によりこれを除外します。

この録音では、話者は MATLAB® という言葉を話しています。データとサンプルレート を読み込みます。

load mtlb

% To hear, type soundsc(mtlb,Fs)

サンプル遅延させ、既知の係数 : により減衰させた信号のコピーを録音に追加して、エコーをモデル化します。タイム ラグを 0.23 秒に、減衰係数を 0.5 に指定します。

timelag = 0.23;
delta = round(Fs*timelag);
alpha = 0.5;

orig = [mtlb;zeros(delta,1)];
echo = [zeros(delta,1);mtlb]*alpha;

mtEcho = orig + echo;

元の信号、エコー、結果の信号をプロットします。

t = (0:length(mtEcho)-1)/Fs;

subplot(2,1,1)
plot(t,[orig echo])
legend('Original','Echo')

subplot(2,1,2)
plot(t,mtEcho)
legend('Total')
xlabel('Time (s)')

% To hear, type soundsc(mtEcho,Fs)

信号の自己相関の不偏推定を計算します。ゼロより大きいラグに対応するセクションを選択してプロットします。

[Rmm,lags] = xcorr(mtEcho,'unbiased');

Rmm = Rmm(lags>0);
lags = lags(lags>0);

figure
plot(lags/Fs,Rmm)
xlabel('Lag (s)')

自己相関は、エコーが到着した時点のラグで鋭いピークになります。出力 に従う IIR システムをとおして信号をフィルター処理することにより、エコーをキャンセルします。

[~,dl] = findpeaks(Rmm,lags,'MinPeakHeight',0.22);

mtNew = filter(1,[1 zeros(1,dl-1) alpha],mtEcho);

フィルター処理した信号をプロットして、元の信号と比較します。

subplot(2,1,1)
plot(t,orig)
legend('Original')

subplot(2,1,2)
plot(t,mtNew)
legend('Filtered')
xlabel('Time (s)')

% To hear, type soundsc(mtNew,Fs)

および () で得られる 11 サンプルの指数シーケンスを 3 つ作成します。stem3 を使用してシーケンスを並べてプロットします。

N = 11;
n = (0:N-1)';

a = 0.4;
b = 0.7;
c = 0.999;

xabc = [a.^n b.^n c.^n];

stem3(n,1:3,xabc','filled')
ax = gca;
ax.YTick = 1:3;
view(37.5,30)

シーケンスの自己相関と相互相関を計算します。継続して追跡する必要がないようにラグを出力します。ゼロ ラグで自己相関が単位値をもつように結果を正規化します。

[cr,lgs] = xcorr(xabc,'coeff');

for row = 1:3
    for col = 1:3
        nm = 3*(row-1)+col;
        subplot(3,3,nm)
        stem(lgs,cr(:,nm),'.')
        title(sprintf('c_{%d%d}',row,col))
        ylim([0 1])
    end
end

ラグの計算範囲を に制限します。

[cr,lgs] = xcorr(xabc,5,'coeff');

for row = 1:3
    for col = 1:3
        nm = 3*(row-1)+col;
        subplot(3,3,nm)
        stem(lgs,cr(:,nm),'.')
        title(sprintf('c_{%d%d}',row,col))
        ylim([0 1])
    end
end

自己相関と相互相関の不偏推定を計算します。既定では、ラグは の間で実行されます。

cu = xcorr(xabc,'unbiased');

for row = 1:3
    for col = 1:3
        nm = 3*(row-1)+col;
        subplot(3,3,nm)
        stem(-(N-1):(N-1),cu(:,nm),'.')
        title(sprintf('c_{%d%d}',row,col))
    end
end

について、28 サンプルの指数シーケンス の自己相関関数を計算します。

a = 0.95;

N = 28;
n = 0:N-1;
lags = -(N-1):(N-1);

x = a.^n;
c = xcorr(x);

を分析的に判断して結果が正しいかを確認します。より高いサンプル レートを使用して連続状況をシミュレートします。 において の場合、シーケンス の自己相関関数は次になります。

fs = 10;
nn = -(N-1):1/fs:(N-1);

dd = (1-a.^(2*(N-abs(nn))))/(1-a^2).*a.^abs(nn);

シーケンスを同じ図にプロットします。

stem(lags,c);
hold on
plot(nn,dd)
xlabel('Lag')
legend('xcorr','Analytic')
hold off

計算を繰り返しますが、今度は自己相関の "不偏推定" を求めます。不偏推定が で求められることを検証します。

cu = xcorr(x,'unbiased');

du = dd./(N-abs(nn));

stem(lags,cu);
hold on
plot(nn,du)
xlabel('Lag')
legend('xcorr','Analytic')
hold off

計算を繰り返しますが、今度は自己相関の "バイアス推定"を求めます。バイアス推定が で求められることを検証します。

cb = xcorr(x,'biased');

db = dd/N;

stem(lags,cb);
hold on
plot(nn,db)
xlabel('Lag')
legend('xcorr','Analytic')
hold off

ゼロ ラグの値が単位元である自己相関の推定を求めます。

cz = xcorr(x,'coeff');

dz = dd/max(dd);

stem(lags,cz);
hold on
plot(nn,dz)
xlabel('Lag')
legend('xcorr','Analytic')
hold off

2 つの 16 サンプルの指数シーケンス および () の相互相関を計算してプロットします。

N = 16;
n = 0:N-1;

a = 0.84;
b = 0.92;

xa = a.^n;
xb = b.^n;

r = xcorr(xa,xb);

stem(-(N-1):(N-1),r)

を分析的に判断して結果が正しいかを確認します。より高いサンプル レートを使用して連続状況をシミュレートします。 において の場合、シーケンス および の相互相関関数は次になります。

fs = 10;
nn = -(N-1):1/fs:(N-1);

cn = (1 - (a*b).^(N-abs(nn)))/(1 - a*b) .* ...
     (a.^nn.*(nn>0) + (nn==0) + b.^-(nn).*(nn<0));

シーケンスを同じ図にプロットします。

hold on
plot(nn,cn)

xlabel('Lag')
legend('xcorr','Analytic')

オペランドの順序を切り替えるとシーケンスが逆になることを確認します。

figure

stem(-(N-1):(N-1),xcorr(xb,xa))

hold on
stem(-(N-1):(N-1),fliplr(r),'--*')

xlabel('Lag')
legend('xcorr(x_b,x_a)','fliplr(xcorr(x_a,x_b))')

20 サンプルの指数シーケンス を生成します。 および の相互相関を計算し、プロットします。ラグを出力してプロットを簡単にします。xcorr では、短いシーケンスの最後にゼロを追加して長いシーケンスと長さを合わせます。

xc = 0.77.^(0:20-1);

[xca,la] = xcorr(xa,xc);
[xcb,lb] = xcorr(xb,xc);

figure

subplot(2,1,1)
stem(la,xca)
subplot(2,1,2)
stem(lb,xcb)
xlabel('Lag')

この例では、Parallel Computing Toolbox™ ソフトウェアおよび CUDA 対応 NVIDIA GPU (Compute Capability 1.3 以上) が必要です。詳細は、GPU System Requirementsを参照してください。

1 kHz でサンプリングされる、加法性ノイズをもつ 10 Hz の正弦波で構成される信号を作成します。gpuArray を使用して、コンピューターの GPU に保存される gpuArray オブジェクトを作成します。

t = 0:0.001:10-0.001;
x = cos(2*pi*10*t) + randn(size(t));
X = gpuArray(x);

ラグ 200 まで正規化された自己相関列を計算します。

[xc,lags] = xcorr(X,200,'coeff');

出力の xcgpuArray オブジェクトです。

gather を使用して、GPU から MATLAB® ワークスペースに倍精度ベクトルとしてデータを転送します。

xc = gather(xc);

入力引数

すべて折りたたむ

入力配列。ベクトル、行列または gpuArray オブジェクトとして指定します。

xcorrgpuArray オブジェクトとの使用の詳細は、GPU 計算 (Parallel Computing Toolbox)および GPU System Requirements を参照してください。

例: sin(2*pi*(0:9)/10) + randn([1 10])/10 では、行ベクトルとしてノイズを含んだ正弦波を指定します。

例: sin(2*pi*[0.1;0.3]*(0:39))' + randn([40 2])/10 では 2 チャネルのノイズを含んだ正弦波を指定します。

例: gpuArray(sin(2*pi*(0:9)/10) + randn([1 10])/10) では、gpuArray オブジェクトとしてノイズを含んだ正弦波を指定します。

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

入力配列。ベクトルまたは gpuArray オブジェクトとして指定します。

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

最大ラグ。整数スカラーとして指定します。maxlag を指定した場合、返される相互相関列の範囲は –maxlagmaxlag です。maxlag を指定しない場合、ラグの範囲は 2N – 1 となります。ここで N は x および y の長さよりも大きくなります。

データ型: single | double

正規化オプション。次のいずれかで指定します。

  • 'none' — スケーリングされていない生データの相互相関。これは、xy の長さが異なる場合に唯一可能なオプションです。

  • 'biased' — バイアスされた相互相関推定。

    R^xy,biased(m)=1NR^xy(m).

  • 'unbiased' — バイアスされていない相互相関推定。

    R^xy,unbiased(m)=1N|m|R^xy(m).

  • 'coeff' — ゼロ ラグでの自己相関が 1 になるようにシーケンスを正規化。

    R^xy,coeff(m)=1R^xx(0)R^yy(0)R^xy(m).

データ型: char

出力引数

すべて折りたたむ

相互相関列。ベクトル、行列または gpuArray オブジェクトを返します。

x がその列の N 個のチャネルを表す M × N 列の信号行列の場合、xcorr(x) では、x のチャネルの自己相関と相互相関をもつ (2M – 1) × N2 列の行列を返します。maxlag を指定した場合、r のサイズは (2 × maxlag – 1) × N2 となります。

たとえば、S が 3 チャネル信号 S=(x1x2x3) の場合、R = xcorr(S) の結果は次のようになります。

R=(Rx1x1Rx1x2Rx1x3Rx2x1Rx2x2Rx2x3Rx3x1Rx3x2Rx3x3).

ラグのインデックス。ベクトルとして返されます。

詳細

すべて折りたたむ

相互相関および自己相関

xcorr の結果は、2 つのランダム シーケンス間の相関の推定または 2 つの確定 的な信号間の確定的相関として解釈できます。

共に定常的なランダム過程の xn と yn 真の相互相関列は、次のようにして求められます。

Rxy(m)=E{xn+myn*}=E{xnynm*},

ここで、−∞ < n < ∞ で、アスタリスク (*) は複素共役を表し、E は期待値演算子です。実際には、無限長のランダム過程の 1 つの実現のある有限部分しか使用できないため、xcorr はシーケンスのみを推定します。

既定の設定では、xcorr では正規化を行わずに生データの相関が計算されます。

R^xy(m)={n=0Nm1xn+myn,m0,R^yx*(m),m<0.

出力ベクトルの c は、次の式で与えられる要素をもちます。

c(m)=R^xy(mN),m=1,2,,2N1.

一般に、相関関数は関数を正しく推定するために正規化を必要とします。この相関の正規化は、入力引数 scaleopt を使用して制御できます。

参照

[1] Buck, John R., Michael M. Daniel, and Andrew C. Singer. Computer Explorations in Signals and Systems Using MATLAB. 2nd Edition. Upper Saddle River, NJ: Prentice Hall, 2002.

[2] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.

拡張機能

R2006a より前に導入