wsst

ウェーブレット シンクロスクイーズド変換

構文

sst = wsst(x)

[sst,f] = wsst(x)

[sst,f,fbparam] = wsst(x)

[___] = wsst(x,fs)

[___] = wsst(x,ts)

[___] = wsst(___,wav)

wsst(___)

[___] = wsst(___,Name=Value)

説明

sst = wsst(x) は、x のウェーブレット シンクロスクイーズド変換を返します。関数 wsst は、解析 Morlet ウェーブレットを使用してシンクロスクイーズド変換を計算します。

シンクロスクイーズド変換では、連続ウェーブレット変換 (CWT) と比較して、エネルギーの不鮮明さが低減します。関数 wsst は、解析ウェーブレットを正規化して L1 ノルムを保持します。詳細については、アルゴリズムを参照してください。

[sst,f] = wsst(x) は、周波数ベクトル f をサンプルあたりのサイクル数で返します。

[sst,f,fbparam] = wsst(x) は、シンクロスクイーズド変換で使用される CWT フィルター バンク パラメーターを返します。

[___] = wsst(x,fs) は、指定されたサンプリング周波数 fs を使用して、シンクロスクイーズド変換を返します。

[___] = wsst(x,ts) は、指定されたサンプリング間隔 ts を使用して、シンクロスクイーズド変換を返します。

[___] = wsst(___,wav) は、wav で指定された解析ウェーブレットを使用して、シンクロスクイーズド変換を計算します。

出力引数なしの wsst(___) は、シンクロスクイーズド変換の大きさを時間と周波数の関数としてプロットします。サンプリング周波数またはサンプリング間隔を指定しない場合、シンクロスクイーズド変換はサンプルあたりのサイクル数でプロットされます。サンプリング周波数を指定した場合、シンクロスクイーズド変換がヘルツ単位でプロットされます。duration を使用してサンプリング間隔を指定した場合、プロットは単位時間あたりのサイクルになります。時間単位は duration から導出されます。

[___] = wsst(___,Name=Value) は、1 つ以上の名前と値の引数で指定される追加オプションを使用して、シンクロスクイーズド変換を返します。

例

音声信号のシンクロスクイーズド変換

既定値を使用して、音声サンプルのウェーブレット シンクロスクイーズド変換を取得します。

load mtlb
sst = wsst(mtlb);

音声信号のシンクロスクイーズド変換と再構成

音声信号を読み込みます。信号は $$F_s=7418$$ Hz でサンプリングされています。シンクロスクイーズド変換を求めます。

load mtlb
% To hear, type sound(mtlb,Fs)
dt = 1/Fs;
t = 0:dt:numel(mtlb)*dt-dt;
[sst,f] = wsst(mtlb,Fs);

シンクロスクイーズド変換の大きさをプロットします。

pcolor(t,f,abs(sst))
shading interp
xlabel("Seconds")
ylabel("Frequency (Hz)")
title("Synchrosqueezed Transform")

逆シンクロスクイーズド変換を求めます。元の信号と再構成したものをプロットします。

xrec = iwsst(sst);
% To hear, type soundsc(xrec,Fs)
plot(t,mtlb)
hold on
plot(t,xrec)
hold off
xlabel("Seconds")
ylabel("Amplitude")
legend("Original","Reconstruction")

シンクロスクイーズド変換と CWT の比較

2 成分信号のシンクロスクイーズド変換を取得し、同じ信号の CWT と比較します。入力は、振幅変調信号と周波数変調信号の組み合わせです。

2 成分信号を作成します。信号は振幅変調信号と周波数変調信号の組み合わせです。サンプル レート 1000 Hz で 10 分の 1 秒間信号をサンプリングします。

Fs = 1e3;
t = 0:1/Fs:0.1;
x1 = (2+0.5*cos(2*pi*10*t)).*cos(2*pi*200*t+10*sin(2*pi*5*t));
x2 = cos(2*pi*75*t);
sig = x1+x2;

信号のシンクロスクイーズド変換を取得します。境界の影響を最小限に抑えるために、信号を対称的に拡張します。また、基底 CWT フィルター バンクのパラメーターも取得します。シンクロスクイーズド変換の大きさをプロットします。

[sst,sstF,fbparams] = wsst(sig,Fs,ExtendSignal=true);
pcolor(t,sstF,abs(sst))
shading interp
colorbar
title("WSST Magnitude")
xlabel("Time (s)")
ylabel("Frequency (Hz)")

基底フィルター バンク パラメーターを使用して CWT フィルター バンクを作成します。次に、フィルター バンクを使用して、信号の CWT を取得します。CWT のスカログラムをプロットし、シンクロスクイーズド変換の結果と比較します。

cfb = cwtfilterbank(fbparams);
[cfs,cfsF] = cwt(sig,FilterBank=cfb);
pcolor(t,cfsF,abs(cfs))
colorbar
shading interp
title("Scalogram")
xlabel("Time (s)")
ylabel("Frequency (Hz)")

2 次チャープのシンクロスクイーズド変換

2 次チャープのウェーブレット シンクロスクイーズド変換を求めます。チャープは 1000 Hz でサンプリングされます。

load quadchirp
[sst,f] = wsst(quadchirp,1000);

変換の大きさをプロットします。

hp = pcolor(tquad,f,abs(sst));
hp.EdgeColor = "none";
title("Wavelet Synchrosqueezed Transform")
xlabel("Time (s)")
ylabel("Frequency (Hz)")

太陽黒点データのシンクロスクイーズド変換

解析 Morse (3,30) ウェーブレットを使用し、太陽黒点データのウェーブレット シンクロスクイーズド変換をプロットします。サンプリング間隔を 1 年に指定します。

load sunspot
wsst(sunspot(:,2),years(1),"Morse",WaveletParameters=[3 30])

Bump ウェーブレットを使用した太陽黒点データのシンクロスクイーズド変換

Bump ウェーブレットを使用し、太陽黒点データのウェーブレット シンクロスクイーズド変換をプロットします。サンプリング間隔を 1 年に指定します。

load sunspot
wsst(sunspot(:,2),years(1),"bump")

入力引数

x — 入力信号

ベクトル

入力信号。ベクトルとして指定します。x は少なくとも 4 つのサンプルをもたなければなりません。

データ型: single | double

fs — サンプリング周波数

正のスカラー

ヘルツ単位のサンプリング周波数。正のスカラーとして指定します。

ts — サンプリング間隔

duration スカラー

サンプリング間隔 (サンプリング周期とも呼ばれる)。duration スカラーとして指定します。ts は、x のサンプル間の時間です。有効な duration は、years、days、hours、minutes、および seconds です。ts と fs の両方を指定することはできません。

例: seconds(1) は、連続する信号サンプル間の 1 秒間の時間差を表す duration スカラーです。

wav — 解析ウェーブレット

"amor" (既定値) | "Morse" | "bump"

シンクロスクイーズド変換の計算に使用される解析ウェーブレット。次のいずれかとして指定します。

• "amor" — 解析 Morlet ウェーブレット

• "Morse" — 解析 Morse ウェーブレット (R2024a 以降)

• "bump" — Bump ウェーブレット

逆シンクロスクイーズド変換で使用する場合、通常、解析 Morlet または解析 Morse (3,30) ウェーブレットが最良の結果を提供します。

名前と値の引数

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

例: VoicesPerOctave=26

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

例: 'VoicesPerOctave',26

VoicesPerOctave — オクターブあたりの音の数

32 (既定値) | 10 ～ 48 の偶数の整数

シンクロスクイーズド変換で使用するオクターブあたりの音の数。10 ～ 48 の偶数の整数として指定します。オクターブあたりの音の数とオクターブ数の積が音階の数になります。

WaveletParameters — Morse ウェーブレット パラメーター

[3,60] (既定値) | スカラーの 2 要素ベクトル

R2024a 以降

Morse ウェーブレット パラメーター。2 要素ベクトルとして指定します。最初の要素は対称性パラメーター (γ) で、1 以上でなければなりません。2 番目の要素は時間-帯域積で、ガンマ以上でなければなりません。γ に対する時間-帯域積の比が 40 を超えることはできません。時間-帯域幅パラメーターとガンマの比は、Morse ウェーブレットに関する文献において "ベータ パラメーター" と呼ばれることもあります。

γ が 3 に等しい場合、Morse ウェーブレットは周波数領域において完全に対称です。この歪度は 0 と等しくなります。ガンマの値が 3 より大きい場合は歪度が正になり、ガンマの値が 3 より小さい場合は歪度が負になります。

モード再構成にシンクロスクイージングを使用している場合、通常、Morse (3,30) ウェーブレットが良好な逆シンクロスクイージングの結果を提供します。

WaveletParameters は、解析 Morse ウェーブレットを指定した場合にのみ有効です。

詳細については、Morse ウェーブレットを参照してください。

ExtendSignal — 信号の対称的な拡張

false または 0 (既定値) | true または 1

信号の対称な拡張の logical 値。数値または logical の 0 (false) または 1 (true) として指定します。true として指定した場合、関数 wsst は信号を対称的に拡張します。信号を対称的に拡張すると境界の影響を軽減できます。

出力引数

sst — シンクロスクイーズド変換

行列

シンクロスクイーズド変換。Na 行 N 列の行列として返されます。ここで、Na はスケールの数、N は x のサンプル数です。

スケールの数は、最小周波数からナイキスト周波数 (サンプル レートの半分) までの周波数範囲が使用され、基底 CWT フィルター バンク オブジェクトによって決定されます。周波数ベクトル f とフィルター バンク パラメーター fbparam を使用して、シンクロスクイーズド変換で使用されるスケールを取得できます。

データ型: single | double

f — 周波数

ベクトル

シンクロスクイーズド変換の周波数。ベクトルとして返されます。周波数は sst の行に対応します。f の長さと sst の行数は等しくなります。

• サンプリング周波数 fs またはサンプリング間隔 ts を指定しない場合、周波数はサンプルあたりのサイクル数になります。

• サンプリング周波数を指定した場合、周波数の単位はヘルツです。

• サンプリング間隔を指定した場合、周波数は単位時間あたりのサイクル数で表され、時間単位は ts と同じ時間単位になります。関数は ts を使用して、f のスケールから周波数への変換を計算します。

fbparam — パラメーター

構造体配列

R2024a 以降

シンクロスクイーズド変換で使用される基底 CWT フィルター バンクのパラメーター。構造体配列として返されます。フィールド名は、シンクロスクイーズド変換で使用される CWT フィルター バンクのプロパティを示します。各フィールドの値はプロパティ値です。

f と fbparam を使用して、シンクロスクイーズド変換で使用されるスケールを取得できます。

アルゴリズム

関数 wsst は、解析ウェーブレットを正規化して L1 ノルムを保持します。これを言い換えると、wsst はウェーブレット バンドパス フィルターのフーリエ変換にスケールの平方根を乗算しないということです。スケールの平方根を乗算すると、各バンドパスの寄与の重み付けが均等ではなくなります。

L1 正規化では、振幅の等しい振動成分がスケールの異なるデータ内にある場合、そのような振動成分は CWT で等しい大きさをもちます。関数 cwt も L1 正規化を使用します。詳細については、CWT の L1 ノルムを参照してください。

参照

[1] Daubechies, Ingrid, Jianfeng Lu, and Hau-Tieng Wu. “Synchrosqueezed Wavelet Transforms: An Empirical Mode Decomposition-like Tool.” Applied and Computational Harmonic Analysis 30, no. 2 (March 2011): 243–61. https://doi.org/10.1016/j.acha.2010.08.002.

[2] Thakur, Gaurav, Eugene Brevdo, Neven S. Fučkar, and Hau-Tieng Wu. “The Synchrosqueezing Algorithm for Time-Varying Spectral Analysis: Robustness Properties and New Paleoclimate Applications.” Signal Processing 93, no. 5 (May 2013): 1079–94. https://doi.org/10.1016/j.sigpro.2012.11.029.

拡張機能

C/C++ コード生成

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

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

• 入力ウェーブレット名 wav はコンパイル時に定数でなければなりません。coder.Constant (MATLAB Coder) を使用してください。

• ウェーブレット シンクロスクイーズド変換アルゴリズムは単精度データに対して敏感です。このため、MATLAB^® コードを使用して取得される結果と生成した C/C++ コードを使用して取得される結果が一致しない場合があります。この不一致は、周波数の再割り当てにおける off-by-one エラーに起因します。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

wsst 関数は GPU 配列入力をサポートしますが、次の使用上の注意および制限があります。

• ウェーブレット シンクロスクイーズド変換アルゴリズムは単精度データに対して敏感です。このため、CPU で取得される結果と GPU で取得される結果が一致しない場合があります。この不一致は、周波数の再割り当てにおける off-by-one エラーに起因します。

詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2016a で導入

R2024a: CWT フィルター バンク パラメーターの取得

関数 wsst がシンクロスクイージング変換で使用する基底 CWT フィルター バンクのパラメーターを取得できます。

R2024a: wsst による解析 Morse ウェーブレットのサポート

関数 wsst は解析 Morse ウェーブレットをサポートします。

R2024a: wsst は C/C++ コード生成と gpuArray オブジェクトをサポートする

関数 wsst は以下をサポートします。

• C/C++ コード生成。C/C++ コードを生成するには、MATLAB Coder™ が必要です。

• gpuArray オブジェクト。この機能を使用するには、Parallel Computing Toolbox™ が必要です。

参考

関数

• iwsst | wsstridge | cwt

オブジェクト

• cwtfilterbank

トピック

• Practical Introduction to Time-Frequency Analysis Using the Continuous Wavelet Transform
• Time-Frequency Reassignment and Mode Extraction with Synchrosqueezing
• ウェーブレット シンクロスクイージング
• Morse ウェーブレット
• 時間-周波数ギャラリー