このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

フーリエ シンクロスクイーズド変換による近接した正弦波の検出

ここでは、フーリエ シンクロスクイーズド変換が鋭いリッジを生成しながらも、任意の間隔で配置された正弦波を分解できない例について示します。2 つの正弦波の近接の度合いがどの程度までであればこれらが別個のものとして検出されるかは、変換で使用するウィンドウの幅によって決まります。

ガウス ウィンドウ $$g(t)=e^{-\pi t^2}$$ でウィンドウ処理された正弦波 $$f(x)=e^{j2\pi \nu x}$$ を考えます。短時間変換は次のようになります。

$$V_{g}f(t,\eta )=e^{j2\pi \nu t}{\int }_{-\infty }^{\infty }{e}^{-\pi (x-t{)}^2}{e}^{-j2\pi (x-t)(\eta -\nu )}dx=e^{-\pi (\eta -\nu {)}^2}{e}^{j2\pi \nu t}.$$

周波数関数として表示した場合、変換は $$\nu$$ における (時間的に) 一定の振動と、そこから減衰するガウスを結合します。瞬時周波数のシンクロスクイーズ推定

$${\Omega }_{g}f(t,\eta )=\frac{1}{j2\pi }\frac{e^{-\pi (\eta -\nu {)}^2}{\displaystyle \frac{\partial }{{\displaystyle \partial t}}}{e}^{j2\pi \nu t}}{e^{-\pi (\eta -\nu {)}^2}{e}^{j2\pi \nu t}}=\nu ,$$

これは、次の標準的な定義を使用して取得される値と等しくなります。

$${\nu }_{inst}=\frac{1}{2\pi }\frac{d}{dt}\arg f(t)=\frac{1}{2\pi }\frac{d}{dt}2\pi \nu t.$$

正弦波の重ね合わせ

$$f(x)=\sum _{k=1}^K{A}_k{e}^{j2\pi {\nu }_{k}x},$$

の場合、短時間変換は次のようになります。

$$V_{g}f(t,\eta )=\sum _{k=1}^K{A}_k{e}^{-\pi (\eta -{\nu }_k{)}^2}{e}^{j2\pi {\nu }_{k}t}.$$

2 つの正弦波から構成される信号のサンプルを 1024 個作成します。1 つの正弦波には $${\omega }_0=\pi /5$$ ラジアン/サンプルの正規化周波数があります。もう 1 つの正弦波の周波数と振幅は 3 倍になります。

N = 1024;
n = 0:N-1;

w0 = pi/5;
x = exp(1j*w0*n)+3*exp(1j*3*w0*n);

関数spectrogramを使用して信号の短時間フーリエ変換を計算します。256 個のサンプルのガウス ウィンドウを、$$\alpha =20$$、隣接するセクション間にある 255 個のオーバーラップのサンプル、および 1024 個の DFT 点と共に使用します。変換の絶対値をプロットします。

Nw = 256;
nfft = 1024;
alpha = 20;

[s,w,t] = spectrogram(x,gausswin(Nw,alpha),Nw-1,nfft,"centered");

surf(t,w/pi,abs(s),EdgeColor="none")
view(2)
axis tight
xlabel("Samples")
ylabel("Normalized Frequency (\times\pi rad/sample)")

関数fsstに実装されているフーリエ シンクロスクイーズド変換の結果は、より鋭くより限定されたスペクトルの推定になります。

[ss,sw,st] = fsst(x,[],gausswin(Nw,alpha));

fsst(x,"yaxis")

正弦波は予測された周波数値での一定の振動として表示されます。リッジからの減衰がガウスであることを確認するには、変換の瞬間的な値をプロットして、ガウスの 2 つのインスタンスを重ね合わせます。$$\alpha$$ およびウィンドウの長さに関して、ガウス振幅と標準偏差を表します。周波数領域ウィンドウの標準偏差は、時間領域ウィンドウの標準偏差の逆数であることを思い出してください。

rstdev = (Nw-1)/(2*alpha);
amp = rstdev*sqrt(2*pi);

instransf = abs(s(:,128));

plot(w/pi,instransf)
hold on
plot(w/pi,[1 3]*amp.*exp(-rstdev^2/2*(w-[1 3]*w0).^2),"--")
hold off
xlabel("Normalized Frequency (\times\pi rad/sample)")
legend(["Transform" "First sinusoid" "Second sinusoid"],Location="best")

フーリエ シンクロスクイーズド変換は、推定された瞬時周波数での信号のエネルギー量を集中させます。

stem(sw/pi,abs(ss(:,128)))
xlabel("Normalized Frequency (\times\pi rad/sample)")
title("Synchrosqueezed Transform")

瞬時周波数のシンクロスクイーズド推定は、正弦波が $$2\Delta$$ より離れている場合にのみ有効です。ここで、ガウス ウィンドウの

$$\Delta ={\displaystyle \frac{1}{{\displaystyle \sigma }}}\sqrt{2\log 2}$$

および $$\sigma$$ は標準偏差です。

前の計算を繰り返しますが、今度は 2 番目の正弦波に $${\omega }_0+1.9\Delta$$ ラジアン/サンプルの正規化周波数があることを指定します。

D = sqrt(2*log(2))/rstdev;

w1 = w0+1.9*D;

x = exp(1j*w0*n)+3*exp(1j*w1*n);

[s,w,t] = spectrogram(x,gausswin(Nw,alpha),Nw-1,nfft,"centered");
instransf = abs(s(:,20));

plot(w/pi,instransf)
hold on
plot(w/pi,[1 3]*amp.*exp(-rstdev^2/2*(w-[w0 w1]).^2),"--")
hold off
xlabel("Normalized Frequency (\times\pi rad/sample)")
legend(["Transform" "First sinusoid" "Second sinusoid"],Location="best")

フーリエ シンクロスクイーズド変換は、$$|{\omega }_1-{\omega }_0|<2\Delta$$ であるため、正弦波をうまく分解できません。スペクトル推定で値が大幅に減ります。

[ss,sw,st] = fsst(x,[],gausswin(Nw,alpha));

stem(sw/pi,abs(ss(:,128)))
xlabel("Normalized Frequency (\times\pi rad/sample)")
title("Synchrosqueezed Transform")

参考

fsst | gausswin | ifsst | spectrogram

関連するトピック

• ヒルベルト変換および瞬時周波数
• 複素チャープの瞬時周波数
• 時間-周波数解析の実践的基礎