時間-周波数解析と連続ウェーブレット変換

この例では次を使用します。

この例では、連続ウェーブレット変換の時間-周波数分解能が可変であることによって、シャープな時間-周波数表現が得られることを示します。

連続ウェーブレット変換 (CWT) は時間-周波数変換であり、非定常信号の解析に最適です。信号が非定常であるとは、その周波数領域表現が時間とともに変化することを意味します。心電図、オーディオ信号、地震データ、気候データなど、多くの信号が非定常です。

双曲線チャープの読み込み

2 つの双曲線チャープで構成される信号を読み込みます。データは 2048 Hz でサンプリングされています。1 つ目のチャープは 0.1 ～ 0.68 秒でアクティブです。2 つ目のチャープは 0.1 ～ 0.75 秒でアクティブです。時間 $t$ における最初のチャープの瞬時周波数 (Hz) は、 $\frac{15 π}{{(0.8 - t)}^{2}} / 2 π$ です。時間 $t$ における 2 番目のチャープの瞬時周波数は、 $\frac{5 π}{{(0.8 - t)}^{2}} / 2 π$ です。信号をプロットします。

load hychirp
plot(t,hychirp)
grid on
title("Signal")
axis tight
xlabel("Time (s)")
ylabel('Amplitude')

時間-周波数解析: フーリエ変換

フーリエ変換 (FT) は、信号内に存在する周波数成分の特定に非常に適しています。しかし、FT はその周波数成分が発生する時間を特定することはできません。

信号の振幅スペクトルをプロットします。0 ～ 200 Hz の領域を拡大します。

sigLen = numel(hychirp);
fchirp = fft(hychirp);
fr = Fs*(0:1/Fs:1-1/Fs);
plot(fr(1:sigLen/2),abs(fchirp(1:sigLen/2)),"x-")
xlabel("Frequency (Hz)")
ylabel("Amplitude")
axis tight
grid on
xlim([0 200])

時間-周波数解析: 短時間フーリエ変換

フーリエ変換は時間情報を提供しません。周波数の変化が発生する時間を特定するため、短時間フーリエ変換 (STFT) 手法では、信号を異なるチャンクに分割し、各チャンクに対して FT を実行します。時間-周波数平面の STFT タイリングを以下に示します。

STFT は、信号イベントが発生するタイミングと周波数の両方に関する情報を提供します。しかし、ウィンドウ (セグメント) サイズの選択が重要です。STFT を使用する時間-周波数解析では、小さいウィンドウサイズを選択すると、時間分解能が向上しますが、周波数分解能は低下します。逆に、大きいウィンドウを選択すると、周波数分解能は向上しますが、時間分解能は低下します。

ウィンドウサイズを選択すると、そのサイズは解析全体にわたって固定されます。信号内に存在すると思われる周波数成分を推定できる場合は、その情報を使用して解析用のウィンドウサイズを選択できます。

最初の時点での 2 つのチャープの瞬時周波数は、約 5 Hz と 15 Hz です。補助関数 helperPlotSpectrogram を使用して、信号のスペクトログラムを 200 ミリ秒の時間ウィンドウサイズでプロットします。helperPlotSpectrogram のソースコードの一覧を付録に示します。この補助関数は、スペクトログラムの上に瞬時周波数を黒の破線の線分としてプロットします。瞬時周波数は、信号の前半では分解されていますが、後半では十分に分解されていません。

helperPlotSpectrogram(hychirp,t,Fs,200)

次に、helperPlotSpectrogram を使用して、スペクトログラムを 50 ミリ秒の時間ウィンドウサイズでプロットします。信号の後半に発生する高い周波数が分解されるようになりましたが、信号の最初にある低い周波数が分解されていません。

helperPlotSpectrogram(hychirp,t,Fs,50)

双曲線チャープのような非定常信号の場合、STFT の使用は問題となります。単一のウィンドウサイズでは、そのような信号の周波数成分全体を分解することはできません。

時間-周波数解析: 連続ウェーブレット変換

連続ウェーブレット変換 (CWT) は、STFT に固有の分解能の問題を解消するために作成されました。時間-周波数平面の CWT タイリングを以下に示します。

平面の CWT タイリングが役に立つ理由は、実際の信号の多くに含まれるゆっくりと振動する成分は長いスケールで発生し、高い周波数イベントは唐突または過渡的という傾向があるからです。ただし、もし高周波数イベントの持続時間が長いほうが自然であるなら、CWT の使用は適切ではなくなります。周波数分解能が低下し、時間分解能も得られません。しかし、ほとんどの場合、これには該当しません。人間の聴覚系もこのように機能します。つまり、低い周波数ほど周波数の位置推定に優れ、高い周波数では時間の位置推定に優れています。

CWT のスカログラムをプロットします。スカログラムは、時間と周波数の関数としてプロットされた CWT の絶対値です。CWT の周波数は対数であるため、プロットでは対数の周波数軸が使用されます。信号内に 2 つの双曲線チャープが存在することは、スカログラムから明らかです。CWT により、信号の持続時間全体にわたって瞬時周波数を正確に推定でき、セグメントの長さを選択する懸念もありません。

cwt(hychirp,Fs)

白の破線は、"円錐状影響圏" として知られるものを示します。円錐状影響圏は、境界の影響を受ける可能性があるスカログラムの領域を示します。詳細については、Boundary Effects and the Cone of Influenceを参照してください。

ウェーブレット係数の大きさが増加する速度を把握するには、補助関数 helperPlotScalogram3d を使用してスカログラムを 3 次元表面としてプロットします。helperPlotScalogram3d のソースコードの一覧を付録に示します。

helperPlotScalogram3d(hychirp,Fs)

補助関数 helperPlotScalogram を使用して、信号のスカログラムと瞬時周波数をプロットします。helperPlotScalogram のソースコードの一覧を付録に示します。瞬時周波数は、スカログラムの特徴とよく一致しています。

helperPlotScalogram(hychirp,Fs)

付録 – 補助関数

helperPlotSpectrogram

function helperPlotSpectrogram(sig,t,Fs,timeres)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.

[px,fx,tx] = pspectrum(sig,Fs,"spectrogram",TimeResolution=timeres/1000);
hp = pcolor(tx,fx,20*log10(abs(px)));
hp.EdgeAlpha = 0;
ylims = hp.Parent.YLim;
yticks = hp.Parent.YTick;
cl = colorbar;
cl.Label.String = "Power (dB)";
axis tight
hold on
title("Time Resolution: "+num2str(timeres)+" ms")
xlabel("Time (s)")
ylabel("Hz");
dt  = 1/Fs;
idxbegin = round(0.1/dt);
idxend1 = round(0.68/dt);
idxend2 = round(0.75/dt);
instfreq1 = abs((15*pi)./(0.8-t).^2)./(2*pi);
instfreq2 = abs((5*pi)./(0.8-t).^2)./(2*pi);
plot(t(idxbegin:idxend1),(instfreq1(idxbegin:idxend1)),'k--');
hold on;
plot(t(idxbegin:idxend2),(instfreq2(idxbegin:idxend2)),'k--');
ylim(ylims);
hp.Parent.YTick = yticks;
hp.Parent.YTickLabels = yticks;
hold off
end

helperPlotScalogram

function helperPlotScalogram(sig,Fs)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.
[cfs,f] = cwt(sig,Fs);

sigLen = numel(sig);
t = (0:sigLen-1)/Fs;

hp = pcolor(t,log2(f),abs(cfs));
hp.EdgeAlpha = 0;
ylims = hp.Parent.YLim;
yticks = hp.Parent.YTick;
cl = colorbar;
cl.Label.String = "Magnitude";
axis tight
hold on
title("Scalogram and Instantaneous Frequencies")
xlabel("Seconds");
ylabel("Hz");
dt  = 1/2048;
idxbegin = round(0.1/dt);
idxend1 = round(0.68/dt);
idxend2 = round(0.75/dt);
instfreq1 = abs((15*pi)./(0.8-t).^2)./(2*pi);
instfreq2 = abs((5*pi)./(0.8-t).^2)./(2*pi);
plot(t(idxbegin:idxend1),log2(instfreq1(idxbegin:idxend1)),"k--");
hold on;
plot(t(idxbegin:idxend2),log2(instfreq2(idxbegin:idxend2)),"k--");
ylim(ylims);
hp.Parent.YTick = yticks;
hp.Parent.YTickLabels = 2.^yticks;
end

helperPlotScalogram3d

function helperPlotScalogram3d(sig,Fs)
% This function is only intended to support this wavelet example.
% It may change or be removed in a future release.
figure
[cfs,f] = cwt(sig,Fs);

sigLen = numel(sig);
t = (0:sigLen-1)/Fs;
surface(t,f,abs(cfs));
xlabel("Time (s)")
ylabel("Frequency (Hz)")
zlabel("Magnitude")
title("Scalogram In 3-D")
set(gca,Yscale="log")
shading interp
view([-40 30])
end