Fourier synchrosqueezed transform
[___] = fsst(___, uses
divide the signal into segments and perform windowing. You can use
any combination of input arguments from previous syntaxes to obtain
the corresponding output arguments.
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 rad/sample and 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.
Compare the synchrosqueezed transform with the short-time Fourier transform (STFT). Compute the STFT using the
spectrogram function. Specify the default parameters used by
A 256-point Kaiser window with β = 10 to window the signal
An overlap of 255 samples between adjoining windowed segments
An FFT length of 256
[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
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
Load a speech signal sampled at . 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.
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)
x— Input signal
Input signal, specified as a vector.
a sinusoid embedded in white Gaussian noise.
Complex Number Support: Yes
fs— Sample rate
Sample rate, specified as a positive scalar. The sample rate is the number of samples per unit time. If the unit of time is seconds, then the sample rate is in Hz.
window— Window used to divide the signal into segments
kaiser(256,10)(default) | integer | vector |
Window used to divide the signal into segments, specified as an integer or as a row or column vector.
window is an integer, then
segments of length
window and windows each segment
with a Kaiser window of that length and β = 10. The overlap
between adjacent segments is
window – 1.
window is a vector, then
segments of the same length as the vector and windows each segment
window. The overlap between adjacent segments
length(window) – 1.
window is not specified, then
segments of length 256 and windows each segment with a 256-sample
Kaiser window with β = 10. The overlap
between adjacent segments is 255. If
x has fewer
than 256 samples, then the function uses a single Kaiser window with
the same length as
x and β = 10.
For a list of available windows, see Windows.
specify a Hann window of length
N + 1.
freqloc— Frequency display axis
Frequency display axis, specified as
'xaxis' — Displays frequency
on the x-axis and time on the y-axis.
'yaxis' — Displays frequency
on the y-axis and time on the x-axis.
This argument is ignored if you call
s— Fourier synchrosqueezed transform
Fourier synchrosqueezed transform, returned as a matrix. Time
increases across the columns of
s and frequency
increases down the rows of
s, starting from zero.
x is real, then its synchrosqueezed spectrum
is one-sided. If
x is complex, then its synchrosqueezed
spectrum is two-sided and centered.
w— Normalized frequencies
Normalized frequencies, returned as a vector. The length
w equals the number of rows in
f— Cyclical frequencies
Cyclical frequencies, returned as a vector. The length of
the number of rows in
Many real-world signals such as speech waveforms, machine vibrations, and physiologic signals can be expressed as a superposition of amplitude-modulated and frequency-modulated modes. For time-frequency analysis, it is convenient to express such signals as sums of analytic signals through
The phases ϕk(t) have time derivatives dϕk(t)/dt that correspond to instantaneous frequencies. When the exact phases are unknown, you can use the Fourier synchrosqueezed transform to estimate them.
The Fourier synchrosqueezed transform is based on the short-time
Fourier transform implemented in the
For certain kinds of nonstationary signals, the synchrosqueezed transform
resembles the reassigned spectrogram because it generates sharper
time-frequency estimates than the conventional transform. The
determines the short-time Fourier transform of a function, f,
using a spectral window, g, and computing
Unlike the conventional definition, this definition has an extra factor of ej2πηt. The transform values are then “squeezed” so that they concentrate around curves of instantaneous frequency in the time-frequency plane. The resulting synchrosqueezed transform is of the form
where the instantaneous frequencies are estimated with the “phase transform”
The transform in the
denominator decreases the influence of the window. To see a simple
example, refer to Detect Closely Spaced Sinusoids. The
definition of Tgf(t,ω)
differs by a factor of 1/g(0) from
other expressions found in the literature.
the factor in the mode-reconstruction step.
Unlike the reassigned spectrogram, the synchrosqueezed transform is invertible and thus can reconstruct the individual modes that compose the signal. Invertibility imposes some constraints on the computation of the short-time Fourier transform:
The number of DFT points is equal to the length of the specified window.
The overlap between adjoining windowed segments is one less than the window length.
The reassignment is performed only in frequency.
To find the modes, integrate the synchrosqueezed transform over a small frequency interval around Ωgf(t,η):
where ɛ is a small number.
The synchrosqueezed transform produces narrow ridges compared to the windowed short-time Fourier transform. However, the width of the short-time transform still affects the ability of the synchrosqueezed transform to separate modes. To be resolvable, the modes must obey these conditions:
For each mode, the frequency must be strictly greater than the rate of change of the amplitude: for all k.
Distinct modes must be separated by at least the frequency bandwidth of the window. If the support of the window is the interval [–Δ,Δ], then for k ≠ m.
For an illustration, refer to Detect Closely Spaced Sinusoids.
 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.
 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.
 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.
Usage notes and limitations:
The window must be double precision.
Duration arrays are not supported for code generation.
Usage notes and limitations:
The length of the window must be smaller than or equal to the length of the input signal.
The syntax with no output arguments is not supported.
This function fully supports GPU arrays. For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).