%% SPLS/2e MATLAB
% Chapter 12
clear; close all;
%% Ex. 12.1
Omega = linspace(-pi,pi,401); H = @(z) z./(z-0.8);
subplot(1,2,1); plot(Omega,abs(H(exp(1j*Omega))),'k'); axis tight;
xlabel('\Omega'); ylabel('|H[e^{j \Omega}]|');
subplot(1,2,2); plot(Omega,angle(H(exp(1j*Omega)))*180/pi,'k'); 
axis tight;
xlabel('\Omega'); ylabel('\angle H[e^{j \Omega}] [deg]');
%% Ex. 12.2
Omega = linspace(0,pi,400); 
H = @(z,gamma_m) (z.^2-1)./(z.^2-sqrt(2)*gamma_m*z+gamma_m^2);
plot(Omega,abs(H(exp(1j*Omega),0.83)),...
    Omega,abs(H(exp(1j*Omega),0.96)),...
    Omega,abs(H(exp(1j*Omega),0.99))); 
text(.27*pi,35,'|\gamma|=1'); 
text(.28*pi,25.5,'|\gamma|=0.96');
text(.35*pi,6.41,'|\gamma|=0.83'); 
set(gca,'xtick',0:pi/4:pi,'ytick',[0 6.41,25.5]);
axis([0 pi 0 40]); xlabel('\Omega'); ylabel('|H[e^{j \Omega}]|');
%% Ex. 12.4
omegac = 10^5; Ba = [omegac]; Aa = [1 omegac]; Fs = 10^6/pi;
[B,A] = impinvar(Ba,Aa,Fs)
%% Ex. 12.5
omegac = 10^5; Ba = [omegac]; Aa = [1 omegac]; Fs = 10^6/pi;
[B,A] = bilinear(Ba,Aa,Fs)
%% Ex. 12.6
omegap = 8/35; omegas = 15/35; 
Ghatp = -2; Ghats = -11; T = pi/35;
[N,omegac] = buttord(omegap,omegas,-Ghatp,-Ghats);
[B,A] = butter(N,omegac)
%% Ex. 12.8
omegap = 15/100; omegas = 10/100; 
Ghatp = -1; Ghats = -6.3; T = pi/100;
[N,omegap] = cheb1ord(omegap,omegas,-Ghatp,-Ghats);
[B,A] = cheby1(N,-Ghatp,omegap,'high')
%% Ex. 12.9
omegap = [1000 2000]/10^4; omegas = [450 4000]/10^4;
Ghatp = -2.1; Ghats = -20; T = pi/10000;
[N,omegac] = buttord(omegap,omegas,-Ghatp,-Ghats);
[B,A] = butter(N,omegac)
%% Drill 12.6
omegas = [1000 2000]/10^4; omegap = [450 4000]/10^4;
Ghatp = -2.1; Ghats = -20; T = pi/10000;
[N,omegap] = cheb1ord(omegap,omegas,-Ghatp,-Ghats);
[B,A] = cheby1(N,-Ghatp,omegap,'stop')
%% Ex. 12.11, Fig. 12.20 (rect. window)
mysinc = @(x) sinc(x/pi); Omega = linspace(-pi,pi,1001);
N = 6; fc = 20000; T = 1/(4*fc); n = [0:N]; 
ha = @(t) 1/(2*T)*mysinc(pi*t/(2*T)); 
wrec = @(n) 1.0*((n>=0)&(n<=N));
hrec = T*ha(n*T-N/2*T).*wrec(n)
Hrec = polyval(hrec,exp(1j*Omega)).*exp(-1j*N*Omega); 
clf; plot(Omega,abs(Hrec),'k');
%% Ex. 12.11, Fig. 12.20 (ham. window)
wham = @(n) (0.54-0.46*cos(2*pi*n/N)).*((n>=0)&(n<=N));
hham = T*ha(n*T-N/2*T).*wham(n)
Hham = polyval(hham,exp(1j*Omega)).*exp(-1j*N*Omega); 
plot(Omega,abs(Hham),'k');
%% Ex. 12.12
N = 10; n = 0:N; Omega = linspace(-pi,pi,1001);
wrec = @(n) 1.0*((n>=0)&(n<=N));
wham = @(n) (0.54-0.46*cos(2*pi*n/N)).*((n>=0)&(n<=N));
Threc = cos(pi*(n-N/2))./(n-N/2).*wrec(n) 
Thham = cos(pi*(n-N/2))./(n-N/2).*wham(n)
Threc(n==N/2) = 0; Thham(n==N/2) = 0;
subplot(221); stem(n,Threc,'k.'); subplot(222); stem(n,Thham,'k.');
THrec = polyval(Threc,exp(1j*Omega)).*exp(-1j*N*Omega);
THham = polyval(Thham,exp(1j*Omega)).*exp(-1j*N*Omega);
subplot(212); plot(Omega,abs(THrec),Omega,abs(THham),'--k');
%% Ex. 12.13
N = 6; N0 = N+1; r = 0:N; Omega0 = 2*pi/N0;
Hr = 1.0*(r*Omega0<=pi/2).*exp(-1j*r*pi*(N0-1)/N0);
Hr(1,fix(N0/2)+2:end) = conj(Hr(1,round(N0/2):-1:2));
h = ifft(Hr)
Omega = linspace(0,2*pi,1001);
H = polyval(h,exp(1j*Omega)).*exp(-1j*N*Omega);
clf; stem(r*Omega0,abs(Hr),'k.');
line(Omega,abs(H),'color',[0 0 0]);
xlabel('\Omega'); ylabel('|H[e^{j\Omega}]|');
%% Sec. 12.9.1, Ex. 12.14
omegac = 1; Omegac = pi/3; N = 20; M = 0; 
k = 1:N; pk = exp(1j*pi*(2*k+N-1)/(2*N));
T = 2/omegac*tan(Omegac/2);
gain = 1/prod(2/T-pk); zeros = -1*ones(1,N);
poles = (1+pk*T/2)./(1-pk*T/2);
B = gain*poly(zeros); A = poly(poles);
Omega = linspace(0,pi,1001);
H = polyval(B,exp(1j*Omega))./polyval(A,exp(1j*Omega));
subplot(121); circle = exp(1j*(0:.01:2*pi));
plot(real(zeros),imag(zeros),'ko',real(poles),imag(poles),'kx');
theta = 0:.01:2*pi; circle = exp(1j*theta); 
line(real(circle),imag(circle),'color',[0 0 0],'linestyle',':');
axis equal; axis([-1.1 1.1 -1.1 1.1]);
subplot(122);
plot(Omega,abs(H),'k-'); axis([0 pi 0 1.1]);
xlabel('\Omega'); ylabel('|H[e^{j\Omega}]|');
set(gca,'xtick',0:pi/3:pi); grid on;
%% Sec. 12.9.1, Ex. 12.15, Fig. 12.24
N = 30; N0 = N+1; r = 0:N; Omega0 = 2*pi/N0;
Hr = 1.0*(r*Omega0<=pi/2);
Hr(1,fix(N0/2)+2:end) = conj(Hr(1,round(N0/2):-1:2));
h = ifft(Hr);
Omega = linspace(0,2*pi,1001);
H = polyval(h,exp(1j*Omega)).*exp(-1j*N*Omega);
subplot(211); stem(0:N,h,'k.'); axis([0 N -.2 .6]);
xlabel('n'); ylabel('h[n]');
subplot(212); stem(r*Omega0,abs(Hr),'k.');
line(Omega,abs(H),'color',[0 0 0]);
xlabel('\Omega'); ylabel('|H[e^{j\Omega}]|');
axis([0 2*pi 0 1.7]); set(gca,'xtick',0:pi/2:2*pi);
%% Sec. 12.9.1, Ex. 12.15, Fig. 12.25
N = 30; N0 = N+1; r = 0:N; Omega0 = 2*pi/N0;
Hr = 1.0*(r*Omega0<=pi/2).*exp(-1j*r*pi*(N0-1)/N0);
Hr(1,fix(N0/2)+2:end) = conj(Hr(1,round(N0/2):-1:2));
h = ifft(Hr);
Omega = linspace(0,2*pi,1001);
H = polyval(h,exp(1j*Omega)).*exp(-1j*N*Omega);
subplot(211); stem(0:N,h,'k.'); axis([0 N -.2 .6]);
xlabel('n'); ylabel('h[n]');
subplot(212); stem(r*Omega0,abs(Hr),'k.');
line(Omega,abs(H),'color',[0 0 0]);
xlabel('\Omega'); ylabel('|H[e^{j\Omega}]|');
axis([0 2*pi 0 1.7]); set(gca,'xtick',0:pi/2:2*pi);
%% Ex. 12.16
N = 50; N0 = N+1; r = 0:N; Omega0 = 2*pi/N0;
Hr = (2-r*Omega0).*(r*Omega0<=2).*exp(-1j*r*pi*(N0-1)/N0);
Hr(1,fix(N0/2)+2:end) = conj(Hr(1,round(N0/2):-1:2));
h = ifft(Hr);
Omega = linspace(0,2*pi,1001);
H = polyval(h,exp(1j*Omega)).*exp(-1j*N*Omega);
subplot(211); stem(0:N,h,'k.'); axis([0 N -.2 .8]);
xlabel('n'); ylabel('h[n]');
subplot(212); stem(r*Omega0,abs(Hr),'k.');
line(Omega,abs(H),'color',[0 0 0]);
xlabel('\Omega'); ylabel('|H[e^{j\Omega}]|');
axis([0 2*pi 0 2.1]); set(gca,'xtick',[0 2 pi 2*pi-2 2*pi]);