%% Fourier Synchrosqueezed Transform and Its Inverse
% Generate a signal sampled at 1024 Hz for 2 seconds.

nSamp = 2048;
Fs = 1024;
t = (0:nSamp-1)'/Fs;

%
% During the first second, the signal consists of a 400 Hz sinusoid and a concave quadratic chirp. Specify a chirp that is symmetric about the interval midpoint, starts and ends at a frequency of 250 Hz, and attains a minimum of 150 Hz.

t1 = t(1:nSamp/2);

x11 = sin(2*pi*400*t1);
x12 = chirp(t1-t1(nSamp/4),150,nSamp/Fs,1750,'quadratic');
x1 = x11+x12;

%
% The rest of the signal consists of two linear chirps of decreasing frequency. One chirp has an initial frequency of 250 Hz that decreases to 100 Hz. The other chirp has an initial frequency of 400 Hz that decreases to 250 Hz.

t2 = t(nSamp/2+1:nSamp);

x21 = chirp(t2,400,nSamp/Fs,100);
x22 = chirp(t2,550,nSamp/Fs,250);
x2 = x21+x22;

%
% Compute the Fourier synchrosqueezed transform of the signal. Specify a 256-sample Kaiser window with a shape parameter β = 100. Use the plotting functionality of fsst to display the result.

sig = [x1;x2];
wind = kaiser(256,120);

[sigtr,ftr,ttr] = fsst(sig,Fs,wind);

fsst(sig,Fs,wind,'yaxis')

%
% Invert the transform to reconstruct the function. Plot the original and inverted signals and the difference between them.

x = ifsst(sigtr,wind);

plot(t,sig,t,x,t,x-sig)
legend('Original','Reconstructed','Difference')

diffnorm = norm(x-sig)

% Copyright 2012 The MathWorks, Inc.