clear all
clc
n=0:127; N=length(n); % When N become bigger,better precision 
n1=0:255; N1=length(n1);
M=(N1-N)/2;
w=2*pi*[0:1/(N1-1):1];
wc1=0.2*pi; wc2=0.4*pi; % cutoff frequency
alfa=(N-1)/2;
w_1=boxcar(2*M-1); % rectangular window function
w1=[w_1',zeros(1,N1-2*M+1)];
hd=sin(wc2*(n-alfa))./(pi*(n-alfa))-sin(wc1*(n-alfa))./(pi*(n-alfa)); % bandpass filter
h_d=[hd,zeros(1,N1-N)];
hn=h_d.*w1;  
[h,W_1]=freqz(hn,1,N1/2);
h=[abs(h'),zeros(1,N1/2)]; 
subplot(221)
plot(n1,hn);grid on;
xlabel('n');ylabel('h(n)');title('Designed bandpass filter in time domain');
subplot(222)
plot(w,h);grid on;
xlabel('w');ylabel('h(w)');title('Designed bandpass filter in frequency domain');
axis([0 pi 0 1.2]);

% Input signal x(n)
noise_mu=0;  % mean=0
noise_var=0.01;  %variance 
noise=randn(size(n)).*sqrt(noise_var) + noise_mu; 
xn=noise+cos(0.1*pi*n)+0.75*cos(0.3*pi*n)+cos(0.5*pi*n);
w__a=[xn,zeros(1,N1-N)];
X=fft(w__a);
subplot(223)
plot(w,abs(X));
xlabel('w');ylabel('x(w)');title('Input signal in frequency domain');
axis([0 pi 0 70]);

% After fir filtering
y=abs(X).*h;
subplot(224)
plot(w,y);
xlabel('w');ylabel('y(w)');title('Output signal in frequency domain');
axis([0 pi 0 70]);