%% Spectrogram Reassignment and Thresholding
% 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 the chirp so that it is symmetric about the interval midpoint, starting and ending at a frequency of 250 Hz and attaining 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;

%
% Add white Gaussian noise to the signal. Specify a signal-to-noise ratio of 20 dB. Reset the random number generator for reproducible results.

SNR = 20;
rng('default')

sig = [x1;x2];
sig = sig + randn(size(sig))*std(sig)/db2mag(SNR);

%
% Compute and plot the spectrogram of the signal. Specify a Kaiser window of length 63 with a shape parameter \beta=17, 10 fewer samples of overlap between adjoining sections, and an FFT length of 256.

nwin = 63;
wind = kaiser(nwin,17);
nlap = nwin-10;
nfft = 256;

spectrogram(sig,wind,nlap,nfft,Fs,'yaxis')

%
% Threshold the spectrogram so that any elements with values smaller than the SNR are set to zero.

spectrogram(sig,wind,nlap,nfft,Fs,'MinThreshold',-SNR,'yaxis')

%
% Reassign each PSD estimate to the location of its center of energy.

spectrogram(sig,wind,nlap,nfft,Fs,'reassign','yaxis')

%
% Threshold the reassigned spectrogram so that any elements with values smaller than the SNR are set to zero.

spectrogram(sig,wind,nlap,nfft,Fs,'reassign','MinThreshold',-SNR,'yaxis')

% Copyright 2012 The MathWorks, Inc.