%% Truncated sinc-filter for low pass filtering

%% Finding the LP filter and its zeros
fc=0.25; % Cutoff according to plot on lecture notes
wc=2*pi*fc;
N=31; % length of filter 
n=-(N-1)/2:(N-1)/2; %time axis 
h_LP=(wc/pi)*sin(pi/2.*n)./(pi/2*n); % not shifted in time, but this is not relevant for magnitude of H(w)
h_LP((N+1)/2)=(wc/pi); % manually adjust the center tap since it is of the form sin(0)/0
zeros_LP=roots(h_LP); % find the zeros of the FIR filter

%% Plotting the zeros
figure(1);
plot(real(zeros_LP),imag(zeros_LP),'o') % Plot them
hold on
plot(0,0,'bx')  % plot the poles (all of them are at z=0
axis([-1.5 1.5 -1.5 1.5]);
axis square

%% Plotting the unit circle
x=linspace(-1,1,1000); 
plot(x,sqrt(1-x.^2),'r-')
plot(x,-sqrt(1-x.^2),'r-')
hold off

%% Finding and plotting the Fourier transform
X=exp(-1i*n'*linspace(-pi,pi,1000)); % Fourier transformation matrix
h_LP=h_LP./sqrt(h_LP*h_LP'); % Normalize to unit energy for later usage 
H=h_LP*X; % Fourier transform
figure(2);
plot(linspace(-pi,pi,1000),abs(H)) % plotting the magnitude of H(w)

%% Removing the zeros in the right half plane (for comparison purposes - what would happen if we only have the zeros at the places we would like to attenuate)
pos=find(real(zeros_LP)<0);  % finding them using hte find-command
zeros_LP_mod=zeros_LP(pos);
h_LP_mod=poly(zeros_LP_mod); % make a filter from the left half plane zeros
h_LP_mod=h_LP_mod./sqrt(h_LP_mod*h_LP_mod'); % normalizing it to the same energy as the original filter.
ste=[-(length(h_LP_mod)-1)/2:(length(h_LP_mod)-1)/2];
X=exp(-1i*ste'*linspace(-pi,pi,1000));% Fourier transformation matrix
H=h_LP_mod*X; % Fourier transform
figure(2);
hold on
plot(linspace(-pi,pi,1000),abs(H),'r') % plotting the magnitude of the modified H(w) in red
hold off