Generate 1024 samples of a signal that consists of a sum of sinusoids embedded in white Gaussian noise. The normalized frequencies of the sinusoids are $$2\pi /5$$ rad/sample and $$4\pi /5$$ rad/sample. The higher frequency sinusoid has 3 times the amplitude of the other sinusoid.

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

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

Compute the Fourier synchrosqueezed transform of the signal. Plot the result.

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

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

axis tight
view(2)
colorbar

Compute the conventional short-time Fourier transform of the signal for comparison. Use the default values of spectrogram. Plot the result.

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

axis tight
view(2)
colorbar

Generate a signal that consists of two chirps. The signal is sampled at 3 kHz for one second. The first chirp has an initial frequency of 400 Hz and reaches 800 Hz at the end of the sampling. The second chirp starts at 500 Hz and reaches 1000 Hz at the end. The second chirp has twice the amplitude of the first chirp.

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

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

Compute and plot the Fourier synchrosqueezed transform of the signal.

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

Compare the synchrosqueezed transform with the short-time Fourier transform (STFT). Compute the STFT using the spectrogram function. Specify the default parameters used by fsst:

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

Plot the absolute value of the STFT.

mesh(t,f,abs(stft))

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

Compute and display the Fourier synchrosqueezed transform of a quadratic chirp that starts at 100 Hz and crosses 200 Hz at t = 1 s. Specify a sample rate of 1 kHz. Express the sample time as a duration scalar.

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

The synchrosqueezing algorithm works under the assumption that the frequency of the signal varies slowly. Thus the spectrum is better concentrated at early times, where the rate of change of frequency is smaller.

Compute and display the Fourier synchrosqueezed transform of a linear chirp that starts at DC and crosses 150 Hz at t = 1 s. Use a 256-sample Hamming window.

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

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

title('Linear Chirp')

Compute and display the Fourier synchrosqueezed transform of a logarithmic chirp. The chirp is sampled at 1 kHz, starts at 20 Hz, and crosses 60 Hz at t = 1 s. Use a 256-sample Kaiser window with β = 20.

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)

Use a logarithmic scale for the frequency axis. The transform becomes a straight line.

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

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)=e^{j2\pi \nu x}$$, windowed with a Gaussian window, $$g(t)=e^{-\pi t^2}$$. The short-time transform is

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

equals the value obtained by using the standard definition,

For a superposition of sinusoids,

the short-time transform becomes

Create 1024 samples of a signal consisting of two sinusoids. One sinusoid has a normalized frequency of $${\omega }_0=\pi /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 $$\alpha =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 $$\alpha$$ 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\Delta$$, where

for a Gaussian window and $$\sigma$$ is the standard deviation.

Repeat the previous calculation, but now specify that the second sinusoid has a normalized frequency of $${\omega }_0+1.9\Delta$$ 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 $$|{\omega }_1-{\omega }_0|<2\Delta$$. 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")

Load a speech signal sampled at $$F_s=7418Hz$$. The file contains a recording of a female voice saying the word "MATLAB®."

load mtlb

% To hear, type sound(mtlb,Fs)

Compute the synchrosqueezed transform of the signal. Use a Hann window of length 256. Display the time on the x-axis and the frequency on the y-axis.

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

Use ifsst to invert the transform. Compare the original and reconstructed signals.

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)