cwt

連続 1 次元ウェーブレット変換

構文

wt = cwt(x)

wt = cwt(x,wname)

[wt,f] = cwt(___,fs)

[wt,period] = cwt(___,ts)

[wt,f,coi] = cwt(___)

[wt,period,coi] = cwt(___,ts)

[___,coi,fb] = cwt(___)

[___,fb,scalingcfs] = cwt(___)

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

cwt(___)

説明

wt = cwt(x) は、x の連続ウェーブレット変換 (CWT) を返します。CWT は、対称性パラメーターであるガンマ ($$\gamma$$) が 3 に等しく時間-帯域積が 60 に等しい解析 Morse ウェーブレットを使用して得られます。cwt で使用されるオクターブあたりの音の数は 10 です。最小スケールと最大スケールは、ウェーブレットの周波数と時間におけるエネルギーの広がりに基づいて自動的に決まります。

関数 cwt では L1 正規化を使用します。L1 正規化では、振幅の等しい振動成分がスケールの異なるデータ内にある場合、そのような振動成分は CWT で等しい大きさをもちます。L1 正規化を使用すると信号の表現の精度が高くなります。CWT の L1 ノルムおよび2 つの複素指数の連続ウェーブレット変換を参照してください。

wt = cwt(x,wname) は、wname で指定された解析ウェーブレットを使用して CWT を計算します。

[wt,f] = cwt(___,fs) は、サンプリング周波数 fs をヘルツ単位で指定し、スケールから周波数への変換 f をヘルツ単位で返します。サンプリング周波数を指定しない場合、cwt は f をサンプルあたりのサイクル数で返します。

[wt,period] = cwt(___,ts) は、サンプリング周期 ts を正の duration スカラーとして指定します。cwt は、ts を使用して、スケールから周期への変換である period を計算します。period は、ts と同じ書式プロパティをもつ duration の配列です。

[wt,f,coi] = cwt(___) は、円錐状影響圏 coi をサンプルあたりのサイクル数で返します。円錐状影響圏をヘルツの単位で返すには、サンプリング周波数 fs をヘルツ単位で指定します。

[wt,period,coi] = cwt(___,ts) は、円錐状影響圏 coi を、ts と同じ書式プロパティをもつ duration の配列として返します。

[___,coi,fb] = cwt(___) は、CWT で使用されるフィルター バンクを返します。cwtfilterbank を参照してください。

[___,fb,scalingcfs] = cwt(___) は、ウェーブレット変換のスケーリング係数を返します。

[___] = cwt(___,Name=Value) は、名前と値の引数を 1 つ以上追加で指定します。たとえば、wt = cwt(x,TimeBandwidth=40,VoicesPerOctave=20) は、時間-帯域積を 40、オクターブあたりの音の数を 20 と指定します。

出力引数なしで cwt(___) を使用すると、CWT のスカログラムがプロットされます。スカログラムは、時間と周波数の関数としてプロットされた CWT の絶対値です。対数スケールで周波数がプロットされます。エッジの影響が顕著になる円錐状影響圏もプロットされます。白い破線の外側のグレーの領域がエッジの影響が顕著になる領域です。入力信号が複素数値の場合、正 (反時計回り) と負 (時計回り) の成分が別々のスカログラムにプロットされます。

サンプリング周波数またはサンプリング周期を指定しない場合、周波数はサンプルあたりのサイクル数としてプロットされます。サンプリング周波数を指定した場合、周波数の単位はヘルツです。サンプリング周期を指定した場合は、スカログラムが時間と周期の関数としてプロットされます。入力信号が timetable の場合は、スカログラムが時間と周波数の関数として Hz 単位でプロットされ、RowTimes が時間軸の基底として使用されます。

スカログラムの点の時間、周波数、および振幅を表示するには、Figure の座標軸のツール バーでデータ ヒントを有効にし、スカログラムで目的の点をクリックします。

メモ

cwt は、プロットの前に現在の Figure をクリアします (clf)。スカログラムをサブプロットにプロットするには、プロット関数を使用します。CWT のスカログラムのサブプロットへのプロットを参照してください。

例

既定値を使用した連続ウェーブレット変換

既定値を使用して、音声サンプルの連続ウェーブレット変換を取得します。

load mtlb;
w = cwt(mtlb);

指定のウェーブレットを使用した連続ウェーブレット変換

音声サンプルを読み込みます。

load mtlb

mtlb.mat ファイルを読み込むと、音声信号 mtlb とサンプル レート Fs がワークスペースに取り込まれます。Bump ウェーブレットを使用して取得した音声サンプルのスカログラムを表示します。

load mtlb
cwt(mtlb,"bump",Fs)

既定の Morse ウェーブレットを使用して取得したスカログラムと比較します。

cwt(mtlb,Fs)

地震データの連続ウェーブレット変換

阪神淡路大震災のデータの CWT を取得します。このデータは、オーストラリアのホバートにあるタスマニア大学の地震計で 1995 年 1 月 16 日 20:56:51 (GMT) から 51 分間にわたって記録された測定値 (垂直加速度、nm/sq.sec) です。サンプリング周波数は 1 Hz です。

load kobe

地震データをプロットします。

plot((1:numel(kobe))./60,kobe)
xlabel("Time (mins)")
ylabel("Vertical Acceleration (nm/s^2)")
title("Kobe Earthquake Data")
grid on
axis tight

CWT、周波数、および円錐状影響圏を求めます。

[wt,f,coi] = cwt(kobe,1);

円錐状影響圏を含むスカログラムを表示します。

cwt(kobe,1)

サンプリング周波数の代わりにサンプリング周期を指定して、CWT、期間、および円錐状影響圏を求めます。

[wt,periods,coi] = cwt(kobe,minutes(1/60));

サンプリング周期を指定するときに生成されるスカログラムを表示します。

cwt(kobe,minutes(1/60))

2 つの複素指数の連続ウェーブレット変換

周波数が 32 Hz と 64 Hz の振幅が異なる 2 つの複素指数を作成します。データは 1000 Hz でサンプリングされます。2 つの複素指数の時間のサポートは互いに素です。

Fs = 1e3;
t = 0:1/Fs:1;
z = exp(1i*2*pi*32*t).*(t>=0.1 & t<0.3)+2*exp(-1i*2*pi*64*t).*(t>0.7);

標準偏差 0.05 の複素ホワイト ガウス ノイズを付加します。

wgnNoise = 0.05/sqrt(2)*randn(size(t))+1i*0.05/sqrt(2)*randn(size(t));
z = z+wgnNoise;

Morse ウェーブレットを使用して、CWT を求めてプロットします。

cwt(z,Fs)

カラー バーの複素指数成分の大きさは、スケールが異なっていても基本的にそれらの振幅になることに注意してください。これは L1 正規化による直接的な結果です。これは、このスクリプトを実行し、データ カーソルで各サブプロットを調べると確認できます。

正弦波とウェーブレット係数の振幅

この例では、信号の振動成分の振幅が対応するウェーブレット係数の振幅と一致することを示します。

時間のサポートが互いに素である 2 つの正弦波で構成される信号を作成します。一方の正弦波は周波数が 32 Hz で振幅は 1 です。もう一方の正弦波は周波数が 64 Hz で振幅は 2 です。この信号を 1000 Hz で 1 秒間サンプリングします。信号をプロットします。

frq1 = 32;
amp1 = 1;
frq2 = 64;
amp2 = 2;

Fs = 1e3;
t = 0:1/Fs:1;
x = amp1*sin(2*pi*frq1*t).*(t>=0.1 & t<0.3)+... 
    amp2*sin(2*pi*frq2*t).*(t>0.6 & t<0.9);

plot(t,x)
grid on
xlabel("Time (sec)")
ylabel("Amplitude")
title("Signal")

信号に適用できる CWT フィルター バンクを作成します。信号の成分の周波数がわかっているため、フィルター バンクの周波数範囲をそれらの周波数を含む狭い範囲に設定します。範囲を確認するには、フィルター バンクの振幅周波数応答をプロットします。

fb = cwtfilterbank(SignalLength=numel(x),SamplingFrequency=Fs,...
    FrequencyLimits=[20 100]);
freqz(fb)

cwt とフィルター バンクを使用して信号のスカログラムをプロットします。

cwt(x,FilterBank=fb)

データ カーソルを使用して、ウェーブレット係数の振幅が正弦波成分の振幅と基本的に等しいことを確認します。結果は次の Figure のようになるはずです。

複数の時系列での CWT フィルター バンクの使用

この例では、CWT フィルター バンクを使用することで、複数の時系列の CWT を取得する場合の計算効率がどのように高まる可能性があるかを示します。

100 行 1024 列の行列 x を作成します。1024 個のサンプルをもつ信号に適した CWT フィルター バンクを作成します。

x = randn(100,1024);
fb = cwtfilterbank;

既定の設定で cwt を使用して、1024 個のサンプルをもつ信号の CWT を取得します。100 個の信号の CWT 係数を格納できる 3 次元配列を作成します。各信号には、1,024 個のサンプルがあります。

cfs = cwt(x(1,:));
res = zeros(100,size(cfs,1),size(cfs,2));

関数 cwt を使用して、行列 x の各行の CWT を取得します。経過時間を表示します。

tic
for k=1:100
    res(k,:,:) = cwt(x(k,:));
end
toc

Elapsed time is 0.928160 seconds.

次に、フィルター バンクのオブジェクト関数 wt を使用して、x の各行の CWT を取得します。経過時間を表示します。

tic
for k=1:100
    res(k,:,:) = wt(fb,x(k,:));
end
toc

Elapsed time is 0.393524 seconds.

CWT からの CUDA コード

この例では、MEX ファイルを生成し、生成された CUDA® コードを使用して連続ウェーブレット変換 (CWT) を実行する方法について説明します。

まず、CUDA 対応 GPU と NVCC コンパイラがあることを確認します。GPU 環境のチェックおよびアプリの設定 (GPU Coder)を参照して、適切な構成があることを確認してください。

GPU Coder 構成オブジェクトを作成します。

cfg = coder.gpuConfig("mex");

1,000 Hz で 100,000 サンプルの信号を生成します。この信号は、時間のサポートが互いに素である 2 つの余弦波で構成されています。

t = 0:.001:(1e5*0.001)-0.001;
x = cos(2*pi*32*t).*(t > 10 & t<=50)+ ...
    cos(2*pi*64*t).*(t >= 60 & t < 90)+ ...
    0.2*randn(size(t));

単精度を使用するように信号をキャストします。GPU 計算は、ほとんどの場合、単精度で実行するほうが効率が高くなります。ただし、NVIDIA® GPU が倍精度をサポートしている場合、倍精度のコードを生成することもできます。

x = single(x);

GPU MEX ファイルとコード生成レポートを生成します。MEX ファイルの生成を許可するには、次の 3 つの入力パラメーターのプロパティ (クラス、サイズ、実数/複素数) を指定しなければなりません。

• coder.typeof(single(0),[1 1e5]) は、実数の single 値を含む長さ 100,000 の行ベクトルを指定します。

• coder.typeof('c',[1 inf]) は、任意の長さの文字配列を指定します。

• coder.typeof(0) は、実数の double 値を指定します。

sig = coder.typeof(single(0),[1 1e5]);
wav = coder.typeof('c',[1 inf]);
sfrq = coder.typeof(0);
codegen cwt -config cfg -args {sig,wav,sfrq} -report

Code generation successful: View report

-report フラグはオプションです。-report を使用すると、コード生成レポートが生成されます。レポートの [概要] タブに、[GPU コード メトリクス] リンクが表示され、ここで、生成された CUDA カーネルの数、割り当てられたメモリ量などの詳細情報が提供されます。

データに対して MEX ファイルを実行し、スカログラムをプロットします。プロットが 2 つの互いに素な余弦波と一致することを確認します。

[cfs,f] = cwt_mex(x,'morse',1e3);
image("XData",t,"YData",f,"CData",abs(cfs),"CDataMapping","scaled")
set(gca,"YScale","log")
axis tight
xlabel("Time (Seconds)")
ylabel("Frequency (Hz)")
title("Scalogram of Two-Tone Signal")

_mex を付けずに上記の CWT コマンドを実行します。MATLAB® と GPU MEX のスカログラムが同一であることを確認します。

[cfs2,f2] = cwt(x,'morse',1e3);
max(abs(cfs2(:)-cfs(:)))

ans = single
    7.3380e-07

既定の周波数軸ラベルの変更

この例では、出力引数なしでプロットを取得した場合の CWT の既定の周波数軸ラベルを変更する方法を示します。

周波数が 32 Hz と 64 Hz の 2 つの正弦波を作成します。データは 1000 Hz でサンプリングされます。2 つの正弦波の時間のサポートは互いに素です。標準偏差 0.05 のホワイト ガウス ノイズを付加します。既定の Morse ウェーブレットを使用して、CWT を求めてプロットします。

Fs = 1e3;
t = 0:1/Fs:1;
x = cos(2*pi*32*t).*(t>=0.1 & t<0.3)+sin(2*pi*64*t).*(t>0.7);
wgnNoise = 0.05*randn(size(t));
x = x+wgnNoise;
cwt(x,1000)

CWT の周波数は対数であるため、プロットでは対数の周波数軸が使用されます。MATLAB の対数軸の単位は 10 のべき乗です。cwtfreqboundsを使用すると、特定の信号長、サンプリング周波数、およびウェーブレットについて、ウェーブレット バンドパスの最小周波数と最大周波数を特定できます。

[minf,maxf] = cwtfreqbounds(numel(x),1000);

MATLAB の既定の設定により、最小周波数から最大周波数までの間の 10 のべき乗である 10 と 100 の位置に周波数の目盛りがあることがわかります。周波数軸の目盛りを増やす場合は、最小周波数から最大周波数までの間の対数的に等間隔である一連の周波数を以下を使用して求めることができます。

numfreq = 10;
freq = logspace(log10(minf),log10(maxf),numfreq);

次に、以下を使用して、現在の座標軸のハンドルを取得し、周波数軸の目盛りとラベルを置き換えます。

AX = gca;
AX.YTickLabelMode = "auto";
AX.YTick = freq;

CWT では、周波数が 2 のべき乗で計算されます。2 のべき乗の周波数の目盛りと目盛ラベルを作成するには、以下を実行します。

newplot
cwt(x,1000)
AX = gca;
freq = 2.^(round(log2(minf)):round(log2(maxf)));
AX.YTickLabelMode = "auto";
AX.YTick = freq;

スカログラムのカラーリングの変更

この例では、プロットする各レベルの最大絶対値でスカログラムの値をスケーリングする方法を示します。

信号を読み込み、既定のスカログラムを表示します。カラーマップを pink(240) に変更します。

load noisdopp
cwt(noisdopp)
colormap(pink(240))

信号の CWT を実行し、ウェーブレット係数と周波数を求めます。

[cfs,frq] = cwt(noisdopp);

各周波数 (レベル) の係数の最大値を効率的に見つけるには、最初に係数の絶対値を転置します。各レベルの最小値を求めます。それぞれのレベルについて、そのレベルの最小値を減算します。

tmp1 = abs(cfs);
t1 = size(tmp1,2);
tmp1 = tmp1';
minv = min(tmp1);
tmp1 = (tmp1-minv(ones(1,t1),:));

tmp1 の各レベルの最大値を求めます。それぞれのレベルについて、そのレベルの最大値で各値を除算します。結果をカラーマップの色の数で乗算します。ゼロのエントリをすべて 1 に設定します。結果を転置します。

maxv = max(tmp1);
maxvArray = maxv(ones(1,t1),:);
indx = maxvArray<eps;
tmp1 = 240*(tmp1./maxvArray);
tmp2 = 1+fix(tmp1);
tmp2(indx) = 1;
tmp2 = tmp2';

結果を表示します。スカログラムの値が各レベルの最大絶対値でスケーリングされます。周波数が線形スケールで表示されます。

t = 0:length(noisdopp)-1;
pcolor(t,frq,tmp2)
shading interp
xlabel("Time (Samples)")
ylabel("Normalized Frequency (cycles/sample)")
title("Scalogram Scaled By Level")
colormap(pink(240))
colorbar

時間-帯域積の変更

この例では、Morse ウェーブレットの時間-帯域積 $$P^2$$ が大きくなるとウェーブレットの包絡線下部の振動が増えることを示します。$$P^2$$ が大きくなると、ウェーブレットの周波数が狭まります。

2 つのフィルター バンクを作成します。1 つ目のフィルター バンクの TimeBandwidth の値は既定の 60 です。2 つ目のフィルター バンクの TimeBandwidth の値は 10 です。SignalLength はどちらのフィルター バンクも 4096 サンプルです。

sigLen = 4096;
fb60 = cwtfilterbank(SignalLength=sigLen);
fb10 = cwtfilterbank(SignalLength=sigLen,TimeBandwidth=10);

フィルター バンクの時間領域のウェーブレットを求めます。

[psi60,t] = wavelets(fb60);
[psi10,~] = wavelets(fb10);

関数 scales を使用して各フィルター バンクのマザー ウェーブレットを検出します。

sca60 = scales(fb60);
sca10 = scales(fb10);
[~,idx60] = min(abs(sca60-1));
[~,idx10] = min(abs(sca10-1));
m60 = psi60(idx60,:);
m10 = psi10(idx10,:);

fb60 フィルター バンクの方が時間-帯域積が大きいため、m60 ウェーブレットの方が m10 ウェーブレットよりも包絡線の下の振動が多いことを確認します。

subplot(2,1,1)
plot(t,abs(m60))
grid on
hold on
plot(t,real(m60))
plot(t,imag(m60))
hold off
xlim([-30 30])
legend("abs(m60)","real(m60)","imag(m60)")
title("TimeBandwidth = 60")
subplot(2,1,2)
plot(t,abs(m10))
grid on
hold on
plot(t,real(m10))
plot(t,imag(m10))
hold off
xlim([-30 30])
legend("abs(m10)","real(m10)","imag(m10)")
title("TimeBandwidth = 10")

m60 と m10 の振幅周波数応答のピークを揃えます。m60 ウェーブレットの周波数応答が m10 ウェーブレットの周波数応答よりも狭いことを確認します。

cf60 = centerFrequencies(fb60);
cf10 = centerFrequencies(fb10);

m60cFreq = cf60(idx60);
m10cFreq = cf10(idx10);

freqShift = 2*pi*(m60cFreq-m10cFreq);
x10 = m10.*exp(1j*freqShift*(-sigLen/2:sigLen/2-1));

figure
plot([abs(fft(m60)).' abs(fft(x10)).'])
grid on
legend("Time-bandwidth = 60","Time-bandwidth = 10")
title("Magnitude Frequency Responses")

CWT のスカログラムのサブプロットへのプロット

この例では、CWT のスカログラムを Figure のサブプロットにプロットする方法を示します。

音声サンプルを読み込みます。データは 7418 Hz でサンプリングされます。CWT の既定のスカログラムをプロットします。

load mtlb
cwt(mtlb,Fs)

信号の連続ウェーブレット変換と CWT の周波数を求めます。

[cfs,frq] = cwt(mtlb,Fs);

関数 cwt により、時間と周波数の座標軸がスカログラムで設定されます。サンプル時間を表すベクトルを作成します。

tms = (0:numel(mtlb)-1)/Fs;

新しい Figure で、元の信号を上のサブプロットにプロットし、スカログラムを下のサブプロットにプロットします。対数スケールで周波数をプロットします。

figure
subplot(2,1,1)
plot(tms,mtlb)
axis tight
title("Signal and Scalogram")
xlabel("Time (s)")
ylabel("Amplitude")
subplot(2,1,2)
surface(tms,frq,abs(cfs))
axis tight
shading flat
xlabel("Time (s)")
ylabel("Frequency (Hz)")
set(gca,"yscale","log")

入力引数

x — 入力信号

実数値または複素数値のベクトル | timetable | gpuArray

入力信号。実数値または複素数値のベクトル、あるいは単一変数の一定間隔でサンプリングされる timetable として指定します。入力 x は少なくとも 4 つのサンプルをもたなければなりません。

関数 cwt は GPU 配列入力も受け入れます。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

データ型: single | double
複素数のサポート: あり

wname — 解析ウェーブレット

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

CWT の計算に使用する解析ウェーブレット。wname の有効なオプションは、"morse"、"amor"、および "bump" で、それぞれ Morse、Morlet (ガボール)、Bump のウェーブレットを示します。

既定では、対称性パラメーターであるガンマ ($$\gamma$$) が 3 に等しく時間-帯域積が 60 に等しい Morse ウェーブレットが使用されます。

データ型: char | string

fs — サンプリング周波数

正のスカラー

ヘルツ単位のサンプリング周波数。正のスカラーとして指定します。fs を指定する場合、ts は指定できません。x が timetable の場合は fs は指定できません。fs は timetable の RowTimes から決定されます。

データ型: single | double

ts — サンプリング周期

スカラーの duration

サンプリング周期 (期間とも呼ばれます)。スカラーの duration として指定します。有効な duration は、years、days、hours、minutes、および seconds です。カレンダー期間は使用できません。ts を指定する場合、fs は指定できません。x が timetable の場合は ts は指定できません。ts は、PeriodLimits の名前と値の引数を設定したときに timetable の RowTimes から決定されます。

例: wt = cwt(x,hours(12))

データ型: duration

名前と値の引数

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

例: wt = cwt(x,"bump",VoicesPerOctave=10) は、Bump ウェーブレットとオクターブあたり 10 の音の数を使用して、x の CWT を返します。

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

例: wt = cwt(x,"ExtendedSignal",true,"FrequencyLimits",[0.1 0.2]) は、入力信号を対称的に拡張し、サイクルあたり 0.1 ～ 0.2 サンプルで周波数範囲を指定します。

ExtendSignal — 入力信号を対称的に拡張

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

入力信号を反射で対称的に拡張するオプション。次のいずれかとして指定します。

• 1 (true) — 対称的に拡張する

• 0 (false) — 対称的に拡張しない

ExtendSignal が false の場合、信号は周期的に拡張されます。信号を対称的に拡張すると境界の影響を軽減できます。

メモ

スケーリング係数と近似合成フィルターをもつ icwt を使用して CWT を反転する場合は、ExtendSignal を false に設定します。

データ型: logical

FrequencyLimits — 周波数の範囲

2 要素スカラー ベクトル

CWT で使用する周波数の範囲。厳密に増加する正のエントリをもつ 2 要素ベクトルとして指定します。

• 最初の要素は最も低いピーク通過帯域周波数を示し、ウェーブレットのピーク周波数 (Hz 単位) と 2 つの時間の標準偏差の積を信号長で除算した値以上でなければなりません。

• 2 番目の要素は最も高いピーク通過帯域周波数を示し、ナイキスト周波数以下でなければなりません。

• 周波数の下限 freqMin に対する周波数の上限 freqMax の比の基底 2 の対数は 1/NV 以上でなければなりません。ここで NV はオクターブあたりの音の数です。

  log2(freqMax/freqMin) ≥ 1/NV.

複素数値の信号の場合、反解析的な部分には (-1) × flimits が使用されます。ここで flimits は FrequencyLimits で指定されたベクトルです。

メモ

許容される範囲外の周波数範囲を指定した場合、cwt では有効な最小値と最大値まで範囲が切り詰められます。CWT の各種のパラメーター表現の周波数範囲を特定するには、cwtfreqbounds を使用します。

例: wt = cwt(x,1000,VoicesPerOctave=10,FrequencyLimits=[80 90])

データ型: double

PeriodLimits — 周期の範囲

2 要素 duration 配列

CWT で使用する周期の範囲。厳密に増加する正のエントリをもつ 2 要素 duration 配列として指定します。

• 最初の要素は 2 × ts 以上でなければなりません。ここで ts はサンプリング周期です。

• 最大周期は、ウェーブレットの 2 つの時間の標準偏差とウェーブレットのピーク周波数の積で信号長を除算した値を超えることはできません。

• 最大周期 maxP に対する最小周期 minP の比の基底 2 の対数は -1/NV 以下でなければなりません。ここで NV はオクターブあたりの音の数です。

  log2(pMin/pMax) ≤ -1/NV.

複素数値の信号の場合、反解析的な部分には (-1) × plimits が使用されます。ここで plimits は PeriodLimits で指定されたベクトルです。

メモ

許容される範囲外の周期範囲を指定した場合、cwt では有効な最小値と最大値まで範囲が切り詰められます。ウェーブレット変換の各種のパラメーター表現の周期範囲を求めるには、cwtfreqbounds を使用します。

例: wt = cwt(x,seconds(0.1),VoicesPerOctave=10,PeriodLimits=[seconds(0.2) seconds(3)])

データ型: duration

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

10 (既定値) | 1 から 48 までの整数

CWT に使用するオクターブあたりの音の数。1 ～ 48 までの整数として指定します。指定したオクターブあたりの音の数を使用して CWT のスケールが離散化されます。ウェーブレットの周波数と時間のエネルギーの広がりは、最小スケールと最大スケールによって自動的に決まります。

TimeBandwidth — Morse ウェーブレットの時間-帯域積

60 (既定値) | 3 以上 120 以下のスカラー

Morse ウェーブレットの時間-帯域積。3 以上 120 以下のスカラーとして指定します。対称性パラメーターのガンマ ($$\gamma$$) は 3 に固定されます。時間-帯域積が大きいウェーブレットほど、時間の広がりは大きくなり、周波数の広がりは狭くなります。Morse ウェーブレットの時間の標準偏差は約 sqrt(TimeBandwidth/2) です。Morse ウェーブレットの周波数の標準偏差は約 1/2 × sqrt(2/TimeBandwidth) です。

TimeBandwidth を指定する場合、WaveletParameters は指定できません。対称性と時間-帯域積の両方を指定するには、代わりに WaveletParameters を使用します。

Morse ウェーブレット の表記法では、TimeBandwidth は P² です。

WaveletParameters — Morse ウェーブレットの対称性と時間-帯域積

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

Morse ウェーブレットの対称性と時間-帯域積。スカラーの 2 要素ベクトルとして指定します。最初の要素は対称性 $$\gamma$$ で、1 以上でなければなりません。2 番目の要素は時間-帯域積で、$$\gamma$$ 以上でなければなりません。$$\gamma$$ に対する時間-帯域積の比が 40 を超えることはできません。

$$\gamma$$ が 3 に等しい場合、Morse ウェーブレットは周波数領域において完全に対称であり、歪度は 0 です。$$\gamma$$ が 3 よりも大きい場合、歪度は正になります。$$\gamma$$ が 3 よりも小さい場合、歪度は負になります。

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

WaveletParameters を指定する場合、TimeBandwidth は指定できません。

FilterBank — CWT フィルター バンク

cwtfilterbank オブジェクト

CWT の計算に使用する CWT フィルター バンク。cwtfilterbank オブジェクトとして指定します。FilterBank を設定する場合、他のオプションは指定できません。CWT の計算についてのすべてのオプションが、フィルター バンクのプロパティとして定義されます。詳細については、cwtfilterbank を参照してください。

x が timetable の場合、fb のサンプリング周波数またはサンプリング周期は timetable の RowTimes で決まるサンプリング周波数またはサンプリング周期と一致していなければなりません。

例: wt = cwt(x,FilterBank=cfb)

出力引数

wt — 連続ウェーブレット変換

行列

連続ウェーブレット変換。複素数値の行列として返されます。既定では、cwt は解析 Morse (3,60) ウェーブレットを使用します。3 は対称性で、60 は時間-帯域積です。cwt で使用されるオクターブあたりの音の数は 10 です。

• x が実数値の場合、wt は Na 行 N 列の行列になります。Na はスケールの数で、N は x に含まれるサンプルの数です。

• x が複素数値の場合、wt は 3 次元行列です。1 ページ目が正のスケール (解析的な部分または反時計回りの成分) の CWT で、2 ページ目が負のスケール (反解析的な部分または時計回りの成分) の CWT になります。

最小スケールと最大スケールは、ウェーブレットの周波数と時間におけるエネルギーの広がりに基づいて自動的に決まります。スケールが決定される方法の詳細については、アルゴリズムを参照してください。

データ型: single | double

f — スケールから周波数への変換

ベクトル

CWT のスケールから周波数への変換。ベクトルとして返されます。サンプリング周波数 fs を指定した場合、f の単位はヘルツです。fs を指定しない場合、cwt は f をサンプルあたりのサイクル数で返します。入力 x が複素数の場合、wt の両方のページにスケールから周波数への変換が適用されます。

period — スケールから周期への変換

配列

スケールから周期への変換。ts と同じ書式プロパティをもつ duration の配列として返されます。各行が周期に対応しています。入力 x が複素数の場合、wt の両方のページにスケールから周期への変換が適用されます。

coi — 円錐状影響圏

実数の配列 | duration の配列

CWT の円錐状影響圏。サンプリング周波数 fs を指定した場合、円錐状影響圏の単位はヘルツです。スカラーの duration である ts を指定すると、円錐状影響圏は ts と同じ書式プロパティをもつ duration の配列になります。入力 x が複素数の場合、wt の両方のページに円錐状影響圏が適用されます。

円錐状影響圏は、CWT でエッジの影響が現れる箇所を示します。エッジの影響により、円錐状影響圏の外側やそれと重なる領域は信頼性が低くなります。詳細については、Boundary Effects and the Cone of Influenceを参照してください。

fb — CWT フィルター バンク

cwtfilterbank オブジェクト

CWT で使用される CWT フィルター バンク。cwtfilterbank オブジェクトとして返されます。cwtfilterbank を参照してください。

scalingcfs — スケーリング係数

実数値または複素数値のベクトル

CWT のスケーリング係数。実数値または複素数値のベクトルとして返されます。scalingcfs の長さは、入力 x の長さと等しくなります。

詳細

解析ウェーブレット

解析ウェーブレットは、負の周波数に対してはフーリエ変換がゼロになる複素数値のウェーブレットです。解析ウェーブレットは、CWT を使用して時間-周波数解析を行う場合に適しています。ウェーブレット係数は複素数値であるため、解析する信号の位相と振幅の情報を提供します。解析ウェーブレットは、実際の非定常信号の周波数成分が時間の関数に応じてどのように変化するかを調べるのに非常に適しています。

解析ウェーブレットは、ほぼ例外なく急減少関数に基づきます。$$\psi (t)$$ が解析的な時間の急減少関数である場合、そのフーリエ変換 $$\widehat{\psi }(\omega )$$ は周波数の急減少関数であり、一部の区間 $$\alpha <\omega <\beta$$ ($$0<\alpha <\beta$$) を除いて小さくなります。直交ウェーブレットと双直交ウェーブレットは、一般に時間のサポートがコンパクトになるように設計されています。時間のサポートがコンパクトなウェーブレットは、時間が急速に減少するウェーブレットに比べると、周波数のエネルギー集中度が比較的低くなります。直交ウェーブレットと双直交ウェーブレットのほとんどは、フーリエ領域において対称ではありません。

信号の結合時間-周波数表現を求めることが目的の場合は、cwt または cwtfilterbank を使用することをお勧めします。どちらの関数も次の解析ウェーブレットをサポートします。

• Morse ウェーブレット ファミリ (既定)

• 解析的な Morlet (ガボール) ウェーブレット

• Bump ウェーブレット

Morse ウェーブレットの詳細については、Morse ウェーブレットを参照してください。フーリエ領域では、角周波数の観点から次のようになります。

• 解析的 Morlet は次のように定義されます。

  

  ここで、 は区間 [0,∞) の指示関数です。

• Bump ウェーブレットは次のように定義されます。

  

  ここで、ϵ = 2.2204×10^-16 です。

直交ウェーブレットまたは双直交ウェーブレットを使用して時間-周波数解析を実行する場合は、modwpt をお勧めします。

時間-周波数解析にウェーブレットを使用する場合、通常はスケールを周波数または周期に変換して結果を解釈します。この変換は cwt および cwtfilterbank で実行されます。関連付けられている対応するスケールは、cwt のオプションの出力引数 fb で scales を使用して取得できます。

用途に適したウェーブレットを選択するためのガイダンスについては、ウェーブレットの選択を参照してください。

ヒント

• 古い関数 cwt 用の構文は引き続き機能しますが、非推奨になりました。現在のバージョンの cwt を使用してください。古いバージョンと現在のバージョンの両方で同じ関数名が使用されています。関数への入力によって、使用されるバージョンが自動的に決定されます。関数 cwt の構文は変更済みを参照してください。

• 複数の CWT を実行する場合 (for ループ内など)、まず cwtfilterbank オブジェクトを作成してから、オブジェクト関数 wt を使用するワークフローを推奨します。このワークフローにより、オーバーヘッドを最小限に抑え、パフォーマンスを最大限に向上することができます。複数の時系列での CWT フィルター バンクの使用を参照してください。

アルゴリズム

最小スケール

最小スケールを特定するには、基底ウェーブレットのピーク周波数 $${\omega }_x$$ を求めます。Morse ウェーブレットの場合、ウェーブレットの $$\pi$$ ラジアンにおけるフーリエ変換がピーク周波数の 10% に等しくなるようにウェーブレットを膨張させます。最小スケールは最大周波数で発生します。

$$s_0=\frac{{\omega }_x^{\text{'}}}{\pi }$$

結果として、最小スケールは (2, $$s_0$$) の最小値になります。Morse ウェーブレットの場合、最小スケールは通常は $$s_0$$ です。Morlet ウェーブレットの場合、最小スケールは通常は 2 です。

最大スケール

CWT の最小スケールと最大スケールは、どちらもウェーブレットの周波数と時間のエネルギーの広がりに基づいて自動的に決まります。最大スケールを特定するために、CWT は次のアルゴリズムを使用します。

Morse ウェーブレットの時間の標準偏差 $${\sigma }_t$$ は約 $$\sqrt{\frac{P^2}{2}}$$ です。ここで $$P^2$$ は時間-帯域積です。周波数の標準偏差 $${\sigma }_f$$ は約 $$\frac{1}{2}\sqrt{\frac{2}{P^2}}$$ です。ウェーブレットを $$s>1$$ のいずれかでスケーリングすると、期間が $$2s{\sigma }_t=N$$ に変わり、ウェーブレットが入力の全長 (N サンプル) に等しくなるまで引き伸ばされます。このウェーブレットを変換したりラップなしでさらに引き伸ばしたりすることはできないため、最大スケールは $$floor\left(\frac{N}{2{\sigma }_t}\right)$$ になります。

ウェーブレット変換のスケールは 2 のべき乗であり、$$s_0{\left(2^{\frac{1}{NV}}\right)}^j$$ で示されます。NV はオクターブあたりの音の数で、j の範囲は 0 から最大スケールまでです。特定の小さいスケールについて、$$s_0$$ は次のようになります。

$$s_0{\left(2^{\frac{1}{NV}}\right)}^j\le \frac{N}{2{\sigma }_t}$$

log2 に変換します。

$$j{\log }_2\left(2^{\frac{1}{NV}}\right)\le {\log }_2\left(\frac{N}{2{\sigma }_t{s}_0}\right)$$

$$j\le NV{\log }_2\left(\frac{N}{2{\sigma }_t{s}_0}\right)$$

したがって、最大スケールは次のようになります。

$$s_0{\left(2^{\frac{1}{NV}}\right)}^{floor\left(NV{\log }_2\left(\frac{N}{2{\sigma }_t{s}_0}\right)\right)}$$

CWT の L1 ノルム

CWT では、積分形式ではエネルギーが維持されます。ただし、CWT を数値的に実装する場合はエネルギーが維持されません。この場合、使用する正規化に関係なく、CWT は正規直交変換ではありません。関数 cwt では L1 正規化を使用します。

ウェーブレット変換では、一般にウェーブレットの L2 正規化が使用されます。L2 ノルムでは、次のように 1/s (s は 0 より大きい) による信号の膨張が定義されます。

$${\Vert x\left(\frac{t}{s}\right)\Vert }_2^2=s{\Vert x\left(t\right)\Vert }_2^2$$

これにより、エネルギーは元のエネルギーの s 倍になります。フーリエ変換に含まれる場合は、$$\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{s}$}\right.$$ で乗算することで、高い周波数のピークの方が低い周波数のピークよりも縮小幅が大きくなるようにスケールに応じて異なる重みが生成されて適用されます。

多くの用途では、L1 正規化の方が適しています。L1 ノルムの定義には値の平方化はないため、維持する係数は $$\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{s}$}\right.$$ ではなく 1/s です。L1 正規化では、L2 のように高周波数の振幅が縮小されるのではなく、すべての周波数の振幅が同じ値に正規化されます。そのため、L1 ノルムを使用した方が信号の表現の精度が高くなります。2 つの複素指数の連続ウェーブレット変換の例を参照してください。

参照

[1] Lilly, J. M., and S. C. Olhede. “Generalized Morse Wavelets as a Superfamily of Analytic Wavelets.” IEEE Transactions on Signal Processing 60, no. 11 (November 2012): 6036–6041. https://doi.org/10.1109/TSP.2012.2210890.

[2] Lilly, J.M., and S.C. Olhede. “Higher-Order Properties of Analytic Wavelets.” IEEE Transactions on Signal Processing 57, no. 1 (January 2009): 146–160. https://doi.org/10.1109/TSP.2008.2007607.

[3] Lilly, J. M. jLab: A data analysis package for MATLAB^®, version 1.6.2. 2016. http://www.jmlilly.net/jmlsoft.html.

[4] Lilly, Jonathan M. “Element Analysis: A Wavelet-Based Method for Analysing Time-Localized Events in Noisy Time Series.” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473, no. 2200 (April 30, 2017): 20160776. https://doi.org/10.1098/rspa.2016.0776.

拡張機能

GPU コード生成

GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

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

• 単精度および倍精度の入力信号がサポートされています。精度はコンパイル時に設定しなければなりません。

• timetable の入力信号はサポートされていません。

• 解析的な Morse ('morse') ウェーブレットと Morlet ('amor') ウェーブレットのみがサポートされます。

• 次の入力引数はサポートされていません。サンプリング周期 (ts)、および名前と値のペア PeriodLimits、NumOctave、FilterBank。

• スケーリング係数の出力とフィルター バンクの出力はサポートされていません。

• プロットはサポートされていません。

GPU 配列

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

この関数は、GPU 配列を完全にサポートします。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2016b で導入

R2018a: 'NumOctaves' の名前と値の引数は削除予定

NumOctaves の名前と値の引数は、将来のリリースで削除される予定です。次のいずれかを使用してください。

• CWT の周波数範囲を変更する場合は名前と値の引数 FrequencyLimits。

• CWT の周期範囲を変更する場合は名前と値の引数 PeriodLimits。

詳細については、cwtfreqbounds を参照してください。

R2016b: 関数 cwt の構文は変更済み

このリリースでは、連続ウェーブレット変換 cwt の更新バージョンが提供されます。単純化された新しい構文を使用すると、連続ウェーブレット解析、周波数や周期の範囲、およびオクターブあたりの音に対して、最適なウェーブレットを簡単に選択できます。ウェーブレットとスケーリングの既定値が提供されるため、指定する必要はありません。

古い関数 cwt 用の構文は引き続き機能しますが、非推奨になりました。更新されたバージョンの cwt を使用してください。古いバージョンと更新されたバージョンの両方で同じ関数名が使用されています。関数への入力によって、使用されるバージョンが自動的に決定されます。

機能	代替方法	互換性の考慮事項
古い cwt	更新された cwt	更新された cwt 構文を使用するように、cwt のすべてのインスタンスを更新します。

参考

アプリ

• ウェーブレット時間-周波数アナライザー

関数

• cwtfilterbank | icwt | cwtfreqbounds

トピック

• Practical Introduction to Time-Frequency Analysis Using the Continuous Wavelet Transform
• Using Wavelet Time-Frequency Analyzer App
• 連続および離散ウェーブレット変換
• CWT による時間-周波数解析
• Boundary Effects and the Cone of Influence
• Morse ウェーブレット
• 時間-周波数ギャラリー