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

findchangepts

信号の急激な変化の検出

構文

ipt = findchangepts(x)

ipt = findchangepts(x,Name,Value)

[ipt,residual] = findchangepts(___)

findchangepts(___)

説明

ipt = findchangepts(x) は、x の平均値が最も大幅に変化するインデックスを返します。

• x が N 個の要素を持つベクトルである場合、findchangepts は、x を 2 つの領域 x(1:ipt-1) と x(ipt:N) に分割します。これらは、各領域の局所的な平均値からの残差 (二乗) 誤差の和を最小化します。

• x が M 行 N 列の行列である場合、findchangepts は、x を 2 つの領域 x(1:M,1:ipt-1) と x(1:M,ipt:N) に分割し、各領域の局所的な M 次元の平均値からの残差誤差の和を最小化する列インデックスを返します。

ipt = findchangepts(x,Name,Value) は、名前と値のペアの引数を使用して追加オプションを指定します。オプションには、レポートする変化点の数や、平均値の代わりに測定する統計量があります。詳細は、変化点の検出を参照してください。

[ipt,residual] = findchangepts(___) は、前述のいずれかの指定を組み込んでモデル化された変化に対する信号の残差誤差も返します。

出力引数なしの findchangepts(___) は、信号および検出された任意の変化点をプロットします。詳細は、'Statistic'を参照してください。

例

1 次元および 2 次元の変化点

8192 Hz でサンプリングされた列車の警笛の録音を含むファイルを読み込みます。信号の平方根平均二乗レベルが最も大きく変化する点を 10 個検出します。

load train

findchangepts(y,'MaxNumChanges',10,'Statistic','rms')

信号の短時間パワー スペクトル密度を計算します。信号を 128 サンプルのセグメントに分割し、各セグメントにハミング ウィンドウを適用します。隣り合ったセグメント間のオーバーラップを 120 サンプル、DFT 点を 128 に指定します。パワー スペクトル密度が最も大きく変化する点を 10 個検出します。

[s,f,t,pxx] = spectrogram(y,128,120,128,Fs);

findchangepts(pow2db(pxx),'MaxNumChanges',10)

変化点の探索オプション

再現可能な結果が必要な場合は、乱数発生器をリセットします。次の条件でランダムな信号を生成します。

• 平均値は 7 つの各領域内では一定だが、領域から領域にかけて急激に変化する。

• 分散は 5 つの各領域内では一定だが、領域から領域にかけて急激に変化する。

rng('default')

lr = 20;

mns = [0 1 4 -5 2 0 1];
nm = length(mns);

vrs = [1 4 6 1 3];
nv = length(vrs);

v = randn(1,lr*nm*nv)/2;

f = reshape(repmat(mns,lr*nv,1),1,lr*nm*nv);
y = reshape(repmat(vrs,lr*nm,1),1,lr*nm*nv);

t = v.*y+f;

作成の手順を強調しながら信号をプロットします。

subplot(2,2,1)
plot(v)
title('Original')
xlim([0 700])

subplot(2,2,2)
plot([f;v+f]')
title('Means')
xlim([0 700])

subplot(2,2,3)
plot([y;v.*y]')
title('Variances')
xlim([0 700])

subplot(2,2,4)
plot(t)
title('Final')
xlim([0 700])

信号の平均値が最も大きく変化する点を 5 つ検出します。

figure
findchangepts(t,'MaxNumChanges',5)

信号の平方根平均二乗レベルが最も大きく変化する点を 5 つ検出します。

findchangepts(t,'MaxNumChanges',5,'Statistic','rms')

信号の平均値と標準偏差が最も大きく変化する点を検出します。

findchangepts(t,'Statistic','std')

オーディオ ファイルの分割

$$F_s=7418Hz$$ でサンプリングされた音声信号を読み込みます。ファイルには、"MATLAB®" という単語を発声している女性の録音音声が含まれています。

load mtlb

信号の分散が大きく変化している点を検出することによって、単語内の母音と子音を識別します。変化点の数を 5 つに制限します。

numc = 5;

[q,r] = findchangepts(mtlb,'Statistic','rms','MaxNumChanges',numc)

q = 5×1

         132
         778
        1646
        2500
        3454

r = -4.4055e+03

信号をプロットして変化点を表示します。

findchangepts(mtlb,'Statistic','rms','MaxNumChanges',numc)

各セグメントの後に一時停止しながら音声を再生するには、次の行のコメントを解除します。

% soundsc(1:q(1),Fs)
% for k = 1:length(q)-1
%     soundsc(mtlb(q(k):q(k+1)),Fs)
%     pause(1)
% end
% soundsc(q(end):length(mtlb),Fs)

平均、RMS レベル、標準偏差および勾配の変化

可変振幅をもつ 2 つの正弦波と線形トレンドで構成される信号を作成します。

vc = sin(2*pi*(0:201)/17).*sin(2*pi*(0:201)/19).* ...
    [sqrt(0:0.01:1) (1:-0.01:0).^2]+(0:201)/401;

信号の平均値が最も大きく変化している点を検出します。この場合、'Statistic' の名前と値のペアはオプションです。残差誤差の最小改善値に 1 を指定します。

findchangepts(vc,'Statistic','mean','MinThreshold',1)

信号の平方根平均二乗レベルが最も大きく変化する点を検出します。残差誤差の最小改善値に 6 を指定します。

findchangepts(vc,'Statistic','rms','MinThreshold',6)

信号の標準偏差が最も大きく変化する点を検出します。残差誤差の最小改善値に 10 を指定します。

findchangepts(vc,'Statistic','std','MinThreshold',10)

信号の平均値と勾配が最も急激に変化している点を検出します。残差誤差の最小改善値に 0.6 を指定します。

findchangepts(vc,'Statistic','linear','MinThreshold',0.6)

2 次元および 3 次元ベジエ曲線の変化点

20 のランダムな制御点を使用して、1000 のサンプルの 2 次元ベジエ曲線を生成します。ベジエ曲線は次のように定義されます。

$$C(t)={\displaystyle \sum _{k=0}^m\left(\genfrac{}{}{0ex}{}{m}{k}\right)t^k(1-t{)}^{m-k}{P}_k}$$、

ここで、${\mathit{P}}_{\mathit{k}}$ は $\mathit{k}$ 番目の $\mathit{m}$ 制御点、$\mathit{t}$ の範囲は 0 ～ 1、$$\left(\genfrac{}{}{0ex}{}{m}{k}\right)$$ は二項係数です。曲線と制御点をプロットします。

m = 20;
P = randn(m,2);
t = linspace(0,1,1000)';

pol = t.^(0:m-1).*(1-t).^(m-1:-1:0);
bin = gamma(m)./gamma(1:m)./gamma(m:-1:1);
crv = bin.*pol*P;

plot(crv(:,1),crv(:,2),P(:,1),P(:,2),'o:')

各セグメントの点とセグメントの平均値からの距離が最小になるように、曲線を 3 つのセグメントに分割します。

findchangepts(crv','MaxNumChanges',3)

曲線を、直線に当てはめた 20 のセグメントに分割します。

findchangepts(crv','Statistic','linear','MaxNumChanges',19)

20 のランダムな制御点を使用して、3 次元ベジエ曲線を生成してプロットします。

P = rand(m,3);
crv = bin.*pol*P;

plot3(crv(:,1),crv(:,2),crv(:,3),P(:,1),P(:,2),P(:,3),'o:')
xlabel('x')
ylabel('y')

曲線を上から可視化します。

view([0 0 1])

各セグメントの点とセグメントの平均値からの距離が最小になるように、曲線を 3 つのセグメントに分割します。

findchangepts(crv','MaxNumChanges',3)

曲線を、直線に当てはめた 20 のセグメントに分割します。

findchangepts(crv','Statistic','linear','MaxNumChanges',19)

入力引数

x — 入力信号

実数ベクトル

入力信号。実数ベクトルで指定します。

例: reshape(randn(100,3)+[-3 0 3],1,300) は、平均値に 2 回の急激な変化を伴うランダムな信号です。

例: reshape(randn(100,3).*[1 20 5],1,300) は、平方根平均二乗レベルに 2 回の急激な変化を伴うランダムな信号です。

データ型: single | double

名前と値のペアの引数

オプションの引数 Name,Value のコンマ区切りペアを指定します。Name は引数名で、Value は対応する値です。Name は引用符で囲まなければなりません。Name1,Value1,...,NameN,ValueN のように、複数の名前と値のペアの引数を、任意の順番で指定できます。

例: 'MaxNumChanges',3,'Statistic','rms','MinDistance',20 は、平方根平均二乗レベルの変化が最も大きく、かつ少なくとも 20 サンプルで隔てられている点を最大で 3 つ検出します。

'MaxNumChanges' — 返す大きな変化の最大数

1 (既定値) | 整数値スカラー

返す大きな変化の最大数。'MaxNumChanges' と整数スカラーで構成されるコンマ区切りのペアとして指定します。最も大きな変化のある点を検出した後、findchangepts は、指定した最大値を超えない範囲でより多くの変化点を含めるために、徐々にその検索基準を緩和します。何らかの探索設定によってこの最大値を超えて返されている場合、それ以降、関数は何も返しません。'MaxNumChanges' が指定されていない場合、関数は最も大きな変化のある点を返します。'MinThreshold' と 'MaxNumChanges' を同時に指定することはできません。

例: findchangepts([0 1 0]) は 2 番目のサンプルのインデックスを返します。

例: findchangepts([0 1 0],'MaxNumChanges',1) は空行列を返します。

例: findchangepts([0 1 0],'MaxNumChanges',2) は 2 番目の点と 3 番目の点のインデックスを返します。

データ型: single | double

'Statistic' — 検出する変化のタイプ

'mean' (既定値) | 'rms' | 'std' | 'linear'

検出する変化のタイプ。'Statistic' と次の値のいずれかで構成されるコンマ区切りペアとして指定します。

• 'mean' — 平均値の変化を検出。出力引数なしで findchangepts を呼び出すと、関数は、信号、変化点、連続変化点によって囲まれる各セグメントの平均値をプロットします。

• 'rms' — 平方根平均二乗レベルの変化を検出。出力引数なしで findchangepts を呼び出すと、関数は信号および変化点をプロットします。

• 'std' — ガウス対数尤度を使用して、標準偏差の変化を検出。出力引数なしで findchangepts を呼び出すと、関数は、信号、変化点、連続変化点によって囲まれる各セグメントの平均値をプロットします。

• 'linear' — 平均値と勾配の変化を検出。出力引数なしで findchangepts を呼び出すと、関数は、信号、変化点、および連続変化点によって囲まれる信号の各部分に最適なラインをプロットします。

例: findchangepts([0 1 2 1],'Statistic','mean') は 2 番目のサンプルのインデックスを返します。

例: findchangepts([0 1 2 1],'Statistic','rms') は 3 番目のサンプルのインデックスを返します。

'MinDistance' — 変化点間のサンプルの最小数

整数値スカラー

変化点間のサンプルの最小数。'MinDistance' と整数スカラーで構成されるコンマ区切りのペアとして指定します。この数値を指定しない場合、既定の設定は、平均値の変化については 1 に、その他の変化については 2 になります。

例: findchangepts(sin(2*pi*(0:10)/5),'MaxNumChanges',5,'MinDistance',1) は 5 インデックスを返します。

例: findchangepts(sin(2*pi*(0:10)/5),'MaxNumChanges',5,'MinDistance',3) は 2 インデックスを返します。

例: findchangepts(sin(2*pi*(0:10)/5),'MaxNumChanges',5,'MinDistance',5) はインデックスを返しません。

データ型: single | double

'MinThreshold' — 残差誤差の合計の最小改善値

実数スカラー

各変化点に対する残差誤差の合計の最小改善値。'MinThreshold' とペナルティを表す実数スカラーで構成されるコンマ区切りのペアとして指定します。このオプションは、予想される各変化点に追加のペナルティを適用することによって、返される大きな変化の数を制限するように作用します。'MinThreshold' と 'MaxNumChanges' を同時に指定することはできません。

例: findchangepts([0 1 2],'MinThreshold',0) は 2 インデックスを返します。

例: findchangepts([0 1 2],'MinThreshold',1) は 1 インデックスを返します。

例: findchangepts([0 1 2],'MinThreshold',2) はインデックスを返しません。

データ型: single | double

出力引数

ipt — 変化点の位置

ベクトル

変化点の位置。整数インデックスのベクトルとして返されます。

residual — 残差誤差

ベクトル

モデル化された変化に対する信号の残差誤差。ベクトルとして返されます。

詳細

変化点の検出

"変化点" は、信号のなんらかの統計プロパティが突然変化するサンプルまたは時点です。対象プロパティは、信号の平均値、その分散、スペクトル特性などが考えられます。

信号の変化点を検出するために、findchangepts はパラメトリックでグローバルな手法を採用します。関数は以下を実行します。

1. 1 点を選択し、信号を 2 つのセクションに分割します。

2. 各セクションの目的の統計プロパティの経験的推定値を計算します。

3. セクション内の各点で、このプロパティの経験的推定値との偏差を計測します。すべての点での偏差を足し合わせます。

4. セクション間の偏差を足し合わせて、残差誤差の合計を求めます。

5. 残差誤差の合計が最小になるまで、分割点の位置を変化させます。

この手続きは、選択した統計が平均値であるときに最もはっきりします。その場合、findchangepts は各セクションの "ベストな" 水平レベルからの残差誤差の合計を最小化します。信号 x₁, x₂, …, x_N が与えられると、関数は

$$\begin{array}{c}J={\displaystyle \sum _{i=1}^{k-1}{\left(x_i-{\langle x\rangle }_1^{k-1}\right)}^2}+{\displaystyle \sum _{i=k}^N{\left(x_i-{\langle x\rangle }_k^N\right)}^2}\\ ={\displaystyle \sum _{i=1}^{k-1}{\left(x_i-\frac{1}{k-1}{\displaystyle \sum _{r=1}^{k-1}{x}_r}\right)}^2}+{\displaystyle \sum _{i=k}^N{\left(x_i-\frac{1}{N-k+1}{\displaystyle \sum _{r=k}^N{x}_r}\right)}^2}\\ =\left(k-1\right)\mathrm{var}\left(\left[\begin{array}{ccc}{x}_1& \cdots & x_{k-1}\end{array}\right]\right)+\left(N-k+1\right)\mathrm{var}\left(\left[\begin{array}{ccc}{x}_k& \cdots & x_N\end{array}\right]\right),\end{array}$$

が最小になるような k を求めます。平均値と分散は通常の定義をもちます(分散は、1 を差し引かないサンプル数によって正規化されます)。

結果は、一般化して他の統計を組み込むことができます。findchangepts は、セクションの経験的推定 χ と偏差測定 Δ が与えられる場合、

$$J\left(k\right)={\displaystyle \sum _{i=1}^{k-1}\Delta \left(x_i;\chi \left(\left[\begin{array}{ccc}{x}_1& \cdots & x_{k-1}\end{array}\right]\right)\right)}+{\displaystyle \sum _{i=k}^N\Delta \left(x_i;\chi \left(\left[\begin{array}{ccc}{x}_k& \cdots & x_N\end{array}\right]\right)\right)}$$

が最小となるような k を求めます。残差誤差を最小化することは対数尤度を最大化することと等価です。平均値が μ、分散が σ² の標準偏差が与えられると、N 個の独立した観測値の対数尤度は

$$\log {\displaystyle \prod _{i=1}^N\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-{\left(x_i-\mu \right)}^2/2{\sigma }^2}}=-\frac{N}{2}\left(\log 2\pi +\log {\sigma }^2\right)-\frac{1}{2{\sigma }^2}{\displaystyle \sum _{i=1}^N{\left(x_i-\mu \right)}^2}.$$

です。

• 'Statistic' を 'mean' に指定した場合、分散は固定になり、関数は直前に得たとおり

  $$\begin{array}{c}{\displaystyle \sum _{i=m}^n\Delta \left(x_i;\chi \left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)|\right)}={\displaystyle \sum _{i=m}^n{\left(x_i-\mu \left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\right)}^2}\\ =(n-m+1)\mathrm{var}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right),\end{array}$$

  を使用します。

• 'Statistic' を 'std' に指定した場合、平均値は固定になり、関数は

  $$\begin{array}{c}{\displaystyle \sum _{i=m}^n\Delta \left(x_i;\chi \left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\right)}=\left(n-m+1\right)\log {\displaystyle \sum _{i=m}^n{\sigma }^2\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)}\\ =\left(n-m+1\right)\log \left(\frac{1}{n-m+1}{\displaystyle \sum _{i=m}^n{\left(x_i-\frac{1}{n-m+1}{\displaystyle \sum _{r=m}^n{x}_r}\right)}^2}\right)\\ =\left(n-m+1\right)\log \mathrm{var}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right).\end{array}$$

  を使用します。

• 'Statistic' を 'rms' に指定した場合、総偏差は 'std' と同じになりますが、平均値はゼロに設定されます。

  $${\displaystyle \sum _{i=m}^n\Delta \left(x_i;\chi \left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\right)}=\left(n-m+1\right)\log \left(\frac{1}{n-m+1}{\displaystyle \sum _{r=n}^m{x}_r^2}\right).$$

• 'Statistic' を 'linear' に指定した場合、関数は

  $$\begin{array}{c}{\displaystyle \sum _{i=m}^n\Delta \left(x_i;\chi \left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\right)}=\left(n-m+1\right)\mathrm{var}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\\ -\frac{{\left({\displaystyle \sum _{i=m}^n\left(x_i-\mu \left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\right)\left(i-\mu \left(\left[\begin{array}{cccc}m& m+1& \cdots & n\end{array}\right]\right)\right)}\right)}^2}{\left(n-m+1\right)\mathrm{var}\left(\left[\begin{array}{cccc}m& m+1& \cdots & n\end{array}\right]\right)}\\ \equiv S_{xx}-\frac{S_{xt}^2}{S_{tt}}.\end{array}$$

  を使用します。

多くの場合、対象信号は複数の変化点を持ちます。変化点の数が既知の場合、手続きの一般化は簡単です。変化点の数が不明の場合は、残差誤差にペナルティ項を追加しなければなりません。変化点が増えると常に残差誤差が少なくなり過適合になるからです。極端なケースでは、すべての点が変化点になり、残差誤差がゼロになります。findchangepts は変化点の数によって線形に増大するペナルティ項を使用します。K 個の変化点が見つかった場合、関数は以下のように最小化します。

$$J\left(K\right)={\displaystyle \sum _{r=0}^{K-1}{\displaystyle \sum _{i=k_r}^{k_{r+1}-1}\Delta \left(x_i;\chi \left(\left[\begin{array}{ccc}{x}_{k_r}& \cdots & x_{k_{r+1}-1}\end{array}\right]\right)\right)}+\beta K,}$$

ここで k₀ と k_K はそれぞれ信号の最初と最後のサンプルです。

• β と表記され 'MinThreshold' で指定される比例定数は、各変化点で追加された固定ペナルティに対応します。findchangepts は、残差誤差の減少がしきい値を満たさない場合は、さらなる変化点の追加を除外します。考えられるすべての変化を返すには、'MinThreshold' をゼロに設定します。

• 使用されているしきい値がわからない場合、または信号の変化点数が大まかにしかわからない場合は、代わりに 'MaxNumChanges' を指定します。このオプションは、関数が指定された値より少ない変化を検出するようになるまで、しきい値を徐々に増やします。

最小化そのものを実行するために、findchangepts は、早い段階で放棄する動的プログラミングに基づく網羅的なアルゴリズムを使用します。

参照

[1] Lavielle, Marc. “Using penalized contrasts for the change-point problem.” Signal Processing. Vol. 85, August 2005, pp. 1501–1510.

[2] Killick, Rebecca, Paul Fearnhead, and Idris A. Eckley. “Optimal detection of changepoints with a linear computational cost.” Journal of the American Statistical Association. Vol. 107, No. 500, 2012, pp. 1590–1598.

参考

cusum

R2016a で導入

×

Open Example (ja_JP)

A modified version of this example exists on your system. Do you want to open this version instead? (ja_JP)

No, overwrite the modified version (ja_JP) Yes (ja_JP)