fc=0.125;
wc=2*pi*fc;
N=31;
n=-(N-1)/2:(N-1)/2;
h_LP=(wc/pi)*sin(wc.*n)./(wc*n);
h_LP((N+1)/2)=(wc/pi);
h_BP=h_LP.*cos(2*pi*.25*n); % moving the basic filter to a new position. see Fourier properties for this one
zeros_BP=roots(h_BP);
figure(1);
plot(real(zeros_BP),imag(zeros_BP),'o')
hold on
plot(0,0,'bx')
axis([-1.5 1.5 -1.5 1.5]);
axis square

x=linspace(-1,1,1000);
plot(x,sqrt(1-x.^2),'r-')
plot(x,-sqrt(1-x.^2),'r-')
hold off

X=exp(-1i*n'*linspace(-pi,pi,1000));
h_BP=h_BP./sqrt(h_BP*h_BP');
H=h_BP*X;
figure(2);
plot(linspace(-pi,pi,1000)/2/pi,abs(H))

%% Removing the zeros not at the unit circle (for comparison purposes - what would happen if we only have the zeros at the places we would like to attenuate)
pos=find(abs(abs(zeros_BP)-1)<.05);
zeros_BP_mod=zeros_BP(pos);
h_BP_mod=poly(zeros_BP_mod);
h_BP_mod=h_BP_mod./sqrt(h_BP_mod*h_BP_mod');
ste=[-(length(h_BP_mod)-1)/2:(length(h_BP_mod)-1)/2];
X=exp(-1i*ste'*linspace(-pi,pi,1000));
H=h_BP_mod*X;
figure(2);
hold on
plot(linspace(-pi,pi,1000)/2/pi,abs(H),'r')
hold off