% Example of time-frequency analysis, modified from "A better way to define
% and describe Morlet wavelets for time-frequency analysis" by Michael X
% Cohen. bioRxiv preprint doi: https://doi.org/10.1101/397182. August 2018.

%%% generate simulated signal
% The simulated signal consists of three components: a pair of weak
% Gaussian peaks at time = 500 and 1000, a strong variable-frequency sine
% wave, and random white noise. The Gaussian peaks are invisible in the raw
% signal but should appear as yellow blobs at time = 500 and 1000 at the
% bottom (low frequency) end of the segmented spectrum, clearly distinct
% from the much stronger diagonal yellow stripe from the swept sine wave.
% User-adjustable parameters
h=.15; % Height of the Gaussian peak
p1=5000;  % x-axis position of the Gaussian peak 1
p2=10000;  % x-axis position of the Gaussian peak 2
w=2500; % width of the Gaussian peak
noise=0.30;
time=1:20000;
f=.0000002.*time;
srate=6000;
z=sin(2*pi.*(max(f)-f./2).*time); % Compute swept sine wave
g=h.*gaussian(time,p1,w)+h.*gaussian(time,p2,w);
signal=g+z+noise.*randn(size(g));
pnts = length(time);
hz = linspace(0,srate/2,floor(pnts/2)+1);
powr = (2*abs(fft(signal)/pnts)).^2;
frex=30;

% call mwavelet function to create time-frequency matrix tf
[frex,tf]=mwavelet(time,signal,srate,frex);

%%% plotting
figure(1)
clf
% Plot time-frequency matrix as a contour plot
contourf(time,frex,log10(tf))
xlabel('Data points'), ylabel('Frequency (Hz)')
title('Time-frequency domain, log mode')


