clear all;
close all;
clc;
Fs=4660;

Rp = 0.5;
Rs = 30;
Ws=[760 920];
Wp=[460 1040];
[N,Wn]=ellipord(Wp,Ws,Rp,Rs,'s')
[Ns,Ds]=ellip(N,Rp,Rs,Wn,'stop','s')
[Nz,Dz]=bilinear(Ns,Ds,Fs)

%Filtrirane na sin signali
n=(0:99)/Fs;
x1=sin(2*pi*800*n);   %propuska se
x2=sin(2*pi*1300*n);    %ne se propuska

figure(1);
plot(n,x1);
title('vremeva oblast - vhod signal X1');
xlabel('otcheti na vremeto');
ylabel('amplituda');

figure(2);
y1=filter(Nz,Dz,x1);
plot(n,y1);
title('vremeva oblast-izhod signal X1');
xlabel('otcheti na vremeto');
ylabel('amplituda');

figure(3);
plot(n,x2);
title('vremeva oblast - vhod signal X2');
xlabel('otcheti na vremeto');
ylabel('amplituda');

figure(4);
y2=filter(Nz,Dz,x2);
plot(n,y2);
title('vremeva oblast-izhod signal X2');
xlabel('otcheti na vremeto');
ylabel('amplituda');


%Chestotna oblast
K=1024;
n=0:(1/Fs):0.025;
[H,w]=freqz(Nz,Dz,K/2);

y2=filter(Nz,Dz,x2);
y1=filter(Nz,Dz,x1);

Px1=fft(x1,K);
px1=abs(Px1(1:(K/2)));
Py1=fft(y1,K);
py1=abs(Py1(1:(K/2)));


figure(5);
plot(w*(Fs/2)/pi,px1);
title('chestotna oblast - vhod signal X1');
xlabel('chesota, Hz');
ylabel('amplituda');

figure(6);
plot(w*(Fs/2)/pi,py1);
title('chestotna oblast - izhod signal X1');
xlabel('chesota, Hz');
ylabel('amplituda');



y2=filter(Nz,Dz,x2);
Px2=fft(x2,K);
px2=abs(Px2(1:(K/2)));
Py2=fft(y2,K);
py2=abs(Py2(1:(K/2)));

figure(7);
plot(w*(Fs/2)/pi,px2);
title('chestotna oblast - vhod signal X2');
xlabel('chesota, Hz');
ylabel('amplituda');


figure(8);
plot(w*(Fs/2)/pi,py2);
title('chestotna oblast - izhod signal X2');
xlabel('chesota, Hz');
ylabel('amplituda');