fsst

フーリエ シンクロスクイーズド変換

構文

s = fsst(x)

[s,w,n] = fsst(x)

[s,f,t] = fsst(x,fs)

[s,f,t] = fsst(x,ts)

[___] = fsst(___,window)

fsst(___)

fsst(___,freqloc)

説明

s = fsst(x) は、入力信号 x のフーリエ シンクロスクイーズド変換を返します。s の各列には、x のウィンドウが適用されたセグメントのシンクロスクイーズド スペクトルが格納されます。

[s,w,n] = fsst(x) は、正規化周波数 w のベクトルとフーリエ シンクロスクイーズド変換を計算した時点のサンプル数 n のベクトルを返します。w は s の列に、n は s の行に対応します。

[s,f,t] = fsst(x,fs) は、巡回周波数 f のベクトルおよび時点のベクトル t をサンプル レート fs で表して返します。

[s,f,t] = fsst(x,ts) は、サンプル時間 ts を duration スカラーとして指定します。t の単位は ts と同じです。f の単位は、ts の単位の逆数です。

[___] = fsst(___,window) は、window を使用して信号をセグメントに分割し、ウィンドウ処理を実行します。前の構文の入力引数を任意に組み合わせて使用し、対応する出力引数を取得することができます。

出力引数を設定せずに fsst(___) を使用すると、現在の Figure ウィンドウにシンクロスクイーズド変換がプロットされます。

fsst(___,freqloc) では、周波数をプロットする軸を指定します。

例

正弦波信号のフーリエ シンクロスクイーズド変換

ホワイト ガウス ノイズに含まれた正弦波の和で構成される 1024 サンプルの信号を生成します。正弦波の正規化周波数は、$$2\pi /5$$ ラジアン/サンプルおよび $$4\pi /5$$ ラジアン/サンプルです。高周波数の正弦波の振幅は、他の正弦波の振幅の 3 倍です。

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

w0 = 2*pi/5;
x = sin(w0*n)+3*sin(2*w0*n);

信号のフーリエ シンクロスクイーズド変換を計算します。結果をプロットします。

[s,w,n] = fsst(x);

mesh(n,w/pi,abs(s))

axis tight
view(2)
colorbar

比較のために信号の通常の短時間フーリエ変換を計算します。spectrogram の既定の値を使用します。結果をプロットします。

[s,w,n] = spectrogram(x);
 
surf(n,w/pi,abs(s),'EdgeColor','none')

axis tight
view(2)
colorbar

線形チャープのシンクロスクイーズド変換

2 つのチャープで構成された信号を生成します。信号は、3 kHz で 1 秒間サンプリングされます。最初のチャープの初期周波数は 400 Hz で、サンプリングの最後には 800 Hz に到達します。2 番目のチャープは 500 Hz から開始し、最後には 1000 Hz に到達します。2 番目のチャープの振幅は最初のチャープ 2 倍です。

fs = 3000;
t = 0:1/fs:1-1/fs;

x1 = chirp(t,400,t(end),800);
x2 = 2*chirp(t,500,t(end),1000);

信号のフーリエ シンクロスクイーズド変換を計算し、プロットします。

fsst(x1+x2,fs,'yaxis')

シンクロスクイーズド変換と短時間フーリエ変換 (STFT) を比較します。STFT を計算するには、関数 spectrogram を使用します。fsst で使用する既定のパラメーターを指定します。

• β = 10 の 256 点カイザー ウィンドウで信号にウィンドウを適用

• 隣り合ったウィンドウ セグメント間の 255 サンプルのオーバーラップ

• 256 の FFT 長

[stft,f,t] = spectrogram(x1+x2,kaiser(256,10),255,256,fs);

STFT の絶対値をプロットします。

mesh(t,f,abs(stft))

xlabel('Time (s)') 
ylabel('Frequency (Hz)')
title('Short-Time Fourier Transform')
axis tight
view(2)

チャープのフーリエ シンクロスクイーズド変換

100 Hz で開始し、t = 1 秒で 200 Hz になる二次チャープのフーリエ シンクロスクイーズド変換を計算して表示します。サンプル レートを 1 kHz に指定します。サンプル時間を duration スカラーとして表します。

fs = 1000;
t = 0:1/fs:2;
ts = duration(0,0,1/fs);

x = chirp(t,100,1,200,'quadratic');

fsst(x,ts,'yaxis')

title('Quadratic Chirp')

シンクロスクイーズ アルゴリズムは、信号の周波数が緩やかに変化するという仮定の下で動作します。したがって、周波数の変化率は小さいため、スペクトルが早い時点で良好に集中します。

DC で開始し、t = 1 秒で 150 Hz になる線形チャープのフーリエ シンクロスクイーズド変換を計算して表示します。256 サンプルのハミング ウィンドウを使用します。

x = chirp(t,0,1,150);

fsst(x,ts,hamming(256),'yaxis')

title('Linear Chirp')

対数チャープのフーリエ シンクロスクイーズド変換を計算し、表示します。このチャープは 1 kHz でサンプリングされ、20 Hz で開始し、t = 1 秒で 60 Hz に達します。β = 20 の 256 サンプル カイザー ウィンドウを使用します。

x = chirp(t,20,1,60,'logarithmic');

[s,f,t] = fsst(x,fs,kaiser(256,20));

clf
mesh(t,f,(abs(s)))

title('Logarithmic Chirp') 
xlabel('Time (s)')
ylabel('Frequency (Hz)')
view(2)

周波数軸に対数スケールを使用します。この変換は直線になります。

ax = gca;
ax.YScale = 'log';
axis tight

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

ここでは、フーリエ シンクロスクイーズド変換が鋭いリッジを生成しながらも、任意の間隔で配置された正弦波を分解できない例について示します。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")

音声信号のフーリエ シンクロスクイーズド変換

$$F_s=7418Hz$$ でサンプリングされた音声信号を読み込みます。ファイルには、"MATLAB®" という単語を発声している女性の録音音声が含まれています。

load mtlb

% To hear, type sound(mtlb,Fs)

信号のシンクロスクイーズド変換を計算します。長さ 256 のハン ウィンドウを使用します。時間が x 軸に、周波数が y 軸に表示されます。

fsst(mtlb,Fs,hann(256),'yaxis')

ifsst を使用して逆変換を行います。元の信号と再構成後の信号を比較します。

sst = fsst(mtlb,Fs,hann(256));

xrc = ifsst(sst,hann(256));

plot((0:length(mtlb)-1)/Fs,[mtlb xrc xrc-mtlb])
legend('Original','Reconstructed','Difference')

% To hear, type sound(xrc-mtlb,Fs)

入力引数

x — 入力信号

ベクトル

入力信号。ベクトルで指定します。

例: cos(pi/4*(0:159))+randn(1,160) は、ホワイト ガウス ノイズに含まれる正弦波を指定します。

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

fs — サンプル レート

1 Hz (既定値) | 正のスカラー

サンプル レート。正のスカラーで指定します。サンプル レートは単位時間あたりのサンプル数です。時間の単位が秒の場合、サンプル レートの単位は Hz です。

データ型: double | single

ts — サンプル時間

duration スカラー

サンプル時間。duration スカラーで指定します。サンプル時間は、x の連続するサンプル間で経過する時間です。

データ型: duration

window — 信号をセグメントに分割するために使用するウィンドウ

kaiser(256,10) (既定値) | 整数 | ベクトル | []

信号をセグメントに分割するために使用するウィンドウ。整数、あるいは行ベクトルまたは列ベクトルで指定します。

• window が整数の場合、fsst は x を長さが window のセグメントに分割し、各セグメントに、その長さおよび β = 10 のカイザー ウィンドウを適用します。隣接するセグメント間のオーバーラップは、window – 1 です。

• window がベクトルの場合、fsst は x をベクトルと同じ長さのセグメントに分割し、window を使用して各セグメントにウィンドウを適用します。隣接するセグメント間のオーバーラップは、length(window) – 1 です。

• window が指定されない場合、fsst は x を長さ 256 のセグメントに分割し、各セグメントに 256 サンプルのカイザー ウィンドウを β = 10 で適用します。隣接するセグメント間のオーバーラップは 255 です。x のサンプル数が 256 未満の場合、関数は、x と同じ長さで β = 10 の単一のカイザー ウィンドウを使用します。

利用可能なウィンドウのリストについては、ウィンドウを参照してください。

例: hann(N+1) と (1-cos(2*pi*(0:N)'/N))/2 は、いずれも長さ N + 1 のハン ウィンドウを指定します。

データ型: double | single

freqloc — 周波数の表示軸

'xaxis' (既定値) | 'yaxis'

周波数の表示軸。'xaxis' または 'yaxis' で指定します。

• 'xaxis' — 周波数が x 軸に、時間が y 軸に表示されます。

• 'yaxis' — 周波数が y 軸に、時間が x 軸に表示されます。

この引数は出力引数で fsst を呼び出している場合に無視されます。

出力引数

s — フーリエ シンクロスクイーズド変換

行列

フーリエ シンクロスクイーズド変換。行列として返されます。s の列数は x の長さと等しくなります。時間は s の列方向に沿って増加し、周波数は s 行方向に沿って 0 から増加します。x が実数の場合、シンクロスクイーズド スペクトルは片側になります。x が複素数の場合、シンクロスクイーズド スペクトルは両側で中央揃えになります。

w — 正規化周波数

ベクトル

ベクトルとして返される正規化周波数。w の長さは s の列数と等しくなります。

n — サンプル数

ベクトル

サンプル数。ベクトルとして返されます。n の長さは s の列数と等しくなります。n の各サンプル番号はウィンドウが適用されたセグメント x の中点になります。

f — 巡回周波数

ベクトル

ベクトルとして返される巡回周波数。f の長さは s の列数と等しくなります。

t — 時点

ベクトル

時点。ベクトルとして返されます。t の長さは s の列数と等しくなります。t の各時間の値はウィンドウが適用されたセグメント x の中点になります。

詳細

フーリエ シンクロスクイーズド変換

音声波形、機械の振動、生体信号などの多くの現実世界の信号は、振幅変調および周波数変調のモードの重ね合わせで表現できます。時間-周波数解析の場合、信号を解析信号の和として次のように表現できるため便利です。

$$f(t)={\displaystyle \sum _{k=1}^K{f}_k(t)}={\displaystyle \sum _{k=1}^K{A}_k(t)e^{j2\pi {\varphi }_k(t)}}.$$

位相 ϕ_k(t) には時間導関数 dϕ_k(t)/dt があり、これは瞬時周波数に対応します。正確な位相がわからない場合は、フーリエ シンクロスクイーズド変換を使用して推定できます。

フーリエ シンクロスクイーズド変換は、関数 spectrogram に実装された短時間フーリエ変換をベースにしています。特定の非定常信号において、シンクロスクイーズド変換では、通常の変換よりもシャープな時間-周波数推定が生成されるため、再代入されたスペクトログラムとの類似性が高くなります。関数 fsst は、スペクトル ウィンドウ g を使用し、次の計算を実行して関数 f の短時間フーリエ変換を求めます。

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

通常の定義と違い、この定義には追加の係数の e^j2πηt があります。次に、変換値は時間-周波数平面で、瞬時周波数の曲線の周囲に集中するように "押し込まれ" ます。結果のシンクロスクイーズド変換は次の形式になります。

$$T_{g}f(t,\omega )={\displaystyle {\int }_{-\infty }^{\infty }{V}_{g}f(t,\eta )\delta (\omega -{\Omega }_{g}f(t,\eta ))d\eta },$$

ここで、瞬時周波数は、次の "位相変換" によって推定されます。

$${\Omega }_{g}f(t,\eta )=\frac{1}{j2\pi }\frac{\frac{\partial }{\partial t}{V}_{g}f(t,\eta )}{V_{g}f(t,\eta )}=\eta -\frac{1}{j2\pi }\frac{V_{\partial g/\partial t}f(t,\eta )}{V_{g}f(t,\eta )}.$$

分母での変換は、ウィンドウの影響を減少させます。簡単な例については、フーリエ シンクロスクイーズド変換による近接した正弦波の検出を参照してください。T_gf(t,ω) の定義は、文献に記載されている他の式とは 1/g(0) の係数が異なっています。fsst は、モードの再構成のステップでこの係数を組み込みます。

再代入されたスペクトログラムと違い、シンクロスクイーズド変換は可逆性があり、このため信号を構成する個々のモードを再構成できます。短時間フーリエ変換の計算では、可逆性によりいくつかの制約が加わります。

• DFT 点の数は指定したウィンドウの長さに等しくなります。

• 隣り合ったウィンドウ セグメント間のオーバーラップは、ウィンドウの長さよりも 1 短くなります。

• 再割り当ては、周波数単位でのみ実行します。

モードを検出するには、Ω_gf(t,η) 付近の小さな周波数範囲でシンクロスクイーズド変換を積分します。

$$f_k(t)\approx \frac{1}{g(0)}{\displaystyle {\int }_{\left|\omega -{\Omega }_k(t)\right|<\epsilon }{T}_{g}f(t,\omega )d\omega },$$

ここで ɛ は小さい数です。

シンクロスクイーズド変換は、ウィンドウを適用した短時間フーリエ変換と比べて狭いリッジを生成します。ただし、短時間変換の幅は、依然として、シンクロスクイーズド変換のモードを分離する機能に影響します。解決可能にするには、モードが次の条件に従う必要があります。

1. それぞれのモードで、周波数は、振幅の変化率よりも厳密に大きくなる必要があります。すなわち、すべての k について $$\frac{d{\varphi }_k(t)}{dt}>\frac{d{A}_k(t)}{dt}$$ です。

2. 個々のモードは、最低でもウィンドウの周波数帯域で分離する必要があります。ウィンドウのサポートが [–Δ,Δ] の間隔の場合、k ≠ m のときに $$\left|\frac{d{\varphi }_k(t)}{dt}-\frac{d{\varphi }_m(t)}{dt}\right|>2\Delta$$ です。

詳細については、フーリエ シンクロスクイーズド変換による近接した正弦波の検出を参照してください。

参照

[1] Auger, François, Patrick Flandrin, Yu-Ting Lin, Stephen McLaughlin, Sylvain Meignen, Thomas Oberlin, and Hau-Tieng Wu. "Time-Frequency Reassignment and Synchrosqueezing: An Overview." IEEE^® Signal Processing Magazine. Vol. 30, November 2013, pp. 32–41.

[2] Oberlin, Thomas, Sylvain Meignen, and Valérie Perrier. "The Fourier-based Synchrosqueezing Transform." Proceedings of the 2014 IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), pp. 315–319.

[3] Thakur, Gaurav, and Hau-Tieng Wu. "Synchrosqueezing-based Recovery of Instantaneous Frequency from Nonuniform Samples." SIAM Journal of Mathematical Analysis. Vol. 43, 2011, pp. 2078–2095.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

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

• ウィンドウは倍精度でなければなりません。

• コード生成では、duration 配列はサポートされていません。

GPU コード生成

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

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

• ウィンドウの長さは入力信号の長さ以下でなければなりません。

• 出力引数のない構文はサポートされません。

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

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

• 出力引数のない構文はサポートされません。

詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列

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

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

バージョン履歴

R2016b で導入

参考

アプリ

• 信号アナライザー

関数

• ifsst | pspectrum | spectrogram | tfridge

トピック

• ヒルベルト変換および瞬時周波数
• 時間-周波数解析の実践的基礎
• 時間-周波数ギャラリー