# Exemplo 6.1
N=128;
k = 0:N-1;
x = 0.5*sin(2*pi*3*k/N) + 0.5*sin(2*pi*49*k/N);
figure(1);
plot(k,x);
ylim([-1 1]);
title("Sinal original");
# resposta impulsiva do primeiro filtro de Haar (passa-baixas)
l = [ 0.5 0.5 ];
# recorta a convolução linear para manter o tamanho do vetor
xl = conv(x,l)(1:N);
figure(2);
plot(k,xl);
ylim([-1 1]);
title("Sinal filtrado (frequencias baixas)");
# resposta impulsiva do segundo filtro de Haar (passa-altas)
h = [ 0.5 -0.5 ];
xh = conv(x,h)(1:N);
figure(3);
plot(k,xh);
ylim([-1 1]);
title("Sinal filtrado (frequencias altas)");
# Observe que x = xl+xh:
printf("Observe que ||xl+xh-x|| = %g\n",norm(xh+xl-x));



# Continuação do exemplo anterior:
# construção dos vetores de aproximação de detalhe
figure(4);
plot(xl(1:2:N));
ylim([-1 1]);
title("Coeficientes de aproximacao");
figure(5)
plot(xh(1:2:N));
ylim([-1 1]);
title("Coeficientes de detalhe");



# Exemplo 6.2
close all;
N=1024;
k=0:N;
# sinal com vários trechos independentes
x = [ sin(2*pi*12*k(1:N/4)/N) 0.5*ones(1,N/4) 0.1*ones(1,N/4) -0.2*ones(1,N/4) ];
# sentinela para permitir a reconstrução perfeita
x(N+1)=x(N);
figure(1);
plot(k,x);
ylim([-1 1]);
title("Sinal original");
# constroi a componente de frequencias baixas
l = [ 0.5 0.5 ];
xl = conv(x,l);
figure(2);
plot(xl(1:2:N));
ylim([-1 1]);
title("Coeficientes de aproximacao");
# constroi a componente de frequencias altas
h = [ 0.5 -0.5 ];
xh = conv(x,h);
figure(3);
plot(xh(1:2:N));
ylim([-1 1]);
title("Coeficientes de detalhe");



# Continuação do exemplo 6.2: compressão
# Reconstrução usando transformada de Haar, zerando o vetor xh:
close all;
Uxl = zeros(1,N+1);
Uxl(1:2:N+1) = xl(1:2:N+1);
vl = conv(Uxl,[ 1 1 ])(2:N+1);
figure(1);
plot(vl);
axis([0 N -1 1]);
title("Compactacao usando a transformada de Haar");
# Calcula erro relativo da reconstrução
printf("Erro da compactação por Haar (50%%) = %g%%\n",100*norm(x(1:N)-vl)/norm(x(1:N)));
# Reconstrução usando DFT, mantendo 50% dos coeficientes (mais altos)
X = fft(x(1:N));
limiar = median(abs(X));
Xcompactado = X.*(abs(X)>limiar);
xnovo = real(ifft(Xcompactado));
figure(2);
plot(xnovo);
axis([0 N -1 1]);
title("Compactacao usando DFT");
# Calcula erro relativo da reconstrução
printf("Erro da compactação por Fourier (50%%) = %g%%\n",100*norm(x(1:N)-xnovo)/norm(x(1:N)));



# Outra tentativa: ao invés de anular o vetor xh,
# vamos manter 50% dos coeficientes de xh (os mais altos).
# A representação terá no total 75% do tamanho de x
# (50% correspondente a xl e 25% correspondente a xh).
xhnovo = xh.*(abs(xh)>=median(abs(xh)));
Uxh = zeros(1,N+1);
Uxh(1:2:N+1) = xhnovo(1:2:N+1);
vh = conv(Uxh,[ -1 1 ])(2:N+1);
figure(1);
plot(vl+vh);
axis([0 N -1 1]);
title("Compactacao usando a transformada de Haar");
# Calcula erro relativo da reconstrução
printf("Erro da compactação por Haar (75%%) = %g%%\n",100*norm(x(1:N)-vl-vh)/norm(x(1:N)));
# Reconstrução usando DFT, mantendo 75%
# dos coeficientes (mais altos)
# A expressao statistics(abs(X))(2) corresponde ao
# primeiro quartil da distribuição de abs(X).
limiar = statistics(abs(X)')(2);
Xcompactado = X.*(abs(X)>=limiar);
xnovo = real(ifft(Xcompactado));
figure(2);
plot(xnovo);
axis([0 N -1 1]);
title("Compactacao usando DFT");
# Calcula erro relativo da reconstrução
printf("Erro da compactação por Fourier (75%%) = %g%%\n",100*norm(x(1:N)-xnovo)/norm(x(1:N)));


# Exemplo 6.6: Figuras 6.10, 6.11 e 6.12
close all;
N=1024;
k=0:N-1;
# sinal do exemplo 6.2
x = [ sin(2*pi*12*k(1:N/4)/N) 0.5*ones(1,N/4) 0.1*ones(1,N/4) -0.2*ones(1,N/4) ];
figure(1);
plot(k,x);
title("Sinal original");
for j=1:4
  A=figure(3*j-1);
  y=dwt(x,j);
  plot(k(1:N/2^j),y(1:N/2^j));
  title(sprintf("Coeficientes de aproximacao de %da etapa(s)",j));
  B=figure(3*j);
  plot(k(N/2^j+1:N/2^(j-1)),y(N/2^j+1:N/2^(j-1)));
  title(sprintf("Coeficientes de detalhes de %da etapa(s)",j));
  C=figure(3*j+1);
  # zera todos os coeficientes de detalhes de todas as etapas
  y(N/2^j+1:N) *= 0;
  plot(k,idwt(y,j));
  title(sprintf("Reconstrucao a partir da aproximacao da %da etapa(s)",j));
  pause;
  close([A B C]);
endfor




# Exemplo 6.9
close all;
N=1024;
k=0:N-1;
# sinal do exemplo 6.2
x = [ sin(2*pi*12*k(1:N/4)/N) 0.5*ones(1,N/4) 0.1*ones(1,N/4) -0.2*ones(1,N/4) ];
figure(1);
plot(k,x);
title("Sinal original");
for j=1:4
  figure(j+1);
  plot(k,dwt(x,j));
  title(sprintf("DWT (Haar) de %d etapa(s)",j));
endfor



# Figuras 6.14 e 6.15
close all;
N=1024;
k=0:N-1;
# sinal do exemplo 6.2
x = [ sin(2*pi*12*k(1:N/4)/N) 0.5*ones(1,N/4) 0.1*ones(1,N/4) -0.2*ones(1,N/4) ];
for j=1:4
  M = N/(2^j);
  # Experimente também com outros bancos de filtros: Haar, Daubechies...
  #y = idwt((dwt(x,j).*[ones(1,M) zeros(1,N-M)]),j);
  #y = idwt((dwt(x,j,"Le Gall 5/3").*[ones(1,M) zeros(1,N-M)]),j,"Le Gall 5/3");
  y = idwt((dwt(x,j,"Daubechies 4-tap").*[ones(1,M) zeros(1,N-M)]),j,"Daubechies 4-tap");
  figure(j);
  plot(k,y);
  title(sprintf("Recontrucao a partir da etapa %d: erro = %2f%%",j,100*norm(y-x)/norm(x)));
endfor



# Exemplo 6.10, figuras 6.16 e 6.17
close all;
M=rgb2gray(imread("double_ferris.jpg"))(1:896,1:672);
figure(1);
imshow(M);
title("Imagem original");
for j=1:3
  N=dwt2(M,j,"Le Gall 5/3");
  Nlog = log(1+abs(N));
  Nlog /= max(max(Nlog));
  figure(j+1);
  imshow(Nlog);
  title(sprintf("DWT (Le Gall 5/3) de %d etapa(s)",j));
endfor



# Exemplo 6.10, figuras 6.16 e 6.17
close all;
M=rgb2gray(imread("double_ferris.jpg"))(1:896,1:672);
s = size(M);
figure(1);
imshow(M);
title("Imagem original");
for j=1:3
  N=dwt2(M,j,"Le Gall 5/3");
  Mask = 0*N;
  Mask(1:s(1)/2,1:s(2)/2) += 1;
  O = idwt2(N.*Mask,j,"Le Gall 5/3");
  figure(j+1);
  imshow(uint8(O));
  title(sprintf("Reconstrucao a partir da etapa %d: erro = %2f%%",j,100*norm(O-double(M))/norm(double(M))));
  s /= 2;
endfor