%% SPLS/2e MATLAB
% Chapter 10
clear; close all;
%% Ex. 10.2
N_0 = 32; n = (0:N_0-1); Omega_0 = 2*pi/N_0;
x_n = [ones(1,5) zeros(1,23) ones(1,4)];
for r = 0:N_0-1,
    X_r(r+1) = sum(x_n.*exp(-j*r*Omega_0*n))/N_0;
end
r = n; stem(r,real(X_r),'k.');
xlabel('r'); ylabel('X_r'); axis([0 31 -.1 0.3]);
X_r = fft(x_n)/N_0; stem(r,real(X_r),'k.');
xlabel('r'); ylabel('X_r'); axis([0 31 -.1 0.3]);
%% Ex. 10.5
Omega = linspace(0,2*pi,1000); 
X = sin(4.5*Omega)./sin(0.5*Omega); 
X(mod(Omega,2*pi)==0) = 4.5/0.5; N_0 = 64; M = 9; 
x = [ones(1,(M+1)/2) zeros(1,N_0-M) ones(1,(M-1)/2)];
Xr = fft(x); Omega_0 = 2*pi/N_0; r = 0:N_0-1;
plot(Omega,abs(X),'k-',Omega_0*r,abs(Xr),'k.'); axis([0 2*pi 0 9.5]);
xlabel('\Omega'); ylabel('|X(\Omega)|');
%% Ex. 10.16
N0 = 32; x = [3 2 3 zeros(1,N0-3)];
Omega0 = 2*pi/N0; r = 0:N0-1; Xr = fft(x); 
Omega = linspace(0,2*pi,1001);
X = @(Omega) exp(-1j*Omega).*(2+6*cos(Omega));
subplot(121); stem(r,abs(Xr),'k.');
xlabel('r'); ylabel('|X_r|'); axis([0 N0 0 8.5]);
line(Omega/Omega0,abs(X(Omega)),'color',[0 0 0]);
subplot(122); stem(r,angle(Xr),'k.');
xlabel('r'); ylabel('\angle X_r'); axis([0 N0 -pi pi]);
line(Omega/Omega0,angle(X(Omega)),'color',[0 0 0]);
%% Ex. 10.19
x = [3 2 3]; h = [2 2]; N0 = length(x)+length(h)-1; n = 0:N0-1;
y = ifft(fft(x,N0).*fft(h,N0));
subplot(121); stem(n,y,'k.'); xlabel('n'); ylabel('y[n]');
axis([-.5 N0-.5 0 12]); set(gca,'ytick',[0,6,10]);
y = conv(x,h);
subplot(122); stem(n,y,'k.'); xlabel('n'); ylabel('y[n]');
axis([-.5 N0-.5 0 12]); set(gca,'ytick',[0,6,10]);
%% Ex. 10.20
N0 = 32; n = 0:N0-1; x = (0.8).^(n); h = (0.5).^(n);
y = ifft(fft(x).*fft(h));
subplot(121); stem(n,y,'k.'); xlabel('n'); ylabel('y[n]');
axis([-.5 N0-.5 0 1.5]); 
y = ((0.8).^(n+1)-(0.5).^(n+1))/(0.8-0.5);
subplot(122); stem(n,y,'k.'); xlabel('n'); ylabel('y[n]');
axis([-.5 N0-.5 0 1.5]); 
%% Sec. 10.8.1, Fig. 10.32
T = 1/1000; N0 = 100; n = (0:N0-1)';
x = cos(2*pi*50*n*T);
X = fft(x)/N0; f = (0:N0-1)/(T*N0);
clf; stem(f-1/(2*T),fftshift(abs(X)),'k.'); 
axis([-500 500 -0.05 0.55]); xlabel('f [Hz]'); ylabel('|X(f)|');
%% Sec. 10.8.1, Fig. 10.33
x = real(ifft(X)*N0); stem(n,x,'k.'); 
axis([0 99 -1.1 1.1]); xlabel('n'); ylabel('x[n]');
%% Sec. 10.8.1, Fig. 10.32 via DFT matrix
W = @(N0) (exp(-1j*2*pi/N0)).^((0:N0-1)'*(0:N0-1));
X = W(N0)*x/N0; stem(f-1/(2*T),fftshift(abs(X)),'k.');
axis([-500 500 -0.05 0.55]); xlabel('f [Hz]'); ylabel('|X(f)|');
%% Sec. 10.8.2
tic; W(N0)*x/N0; toc
tic; for i=1:100, W(N0)*x/N0; end; toc
W100 = W(100); tic; for i=1:100, W100*x/N0; end; toc
tic; for i=1:100, fft(x)/N0; end; toc
x1 = rand(1015,1); tic; for i=1:100; fft(x1)/1015; end; T1 = toc
x2 = [x1;zeros(4,1)]; tic; for i=1:100; fft(x2)/1019; end; T2 = toc
x3 = [x2;zeros(5,1)]; tic; for i=1:100; fft(x3)/1024; end; T3 = toc