clear all
close all
clc

Fs = 6100;
Fn = Fs / 2;
Rp = 3;
Rs = 25;
Wp = [500 2000]/Fn;
Ws = [707 1414]/Fn;
K = 1024;

f1=Fs/pi*tan((pi*500)/Fs)
f2=Fs/pi*tan((pi*2000)/Fs)
f3=Fs/pi*tan((pi*707)/Fs)
f4=Fs/pi*tan((pi*1414)/Fs)

figure(1)
hold on;
plot([0 500],[25 25]);
plot([500 500],[0 25]);
plot([707 1414],[3 3]);
plot([2000 2000],[0 25]);
plot([2000 Fn],[25 25]);
axis([0 Fn 0 45]);
xlabel('f,Hz');
ylabel('a,dB');
title('Габарит на ЛФ');

[N,Wn] = buttord(Wp,Ws,Rp,Rs);
[Nz,Dz] = butter(N,Wn);


figure(2)
hold on;
plot([0 500],[25 25]);
plot([500 500],[0 25]);
plot([707 1414],[3 3]);
plot([2000 2000],[0 25]);
plot([2000 Fn],[25 25]);
axis([0 Fn 0 45]);

[H,w] = freqz(Nz,Dz);
ma = -20*log10(abs(H));
plot(w*Fs/(2*pi),ma,'r-.');
xlabel('f,Hz');
ylabel('a,dB');
title('Габарит на ЛВ и крива на затихване');
legend('Габарит','Затихване');


figure(3)
[z,p,k] = butter(N,Wn);

subplot(211);
zplane(Nz,Dz);
title('ПНД');

fi = angle(H);
fi1 = unwrap(fi);
subplot(212)
plot(w*Fs/(2*pi), fi1*180/pi);
title('ФЧХ');
xlabel('Честота (rad/s)');
ylabel('ФЧХ (rad)');

figure(4)
[h,t] = impz(Nz,Dz);
stem(t,h,'m');
grid;
title('Импулсна характеристика');
xlabel('Отчети на времето');
ylabel('Амплитуда (volts)');


figure(5)

n = (0:199)/Fs;
x = sin(2*pi*1000*n);
y = filter(Nz,Dz,x);
subplot(221);
plot(n,x);
title('Входен сигнал');
xlabel('Амплитуда (volts)');
ylabel('Отчети на времето');

subplot(223);
plot(n,y);
title('Изходен сигнал');
xlabel('Амплитуда (volts)');
ylabel('Отчети на времето');

Px = fft(x,K);
px = abs(Px(1:(K/2)));
Py = fft(y,K);
py = abs(Py(1:(K/2)));

subplot(222);
f = (0:(length(Px) - 1) / 2);
plot(f,px);
xlabel('Честота (Hz)');
title('Спектър на входен сигнал');

subplot(224)
f1 = (0:(length(Py)-1)/2);
plot(f1,py);
xlabel('Честота (Hz)');
title('Спектър на изходен сигнал');

b = [0.0177  0   -0.0887   0   0.1774   0   -0.1774   0   0.0887   0    -0.0177];
a = [1.0000   -3.0507    5.0514   -5.9813    5.7213         -4.3294    2.5774   -1.1966    0.4323   -0.1060    0.0151];
[b0,B,A] = dir2cas(b,a)

figure(6)
hold on;

%WsNew=373.92;
%WpNew=3450.15;
Wp1 = [511 3233];
Ws1 = [740 1731];


[N2,Wn2]=buttord(Wp1,Ws1,Rp,Rs,'s')
[Ns,Ds]=butter(N2,Wn2,'s')

printsys(Ns,Ds)%принтира Ns и Ds



[Ts,ws]=freqs(Ns,Ds);%изчислява аналоговата предавателна функция
ma1=-20*log10(abs(Ts));%затихването в децибели

plot(ws,ma1,'g--'); %крива на затихване на аналогов ЛФ
plot(w*Fs/(2*pi),ma,'r-.');%крива на затихване на цифров ЛФ

ylabel('a,dB');
xlabel('f,Hz');
title('Криви на затихването на аналогов ЛФ и цифров ЛФ');
legend('Крива на затихване на аналогов ЛФ','Крива на затихване на цифров ЛФ');

plot([0 500],[25 25]);
plot([500 500],[0 25]);
plot([707 1414],[3 3]);
plot([2000 2000],[0 25]);
plot([2000 Fn],[25 25]);
axis([0 Fn 0 45]);
xlabel('f,Hz');
ylabel('a,dB');
title('Габарит на ЛФ');