cwt
連続 1 次元ウェーブレット変換
構文
説明
wt = cwt(x)x の連続ウェーブレット変換 (CWT) を返します。CWT は、対称性パラメーターであるガンマ () が 3 に等しく時間-帯域積が 60 に等しい解析 Morse ウェーブレットを使用して得られます。cwt で使用されるオクターブあたりの音の数は 10 です。最小スケールと最大スケールは、ウェーブレットの周波数と時間におけるエネルギーの広がりに基づいて自動的に決まります。
cwt 関数は、L1 正規化を使用します。L1 正規化では、振幅の等しい振動成分がスケールの異なるデータ内にある場合、そのような振動成分は CWT で等しい大きさをもちます。L1 正規化を使用すると信号の表現の精度が高くなります。CWT の L1 ノルムおよび2 つの複素指数の連続ウェーブレット変換を参照してください。
[___, は、CWT で使用されるフィルター バンクを返します。coi,fb] = cwt(___)cwtfilterbank を参照してください。
[___, は、ウェーブレット変換のスケーリング係数を返します。fb,scalingcfs] = cwt(___)
[___] = cwt(___, は、名前と値の引数を 1 つ以上追加で指定します。たとえば、Name=Value)wt = cwt(x,TimeBandwidth=40,VoicesPerOctave=20) は、時間-帯域積を 40、オクターブあたりの音の数を 20 と指定します。
出力引数なしで cwt(___) を使用すると、現在の Figure ウィンドウまたは指定されたターゲットの親コンテナーに CWT スカログラムがプロットされます。スカログラムは、時間と周波数の関数としてプロットされた CWT の絶対値です。対数スケールで周波数がプロットされます。エッジの影響が顕著になる円錐状影響圏もプロットされます。白い破線の外側のグレーの領域がエッジの影響が顕著になる領域です。
入力信号が複素数値の場合、ターゲットの親コンテナーを指定しないと、関数は正 (反時計回り) と負 (時計回り) の成分を、現在の Figure ウィンドウ内の別々のスカログラムにプロットします。それ以外の場合、関数は正と負の成分を連結したものをターゲットの親コンテナーにプロットします。
サンプリング周波数またはサンプリング周期を指定しない場合、周波数はサンプルあたりのサイクル数としてプロットされます。サンプリング周波数を指定した場合、周波数の単位はヘルツです。サンプリング周期を指定した場合は、スカログラムが時間と周期の関数としてプロットされます。入力信号が timetable の場合は、スカログラムが時間と周波数の関数として Hz 単位でプロットされ、RowTimes が時間軸の基底として使用されます。
スカログラムの点の時間、周波数、および振幅を表示するには、Figure の座標軸のツール バーでデータ ヒントを有効にし、スカログラムで目的の点をクリックします。
メモ
cwt 関数は、スカログラムをプロットする前に現在の Figure をクリアします。スカログラムをサブプロットに表示する方法については、CWT のスカログラムのサブプロットへのプロットを参照してください。
例
既定値を使用して、音声サンプルの連続ウェーブレット変換を取得します。
load mtlb;
w = cwt(mtlb);mtlb ファイルを読み込みます。このファイルには、音声サンプル mtlb とサンプル レート Fs が含まれています。
load mtlb解析 Morlet ウェーブレットを使用して得られた音声サンプルのスカログラムを表示します。
cwt(mtlb,"amor",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))
周波数が 32 Hz と 64 Hz の振幅が異なる 2 つの複素指数を作成します。データは 1000 Hz で 1 秒間サンプリングされます。指数の時間のサポートは互いに素です。
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*randn(size(t))); z = z+wgnNoise;
Morse ウェーブレットを使用して、CWT を求めてプロットします。
cwt(z,Fs)
カラー バーの複素指数成分の大きさは、スケールが異なっていても基本的にそれらの振幅になることに注意してください。これは 正規化による直接的な結果です。これは、このスクリプトを実行し、データ カーソルで各サブプロットを調べると確認できます。
R2026a 以降
NPG2006 データセット[5]を読み込みます。複素数値で表されたこのデータは、渦に捕捉された海中フロートの軌跡です。東方向と北方向への変位をプロットします。三角形は初期位置を示します。
load npg2006 plot(npg2006.cx) hold on plot(npg2006.cx(1),"^",MarkerSize=11,Color="r", ... MarkerFaceColor="r") hold off grid on xlabel("Eastward Displacement (km)") ylabel("Northward Displacement (km)")
データの連続ウェーブレット変換を可視化します。サンプリング周期は 4 時間です。データが複素数値であるため、関数は正 (反時計回り) と負 (時計回り) の成分を別々のスカログラムにプロットします。フロートの時計回りの回転は、時計回りの回転のスカログラムで捉えられます。
cwt(npg2006.cx,hours(4))
この例では、信号の振動成分の振幅が対応するウェーブレット係数の振幅と一致することを示します。
時間のサポートが互いに素である 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)
データ カーソルを使用して、ウェーブレット係数の振幅が正弦波成分の振幅と基本的に等しいことを確認します。
R2026a 以降
この例では、さまざまな境界拡張がスカログラムにどのような影響を与えるかを示します。
時間のサポートが互いに素である 2 つの正弦波で構成される信号を生成します。正弦波の周波数はそれぞれ 100 Hz および 400 Hz です。信号は 10 kHz で 1/10 秒間サンプリングされます。
frq1 = 100; frq2 = 400; Fs = 10e3; tm = 0:1/Fs:0.1-1/Fs; len = length(tm); brk = tm(floor(len/2)+1); sig = sin(2*pi*frq1*tm).*(tm>=0 & tm<brk)+ ... sin(2*pi*frq2*tm).*(tm>=brk & tm<=tm(end)); plot(tm,sig) grid on title("Signal") xlabel("Time (s)") ylabel("Amplitude")
サポートされている境界拡張ごとに、その拡張を使用して信号の CWT を求めます。
[wtr,f] = cwt(sig,Fs,Boundary="reflection"); wtp = cwt(sig,Fs,Boundary="periodic"); wtz = cwt(sig,Fs,Boundary="zeropad");
各 CWT のスカログラムをプロットします。境界で異なる周波数が現れるため、周期的な境界拡張 ("periodic") を使用するとラップ アラウンド効果が生じます。対称境界拡張 ("reflection") を使用すると、正弦波を反射するため、アーティファクトが生じます。この信号では、ゼロ パディング ("zeropad") を使用しても境界アーティファクトが生じません。
figure(Position=[0 0 600 600]); t=tiledlayout(3,1); nexttile pcolor(tm,f,abs(wtp),EdgeColor="none") xticks([]) yscale("log") ylabel("Frequency (Hz)") title("Boundary Extension: Periodic") nexttile pcolor(tm,f,abs(wtr),EdgeColor="none") xticks([]) yscale("log") ylabel("Frequency (Hz)") title("Boundary Extension: Reflection") nexttile pcolor(tm,f,abs(wtz),EdgeColor="none") yscale("log") title("Boundary Extension: Zero Padding") xlabel("Time (s)") ylabel("Frequency (Hz)") title(t,"Scalograms")
この例では、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.
この例では、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 ウェーブレットの時間-帯域積 が大きくなるとウェーブレットの包絡線下部の振動が増えることを示します。 が大きくなると、ウェーブレットの周波数が狭まります。
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 ウェーブレットよりも包絡線の下の振動が多いことを確認します。
tiledlayout(2,1) nexttile 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") nexttile 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 のスカログラムを Figure のサブプロットにプロットする方法を示します。
音声サンプルを読み込みます。データは 7418 Hz でサンプリングされます。CWT の既定のスカログラムをプロットします。
load mtlb
cwt(mtlb,Fs)
信号の CWT および CWT のスケールから周波数への変換を取得します。
[cfs,frq] = cwt(mtlb,Fs);
関数 cwt により、時間と周波数の座標軸がスカログラムで設定されます。サンプル時間を表すベクトルを作成します。
tms = (0:numel(mtlb)-1)/Fs;
新しい Figure で、元の信号を上のサブプロットにプロットし、スカログラムを下のサブプロットにプロットします。対数スケールで周波数をプロットします。
figure tiledlayout(2,1) nexttile plot(tms,mtlb) axis tight title("Signal and Scalogram") xlabel("Time (s)") ylabel("Amplitude") nexttile surface(tms,frq,abs(cfs)) axis tight shading flat xlabel("Time (s)") ylabel("Frequency (Hz)") set(gca,"yscale","log")
R2026a 以降
指定したターゲット座標軸とパネル コンテナーに 4 つの信号の CWT スカログラムをプロットします。
ターゲット座標軸での CWT スカログラムのプロット
音声サンプルを読み込みます。サンプル レートは 7418 Hz です。阪神淡路大震災のデータを読み込みます。サンプル レートは 1 Hz です。
load mtlb FsSpeech = 7418; load kobe FsKobe = 1;
新しい Figure ウィンドウの南西隅と北東隅に 2 つの座標軸を作成します。
fig = figure(Position=[100 100 600 600]); ax1 = axes(fig,Position=[0.1 0.1 0.55 0.35]); ax2 = axes(fig,Position=[0.4 0.6 0.55 0.35]);
音声データと阪神淡路大震災のデータの CWT スカログラムを、Figure の南西軸と北東軸にそれぞれプロットします。
cwt(mtlb,FsSpeech,Parent=ax1) title(ax1,"Magnitude Scalogram - Speech") cwt(kobe,FsKobe,Parent=ax2) title(ax2,"Magnitude Scalogram - Kobe")
ターゲット UI 座標軸での CWT スカログラムのプロット
ノイズがある Doppler 信号を読み込みます。
load noisdopp新しい UI Figure ウィンドウの南西隅に座標軸を作成します。
uif = uifigure(Position=[100 100 720 540]); ax3 = uiaxes(uif,Position=[60 70 400 200]);
ノイズを含む Doppler 信号の CWT スカログラムを Figure の座標軸にプロットします。
cwt(noisdopp,Parent=ax3)
title(ax3,"Magnitude Scalogram - Doppler")
ターゲット パネル コンテナーでの CWT スカログラムのプロット
時間のサポートが互いに素である 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);
UI Figure ウィンドウの北東隅にパネル コンテナーを追加します。
p = uipanel(uif,Position=[310 330 350 200], ... Title="CWT Scalogram in Panel Container");
正弦波信号の CWT スカログラムをパネル コンテナーにプロットします。
cwt(x,Fs,Parent=p)
R2026a 以降
1 kHz で 1 秒間サンプリングされた信号を生成します。この信号は、3 つの複素数値の正弦波とホワイト ノイズで構成されます。各正弦波はそれぞれ異なる周波数をもちます。そのうち 2 つの周波数は負の値です。正弦波の時間のサポートは互いに素です。
Fs = 1e3; t = 0:1/Fs:1; z = exp(1i*2*pi*150*t).*(t>=0.1 & t<0.3) + ... exp(-1i*2*pi*300*t).*(t>=0.3 & t<0.7) + ... exp(-1i*2*pi*100*t).*(t>0.7); wgnNoise = 0.05/sqrt(2)*(randn(size(t)) + ... 1i*randn(size(t))); sig = z + wgnNoise;
ターゲットの親コンテナーを指定せずに、信号のスカログラムをプロットします。信号が複素数値であるため、関数は正と負の成分を別々のスカログラムにプロットします。関数は周波数に対して対数スケールを使用します。
cwt(sig,Fs)
新しい Figure ウィンドウに座標軸を作成します。名前と値の引数 Parent の構文を使用して、座標軸にスカログラムをプロットします。また、CWT の係数と周波数を取得します。負の成分は、プロット内の負の周波数に対応します。関数は周波数に対して線形スケールを使用します。
fig = figure; ax = axes(fig,Position=[0.1 0.1 0.5 0.6]); [cfs,f] = cwt(sig,Fs,Parent=ax);
この関数は、CWT 係数を 3 次元配列で返します。1 ページ目は正 (反時計回り) 成分に対応し、2 ページ目は負 (時計回り) 成分に対応します。
係数と周波数を使用して、新しい Figure にターゲット座標軸のプロットを再作成できます。2 つのページを並べ替えて連結します。
flpdp = flip(cfs(:,:,2)); cfs2 = cat(1,cfs(:,:,1),flpdp);
周波数が正の値から負の値に並ぶように、周波数ベクトルを並べ替えて連結します。
flpdv = flip(-f); f2 = cat(1,f,flpdv);
新しい Figure に座標軸 ax2 を追加します。surf 関数を使用して、連結した係数の絶対値を ax2 にプロットします。
ax2 = axes(figure); surf(ax2,t,f2,abs(cfs2),EdgeColor="none") view(ax2,[0,90]) axis(ax2,"tight") set(ax2,YDir="normal") xlabel(ax2,"Time (secs)") ylabel(ax2,"Frequency (Hz)") c = colorbar(ax2); c.Label.String = "Magnitude"; title(ax2,"Recreated Scalogram")
R2026a 以降
NPG2006 データ セットを読み込みます。この複素数値データは、渦に捕捉された海中フロートの軌跡です。サンプル周期は 4 時間です。東方向と北方向への変位をプロットします。三角形は初期位置を示します。
load npg2006 sig = npg2006.cx; plot(sig) hold on plot(npg2006.cx(1),"^",MarkerSize=11,Color="r", ... MarkerFaceColor=[1 0 0 ]) hold off grid on xlabel("Eastward Displacement (km)") ylabel("Northward Displacement (km)")
新しい UI Figure ウィンドウに座標軸を作成します。名前と値の引数 Parent を使用して座標軸にスカログラムをプロットします。また、周期 p と円錐状影響圏 coi を取得します。サンプル周期を時間単位で指定すると、p と coi も時間単位の duration 配列になります。p と coi は duration 配列ですが、ターゲットの親コンテナーと複素数値信号を指定して cwt を呼び出すと、関数が周波数に対してスカログラムをプロットする点に注意してください。
uif = uifigure(Position=[100 100 720 540]); axf = uiaxes(uif,Position=[15 105 575 400]); [cfs,p,coi] = cwt(sig,hours(4),Parent=axf);
スカログラムの時計回り成分と反時計回り成分を周期を基準にして可視化するには、ターゲットの親コンテナーや出力引数を指定せずに cwt を使用します。
figure cwt(sig,hours(4))
入力引数
名前と値の引数
出力引数
詳細
解析ウェーブレットは、負の周波数に対してはフーリエ変換がゼロになる複素数値のウェーブレットです。解析ウェーブレットは、CWT を使用して時間-周波数解析を行う場合に適しています。ウェーブレット係数は複素数値であるため、解析する信号の位相と振幅の情報を提供します。解析ウェーブレットは、実際の非定常信号の周波数成分が時間の関数に応じてどのように変化するかを調べるのに非常に適しています。
解析ウェーブレットは、ほぼ例外なく急減少関数に基づきます。 が解析的な時間の急減少関数である場合、そのフーリエ変換 は周波数の急減少関数であり、一部の区間 () を除いて小さくなります。直交ウェーブレットと双直交ウェーブレットは、一般に時間のサポートがコンパクトになるように設計されています。時間のサポートがコンパクトなウェーブレットは、時間が急速に減少するウェーブレットに比べると、周波数のエネルギー集中度が比較的低くなります。直交ウェーブレットと双直交ウェーブレットのほとんどは、フーリエ領域において対称ではありません。
信号の結合時間-周波数表現を求めることが目的の場合は、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 フィルター バンクの使用を参照してください。
アルゴリズム
参照
[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.
[5] Lilly, J. M., and J.-C. Gascard. “Wavelet Ridge Diagnosis of Time-Varying Elliptical Signals with Application to an Oceanic Eddy.” Nonlinear Processes in Geophysics 13, no. 5 (September 14, 2006): 467–83. https://doi.org/10.5194/npg-13-467-2006.
