# Figura 1.1
# Cria um domínio temporal com 800 amostras/seg
# entre 0 e 4*pi (esta é uma amostragem arbitrária,
# mas que permite visualizar a função
# "como se fosse contínua")
t = 0:(1/800):4*pi;
# Define a função do exemplo da página 4
x = 0.75*sin(3*t) + 0.5*sin(7*t);
# Mostra o resultado. O argumento "." impede os
# pontos adjacentes de serem conectados por
# segmentos de reta.
plot(t,x,".");



# Figura 1.2
# Seleciona pontos no domínio separados de
# 100 amostras (0.125seg)
samples = 1:100:length(t);
ts = t(samples);
# (Re-)Amostra o sinal
xs = x(samples);
# Mostra o sinal amostrado
plot(ts,xs,".");



# transpõe o sinal para a faixa audível e toca
# gera 5 segundos de sinal a 8000 amostras/segundo
T = 0:(1/8000):5;
# multiplica as frequências por 1000 para cairem
# na faixa audível. As frequências originais eram
# 3/(2*pi) = 0.47746 Hz e 7/(2*pi) = 1.1141 Hz,
# as novas serão 477.46 Hz e 1114.1 Hz
X = 0.75*sin(1000*3*T) + 0.5*sin(1000*7*T);
wavwrite(0.5*X',"tmp.wav");
system("play tmp.wav");



# Alternativa à Fig. 1.3: sinal amostrado com
# quantização de 4 bits (16 valores)
# Normaliza o sinal amostrado
lb = min(xs);
ub = max(xs);
xsnorm = (xs-lb)/(ub-lb);
# Quantiza usando 16 valores (0...15)
xqnorm = floor(16*xsnorm);
# Coloca de volta na faixa de valores originais
# [lb...ub], usando como codewords os valores
# médios das 16 faixas entre lb e ub
xq = (ub-lb)*xqnorm/16+lb+0.5*(ub-lb)/16;
plot(t,x,".",ts,xq,"*");

# Plota apenas o erro de quantização
figure(2);
plot(ts,xq-xs,".");
axis([0 14 -1.5 1.5]);
close all

# toca sinal quantizado
# 1. normaliza
lb = min(X);
ub = max(X);
XNORM = (X-lb)/(ub-lb);
# 2. quantiza
XQNORM = floor(16*XNORM);
XQ = (ub-lb)*XQNORM/16+lb+0.5*(ub-lb)/16;
wavwrite(0.5*XQ',"tmp.wav");
system("play tmp.wav");



# toca erro de quantização
wavwrite(0.5*(XQ-X)',"tmp.wav");
system("play tmp.wav");



# Figura 1.3
# Sinal amostrado adicionado de ruído uniforme 
xr = xs + 0.5*(rand(1,length(xs))-0.5);
# Mostra o sinal original e o sinal amostrado
# c/ ruído
plot(t,x,".",ts,xr,"*");



# toca sinal com ruído
R = 0.5*(rand(1,length(T))-0.5);
wavwrite(0.5*(X+R)',"tmp.wav");
system("play tmp.wav");



# toca só ruído
wavwrite(0.5*R',"tmp.wav");
system("play tmp.wav");



# Figura 1.4
# Define um domínio espacial com 201x201 pontos
x = 0:(1/200):1;
# Lembre-se que os índices crescentes da matriz
# refletem valores descrescentes do eixo vertical
y = 1:-(1/200):0;
# Calcula a função do exemplo da página 7
f = 1.5*cos(7*y)'*cos(2*x)+ 0.75*sin(3*y)'*cos(5*x) ...
  -1.3*cos(15*y)'*sin(9*x) + 1.1*sin(11*y)'*sin(13*x);
# Plota o gráfico 3D da função f(x,y)
figure(1);
mesh(x,y,f);
#colormap("gray");
# Normaliza f
lb = min(min(f));
ub = max(max(f));
g = (f-lb)/(ub-lb);
# Exibe (em outra janela) a imagem representada por f
figure(2);
imshow(g);

# quantiza a imagem com 4 bits
gq = floor(16*g)/16+1/32;
figure(3);
imshow(gq);

# quantiza a imagem com 2 bits
gq = floor(4*g)/4+1/2;
figure(3);
imshow(gq);

# quantiza a imagem com 1 bit
gq = floor(2*g)/2+1/4;
figure(3);
imshow(gq);



# Exemplos adicionais
# (Re)Amostra a imagem preservando 1/4
# dos pontos originais
xsamples = 1:4:length(x);
ysamples = 1:4:length(y);
gs = g(ysamples,xsamples);
# Exibe a imagem reamostrada
figure(3);
imshow(gs);
# Quantiza a imagem usando 4 bits
# a expressão da quantização aqui é bem
# mais simples que a anterior, pois os
# extremos de cor [lb=0...ub=1] são
# conhecidos a priori
gq = floor(16*gs)/16+1/32;
figure(4);
imshow(gq);



# Figura 1.5: Imagem original com ruído aditivo
gn = g + 0.1*(rand(length(y),length(x))-0.5);
figure(5);
imshow(gn);



# Seção 1.5



# Exemplo 1.10
close all;
# Domínio temporal de 0 a 2*pi com 100 amostras
t = linspace(0,2*pi,100);
# Define duas funções equivalentes usando senos e cossenos (x1)
# e exponenciais complexas (x2)
x1 = sin(t)+3*sin(-2*t)-2*cos(-5*t);
x2 = (e.^(i*t)-e.^(-i*t))/(2*i)+(-3*e.^(2*i*t)+3*e.^(-2*i*t))/(2*i) \
     -(e.^(5*i*t)+e.^(-5*i*t));
# Plota as funções x1 e x2
figure(1);
plot(t,x1);
figure(2);
plot(t,real(x2));


# Figura 1.6
# Fecha todas as figuras
close all;
# Cria um domínio temporal de 10s com 200 amostras
t = linspace(0,10,200);
# Define uma exponencial complexa com
# frequência 0.5 rad/seg
f = e.^(0.5*i*t);
# Plota as partes real e imaginária da função f
figure(1);
plot(t,real(f),t,imag(f),".");
# Idem para a frequência 2 rad/seg
g = e.^(2*i*t);
figure(2);
plot(t,real(g),t,imag(g),".");

# Alternativa: plota os mesmos gráficos como uma função de R em C
figure(3);
plot3(t,f);
figure(4);
plot3(t,g);


# Ver nesse ponto exemplos em video das exponenciais complexas
# (há um link separado com o código que gera estes vídeos)

# Fecha todas as figuras
close all;
# Cria um domínio temporal de 10s com 1000 amostras
t = linspace(0,10,1000);
e1 = exp(i*2*pi*0.2*t);
e2 = 0.5*exp(i*2*pi*(-1)*t);
e3 = 0.25*exp(i*2*pi*3.8*t);
e4 = e1+e2+e3;
figure(1);
plot(t,real(e1),t,imag(e1),".");
figure(2);
plot(t,real(e2),t,imag(e2),".");
figure(3);
plot(t,real(e3),t,imag(e3),".");
figure(4);
plot(t,real(e4),t,imag(e4),".");


# Seção 1.5.2


# Figura 1.7
close all;
# Define todos os pares (p,q) ilustrados na figura 1.7
PQ = [10 0; 0 10; 30 8; 5 10];
# Percorre os 4 cenários
for j=1:4
  p=PQ(j,1); q=PQ(j,2);
  # Define domínios espaciais (horizontal e vertical)
  # usando as frequências p e q (em Hz)
  x = linspace(0,1,200);
  y = linspace(1,0,200);
  # Define a exponencial complexa no plano xy
  f = e.^(2*pi*i*q*y)'*e.^(2*pi*i*p*x);
  figure(j);
  # Corrige a escala para [0,1]
  imshow(0.5*(real(f)+1));
  title(sprintf("p=%d, q=%d",p,q));
endfor



# Ver nesse ponto exemplos de reconstrução temporal/espacial e espectral (bell*.wav, double_ferris*.avi)
# (há um link separado com o código que gera estes vídeos)


# Seção 1.6.1


# Figura 1.8
close all;
# Exemplo simples de sinal 1D amostrado a 20Hz
t = linspace(0,1,21);
x = sin(44*pi*t);
figure(1);
plot(t,x,"*");


# Figura 1.9
close all;
t = linspace(0,1,21);
x = sin(44*pi*t);
# A mesma função anterior, amostrada a 200Hz (linha contínua),
# ilustrando aliasing temporal
t2 = linspace(0,1,201);
x2 = sin(44*pi*t2);
plot(t,x,"*",t2,x2);


#Seção 1.6.5


# Figura 1.10
# Exemplo de aliasing espacial: a mesma função é amostrada em
# grids de 60x60, 100x100, 300x300 e 1000x1000
close all;
N=[60 100 300 1000];
for j=1:4
  x = linspace(0,1,N(j));
  y = linspace(1,0,N(j));
  f = e.^(2*pi*i*70*y)'*e.^(2*pi*i*50*x);
  figure(j);
  imshow(0.5*(imag(f)+1));
  title(sprintf("amostragem com %dx%d pontos",N(j),N(j)));
  drawnow();
endfor

# Alternativa: mostra a variação do tamanho do grid de 10 em 10
close all;
figure(1);
for n=50:10:400
  x = linspace(0,1,n);
  y = linspace(1,0,n);
  f = e.^(2*pi*i*70*y)'*e.^(2*pi*i*50*x);
  imshow(0.5*(imag(f)+1));
  title(sprintf("amostragem com %dx%d pontos",n,n));
  drawnow();
  pause(0.2);
endfor




# Seção 1.10.3


# Exemplo 1.27 adaptado à função f(t)=t
T=1;
t=-1:0.01:1;
f=t;
# Computa os coeficientes de Fourier calculados
# explicitamente usando integração por partes
alfa0 = 0;
K=1:100;
for k=K
  omega(k) = k*pi/T;
  alfa(k) = -(-1)^k/(i*omega(k));
endfor
# Figura 1.13, adaptada ao exemplo f(t)=t,
# na forma de "animação".
# Faz a ressíntese de f(t) usando
# progressivamente mais harmônicos 
g=0*f+alfa0;
for k=K
  plot(t,f,".",t,g);
  axis([-1 1 -1.2 1.2]);
  title(sprintf("Ressintese de f(t)=t com %d harmonico(s)",k));
  drawnow;
  g += real(alfa(k)*e.^(i*omega(k)*t)
            +conj(alfa(k))*e.^(-i*omega(k)*t));
  pause(2*0.9^k);
endfor



# Exemplo 1.27, com a função original f(t)=t^2
f=t.*t;
# Computa os coeficientes de Fourier calculados
# explicitamente usando integração por partes
beta0 = 2/6;
for k=K
  beta(k) = 2*alfa(k)/(i*omega(k));
endfor
# Figura 1.13, adaptada ao exemplo f(t)=t^2,
# na forma de "animação".
# Faz a ressíntese de f(t) usando
# progressivamente mais harmônicos 
g=0*f+beta0;
for k=K
  plot(t,f,".",t,g);
  axis([-1 1 -0.2 1.2]);
  title(sprintf("Ressintese de f(t)=t^2 com %d harmonico(s)",k));
  drawnow;
  g += real(beta(k)*e.^(i*omega(k)*t)
            +conj(beta(k))*e.^(-i*omega(k)*t));
  pause(4*0.9^k);
endfor



# Exemplo 1.27 adaptado à função f(t)=t^3
f=t.*t.*t;
# Computa os coeficientes de Fourier calculados
# explicitamente usando integração por partes
gama0 = 0;
for k=K
  gama(k) = (3*beta(k)-(-1)^k)/(i*omega(k));
endfor
# Figura 1.13, adaptada ao exemplo f(t)=t^3,
# na forma de "animação".
# Faz a ressíntese de f(t) usando
# progressivamente mais harmônicos 
g=0*f+gama0;
for k=K
  plot(t,f,".",t,g);
  axis([-1 1 -1.2 1.2]);
  title(sprintf("Ressintese de f(t)=t^3 com %d harmonico(s)",k));
  drawnow;
  g += real(gama(k)*e.^(i*omega(k)*t)
            +conj(gama(k))*e.^(-i*omega(k)*t));
  pause(2*0.9^k);
endfor



# Figura 1.14, adaptada ao exemplo f(t)=t^3,
# na forma de "animação".
t=-3:0.01:3;
f=t.*t.*t;
g=0*f+gama0;
for k=K
  plot(t,g);
  axis([-3 3 -1.2 1.2]);
  title(sprintf("Ressintese de f(t)=t^3 com %d harmonico(s)",k));
  drawnow;
  g += real(gama(k)*e.^(i*omega(k)*t)
            +conj(gama(k))*e.^(-i*omega(k)*t));
  pause(2*0.9^k);
endfor


# Figura 1.14, adaptada ao exemplo f(t)=sign(t),
# na forma de "animação".
f=sign(t);
# Computa os coeficientes de Fourier calculados
# explicitamente usando integração por partes
delta0 = 0;
for k=K
  delta(k) = mod(k,2)*2/(i*omega(k));
endfor
g=0*f+delta0;
for k=K
  plot(t,g);
  axis([-3 3 -1.2 1.2]);
  title(sprintf("Ressintese de f(t)=sign(t) com %d harmonico(s)",k));
  drawnow;
  g += real(delta(k)*e.^(i*omega(k)*t)
            +conj(delta(k))*e.^(-i*omega(k)*t));
  pause(2*0.9^k);
endfor