
%% Calculate normalized frequencies corresponding to 2N-1 zeros uniformly on the unit cirlce in the left half plane
N=201;
zero_freqs_nw=linspace(.25,.5,N); %places them in the "northwest" quarter
zero_freqs_sw=linspace(-.25,-.5,N);zero_freqs_sw=zero_freqs_sw(1:end-1);%places them in the "southwest" quarter, and removes the duplicated z=-1
zero_freq=[zero_freqs_sw zero_freqs_nw]; %this contains frequencies corresponding to all zeros


%% This iscompjutes the zeros in the z-domain 
%If you change "1", then you move the zeros away from the unit circle, but still in the left half plane
zros=1*exp(i*2*pi*zero_freq);  

%% We dont have any poles
poles=[];

%% Activate this and see what happens. We will go through this later in the course.
%poles=.95*zros;

%% Calculate the DTFT (square magnitude)
DTFT=[];
for f=linspace(-.5,.5,1001)
    z=exp(i*2*pi*f);
    temp=1;
    for nbr=1:numel(zros)
        temp=temp*(z-zros(nbr));
    end
    for nbr=1:numel(poles)
        temp=temp/(z-poles(nbr));
    end
    DTFT=[DTFT temp];
end
DTFT=abs(DTFT).^2;

%% normalize the DTFT to unit energy
df=1/1000;  % since  f=linspace(-.5,.5,1001)
Energy=sum(DTFT)*df;  %Integral 
DTFT=DTFT./Energy;
%sum(DTFT)*df If activated, you get 1

%% Plot DTFT
figure(1)
plot(linspace(-.5,.5,1001),10*log10(DTFT));
xlabel('f')
ylabel('DTFT (sqaure magnitude)')

%% Plot pole-zero diagram (with the default setting, you need to zoom in to see what happens)
figure(2)
x=linspace(-1,1,1001);
hold off
plot(x,sqrt(1-x.^2),'k');
hold on
plot(x,-sqrt(1-x.^2),'k');
plot([-1.3 1.3],[0 0 ],'r');
plot([0 0 ],[-1.3 1.3],'r');
for n=1:numel(zros)
    plot(real(zros(n)),imag(zros(n)),'bo');
end
for n=1:numel(poles)
    plot(real(poles(n)),imag(poles(n)),'bx');
end
axis square





