function wavelet_demo;

clear all
close all
clc

disp('-----------------------------------------------------------------------------------')
disp('Hi, this is a demo of wavelets which explains important commands')
disp('to compute and display wavelet coefficients.')
disp('Change the parameters and play with the commands & figures to understand them more.')
disp('enjoy !!');
disp('                                  - Aniruddha J. Joshi')
disp('                                 (PhD, CSE, IIT Bombay')
disp('                                 joshi . aniruddha [at] gmail . com')
disp('-----------------------------------------------------------------------------------')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
pause

%% sine
% % t=0:.001:1; s1=cos(2*pi*10*t);
% t=0:.001:1; s=cos(2*pi*10*t)+0.5^0.5*randn(1,length(t));

%% combination -- sine + noisy sine + zero
t=0:.001:0.33; s1=cos(2*pi*10*t);
t=0:.001:0.33; s=cos(2*pi*10*t)+0.5^0.5*randn(1,length(t));
s2= 0.005*rand(339,1); s2(160)= 1; 
t=0:0.1:5; y=exp(-5*t); s2(161:180)= y(1:20); s2= s2';
s= [s1 s s2];

%% sudden spike
% s= 0.005*rand(1001,1); s(400)= 1; 
% t=0:0.1:5; y=exp(-5*t); s(401:420)= y(1:20); s= s';
 
%% issues with Fourier

disp('Issues with Fourier');
disp('-------------------');
disp('Even though the two signals are different, their Fourier transforms look');
disp('the same. The inverse transform can give only signal 1 back, as it assumes');
disp('stationarity of the signal.');
disp(' ');
disp(' ');
pause;

t=0:.001:1;
y1=cos(2*pi*10*t)+cos(2*pi*25*t)+cos(2*pi*50*t)+cos(2*pi*100*t);
y1= y1';
Fs = 1000; T = 1/Fs; L = size(y1,1); NFFT = 2^nextpow2(L); 
Y1 = fft(y1,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2);

h= figure;
set(h,'position',[50,50,1200,600])
subplot(3,5,1:2); plot(y1,'k','LineWidth',1.2); axis([1,1001,-3.2,4]);
title('cos(10)+cos(25)+cos(50)+cos(100)','FontSize',14);
subplot(3,5,[6:7 11:12]); plot(f(1:120),2*abs(Y1(1:120)),'k');%plot(f,2*abs(Y1(1:NFFT/2)));

t=0:.001:0.25;
y2= [cos(2*pi*10*t) cos(2*pi*25*t) cos(2*pi*50*t) cos(2*pi*100*t)];
y2= y2';
Fs = 1000; T = 1/Fs; L = size(y2,1); NFFT = 2^nextpow2(L); 
Y2 = fft(y2,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2);

subplot(3,5,4:5); plot(y2); axis([1,1001,-1.2,1.2]); 
title('[cos(10) cos(25) cos(50) cos(100)]','FontSize',14);
subplot(3,5,[9:10 14:15]); plot(f(1:120),2*abs(Y2(1:120)));

pause; close all

%% 

disp('Instead ...');
disp('Wavelet transform also keeps the time information.')
disp('This CWT plot shows two different WTs for the same signals.')
disp('It can be observed that the fluctuations are at different time')
disp('locations and also at different scales.');
disp(' ');
disp(' ');
pause
h= figure;
set(h,'position',[50,50,1200,600])
subplot(3,5,1:2); plot(y1,'k','LineWidth',1.2); axis([1,1001,-3.2,4]);
title('cos(10)+cos(25)+cos(50)+cos(100)','FontSize',14);
subplot(3,5,[6:7 11:12]); C1= cwt(y1,1:64,'db6','plot');

subplot(3,5,4:5); plot(y2); axis([1,1001,-1.2,1.2]); 
title('[cos(10) cos(25) cos(50) cos(100)]','FontSize',14);
subplot(3,5,[9:10 14:15]); C2= cwt(y2,1:64,'db6','plot3');

pause; close all

%%
disp('It becomes much clear in 3D CWT plot.')
disp('You can observe the transform from different angles by rotating the plot.')
disp(' ');
disp(' ');
pause
h= figure;
set(h,'position',[50,50,1200,600])
subplot(3,5,1:2); plot(y1,'k','LineWidth',1); axis([1,1001,-3.2,4]);
title('cos(10)+cos(25)+cos(50)+cos(100)','FontSize',14);
subplot(3,5,[6:7 11:12]); mesh(C1);

subplot(3,5,4:5); plot(y2); axis([1,1001,-1.2,1.2]); 
title('[cos(10) cos(25) cos(50) cos(100)]','FontSize',14);
subplot(3,5,[9:10 14:15]); mesh(C2);

pause; close all

%% plot different approximations and details using wavedec

disp('simple sine wave and wavedec + wrcoef')
disp('                     ----------------');
disp('Observe the approximations and details at 5 levels.')
disp('As the input sine wave was noise-free, the details are absent and')
disp('the approximations are the sine wave itself.');
disp(' ');
disp(' ');
pause
h= figure;
set(h,'position',[50,50,1200,600])
%t=0:.001:1; s1=cos(2*pi*10*t);
subplot(6,2,1); plot(s1,'r'); %axis([1,1001,-1.2,1.2]);
title('Orig. signal and approx. 1 to 5.');
subplot(6,2,2); plot(s1,'r'); %axis([1,1001,-1.2,1.2]);
title('Orig. signal and details 1 to 5.');
[C,L] = wavedec(s1,5,'coif3');
for i = 1:5
    A(i,:) = wrcoef('a',C,L,'coif3',i);
    D(i,:) = wrcoef('d',C,L,'coif3',i);
end
for i = 1:5
    subplot(6,2,2*i+1); plot(A(i,:),'b'); ylabel(['a' num2str(i)]); %axis([1,1001,-1.4,1.4]); 
    subplot(6,2,2*i+2); plot(D(i,:),'g'); ylabel(['d' num2str(i)]); %axis([1,1001,-1.4,1.4]); 
end
pause; close all
clear A D

%%
disp('sine wave with noise and wavedec + wrcoef')
disp('               -----                     ');
disp('Here, the approximations and details are different at different')
disp('resolutions, due to added noise.')
disp('In this case, a3 to a6 provide the average signal characteristics.')
disp('a7 and a8 are the baseline wander in the signal.');
disp('The detail coefficients captuire the high frequency components')
disp('(sometimes noise) in the signal');
disp(' ');
disp(' ');
pause
h= figure;
set(h,'position',[50,50,1200,600])
for NN= 1:8
    WAVELET= 'coif3';
    %t=0:.001:1; s=cos(2*pi*10*t)+0.5^0.5*randn(1,length(t));
    subplot(NN+1,2,1); plot(s,'r');  axis([1,1001,-1.2,1.2]);
    title(['Orig. signal+noise and approx. 1 to ' num2str(NN)]);
    subplot(NN+1,2,2); plot(s,'r');  axis([1,1001,-1.2,1.2]);
    title(['Orig. signal+noise and details 1 to ' num2str(NN)]);
    [C,L] = wavedec(s,NN,WAVELET);
    for i = 1:NN
        i
        A(i,:) = wrcoef('a',C,L,WAVELET,i); 
        D(i,:) = wrcoef('d',C,L,WAVELET,i);
    end
    for i = 1:NN
        subplot(NN+1,2,2*i+1); plot(A(i,:),'b'); ylabel(['a' num2str(i)]); %axis([1,1001,-2,2]); 
        subplot(NN+1,2,2*i+2); plot(D(i,:),'g'); ylabel(['d' num2str(i)]); %axis([1,1001,-2,2]); 
    end
    pause; 
end
pause; close all
clear A D

%% plot different approximations and details with discrete wavelet transform 

disp('sine wave with noise and three ways of computing dwt')
disp('                         ---------------------------');
disp('They are: wavedec with wrcoef, dwt and wavedec')
disp('Observe the differences at each step.')
disp('wavedec with wrcoef: does not resample.')
disp('wavedec (3rd column) collects the dwt coefficients from dwt (2nd column).')
disp('Follow the colors to understand them.')
disp(' ');
disp(' ');
pause
h1= figure; set(h1,'position',[20,50,500,600])
h2= figure; set(h2,'position',[520,50,500,600])
h3= figure; set(h3,'position',[1020,50,250,600])
COLOR= [{'r'} {'g'} {'c'} {'m'} {'y'}];

for NN= 1:8
    WAVELET= 'coif3';
    
    %% wavedec with wrcoef

    %t=0:.001:1; s=cos(2*pi*10*t)+0.5^0.5*randn(1,length(t));
    figure(h1);subplot(NN+1,2,1); plot(s,'k');  %axis([1,1001,-3.2,3.2]);
    title(['wavedec-wrcoef and approx. 1 to ' num2str(NN)]);
    figure(h1);subplot(NN+1,2,2); plot(s,'k');  %axis([1,1001,-3.2,3.2]);
    title(['wavedec-wrcoef and details 1 to ' num2str(NN)]);
    [C,L] = wavedec(s,NN,WAVELET);
    for i = 1:NN
        A(i,:) = wrcoef('a',C,L,WAVELET,i);
        D(i,:) = wrcoef('d',C,L,WAVELET,i);
    end
    for i = 1:NN
        subplot(NN+1,2,2*i+1); plot(A(i,:),COLOR{mod(i,5)+1}); ylabel(['a' num2str(i)]); %axis([1,1001,-3.2,3.2]); 
        subplot(NN+1,2,2*i+2); plot(D(i,:),COLOR{mod(i,5)+1}); ylabel(['d' num2str(i)]); %axis([1,1001,-3.2,3.2]); 
    end
    
    %% dwt

    %t=0:.001:1; s=cos(2*pi*10*t)+0.5^0.5*randn(1,length(t));
    figure(h2);subplot(NN+1,2,1); plot(s,'k');  %axis([1,1001,-3.2,3.2]);
    title(['dwt and approx. 1 to ' num2str(NN)]);
    figure(h2);subplot(NN+1,2,2); plot(s,'k');  %axis([1,1001,-3.2,3.2]);
    title(['dwt and details 1 to ' num2str(NN)]);

    if NN==1; 
        [CA{1},CD{1}] = dwt(s,WAVELET); 
    else
        [CA{NN} CD{NN}]= dwt(CA{NN-1},WAVELET);
    end
    for i = 1:NN
        subplot(NN+1,2,2*i+1); plot(CA{i},COLOR{mod(i,5)+1}); ylabel(['a' num2str(i)]); %axis([1,1001,-3.2,3.2]); 
        subplot(NN+1,2,2*i+2); plot(CD{i},COLOR{mod(i,5)+1}); ylabel(['d' num2str(i)]); %axis([1,1001,-3.2,3.2]); 
    end
    
    %% wavedec

    %t=0:.001:1; s=cos(2*pi*10*t)+0.5^0.5*randn(1,length(t));
    figure(h3);subplot(NN+1,1,1); plot(s,'k'); axis([1,1001,-3.2,3.2]);
    title('wavedec');

    for i = 1:NN
        [C L]= wavedec(s,i,WAVELET); L= cumsum([1 L]);
        figure(h3);subplot(NN+1,1,i+1); hold on;
        plot(L(1):L(2)-1,C(L(1):L(2)-1),'b');
        for j= 2:length(L)-2
            plot(L(j):L(j+1)-1,C(L(j):L(j+1)-1),COLOR{mod(i-2-j-1,5)+1}); %axis([1,1001,-3.2,3.2]); 
        end
    end
    
    pause; 
end
pause; close all
clear A D

%% cwt

disp('sine wave with noise and cwt')
disp('                         ---')
disp('CWT is pictorially easy to understand, as ALL the scales')
disp('at each time location are visible')
disp(' ');
disp(' ');
pause
%t=0:.001:1; s=cos(2*pi*10*t)+0.5^0.5*randn(1,length(t));
h1= figure; set(h1,'position',[20,50,400,600])
h2= figure; set(h2,'position',[420,50,400,600])
h3= figure; set(h3,'position',[820,50,400,600])
WAVELET= 'coif3';

figure(h1);
subplot(4,1,1); plot(s,'k');  axis([1,1001,-3.2,3.2]);
subplot(4,1,2:4); C= cwt(s,2.^[-4:0.1:4],WAVELET,'plot');
figure(h2);
subplot(4,1,1); plot(s,'k');  axis([1,1001,-3.2,3.2]);
subplot(4,1,2:4); C= cwt(s,4:64,WAVELET,'plot');
figure(h3);
subplot(4,1,1); plot(s,'k');  axis([1,1001,-3.2,3.2]);
subplot(4,1,2:4); C= cwt(s,4:4:400,WAVELET,'plot');

pause; close all

%% discrete cwt

disp('sine wave with noise and discrete cwt')
disp('                         ------------ (???)')
disp('We have chosen few scales of CWT to understand what happens at')
disp('each scale of the CWT')
disp('It can be observed that each scale captures its relative actitivies.')
disp(' ');
disp(' ');
pause

COLOR= [{'r'} {'g'} {'b'} {'c'} {'m'} {'y'} {'k'}];

%t=0:.001:1; s=cos(2*pi*10*t)+0.5^0.5*randn(1,length(t));
h1= figure; set(h1,'position',[20,50,1200,600])
% h2= figure; set(h2,'position',[420,50,400,600])
% h3= figure; set(h3,'position',[820,50,400,600])
WAVELET= 'coif3';
LEVELS= 2.^(-4:8);
C= cwt(s,LEVELS,WAVELET);
figure(h1);
subplot(4,1,1); plot(s,'k');  %axis([1,1001,-3.2,3.2]);
for i= 1:size(C,1)
    subplot(4,1,2:4); hold on; plot(LEVELS(i)+C(i,:),COLOR{mod(i,7)+1})
end

h2= figure; set(h2,'position',[20,50,600,600]);  C2= cwt(s,2.^(-4:0.1:8),WAVELET); mesh(C2);
h3= figure; set(h3,'position',[620,50,600,600]); C3= cwt(s,2.^(-4:8),WAVELET); mesh(C3);

pause; close all

%% complex wavelet

disp('complex mother wavelet')
disp('----------------------')
disp('It provides both real and imaginary coefficients.')
disp('We used both +ve and -ve information on both to detect important points');
disp(' ');
disp(' ');
pause

NLevels= [1:4];
Wavelet= 'fbsp1-1.5-1';
% CWT= cwt(Data,NLevels,Wavelet,'plot');
CWT= cwt(s,NLevels,Wavelet);
size(CWT)

h2= figure; set(h2,'position',[20,50,1200,600]);
subplot(4,1,1); plot(real(CWT(1,real(CWT(1,:))>=0))); ylabel('real positive'); axis([1,501,0,0.35]); 
subplot(4,1,2); plot(real(CWT(1,real(CWT(1,:))<0))); ylabel('real negative'); axis([1,501,-0.35,0]); 
subplot(4,1,3); plot(imag(CWT(1,imag(CWT(1,:))>=0))); ylabel('imag positive'); axis([1,501,0,1]); 
subplot(4,1,4); plot(imag(CWT(1,imag(CWT(1,:))<0))); ylabel('imag negative'); axis([1,501,-1,0]); 

pause; close all

%% reconstruction

disp('reconstruction')
disp('In above decomposition using wavedec, we assume that levels 3 to 6 are')
disp('the signals characteristics and thus, we just keep them to reconstruct')
disp('back the de-noised signal.')
disp(' ');
disp(' ');

h3= figure; set(h3,'position',[20,50,1200,600])

NN= 8;
WAVELET= 'coif3';
figure(h3);subplot(NN+1,1,1); plot(s,'k'); axis([1,1001,-3.2,3.2]);
title('wavedec');
COLOR= [{'r'} {'g'} {'c'} {'m'} {'y'}];

for i = 1:NN
    [C L]= wavedec(s,i,WAVELET); L1= L;
    L= cumsum([1 L]);
    figure(h3);subplot(NN+1,1,i+1); hold on;
    plot(L(1):L(2)-1,C(L(1):L(2)-1),'b');
    for j= 2:length(L)-2
        plot(L(j):L(j+1)-1,C(L(j):L(j+1)-1),COLOR{mod(i-2-j-1,5)+1}); %axis([1,1001,-3.2,3.2]);
    end
end

C1= zeros(size(C));
for i = 3:6
    C1(L(i):L(i+1)-1)= C(L(i):L(i+1)-1);
end
%figure; plot(C); hold on; plot(C1,'r.');
S= waverec(C1,L1,WAVELET);
figure; plot(s); hold on; plot(S,'r');

pause; close all