境界効果と円錐状影響圏

このトピックでは、円錐状影響圏 (COI) について、およびそれを計算するために Wavelet Toolbox™ が使用する規則について説明します。また、このトピックでは、スカログラムのプロットで COI をどのように解釈するか、またcwtfilterbankとcwtで COI がどのように計算されるかについても説明します。

阪神淡路大震災の地震計信号を読み込みます。阪神淡路大震災の地震計信号のスカログラムをプロットします。このデータは 1 ヘルツでサンプリングされています。

load kobe
cwt(kobe,1)

このプロットには、スカログラムに加え、白い破線、および白い線の端から時間軸と周波数軸までの灰色の網掛け領域も表示されます。サンプル レートの代わりにサンプリング間隔を使用して同じデータをプロットします。スカログラムは周波数ではなく周期で表示されるようになりました。

cwt(kobe,seconds(1))

白い破線の向きは逆になっていますが、線と網掛け領域はまだ存在しています。

この白い線は、"円錐状影響圏" として知られるものを示します。円錐状影響圏は、この線、および線の端から周波数軸 (または周期軸) と時間軸までの網掛け領域で構成されます。円錐状影響圏は、"エッジ効果の副作用" による影響を受ける可能性があるスカログラムの領域を示します。これは、引き伸ばされたウェーブレットが観測間隔のエッジを超えて広がる領域によって、スカログラムが受ける影響のことです。白い線で囲まれた網掛けのない領域内では、スカログラムによって得られる情報が必ずデータの正確な時間-周波数表現となります。白い線の外側にある網掛け領域では、エッジ効果を受ける可能性があるため、スカログラムの情報を疑わしいものとして扱う必要があります。

中央に配置されたインパルスの CWT

円錐状影響圏について理解するため、まず、中央に配置された長さ 1024 サンプルのインパルス信号を作成します。cwtfilterbankを使用して、既定値で CWT フィルター バンクを作成します。wtを使用して、インパルスの CWT 係数と周波数を返します。より見やすくするために、各周波数 (各スケール) の最大絶対値が 1 になるように CWT 係数を正規化します。

x = zeros(1024,1);
x(512) = 1;
fb = cwtfilterbank;
[cfs,f] = wt(fb,x);
cfs = cfs./max(cfs,[],2);

スカログラムに補助関数 helperPlotScalogram を使用します。helperPlotFunction のコードはこの例の最後にあります。インパルスの位置を線でマークします。

ax = helperPlotScalogram(f,cfs);
hl = line(ax,[512 512],[min(f) max(f)], ...
    [max(abs(cfs(:))) max(abs(cfs(:)))]);
title("Scalogram of Centered Impulse")

黒い実線は、時間領域におけるインパルスの位置を示しています。周波数が減少すると、時間領域ではインパルスを中心としたゼロでない CWT 係数の幅が増加することに注意してください。逆に、周波数が増加すると、ゼロでない CWT 係数の幅は減少し、インパルスの中心に集まってきます。低い周波数はより長いスケールのウェーブレットに対応し、高い周波数はより短いスケールのウェーブレットに対応します。ウェーブレットが長くなると、インパルスの影響は長時間持続します。つまり、ウェーブレットが長いほど、信号による影響の持続期間が長くなります。ある時点を中心とするウェーブレットの場合、ウェーブレットを伸縮すると、ウェーブレットが "見る" 信号の量が増減します。これは、ウェーブレットの "円錐状影響圏" と呼ばれます。

境界効果

前のセクションでは、観測期間またはデータ期間の中央にあるインパルスの円錐状影響圏について説明しました。では、ウェーブレットがデータの先頭または末尾の近くにある場合はどうなるでしょうか。ウェーブレット変換では、ウェーブレットを膨張させるだけでなく、時間領域での平行移動も行います。データの先頭または末尾の近くにあるウェーブレットは、必然的に観測期間外のデータを "見る" ことになります。データの先頭および末尾の近くにあるウェーブレットの係数が境界の外側に広がるウェーブレットによって受ける影響を補正するため、さまざまな手法が使用されます。cwtfilterbank関数とcwt関数を使用すると、信号の対称的な折り返しまたは周期的な拡張によって境界を処理することができます。ただし、どの手法を使用する場合でも、ウェーブレット係数は対象信号の範囲外にある値によって影響を受けるため、境界付近のウェーブレットの係数を解釈する際には注意が必要です。さらに、観測期間外のデータによってウェーブレット係数が影響を受ける範囲は、スケール (周波数) によって異なります。スケールが長くなるほど、円錐状影響圏も大きくなります。

先ほどと同じインパルスの例を使用しますが、今度は 2 つのインパルスを使用し、データの先頭と末尾に 1 つずつ配置します。円錐状影響圏も返します。より見やすくするために、各周波数 (各スケール) の最大絶対値が 1 になるように CWT 係数を正規化します。

dirac = zeros(1024,1);
dirac([1 1024]) = 1;
[cfs,f,coi] = wt(fb,dirac);
cfs = cfs./max(cfs,[],2);
helperPlotScalogram(f,cfs)
title("Scalogram of Two-Impulse Signal")

これを見ると、観測期間で最も外側に位置する境界に対する円錐状影響圏が、ウェーブレットのスケールに依存する程度まで観測期間内に広がっていることがはっきり分かります。そのため、観測期間の十分内側にあるウェーブレットの係数は、信号の境界にあるウェーブレットが観測するデータによって影響を受ける可能性があります。また、何らかの方法で信号を拡張した場合は、実際の信号境界の手前でも影響を受ける可能性があります。

前の図からだけでも、cwtfilterbankによって返された円錐状影響圏、つまりcwt関数によってプロットされた円錐状影響圏と、2 つのインパルス信号のスカログラム係数がゼロでない領域との間に、顕著な類似性を見て取ることができます。

ウェーブレット係数の解釈では、このような境界効果を理解することが重要です。しかし、各スケールにおける円錐状影響圏の範囲を決定するための数学的に厳密なルールは存在しません。Nobach ら[2]は、各スケールにおける円錐状影響圏の範囲を、ウェーブレット変換の振幅がピーク値の 2% に減衰する点として定義しています。連続ウェーブレット解析で使用されるウェーブレットの多くは時間の経過と共に指数関数的に減衰するため、Torrence と Compo[3] は、各スケールにおける円錐状影響圏の境界を時定数 $$1/e$$ で規定しています。Lilly[1] は、Morse ウェーブレットに対し、ウェーブレットのエネルギーの約 95% を含む期間である "ウェーブレット フットプリント" の概念を用いています。Lilly は、各スケールにおいて、観測期間の先頭でウェーブレット フットプリントの 1/2 を加算し、観測期間の末尾でウェーブレット フットプリントの 1/2 を減算することにより、COI を規定しています。

cwtfilterbank関数およびcwt関数は、$$1/e$$ ルールの近似値を使用して COI を規定します。この近似では、観測期間の先頭で各スケールに対して 1 つの時間領域標準偏差を加算し、観測期間の末尾で各スケールに対して 1 つの時間領域標準偏差を減算します。この対応関係を示す前に、計算された COI を前のプロットに追加します。

helperPlotScalogram(f,cfs,coi)
title("Scalogram with Cone of Influence")

計算された COI は、信号の先頭と末尾にあるインパルスによって著しく影響を受ける範囲の境界を適切に近似していることがわかります。

cwtfilterbankおよびcwtがこのルールを明示的に計算する方法を示すため、解析的 Morlet ウェーブレットと既定の Morse ウェーブレットの 2 つの例を考えてみましょう。まず、解析的 Morlet ウェーブレットについて見ていきます。この場合、1 つの時間領域標準偏差によるルールが、Torrence と Compo[3] の使用した時間領域折り返しの式と厳密に一致します。

fb = cwtfilterbank(Wavelet="amor");
[~,f,coi] = wt(fb,dirac);

Torrence と Compo は、COI を $$\sqrt{2}s$$ で表しています。ここで、$$s$$ はスケールです。cwtfilterbankおよびcwtによる解析的 Morlet ウェーブレットの場合、これは次のように与えられます。

cf = 6/(2*pi);
predtimes = sqrt(2)*cf./f;

cwtfilterbankによって返される COI を、Torrence および Compo が使用した式と共にプロットします。

plot(1:1024,coi,"k--",LineWidth=2)
hold on
grid on
plot(predtimes,f,"r*")
plot(1024-predtimes,f,"r*")
hold off
set(gca,YScale="log")
axis tight
legend("COI","Predicted COI",Location="best")
xlabel("Samples")
ylabel("Hz")
title("Cone of Influence - Analytic Morlet Wavelet")

最後に、cwtfilterbankおよびcwtによる既定の Morse ウェーブレットで同じことを行う例を示します。既定の Morse ウェーブレットでは、時間領域標準偏差は 5.5008 で、ピーク周波数は 0.2995 サイクル/サンプルです。ウェーブレット バンドパス フィルターの中心周波数および時間領域標準偏差のルールを使用して COI の予測を取得し、cwtfilterbankによって返される値と比較します。

fb = cwtfilterbank;
[~,f,coi] = wt(fb,dirac);
sd = 5.5008;
cf = 0.2995;
predtimes = cf./f*sd;
figure
plot(1:1024,coi,"k--",LineWidth=2)
hold on
grid on
plot(predtimes,f,"r*")
plot(1024-predtimes,f,"r*")
hold off
set(gca,'yscale',"log")
axis tight
legend("COI","Predicted COI",Location="best")
xlabel("Samples")
ylabel("Hz")
title("Cone of Influence - Default Morse Wavelet")

付録

この例では次の補助関数が使用されています。

helperPlotScalogram

function varargout = helperPlotScalogram(f,cfs,coi)
nargoutchk(0,1);
ax = newplot;
surf(ax,1:1024,f,abs(cfs),EdgeColor="none")
ax.YScale = "log";
clim([0.01 1])
colorbar
grid on
ax.YLim = [min(f) max(f)];
ax.XLim = [1 size(cfs,2)];
view(0,90)

xlabel("Time")
ylabel("Cycles/Sample")

if nargin == 3
    hl = line(ax,1:1024,coi,ones(1024,1));
    hl.Color = "k";
    hl.LineWidth = 2;
end

if nargout > 0
    varargout{1} = ax;
end

end

参照

[1] Lilly, J. M. “Element analysis: a wavelet-based method for analysing time-localized events in noisy time series.” Proceedings of the Royal Society A. Volume 473: 20160776, 2017, pp. 1–28. dx.doi.org/10.1098/rspa.2016.0776.

[2] Nobach, H., Tropea, C., Cordier, L., Bonnet, J. P., Delville, J., Lewalle, J., Farge, M., Schneider, K., and R. J. Adrian. "Review of Some Fundamentals of Data Processing." Springer Handbook of Experimental Fluid Mechanics (C. Tropea, A. L. Yarin, and J. F. Foss, eds.). Berlin, Heidelberg: Springer, 2007, pp. 1337–1398.

[3] Torrence, C., and G. Compo. "A Practical Guide to Wavelet Analysis." Bulletin of the American Meteorological Society. Vol. 79, Number 1, 1998, pp. 61–78.

参考

アプリ

• 信号アナライザー (Signal Processing Toolbox)

関数

• cwt | cwtfilterbank | pspectrum (Signal Processing Toolbox)

トピック

• Morse ウェーブレット
• 連続および離散ウェーブレット変換