# Figura 2.2

T=1;
N=128;
t=linspace(0,T,N+1)(1:N);
x = 2*cos(2*pi*5*t) + 0.8*sin(2*pi*12*t) + 0.3*cos(2*pi*47*t);
figure(1);
plot(t,x);

# Figura 2.3

X = fft(x);
c = X/N;
# Mostra coeficientes não-nulos
for k=(1:N)(abs(c)>1e-8)
  printf("c(%g)=%.2f%+.2fi\n",k-1,real(c(k)),imag(c(k)));
endfor
E = N*abs(c).^2;
# Mostra energia em cada componente
for k=(1:N)(abs(c)>1e-8)
  printf("E(%g)=%.2f\n",k-1,E(k));
endfor
figure(2);
# plota o trecho (-N/2,+N/2], lembrando
# que os índices em Octave começam em 1
plot(-N/2+1:N/2,[E(N/2+2:N) E(1:N/2+1)]);

# Figura 2.4

plot(0:N-1,E);

# Figura 2.5

plot(0:N/2,[E(1) 2*E(2:N/2) E(N/2)]);

# Figura 2.6

Y = X;
# seleciona índices entre 40 e N-40
Y((40+1):(N-40+1)) *= 0;
y = ifft(Y);
figure(3);
plot(t,y);

# Figura 2.7

tt = linspace(0,1,N/2+1)(1:N/2);
xx = 2*cos(2*pi*5*tt) + 0.8*sin(2*pi*12*tt) + 0.3*cos(2*pi*47*tt);
figure(1);
plot(tt,xx);
XX = fft(xx);
cc = XX/(N/2);
EE = (N/2)*abs(cc).^2;
figure(2);
plot(0:N/4,EE(1:N/4+1));

# Figuras 2.8-2.10

close all;
N=256;
t = 0:N-1;
R = [50 10 2 1];
for j=1:4
  x = [ones(1,R(j)) zeros(1,N-R(j))];
  X = fft(x);
  c = X/N;
  E = N*abs(c).^2;
  figure(j);
  plot(-N/2+1:N/2,[E(N/2+2:N) E(1:N/2+1)]);
endfor

# Figura 2.11

close all;
N=256;
t=0:N-1;
x = randn(1,N);
figure(1);
plot(t,x);
X = fft(x);
c = X/N;
E = N*abs(c).^2;
figure(2);
plot(-N/2+1:N/2,[E(N/2+2:N) E(1:N/2+1)]);


# Figura 2.12, usando bell.wav

close all;
page_screen_output(0)
[x,SR]=wavread("bell.wav");
N=length(x);
figure(1);
plot(x);
X=fft(x);
figure(2);
plot((0:N/2)*SR/N,abs(X(1:N/2+1)));
figure(3);
plot((0:2999)*SR/N,abs(X(1:3000)));
[S,I]=sort(abs(X),"descend");
sig = sinetone(0,SR,N/SR,0);
for k=1:30
  printf("freq=%.2fHz, abs(X(%f))=%.2f\n",(I(k)-1)*SR/N,I(k)-1,S(k));
  sig += sinetone(min(I(k)-1,N-I(k)+1)*SR/N,SR,N/SR,S(k)/N);
endfor
wavwrite(sig,SR,"tmp.wav");
system("aplay tmp.wav");

# Figura 2.13

close all;
M=imread("double_ferris.jpg");
figure(1);
imshow(M);
N=fft2(M);
figure(2);
imshow(255*abs(N)/max(max(max(abs(N)))));
Ncenter=N;
Ncenter(:,:,1)=fftshift(N(:,:,1));
Ncenter(:,:,2)=fftshift(N(:,:,2));
Ncenter(:,:,3)=fftshift(N(:,:,3));
figure(3);
imshow(255*abs(Ncenter)/max(max(max(abs(Ncenter)))));
Nlog = log(1+abs(Ncenter));
minimo=min(min(min(Nlog)));
maximo=max(max(max(Nlog)));
Nshow=M;
Nshow(:,:,:)=255*(Nlog(:,:,:)-minimo)/(maximo-minimo);
figure(4);
imshow(Nshow);
N=fft2(M);
[S,I]=sort(abs(N));
C=0.001*max(max(max(abs(N))));
N=N.*(abs(N)>C);
printf("Usando %f%% maiores coeficientes de Fourier.\n",100*(sum(sum(sum(abs(N)>C)))/prod(size(N))));
L=round(real(ifft2(N)));
K=M;
K(:,:,:)=L(:,:,:);
figure(5);
imshow(K);


# Figura 2.14

close all;
O=zeros(300,300);
O(100:150,100:110) += 1;
figure(1);
imshow(O);
N=fft2(O);
figure(2);
imshow(abs(N)/max(max(max(abs(N)))));
Ncenter=fftshift(N);
figure(3);
imshow(abs(Ncenter)/max(max(max(abs(Ncenter)))));
Nlog = log(1+abs(Ncenter));
minimo=min(min(min(Nlog)));
maximo=max(max(max(Nlog)));
Nshow=O;
Nshow(:,:,:)=(Nlog(:,:,:)-minimo)/(maximo-minimo);
figure(4);
imshow(Nshow);