Main Content

Detect Closely Spaced Sinusoids with the Fourier Synchrosqueezed Transform

This example shows that the Fourier Synchrosqueezed Transform, despite the sharp ridges it produces, cannot resolve arbitrarily spaced sinusoids. The width of the window used by the transform dictates how closely spaced two sinusoids can be and still be detected as distinct.

Consider a sinusoid, f(x)=ej2πνx, windowed with a Gaussian window, g(t)=e-πt2. The short-time transform is

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

When viewed as a function of frequency, the transform combines a constant (in time) oscillation at ν with Gaussian decay away from it. The synchrosqueezing estimate of the instantaneous frequency,

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

equals the value obtained by using the standard definition,

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

For a superposition of sinusoids,

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

the short-time transform becomes

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

Create 1024 samples of a signal consisting of two sinusoids. One sinusoid has a normalized frequency of ω0=π/5 rad/sample. The other sinusoid has three times the frequency and three times the amplitude.

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

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

Use the spectrogram function to compute the short-time Fourier transform of the signal. Use a 256-sample Gaussian window with α=20, 255 samples of overlap between adjoining sections, and 1024 DFT points. Plot the absolute value of the transform.

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

The Fourier synchrosqueezed transform, implemented in the fsst function, results in a sharper, better localized estimate of the spectrum.

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

fsst(x,"yaxis")

The sinusoids are visible as constant oscillations at the expected frequency values. To see that the decay away from the ridges is Gaussian, plot an instantaneous value of the transform and overlay two instances of a Gaussian. Express the Gaussian amplitude and standard deviation in terms of α and the window length. Recall that the standard deviation of the frequency-domain window is the reciprocal of the standard deviation of the time-domain window.

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

The Fourier synchrosqueezed transform concentrates the energy content of the signal at the estimated instantaneous frequencies.

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

The synchrosqueezed estimates of instantaneous frequency are valid only if the sinusoids are separated by more than 2Δ, where

Δ=1σ2log2

for a Gaussian window and σ is the standard deviation.

Repeat the previous calculation, but now specify that the second sinusoid has a normalized frequency of ω0+1.9Δ rad/sample.

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

The Fourier synchrosqueezed transform cannot resolve the sinusoids well because |ω1-ω0|<2Δ. The spectral estimates decrease significantly in value.

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

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

See Also

| | |

Related Topics