# Exemplo da seção 4.2.1
# Resposta do filtro w(n) = 0.5*x(n)+0.5*x(n-1) às entradas
# da forma x(k) = sen(2*pi*q*k/N), para k=0,...,N-1
N=100; k=0:N-1;
W=[];
figure(1);
for q=0:N/2
  # x(k) é a q-ésima função "básica"
  x = sin(2*pi*q*k/N);
  # Coeficientes A e B da expressão de w(k)
  A = 0.5*(1+cos(2*pi*q/N)); B = 0.5*sin(2*pi*q/N);
  w = A*sin(2*pi*q*k/N) - B*cos(2*pi*q*k/N);
  # Mostra a saída do filtro para cada função básica.
  # Observe como a amplitude diminui e a fase inicial
  # também varia.
  plot(w);
  set(gca(),"ylim",[-2 2]);
  pause(0.5);
  # W contém a resposta em magnitude (ou seja, o fator de
  # escala) do filtro para as funções "básicas". Os "eps"
  # só estão aqui por causa da primeira função, que é
  # x(k) = sen(2*pi*0*k/N) = sen(0) = 0, com norma(x)=0.
  W(q+1) = (norm(w)+eps)/(norm(x)+eps);
endfor
# Mostra o gráfico da resposta em magnitude 
figure(2);
plot(W);



# Exemplo 4.1 e figuras 4.1 e 4.2
# Retoma a função do exemplo na seção 2.3
T=1;
N=128;
t=linspace(0,T,N);
x = 2*cos(2*pi*5*t) + 0.8*sin(2*pi*12*t) + 0.3*cos(2*pi*47*t);
# Calcula a saída do filtro da média. A função shift é usada
# para obter o vetor (x(N),x(1),x(2),...,x(N-1))
w = 0.5*(x+shift(x,1));
# Figura 4.1
figure(1);
plot(t,w,t,x);
# Calcula as FFTs do sinal original e do sinal filtrado
X = fft(x); W = fft(w);
c = X/N; d = W/N;
E = N*abs(c).^2; F = N*abs(d).^2;
# Figura 4.2
figure(2);
plot(0:N/2,E(1:N/2+1),0:N/2,F(1:N/2+1));
# Como é difícil observar as diferenças dos espectros,
# as 3 figuras a seguir fazem um zoom dos 3 picos dos
# espectros de x (sinal original) e w (sinal filtrado)
figure(3);
plot(3:7,E(4:8),3:7,F(4:8));
figure(4);
plot(10:15,E(11:16),10:15,F(11:16));
figure(5);
plot(40:50,E(41:51),40:50,F(41:51));



# Exemplo 4.3 e Figura 4.3:
# O filtro da média w(n) = 0.5*x(n)+0.5*x(n-1) pode ser
# escrito como w = x*l, onde * é a convolução. Os
# coeficientes do filtro são l(0)=1/2 e l(1)=1/2.
N = 128;
l = [ 0.5 0.5 zeros(1,N-2) ];
# Mostra os espectros de magnitude e fase de L (Figura 4.3)  
L = fft(l);
figure(1);
plot(-N/2+1:N/2,fftshift(abs(L)));
figure(2);
plot(-N/2+1:N/2,fftshift(arg(L)));



# Exemplo 4.4: filtro passa-alta projetado a partir
# de uma resposta em frequência ideal.
close all;
N = 128;
# H deverá cortar todas as frequências até k=29 (inclusive)
H = [ zeros(1,30) ones(1,69) zeros(1,29) ] ;
# Calcula a resposta impulsiva do filtro
h = real(ifft(H));
# Figura 4.4: resposta impulsiva
figure(1);
plot(h);
# Resposta em frequência do filtro
figure(2);
plot(0:N-1,H);
set(gca(),"ylim",[-0.5 1.5]);

# Versões computacionalmente mais "baratas" (com menos taps)
# do mesmo filtro passa-alta. Os limiares abaixo são usados
# para "cortar" coeficientes do filtro e com isso diminuir
# o número de termos na equação de convolução
EPSVALS =  [ 0.01 0.05 0.1 ];
for k = 1:length(EPSVALS)
  epsilon = EPSVALS(k);
  # seleciona os coeficientes acima do limiar
  mask = (abs(h)>=epsilon);
  # calcula o número de coeficientes do filtro
  ntaps = sum(mask);
  # heps é o filtro simplificado com limiar epsilon
  heps = h.*mask;
  # Mostra a resposta impulsiva com os coeficientes eliminados
  figure(2*(k-1)+1);
  plot(heps);
  title(sprintf("Epsilon = %g, Filtro com %d taps",epsilon,ntaps));
  # Compara a resposta do filtro simplificado com a resposta ideal
  # do filtro original (pontilhada)
  figure(2*k);
  plot(0:N-1,real(fft(heps)),0:N-1,H,".");
endfor



# Colocando em perspectiva a "perfeição" da resposta em frequência
# do filtro original (com banda de corte nula). Este exemplo (que
# não está no livro) se contrapõe à idéia de que o filtro desenhado
# a partir da resposta em frequência H tenha exatamente o comportamento
# de filtro ideal que o gráfico de H poderia sugerir.
# Para este exemplo, construiremos "na mão" a resposta em magnitude
# do filtro H para funções "básicas" do tipo cos(2*pi*f*n/N),
# inicialmente para frequências do tipo f=0,1,...,N-1, e
# posteriormente para outras frequências intermediárias.
close all;
FREQS = 0:N-1;
R = [];
for k=1:length(FREQS);
  # cria função básica de frequência k
  x = cos(2*pi*FREQS(k)*(0:N-1)/N);
  # filtra por h
  y = real(ifft(fft(x).*H));
  # constrói gráfico de magnitude da resposta
  # (fator de escala da saída do filtro)
  R(k) = norm(y)/norm(x);
endfor
# Mostra o gráfico da resposta de magnitude.
figure(1);
plot(FREQS,R);
set(gca(),"ylim",[-0.5 1.5]);
# Repete a mesma construção com frequências não-inteiras
# Observe que a única diferença do código abaixo é que
# agora as frequências são f=0,0.1,0.2,...,N-1
FREQS = 0:0.1:N-1;
R = [];
for k=1:length(FREQS)
  # cria função básica com a k-ésima frequência da lista
  x = cos(2*pi*FREQS(k)*(0:N-1)/N);
  # filtra por h
  y = real(ifft(fft(x).*H));
  # constrói gráfico de magnitude da resposta
  R(k) = norm(y)/norm(x);
endfor
figure(2);
plot(FREQS,R);
set(gca(),"ylim",[-0.5 1.5]);

# parei aqui (10/10/2013)

# Figura 4.6: Magnitude da resposta em frequência
# da "máscara" de suavização de 9 pontos
M = N = 100;
h = zeros(M,N);
# Cria um bloco 3x3 centralizado em (0,0) 
for j=-1:1
  for k=-1:1
    h(mod(j,M)+1,mod(k,N)+1) = 1/9;
  endfor
endfor
# Duas visualizações alternativas para a DFT.
H = fft2(h);
figure(1);
plot3(-M/2+1:M/2,-N/2+1:N/2,abs(fftshift(H)));
figure(2);
imshow(abs(fftshift(H)));



# Figuras 4.7 e 4.8: Uma imagem "original", uma versão ruidosa,
# e duas versões filtradas
close all;
M=rgb2gray(imread("double_ferris.jpg"));
S=size(M);
# Imagem original
figure(1);
imshow(M);
# Versão com ruído adicionado:
Mnoisy = M + 20*(rand(S)-0.5);
figure(2);
imshow(Mnoisy);
# Filtragem usando média de 9 pontos.
# Define a resposta impulsiva do filtro
d = zeros(S(1),S(2));
for j=-1:1
  for k=-1:1
    d(mod(j,S(1))+1,mod(k,S(2))+1) = 1/9;
  endfor
endfor
# O código abaixo faz a filtragem no domínio da frequência.
# A filtragem usando a convolução direta poderia ser mais
# rápida numa implementação em C (por exemplo), mas em uma
# linguagem interpretada (como Octave ou Matlab) a versão
# espectral é muito mais rápida.
Mdenoised = uint8(real(ifft2(fft2(Mnoisy).*fft2(d))));
figure(3);
imshow(Mdenoised);
# Passa a imagem ruidosa 5 vezes no mesmo filtro. Neste
# caso, mesmo em C o processamento espectral seria
# provavelmente mais rápido, pois basta multiplicar
# pelo espectro elevado à 5a potência.
Mdenoised5 = uint8(real(ifft2(fft2(Mnoisy).*fft2(d).^1000)));
figure(4);
imshow(Mdenoised5);



# Figuras 4.9 e 4.10: detecção de bordas na horizontal e vertical
close all;
figure(1);
imshow(M);
# filtro para detecção de bordas na horizontal
d = zeros(S(1),S(2));
d(1,1)=1; d(S(1),1)=-1;
Mhorizontal = abs(ifft2(fft2(M).*fft2(d)));
maximo = max(max(Mhorizontal));
# mostra a imagem filtrada
figure(2);
imshow(1-Mhorizontal/maximo);
# filtro para detecção de bordas na vertical
d = zeros(S(1),S(2));
d(1,1)=1; d(1,S(2))=-1;
Mvertical = abs(ifft2(fft2(M).*fft2(d)));
maximo = max(max(Mvertical));
# mostra a imagem filtrada
figure(3);
imshow(1-Mvertical/maximo);
# combina os resultados das filtragens horizontal e vertical
Mbordas = sqrt(Mhorizontal.^2+Mvertical.^2);
maximo = max(max(Mbordas));
# bordas combinadas
figure(4);
imshow(1-Mbordas/maximo);