Main Content

fsst

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

説明

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

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

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

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

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

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

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

すべて折りたたむ

ホワイト ガウス ノイズに含まれた正弦波の和で構成される 1024 サンプルの信号を生成します。正弦波の正規化周波数は、2π/5 ラジアン/サンプルおよび 4π/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

Figure contains an axes object. The axes object contains an object of type surface.

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

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

axis tight
view(2)
colorbar

Figure contains an axes object. The axes object contains an object of type surface.

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')

Figure contains an axes object. The axes object with title Fourier Synchrosqueezed Transform, xlabel Time (ms), ylabel Frequency (kHz) contains an object of type image.

シンクロスクイーズド変換と短時間フーリエ変換 (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)

Figure contains an axes object. The axes object with title Short-Time Fourier Transform, xlabel Time (s), ylabel Frequency (Hz) contains an object of type surface.

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')

Figure contains an axes object. The axes object with title Quadratic Chirp, xlabel Time, ylabel Frequency (MHz) contains an object of type surface.

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

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

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

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

title('Linear Chirp')

Figure contains an axes object. The axes object with title Linear Chirp, xlabel Time, ylabel Frequency (MHz) contains an object of type surface.

対数チャープのフーリエ シンクロスクイーズド変換を計算し、表示します。このチャープは 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)

Figure contains an axes object. The axes object with title Logarithmic Chirp, xlabel Time (s), ylabel Frequency (Hz) contains an object of type surface.

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

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

Figure contains an axes object. The axes object with title Logarithmic Chirp, xlabel Time (s), ylabel Frequency (Hz) contains an object of type surface.

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

ガウス ウィンドウ g(t)=e-πt2 でウィンドウ処理された正弦波 f(x)=ej2πνx を考えます。短時間変換は次のようになります。

Vgf(t,η)=ej2πνt-e-π(x-t)2e-j2π(x-t)(η-ν)dx=e-π(η-ν)2ej2πνt.

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

Ωgf(t,η)=1j2πe-π(η-ν)2tej2πνte-π(η-ν)2ej2πνt=ν,

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

νinst=12πddtargf(t)=12πddt2πνt.

正弦波の重ね合わせ

f(x)=k=1KAkej2πνkx,

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

Vgf(t,η)=k=1KAke-π(η-νk)2ej2πνkt.

2 つの正弦波から構成される信号のサンプルを 1024 個作成します。1 つの正弦波には ω0=π/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 個のサンプルのガウス ウィンドウを、α=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)")

Figure contains an axes object. The axes object with xlabel Samples, ylabel Normalized Frequency ( times pi blank rad/sample) contains an object of type surface.

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

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

fsst(x,"yaxis")

Figure contains an axes object. The axes object with title Fourier Synchrosqueezed Transform, xlabel Samples, ylabel Normalized Frequency ( times pi blank radians/sample) contains an object of type image.

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

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")

Figure contains an axes object. The axes object with xlabel Normalized Frequency ( times pi blank rad/sample) contains 3 objects of type line. These objects represent Transform, First sinusoid, Second sinusoid.

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

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

Figure contains an axes object. The axes object with title Synchrosqueezed Transform, xlabel Normalized Frequency ( times pi blank rad/sample) contains an object of type stem.

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

Δ=1σ2log2

および σ は標準偏差です。

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

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")

Figure contains an axes object. The axes object with xlabel Normalized Frequency ( times pi blank rad/sample) contains 3 objects of type line. These objects represent Transform, First sinusoid, Second sinusoid.

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

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

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

Figure contains an axes object. The axes object with title Synchrosqueezed Transform, xlabel Normalized Frequency ( times pi blank rad/sample) contains an object of type stem.

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

load mtlb

% To hear, type sound(mtlb,Fs)

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

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

Figure contains an axes object. The axes object with title Fourier Synchrosqueezed Transform, xlabel Time (ms), ylabel Frequency (kHz) contains an object of type image.

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')

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent Original, Reconstructed, Difference.

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

入力引数

すべて折りたたむ

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

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

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

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

データ型: double | single

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

データ型: duration

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

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

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

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

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

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

データ型: double | single

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

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

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

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

出力引数

すべて折りたたむ

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

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

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

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

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

詳細

すべて折りたたむ

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

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

f(t)=k=1Kfk(t)=k=1KAk(t)ej2πϕk(t).

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

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

Vgf(t,η)=f(x)g(xt)ej2πη(xt)dx.

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

Tgf(t,ω)=Vgf(t,η)δ(ωΩgf(t,η))dη,

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

Ωgf(t,η)=1j2πtVgf(t,η)Vgf(t,η)=η1j2πVg/tf(t,η)Vgf(t,η).

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

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

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

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

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

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

fk(t)1g(0)|ωΩk(t)|<εTgf(t,ω)dω,

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

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

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

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

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

参照

[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.

拡張機能

バージョン履歴

R2016b で導入