% Continuous wavelet transform with Morlet
close all;  clc; clear all;
load data;
N=length(acc); 
sig=acc*986; % convert acceleration unit to in/sec^2 
t=(0:N-1)'*dt;       % time vector
sig=sig(:);   t=t(:);
%%
fs=1/dt;
fmax=fs/2;
nf=floor(N/8);
freqs=linspace(1/N,fmax/fs,nf);
wo=6; sigma=2; % morlet parameters
scales  = wo./(2*pi*freqs);
w=-pi:2*pi/(N-1):pi; 
% form nf identical rows containing Fsig=fft(u)
SIGNAL= repmat(fft(sig'),nf,1); 
% shifted frequencies 
wshifted=kron(scales',ifftshift(w)); 
%FT of scaled and shifted wavelets
WAVELETS=exp(-sigma*(wshifted- wo).^2); 
%% convolution theorem 
cwt = fliplr(ifft((WAVELETS.*SIGNAL)')');
%%
A=abs(cwt);
P=angle(cwt);
dur=(N-1)*dt;  df=1/dur;
f=0:df:fmax;

figure
imagesc(t,f,A);
set(gca,'ydir','normal');
xlabel('Time (s)')  
ylabel('Frequency (Hz)')  
title('Wavelet norm')

