# fsst

Fourier synchrosqueezed transform

## Syntax

## Description

returns the Fourier Synchrosqueezed Transform of the input signal,
`s`

= fsst(`x`

)`x`

. Each column of `s`

contains the
synchrosqueezed spectrum of a windowed segment of `x`

.

## Examples

### Fourier Synchrosqueezed Transform of Sinusoidal Signal

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

### Synchrosqueezed Transform of Linear Chirps

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`

:

A 256-point Kaiser window with

*β*= 10 to window the signalAn 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)

### Fourier Synchrosqueezed Transform of Chirps

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

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

$${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 )}\phantom{\rule{0.16666666666666666em}{0ex}}dx={e}^{-\pi (\eta -\nu {)}^{2}}{e}^{j2\pi \nu t}.$$

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,

$${\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 ,$$

equals the value obtained by using the standard definition,

$${\nu}_{inst}=\frac{1}{2\pi}\frac{d}{dt}\mathrm{arg}f(t)=\frac{1}{2\pi}\frac{d}{dt}2\pi \nu t.$$

For a superposition of sinusoids,

$$f(x)=\sum _{k=1}^{K}{A}_{k}{e}^{j2\pi {\nu}_{k}x},$$

the short-time transform becomes

$${V}_{g}f(t,\eta )=\sum _{k=1}^{K}{A}_{k}{e}^{-\pi (\eta -{\nu}_{k}{)}^{2}}{e}^{j2\pi {\nu}_{k}t}.$$

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

$$\Delta ={\displaystyle \frac{1}{{\displaystyle \sigma}}}\sqrt{2\mathrm{log}2}$$

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

### Fourier Synchrosqueezed Transform of Speech Signal

Load a speech signal sampled at $${F}_{s}=7418\phantom{\rule{0.2777777777777778em}{0ex}}Hz$$. 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)`

## Input Arguments

`x`

— Input signal

vector

Input signal, specified as a vector.

**Example: **`cos(pi/4*(0:159))+randn(1,160)`

specifies
a sinusoid embedded in white Gaussian noise.

**Data Types: **`single`

| `double`

**Complex Number Support: **Yes

`fs`

— Sample rate

1 Hz (default) | positive scalar

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.

**Data Types: **`double`

| `single`

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

If

`window`

is an integer, then`fsst`

divides`x`

into 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.If

`window`

is a vector, then`fsst`

divides`x`

into segments of the same length as the vector and windows each segment using`window`

. The overlap between adjacent segments is`length(window)`

– 1.If

`window`

is not specified, then`fsst`

divides`x`

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

**Example: **`hann(N+1)`

and `(1-cos(2*pi*(0:N)'/N))/2`

both
specify a Hann window of length `N`

+ 1.

**Data Types: **`double`

| `single`

`freqloc`

— Frequency display axis

`'xaxis'`

(default) | `'yaxis'`

Frequency display axis, specified as `'xaxis'`

or `'yaxis'`

.

`'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 `fsst`

with
output arguments.

## Output Arguments

`s`

— Fourier synchrosqueezed transform

matrix

Fourier synchrosqueezed transform, returned as a matrix. The number of columns of
`s`

equals the length of `x`

. Time
increases across the columns of `s`

and frequency
increases down the rows of `s`

, starting from zero. If
`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

vector

Normalized frequencies, returned as a vector. The length
of `w`

equals the number of rows in `s`

.

`f`

— Cyclical frequencies

vector

Cyclical frequencies, returned as a vector. The length of `f`

equals
the number of rows in `s`

.

## More About

### Fourier Synchrosqueezed Transform

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

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

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 `spectrogram`

function. 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 `fsst`

function determines the
short-time Fourier transform of a function, *f*, using a spectral
window, *g*, and computing

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

Unlike the conventional definition, this definition has an extra
factor of *e*^{j2πη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

$${T}_{g}f(t,\omega )={\displaystyle {\int}_{-\infty}^{\infty}{V}_{g}f(t,\eta )\text{\hspace{0.17em}}\delta (\omega -{\Omega}_{g}f(t,\eta ))\text{\hspace{0.17em}}d\eta},$$

where the instantaneous frequencies are estimated with the “phase transform”

$${\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 )}.$$

The transform in the denominator decreases the influence of the
window. To see a simple example, refer to Detect Closely Spaced Sinusoids with the Fourier Synchrosqueezed Transform. The
definition of *T _{g}f*(

*t*,

*ω*) differs by a factor of 1/

*g*(0) from other expressions found in the literature.

`fsst`

incorporates 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 Ω* _{g}f*(

*t*,

*η*):

$${f}_{k}(t)\approx \frac{1}{g(0)}{\displaystyle {\int}_{\left|\omega -{\Omega}_{k}(t)\right|<\epsilon}{T}_{g}f(t,\omega )\text{\hspace{0.17em}}d\omega},$$

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: $$\frac{d{\varphi}_{k}(t)}{dt}>\frac{d{A}_{k}(t)}{dt}$$ 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 $$\left|\frac{d{\varphi}_{k}(t)}{dt}-\frac{d{\varphi}_{m}(t)}{dt}\right|>2\Delta $$ for

*k*≠*m*.

For an illustration, refer to Detect Closely Spaced Sinusoids with the Fourier Synchrosqueezed Transform.

## References

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

## Extended Capabilities

### C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

The window must be double precision.

Duration arrays are not supported for code generation.

### GPU Code Generation

Generate CUDA® code for NVIDIA® GPUs using GPU Coder™.

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.

### Thread-Based Environment

Run code in the background using MATLAB® `backgroundPool`

or accelerate code with Parallel Computing Toolbox™ `ThreadPool`

.

Usage notes and limitations:

The syntax with no output arguments is not supported.

For more information, see Run MATLAB Functions in Thread-Based Environment.

### GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

This function fully supports GPU arrays. For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

## Version History

**Introduced in R2016b**

## Open Example

You have a modified version of this example. Do you want to open this example with your edits?

## MATLAB Command

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)