close all;
echo on; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CS515 DEMO3: Introduction to Wavelet Decompositions 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Two-Scale Filters: orthogonal: haar, biorthogonal: 7/9

haar = [1 1] / sqrt(2);

f7   = [-0.06453888262876  -0.04068941760920   0.41809227322204 ...
         0.78848561640637   0.41809227322204  -0.04068941760920 ...
        -0.06453888262876 ]; 

f9   = [ 0.03782845550726  -0.02384946501956  -0.11062440441844 ...
         0.37740285561283   0.85269867900889   0.37740285561283 ...
        -0.11062440441844  -0.02384946501956   0.03782845550726 ];

% Plot the HAAR system
[t,y] = cascade_ON(haar,9);
echo off;

subplot(2,1,1); 
plot(t(1,:),y(1,:)); 
axis([-.5,1.5,0,1.5]);
title('phi(t) := the HAAR scaling function');
subplot(2,1,2); 
plot(t(2,:),y(2,:));
axis([-.5,1.5,-1.5,1.5]);
title('psi(t) := the HAAR wavelet');
wait;

% Plot the Fourier transforms of the HAAR system

% Freq domain parameters used by HAAR and 7/9
N = pow2(14);
w = linspace(-N/2,N/2,N);
ax = [-200,200,0,1];

z = abs(fft(y',N))';
scale = 1/max(z(1,:));

figure;
subplot(2,1,1); 
plot(w,scale*fftshift(z(1,:))); 
axis(ax);
set(gca,'YTick',[0,1]); 
set(gca,'XTick',[0]);
title('Frequency domain: |phi\^(w)|');
subplot(2,1,2); 
plot(w,scale*fftshift(z(2,:)));
axis(ax);
set(gca,'YTick',[0,1]);
set(gca,'XTick',[0]);
title('Frequency domain: |psi\^(w)|');
wait;

echo on;
% Plot the 7/9 biorthogonal system
[t,y] = cascade_BI(f7,f9,9);

echo off;
close all;
subplot(2,2,1); 
plot(t(1,:),y(1,:));
axis([-4,4,-.5,1.5]);
title('phi7(t) := scaling function');
subplot(2,2,2); 
plot(t(2,:),y(2,:));
axis([-3,5,-1,2]);
title('psi7(t) := wavelet');
subplot(2,2,3); 
plot(t(3,:),y(3,:));
axis([-4,4,-.5,1.5]);
title('phi9(t) := dual scaling function');
subplot(2,2,4); 
plot(t(4,:),y(4,:));
axis([-3,5,-1,2]);
title('psi9(t) := dual wavelet');
wait;

% Plot the Fourier transforms of the 7/9 system
figure;
z = abs(fft(y',N))';
scale = 1/max(z(1,:));

subplot(2,2,1); 
plot(w,scale*fftshift(z(1,:))); 
title('|phi7\^(w)|');
axis(ax);
set(gca,'YTick',[0,1]);
set(gca,'XTick',0);
subplot(2,2,2); 
plot(w,scale*fftshift(z(2,:)));
title('|psi7\^(w)|');
axis(ax);
set(gca,'YTick',[0,1]);
set(gca,'XTick',0);
subplot(2,2,3); 
scale = 1/max(z(3,:));
plot(w,scale*fftshift(z(3,:))); 
title('|phi9\^(w)|');
axis(ax); 
set(gca,'YTick',[0,1]); 
set(gca,'XTick',0);
subplot(2,2,4); 
plot(w,scale*fftshift(z(4,:)));
title('|psi9\^(w)|');
axis(ax); 
set(gca,'YTick',[0,1]); 
set(gca,'XTick',0);
wait;

echo on;
% Example 1: Piecewise Linear function
N = 1024; 
t = linspace(-1,1,N+1); t = t(1:N);
% y = spval(spmak[0:2],1),2*t+1);
y = bsp(2,2*t+1);

echo off;
close all; 
plot(t,y); 
title('Example 1: Piecewise Linear');
legend('f(t) := B2(2t+1)');
set(gca,'Xtick',[-1,-.5,0,.5,1]);
wait;
hold on;
index = find((t==-.5) + (t==0) + (t==.5));
plot(t(index),y(index),'r*');
legend('f(t) := B2(2t+1)','features');
wait;

% Foward wavelet transform
wc1 = FWT_PO(y,0,haar);
wc2 = FWT_PBS(y,0,f9,f7);

% Plot time-freq tiles
plot_tile(wc1,3,'HAAR System'); 
wait;
plot_tile(wc2,3,'7/9 System'); 
wait;

% MRA plot
close all;
plot_ON(haar,y,1,'HAAR System',1); 
plot_BI(f9,f7,y,1,'7/9 System',2); 
wait;

echo on;
% Example 2: Quadratic B-Spline
N = 1024; 
t = linspace(-1,1,N+1); t = t(1:N);
% y = spval(spmak[0:3],1),2*t+1.5);
y = bsp(3,2*t+1.5);

echo off;
close all; 
plot(t,y);
title('Example 2: Piecewise Quadratic');
legend('f(t) := B3(2t+1.5)'); 
wait;
hold on;
index = find((t==-.75) + (t==-.25) + (t==.25) + (t==.75));
plot(t(index),y(index),'r*');
legend('f(t) := B3(2t+1.5)','features'); 
wait;

% Forward wavelet transform
wc1 = FWT_PO(y,0,haar);
wc2 = FWT_PBS(y,0,f9,f7);

% Plot time-frequency tiles
plot_tile(wc1,3,'HAAR System'); 
wait;
plot_tile(wc2,3,'7/9 System');
wait;

% MRA plot
close all;
plot_ON(haar,y,1,'HAAR System',1); 
plot_BI(f9,f7,y,1,'7/9 System',2); 

