close all;
echo on; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CS515 DEMO4: Cascade algorithm and the Unitary Extension Principle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Cascade plot of the D2 system

h0 = daub_filt(2) / sqrt(2);
h1 = fliplr(cos([0:3]*pi).*h0);

echo off;
close;
for j = 1:5
   [t,y] = cascade_ON(h0,j);
   plot(t(1,:),y(1,:),'.-'); 
   axis([-1,4,-.5,1.5]);
   itext = ['iter ',num2str(j)];
   title(['Cascade evaluation of D2, ',itext]);
   wait;
end

[t,y] = cascade_ON(h0,8);
plot(t(1,:),y(1,:)); 
axis([-1,4,-.5,1.5]);
title(['Cascade evaluation of D2']);
wait;

% Show wavelet system and frequency response
z0 = fftshift(fft(h0,1024));
z1 = fftshift(fft(h1,1024));
w = linspace(-pi,pi,1024);

subplot(2,2,1);
plot(t(1,:),y(1,:));
axis([-1,4,-.5,1.5]);
title('phi(t)');

subplot(2,2,2);
plot(w,abs(z0),'r');
title('|H0(w)|');
axis([-pi,pi,0,1]);
xtics = [-pi,0,pi];
set(gca,'XTick',xtics);
set(gca,'XTickLabel','-pi||pi');

subplot(2,2,3);
plot(t(2,:),y(2,:));
axis([-3,2,-1.5,2]);
title('psi(t)');

subplot(2,2,4);
plot(w,abs(z1),'r');
title('|H1(w)|');
axis([-pi,pi,0,1]);
xtics = [-pi,0,pi];
set(gca,'XTick',xtics);
set(gca,'XTickLabel','-pi||pi');

suptitle('The D2 Wavelet System');

wait;

close;
echo on;
% Cascade plot of the tight-frame from B2

[h0,h1] = RSFilt_TF(2);
h0 = h0 / sqrt(2);
h1 = h1 / sqrt(2);

echo off;
close;
for j = 1:4
   [t,y] = cascade_TF(h0,h1,j);
   plot(t(1,:),y(1,:),'.-'); 
   axis([-2,2,-.5,1.5]);
   itext = ['iter ',num2str(j)];
   title(['Cascade evaluation of B2(t+1), ',itext]);
   wait;
end

plot(t(1,:),y(1,:)); 
axis([-2,2,-.5,1.5]);
set(gca,'XTick',[-2:2]);
title(['Cascade evaluation of B2(t+1)']);
wait;

% Show wavelet system and frequency response
z0 = fftshift(fft(h0,1024));
z1 = fftshift(fft(h1(1,:),1024));
z2 = fftshift(fft(h1(2,:),1024));
w = linspace(-pi,pi,1024);

subplot(3,2,1);
plot(t(1,:),y(1,:));
axis([-2,2,-.5,1.5]);
set(gca,'XTick',[-2 0 2]);
set(gca,'XTickLabel','-2||2');
title('phi(t)');

subplot(3,2,2);
plot(w,abs(z0),'r');
title('|H0(w)|');
axis([-pi,pi,0,1]);
xtics = [-pi,0,pi];
set(gca,'XTick',xtics);
set(gca,'XTickLabel','-pi||pi');

subplot(3,2,3);
plot(t(2,:),y(2,:));
axis([-2,2,-1,1]);
set(gca,'XTick',[-2 0 2]);
set(gca,'XTickLabel','-2||2');
title('psi1(t)');

subplot(3,2,4);
plot(w,abs(z1),'r');
title('|H1(w)|');
axis([-pi,pi,0,1]);
xtics = [-pi,0,pi];
set(gca,'XTick',xtics);
set(gca,'XTickLabel','-pi||pi');

subplot(3,2,5);
plot(t(3,:),y(3,:));
axis([-2,2,-1,1.5]);
set(gca,'XTick',[-2 0 2]);
set(gca,'XTickLabel','-2||2');
title('psi2(t)');

subplot(3,2,6);
plot(w,abs(z2),'r');
title('|H2(w)|');
axis([-pi,pi,0,1]);
xtics = [-pi,0,pi];
set(gca,'XTick',xtics);
set(gca,'XTickLabel','-pi||pi');

suptitle('The Tight-Frame from B2');

wait;
close;
echo off;

y = ReadSignal('Greasy')';

n = length(y);
a = max(abs(y))*1.2;
Y = fftshift(abs(fft(y)))/n;
A = max(Y)*1.2;


% D2 "filter bank"

h0 = daub_filt(2) / sqrt(2);
h1 = fliplr(cos([0:3]*pi).*h0);

d = DownDyadHi(y,h0);
c = DownDyadLo(y,h0);

subplot(3,1,1);
plot(y,'k');
axis([1,n,-a,a]);
set(gca,'YTick',[]);
title('The word "Greasy" processed by D2');
ylabel('V0,0');

subplot(3,1,2);
plot(d);
axis([1,n/2,-a,a]);
set(gca,'YTick',[]);
ylabel('V1,1');

subplot(3,1,3);
plot(c);
axis([1,n/2,-a,a]);
set(gca,'YTick',[]);
ylabel('V1,0');

wait;

d = 2*UpDyadHi(d,h0);
c = 2*UpDyadLo(c,h0);

D = fftshift(abs(fft(d)))/n;
C = fftshift(abs(fft(c)))/n;

subplot(3,2,1);
plot(y,'k');
axis([1,length(y),-a,a]);
set(gca,'YTick',[]);
set(gca,'XTick',[]);
title('Time');
ylabel('V0,0');

subplot(3,2,2);
plot(Y,'k');
axis([1,n,0,A]);
set(gca,'YTick',[]);
set(gca,'XTick',n/2); 
set(gca,'XTickLabel','');
title('Frequency (magnitude)');

subplot(3,2,3);
plot(d);
axis([1,length(d),-a,a]);
set(gca,'YTick',[]);
set(gca,'XTick',[]);
ylabel('V1,1');

subplot(3,2,4);
plot(D);
axis([1,n,0,A]);
hold on;
plot(fftshift(abs(fft(h1,n).^2*A)),'r:');
set(gca,'YTick',[]);
set(gca,'XTick',n/2); 
set(gca,'XTickLabel','');

subplot(3,2,5);
plot(c);
axis([1,length(c),-a,a]);
set(gca,'YTick',[]);
set(gca,'XTick',[]);
ylabel('V1,0');
xlabel('Lower levels have been reconstructed to V0');

subplot(3,2,6);
plot(C);
axis([1,n,0,A]);
hold on;
plot(fftshift(abs(fft(h0,n).^2*A)),'r:');
set(gca,'YTick',[]);
set(gca,'XTick',n/2); 
set(gca,'XTickLabel','');

suptitle(['D2 Feature Detection']);

wait;

close;

[h0,h1] = RSFilt_TF(2);
h2 = h1(2,:) / sqrt(2);
h1 = h1(1,:) / sqrt(2);
h0 = h0 / sqrt(2);

d2 = DownDyad_TF(y,h2);
d1 = DownDyad_TF(y,h1);
c = DownDyad_TF(y,h0);

subplot(4,1,1);
plot(y,'k');
axis([1,n,-a,a]);
set(gca,'YTick',[]);
title('The word "Greasy" processed by the Tight-Frame from B2');
ylabel('V0,0');

subplot(4,1,2);
plot(d2);
axis([1,n/2,-a,a]);
set(gca,'YTick',[]);
ylabel('V1,2');

subplot(4,1,3);
plot(d1);
axis([1,n/2,-a,a]);
set(gca,'YTick',[]);
ylabel('V1,1');

subplot(4,1,4);
plot(c);
axis([1,n/2,-a,a]);
set(gca,'YTick',[]);
ylabel('V1,0');

wait;

d2 = 2*UpDyad_TF(d2,h2);
d1 = 2*UpDyad_TF(d1,h1);
c = 2*UpDyad_TF(c,h0);

D2 = fftshift(abs(fft(d2)))/n;
D1 = fftshift(abs(fft(d1)))/n;
C = fftshift(abs(fft(c)))/n;

subplot(4,2,1);
plot(y,'k');
axis([1,n,-a,a]);
set(gca,'YTick',[]);
set(gca,'XTick',[]);
title('Time');
ylabel('V0,0');

subplot(4,2,2);
plot(Y,'k');
axis([1,n,0,A]);
set(gca,'YTick',[]);
set(gca,'XTick',n/2); 
set(gca,'XTickLabel','');
title('Frequency (magnitude)');

subplot(4,2,3);
plot(d2);
axis([1,n,-a,a]);
set(gca,'YTick',[]);
set(gca,'XTick',[]);
ylabel('V1,2');

subplot(4,2,4);
plot(D2);
axis([1,n,0,A]);
hold on;
plot(fftshift(abs(fft(h2,n).^2*A)),'r:');
set(gca,'YTick',[]);
set(gca,'XTick',n/2); 
set(gca,'XTickLabel','');

subplot(4,2,5);
plot(d1);
axis([1,n,-a,a]);
set(gca,'YTick',[]);
set(gca,'XTick',[]);
ylabel('V1,1');

subplot(4,2,6);
plot(D1);
axis([1,n,0,A]);
hold on;
plot(fftshift(abs(fft(h1,n).^2*A)),'r:');
set(gca,'YTick',[]);
set(gca,'XTick',n/2); 
set(gca,'XTickLabel','');

subplot(4,2,7);
plot(c);
axis([1,n,-a,a]);
set(gca,'YTick',[]);
set(gca,'XTick',[]);
ylabel('V1,0');
xlabel('Lower levels have been reconstructed to V0');

subplot(4,2,8);
plot(C);
axis([1,n,0,A]);
hold on;
plot(fftshift(abs(fft(h0,n).^2*A)),'r:');
set(gca,'YTick',[]);
set(gca,'XTick',n/2); 
set(gca,'XTickLabel','');

suptitle(['Tight-Frame Feature Detection']);

wait;

echo off;
close;

h0 = daub_filt(2) / sqrt(2);
h1 = fliplr(cos([0:3]*pi).*h0);

m = 3;
N = length(y);
subplot(m+1,2,1);
plot(y,'k');
axis([1,N,-a,a]);
set(gca,'YTick',[]);
ylabel('V0,0');
title('time');
subplot(m+1,2,2);
hold on;
H0 = abs(fft(h0,n));
H1 = abs(fft(h1,n));
plot(fftshift(H0).^2*A,'b');
plot(fftshift(H1).^2*A,'r');
plot(Y,'k');
axis([1,length(y),0,A]);
set(gca,'YTick',[]);
set(gca,'XTick',length(y)/2); 
set(gca,'XTickLabel','');
title('frequency');

wait;

hold off;

c = y;
for j = 1:m
   d = DownDyadHi(c,h0);
   c = DownDyadLo(c,h0);
   n = length(c);

   % Fit to original scale and periodize
   D = UpSample(d); 
   C = UpSample(c);
   H0 = UpSample(h0);
   H1 = UpSample(h1);
   for k=2:j
      D = UpSample(D);
      C = UpSample(C);
      H0 = UpSample(H0);
      H1 = UpSample(H1);
   end
   D = abs(fft(D))/n;
   C = abs(fft(C))/n;
   H0 = abs(fft(H0,N));
   H1 = abs(fft(H1,N));

   subplot(m+1,2,2*j-1);
   plot(d,'k');
   axis([1,n,-a,a]);
   set(gca,'YTick',[]);
   ylabel(['V',num2str(j),',1']);
   if j == 1,title('time'),end;
   subplot(m+1,2,2*j);
   plot(fftshift(D),'k');
   axis([1,length(D),0,A]);
   set(gca,'YTick',[]);
   set(gca,'XTick',[1,length(D)/2,length(D)]); 
   set(gca,'XTickLabel','-pi|0|pi');
   if j == 1,title('frequency'),end;

   subplot(m+1,2,2*j+1);
   plot(c,'k');
   axis([1,n,-a,a]);
   set(gca,'YTick',[]);
   ylabel(['V',num2str(j),',0']);
   subplot(m+1,2,2*j+2);
   hold on;
   if j < m
      plot(fftshift(H0.^2*A),'b');
      plot(fftshift(H1.^2*A),'r');
   end
   plot(fftshift(C),'k');
   axis([1,length(C),0,A]);
   set(gca,'YTick',[]);
   set(gca,'XTick',[1,length(C)/2,length(C)]); 
   set(gca,'XTickLabel','-pi|0|pi');

   wait;
   hold off;
end

suptitle('Three FWT levels of "Greasy" using D2');