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を参照してください。

メモ

関数 findchangepts は、プロットの前に現在の Figure をクリアします (clf)。信号および検出された変化点をサブプロットにプロットするには、プロット関数を使用します。オーディオ ファイルの分割を参照してください。

例

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);

変化点インデックスに基づいて音声信号の信号マスクを作成します。信号マスクの使用の詳細については、signalMaskを参照してください。

t = (0:length(mtlb)-1)/Fs;
roitable = ([[1;q] [q;length(mtlb)]]);

x = ["M" "A" "T" "L" "A" "B"]';
c = categorical(x,unique(x,"stable"));

msk = signalMask(table(t(roitable),c),SampleRate=Fs,RightShortening=1);
roimask(msk)

ans=6×2 table
           Var1            c
    ___________________    _

          0    0.017525    M
    0.01766     0.10461    A
    0.10475     0.22162    T
    0.22176     0.33675    L
    0.33688     0.46535    A
    0.46549     0.53909    B

音声信号および検出された変化点を関心領域とともに信号マスクからサブプロットにプロットします。

• 上のサブプロットで、関数 plotsigroi を使用して信号マスク領域を可視化します。カラー バーが上に表示されるように設定を調整します。

• 下のサブプロットで、元の音声信号をプロットし、検出された変化点を垂直線として追加します。

subplot(2,1,1)
plotsigroi(msk,mtlb)
colorbar("off")
nc = numel(c)-1;
colormap(gca,lines(nc));
colorbar(TickLabels=categories(c),Ticks=1/2/nc:1/nc:1, ...
    TickLength=0,Location="northoutside")
xlabel("")

subplot(2,1,2)
plot(t,mtlb)
hold on
xline(q/Fs)
hold off
xlim([0 t(end)])
xlabel("Seconds")

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

% for k = 1:length(roitable)
%     intv = roitable(k,1):roitable(k,2);
%     soundsc(mtlb(intv).*hann(length(intv)),Fs)
%     pause(.5)
% end

平均、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{m}$ 個の制御点の $\mathit{k}$ 番目、$\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

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前は、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みました。

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

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

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

返す大きな変化の最大数。整数スカラーとして指定します。最も大きな変化のある点を検出した後、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"

検出する変化のタイプ。次の値のいずれかとして指定します。

• "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 — 変化点間のサンプルの最小数

整数値スカラー

変化点間のサンプルの最小数。整数スカラーとして指定します。この数値を指定しない場合、既定の設定は、平均値の変化については 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 と 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}\mathrm{mean}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)=\frac{1}{n-m+1}{\displaystyle \sum _{r=m}^n{x}_r,}\\ \mathrm{var}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)=\frac{1}{n-m+1}{\displaystyle \sum _{r=m}^n{\left(x_r-\mathrm{mean}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\right)}^2\equiv \frac{{S_{xx}|}_m^n}{n-m+1}},\end{array}$$

ここで、"二乗和" は次のようになります。

$${S_{xy}|}_m^n\equiv {\displaystyle \sum _{r=m}^n\left(x_r-\mathrm{mean}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\right)\left(y_r-\mathrm{mean}\left(\left[\begin{array}{ccc}{y}_m& \cdots & y_n\end{array}\right]\right)\right)},$$

findchangepts は、

$$\begin{array}{c}J={\displaystyle \sum _{i=1}^{k-1}{\left(x_i-\mathrm{mean}\left(\left[\begin{array}{ccc}{x}_1& \cdots & x_{k-1}\end{array}\right]\right)\right)}^2}+{\displaystyle \sum _{i=k}^N{\left(x_i-\mathrm{mean}\left(\left[\begin{array}{ccc}{x}_k& \cdots & x_N\end{array}\right]\right)\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 を求めます。結果は、一般化して他の統計を組み込むことができます。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-\mathrm{mean}\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-\mathrm{mean}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\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=m}^n{x}_r^2}\right).$$

• Statistic を "linear" として指定した場合、関数は、信号値とその値に近似する最小二乗線形予測の差の二乗和を合計偏差として使用します。この量は、"残差二乗和" または "SSE" と呼ばれます。x_m, x_m+1, …, x_n での最適な近似線は次のとおりです。

  $$\widehat{x}\left(t\right)=\frac{{S_{xt}|}_m^n}{{S_{tt}|}_m^n}\left(t-\mathrm{mean}\left(\left[\begin{array}{ccc}{t}_m& \cdots & t_n\end{array}\right]\right)\right)+\mathrm{mean}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)$$

  また、SSE は次のとおりです。

  $$\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-\widehat{x}\left(t_i\right)\right)}^2}\\ ={S_{xx}|}_m^n-\frac{{S_{xt}^2|}_m^n}{{S_{tt}|}_m^n}\\ =\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-\mathrm{mean}\left(\left[\begin{array}{ccc}{x}_m& \cdots & x_n\end{array}\right]\right)\right)\left(i-\mathrm{mean}\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)}.\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] 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.

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

拡張機能

C/C++ コード生成

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

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

• 入力引数 Statistic を指定する場合、それはコンパイル時の定数でなければなりません。

バージョン履歴

R2016a で導入

参考

cusum

×

例を開く

この例の変更されたバージョンがあります。編集された方の例を開きますか?

いいえ、変更されたバージョンを上書きします。 はい