% eda09_13
% complex z-plane analysis of filters, case 4, notch

% i used here as imaginary unit!
clear i;

% standard time/frequency setup
N=32;
Nf=N/2+1;
Dt=1.0;
tmax=Dt*(N-1); 
t=Dt*[0:N-1]'; 
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]'; 
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]'; 

% two roots
F1 = [1, 0.9*i]';
zp1 = -F1(1)/F1(2);
F2 = conj(F1);
zp2 = -F2(1)/F2(2);

% two poles
F3 = [1, 0.8*i]';
zp3 = -F3(1)/F3(2);
F4 = conj(F3);
zp4 = -F4(1)/F4(2);

disp(sprintf('zeros at (%f %f) and (%f %f)',real(zp1),imag(zp1),real(zp2),imag(zp2)));
disp(sprintf('zeros at (%f %f) and (%f %f)',real(zp3),imag(zp3),real(zp4),imag(zp4)));

% complex z-plane
Nz=128;
zmin=-2;
zmax=2;
dz=(zmax-zmin)/(Nz-1); 
zr=zmin+dz*[0:Nz-1]'; 
zi=zmin+dz*[0:Nz-1]';

% evaluate filter polynomial everywhere in complex z-plane
S = zeros(Nz,Nz);
for n = [1:Nz]
    for m = [1:Nz]
        v1 = (F1(1) + F1(2)*(zr(n)+i*zi(m)));
        v2 = (F2(1) + F2(2)*(zr(n)+i*zi(m)));
        v3 = 1/(F3(1) + F3(2)*(zr(n)+i*zi(m)));
        v4 = 1/(F4(1) + F4(2)*(zr(n)+i*zi(m)));
        S(m,n) = log(abs( v1*v2*v3*v4 )^2);
    end
end

% plot amplitude of polynomial
figure(1);
clf;
subplot(1,2,1);
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([zmin, zmax, zmin, zmax]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(S))-min(min(S));
if( range==0 )
    range=1;
end
imagesc( [zmin, zmax], [zmin, zmax], (S-min(min(S)))/range);
hold on;
xlabel('real z');
ylabel('imag z');
title('|S|^2');

% plot Fourier transform points on the unit circle
plot(0,0,'wo','LineWidth',2);
z=exp(-i*w*Dt);
plot(real(z),imag(z),'k+','LineWidth',2);


% plot zeros and zeros of filter
plot(real(zp1),imag(zp1),'w*','LineWidth',2);
plot(real(zp2),imag(zp2),'w*','LineWidth',2);
plot(real(zp3),imag(zp3),'w+','LineWidth',2);
plot(real(zp4),imag(zp4),'w+','LineWidth',2);


% plot frequency response of filter
subplot(1,2,2);
set(gca,'LineWidth',2);
hold on;
Sz = abs( ((F1(1)+F1(2)*z) .* (F2(1)+F2(2)*z)) ./ ((F3(1)+F3(2)*z) .* (F4(1)+F4(2)*z)) ) .^2;
axis( [0, f(Nf), 0, max(Sz)] );
plot(f(1:Nf),Sz(1:Nf),'k-','LineWidth',2);
plot(f(1:Nf),Sz(1:Nf),'ko','Linewidth',2);
xlabel('frequency, Hz');
ylabel('|S|^2');