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 を使用して、シンクロスクイーズド変換で使用されるスケールを取得できます。

制限

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

アルゴリズム

関数 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) を使用してください。

• 制限も参照してください。

GPU 配列

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

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

• 制限を参照してください。

バージョン履歴

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 | cwtfilterbank | cwt

トピック

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