このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

cwt

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

旧バージョンの cwt の詳細については、cwt を参照してください。旧バージョンは非推奨になりました。

構文

wt = cwt(x)

wt = cwt(x,wname)

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

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

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

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

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

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

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

cwt(___)

説明

wt = cwt(x) は、x の連続ウェーブレット変換 (CWT) を返します。入力 x は実数値または複素数値のベクトルか単一変数の一定間隔でサンプリングされる timetable で、少なくとも 4 つのサンプルを含まなければなりません。CWT は、対称性パラメーター (ガンマ) が 3 に等しく時間-帯域積が 60 に等しい解析 Morse ウェーブレットを使用して得られます。cwt で使用されるオクターブあたりの音の数は 10 です。最小スケールと最大スケールは、ウェーブレットの周波数と時間におけるエネルギーの広がりに基づいて自動的に決まります。x が実数値の場合、wt は各行が 1 つのスケールに対応する 2 次元行列です。wt の列サイズは x の長さと等しくなります。x が複素数値の場合、wt は 3 次元行列です。1 ページ目が正のスケール (解析的な部分または反時計回りの成分) の CWT で、2 ページ目が負のスケール (反解析的な部分または時計回りの成分) の CWT になります。

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

wt = cwt(x,wname) は、wname で指定された解析ウェーブレットを使用して CWT を計算します。wname の有効なオプションは、'morse'、'amor'、および 'bump' で、それぞれ Morse、Morlet (ガボール)、Bump のウェーブレットを示します。wname を指定しない場合、wname は既定で 'morse' になります。

[wt,f] = cwt(___,fs) では、サンプリング周波数 fs を Hz 単位で正のスカラーとして指定します。cwt は、fs を使用してスケールから周波数への変換を決定し、周波数 f を Hz 単位で返します。サンプリング周波数を指定しない場合、cwt は f をサンプルあたりのサイクル数で返します。入力 x が複素数の場合、wt の両方のページにスケールから周波数への変換が適用されます。x が timetable の場合は fs は指定できません。fs は timetable の RowTimes から決定されます。

[wt,period] = cwt(___,ts) では、サンプリング周期 ts を正の duration スカラーとして指定します。duration は、年、日、時間、分、または秒で指定できます。cwt は、ts を使用してスケールから周期への変換を計算し、期間を period で返します。period における duration の配列の書式プロパティは ts と同じになります。入力 x が複素数の場合、wt の両方のページにスケールから周期への変換が適用されます。x が timetable の場合は ts は指定できません。ts は、'PeriodLimits' の名前と値のペアを設定したときに timetable の RowTimes から決定されます。

[wt,f,coi] = cwt(___,fs) は、CWT のエッジの影響が顕著になる円錐状影響圏 coi を返します。CWT の円錐状影響圏の単位は Hz です。入力 x が複素数の場合、wt の両方のページに円錐状影響圏が適用されます。

[wt,period,coi] = cwt(___,ts) は、CWT のエッジの影響が顕著になる円錐状影響圏 coi を返します。CWT の円錐状影響圏の単位は周期です。入力 x が複素数の場合、wt の両方のページに円錐状影響圏が適用されます。

[___] = cwt(___,Name,Value) は、1 つ以上のペアの引数 Name,Value で指定された追加オプションを使用して CWT を返します。

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

[___,fb,scalingcfs] = cwt(___) は、解析ウェーブレットが 'morse' または 'amor' の場合にスケーリング係数を返します。スケーリング係数は Bump ウェーブレットではサポートされていません。

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

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

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

メモ

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

例

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

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

load mtlb;
w = cwt(mtlb);

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

既定の Morse ウェーブレットの代わりに Bump ウェーブレットを使用して、音声サンプルの連続ウェーブレット変換を求めます。

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

既定の Morse ウェーブレットを使用して、CWT から得られた結果を比較します。

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('mins')
ylabel('nm/s^2')
grid on
title('Kobe Earthquake Data')

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]);
figure
freqz(fb)

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

figure
cwt(x,'FilterBank',fb)

このスクリプトを実行し、データ カーソルを使用してウェーブレット係数の振幅が正弦波成分の振幅と基本的に等しいことを確認します。

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

この例では、複数の時系列の CWT を行う場合の計算効率を CWT フィルター バンクを使用して高める方法を示します。

1995 年の阪神淡路大震災の発生時に記録された地震計データを読み込みます。このデータは、オーストラリアのホバートにあるタスマニア大学の地震計で 1995 年 1 月 16 日 20:56:51 (GMT) から 51 分間にわたって 1 秒間隔で記録された測定値 (垂直加速度、nm/sq.sec) です。データに適用できる CWT フィルター バンクを作成します。

load kobe
fb = cwtfilterbank('SignalLength',numel(kobe),'SamplingFrequency',1);

関数 cwt を使用してデータの CWT を 250 回実行します。これに要した経過時間を表示します。

num = 250;
tic;
for k=1:num
    cfs = cwt(kobe);
end
toc

Elapsed time is 6.551628 seconds.

次に、フィルター バンクのオブジェクト関数 wt を使用してデータの CWT を実行します。フィルター バンクを使用した方が速いことを確認します。

tic;
for k=1:num
    cfs = wt(fb,kobe);
end
toc

Elapsed time is 3.782376 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 値を指定します。

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

Code generation successful: To view the report, open('codegen/mex/cwt/html/report.mldatx').

-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('Seconds')
ylabel('Hz')
title('CWT of Two-Tone Signal')

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

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

ans =

  single

  6.7583e-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
ylabel('Frequency')
title('Scalogram Scaled By Level')
colormap(pink(240))

時間-帯域積の変更

この例では、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))
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))
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 の計算に使用する解析ウェーブレット。'morse'、'amor'、または 'bump' として指定します。これらの文字ベクトルは、それぞれ Morse、Morlet (ガボール)、および Bump の解析ウェーブレットを示します。

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

fs — サンプリング周波数

正のスカラー

サンプリング周波数 (Hz 単位)。正のスカラーとして指定します。fs を指定する場合、ts は指定できません。

データ型: single | double

ts — サンプリング周期

スカラーの duration

サンプリング周期 (期間とも呼ばれます)。スカラーの duration として指定します。有効な duration は、years、days、hours、minutes、および seconds です。カレンダー期間は使用できません。ts を指定する場合、fs は指定できません。

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

データ型: duration

名前と値の引数

例: 'ExtendSignal',false は、信号が拡張されないことを示します。

オプションの引数 Name,Value のコンマ区切りペアを指定します。Name は引数名で、Value は対応する値です。Name は引用符で囲まなければなりません。Name1,Value1,...,NameN,ValueN のように、複数の名前と値のペアの引数を、任意の順番で指定できます。

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

true (既定値) | false

入力信号を反射で対称的に拡張するオプション。'ExtendSignal' と true または false のいずれかで構成されるコンマ区切りのペアとして指定します。信号を対称的に拡張すると境界の影響を軽減できます。true を指定した場合、信号が拡張されます。false を指定した場合、信号は拡張されません。

FrequencyLimits — 周波数の範囲

2 要素スカラー ベクトル

CWT で使用する周波数の範囲。厳密に増加する正のエントリをもつ 2 要素ベクトルとして指定します。最初の要素は最も低いピーク通過帯域周波数を示し、ウェーブレットのピーク周波数 (Hz 単位) と 2 つの時間の標準偏差の積を信号長で除算した値以上でなければなりません。2 番目の要素は最も高いピーク通過帯域周波数を示し、ナイキスト周波数以下でなければなりません。周波数の下限 freqMin に対する周波数の上限 freqMax の比の基底 2 の対数は 1/NV 以上でなければなりません。ここで NV はオクターブあたりの音の数です。

log2(freqMax/freqMin) ≥ 1/NV.

許容される範囲外の周波数範囲を指定した場合、cwt では有効な最小値と最大値まで範囲が切り詰められます。CWT の各種のパラメーター表現の周波数範囲を特定するには、cwtfreqbounds を使用します。複素数値の信号の場合、反解析的な部分には (-1)*flimits が使用されます。ここで flimits は FrequencyLimits で指定されたベクトルです。

例: 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.

許容される範囲外の周期範囲を指定した場合、cwt では有効な最小値と最大値まで範囲が切り詰められます。ウェーブレット変換の各種のパラメーター表現の周期範囲を求めるには、cwtfreqbounds を使用します。複素数値の信号の場合、反解析的な部分には (-1)*plimits が使用されます。ここで plimits は PeriodLimits で指定されたベクトルです。

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

データ型: duration

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

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

CWT に使用するオクターブあたりの音の数。'VoicesPerOctave' と 4 から 48 までの偶数の整数で構成されるコンマ区切りのペアとして指定します。指定したオクターブあたりの音の数を使用して CWT のスケールが離散化されます。ウェーブレットの周波数と時間のエネルギーの広がりは、最小スケールと最大スケールによって自動的に決まります。

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

60 (既定値) | 3 より大きく 120 以下のスカラー

Morse ウェーブレットの時間-帯域積。'TimeBandwidth' と 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 ウェーブレットの対称性と時間-帯域積。'WaveletParameters' とスカラーの 2 要素ベクトルで構成されるコンマ区切りのペアとして指定します。最初の要素は対称性 $$\gamma$$ で、1 以上でなければなりません。2 番目の要素は時間-帯域積で、$$\gamma$$ よりも厳密に大きくなければなりません。$$\gamma$$ に対する時間-帯域積の比が 40 を超えることはできません。

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

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

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

NumOctaves — オクターブの数

正の整数

オクターブの数。'NumOctaves' と正の整数で構成されるコンマ区切りのペアとして指定します。オクターブの数が log2(fmax/fmin) を超えることはできません。ここで fmax と fmin は CWT の周波数 (または周期) の最大値と最小値です。これらは、信号長、サンプリング周波数、およびウェーブレットで決まります。詳細は、cwtfreqbounds を参照してください。

名前と値のペア 'NumOctaves' は非推奨であり、将来のリリースで削除される予定です。CWT の周波数または周期の範囲を変更するときは、名前と値のペア 'FrequencyLimits' または 'PeriodLimits' を使用することを推奨します。'NumOctaves' と 'FrequencyLimits' または 'PeriodLimits' の両方の名前と値のペアを指定することはできません。

FilterBank — CWT フィルター バンク

CWT フィルター バンク オブジェクト

CWT の計算に使用する CWT フィルター バンク。CWT フィルター バンク オブジェクトとして指定します。名前と値のペア 'FilterBank' を使用する場合、他のオプションは指定できません。CWT の計算についてのすべてのオプションが fb のプロパティで定義されます。

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

例: wt = cwt(x,'FilterBank',fb)

出力引数

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 の単位は Hz です。fs を指定しない場合、cwt は f をサンプルあたりのサイクル数で返します。

period — 周期

配列

周期。duration の配列として返されます。duration の形式は ts と同じになります。各行が周期に対応しています。

coi — 円錐状影響圏

実数の配列 | duration の配列

CWT の円錐状影響圏。実数の配列または duration の配列のいずれかとして返されます。円錐状影響圏は、CWT でエッジの影響が現れる箇所を示します。サンプリング周波数 fs を指定した場合、円錐状影響圏の単位は Hz です。スカラーの duration ts を指定した場合、円錐状影響圏の単位は周期です。エッジの影響により、円錐状影響圏の外側やそれと重なる領域は信頼性が低くなります。

詳細は、Boundary Effects and the Cone of Influenceを参照してください。

fb — CWT フィルター バンク

CWT フィルター バンク オブジェクト

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

scalingcfs — スケーリング係数

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

CWT のスケーリング係数。解析ウェーブレットが 'morse' または 'amor' の場合に実数値または複素数値のベクトルとして返されます。scalingcfs の長さは、入力 x の長さと等しくなります。

スケーリング係数は Bump ウェーブレットではサポートされていません。

詳細

解析ウェーブレット

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

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

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

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

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

• Bump

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

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

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

ヒント

• 複数の 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' は削除予定

R2018a 以降は非推奨

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

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

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

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

参考

cwtfilterbank | icwt | cwtfreqbounds

トピック

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

×

例を開く

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

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