scilab

//BcIb4
//порядок 4
//фильтр баттерворта - каскадная схема из прямых форм
//фильтр Чебышева II рода, прямая форма II (каноническая)
//_______________________________________________________
//Модельный сигнал
N=1024;
dT=0.01;
x=0:dT:N*dT;
sig=0.7*sin(2*%pi*2*x)+0.35*cos(2*%pi*4*x);
y=sig+grand(1, (N+1), "nor", 0, 0.2);
clf(0);
scf(0);
subplot(211);
plot2d(x,y);
xtitle("График модельного сигнала","Время t, с", "Амплитуда, ед.");
//Спектр модельного сигнала
spectr = abs(fft(y));
dW = 1/(dT*N); //Шаг по частота
W=0:dW:1/dT; //Диапазон частот
_W_ = length(W);
_spectr_=length(spectr);
subplot(212);
plot2d(W(1:_W_/2),spectr(1:_spectr_/2)); //Делим для вывода спектра до частоты Найквиста
xtitle('Спектр сигнала','Частота f, Гц','Амплитуда, ед.');
//Аналоговые фильтры
hsButt=analpf(4, 'butt', [0 0], 0.3*2/dT); //строим вид передаточной функции для фильров
hsCheb2=analpf(4,'cheb2',[0.1 0.1],0.3*2/dT);
ls_Butt=syslin('c', hsButt); //задаём линейную систему
ls_Cheb2=syslin('c', hsCheb2);
scf(1);
clf(1);
bode(ls_Butt); //Диаграмма Боде для полученной линейной системы
xtitle("АЧХ и ФЧХ фильтра Баттерворта");
scf(2);
clf(2);
bode(ls_Cheb2);
xtitle("АЧХ и ФЧХ фильтра Чебышёва II рода");
//Обработка спектра сигнала фильтрами
FButt=repfreq(ls_Butt,W);
FCheb2=repfreq(ls_Cheb2,W);
SpectrButt=fft(y).*FButt;//Обработка спектра
SpectrCheb2=fft(y).*FCheb2;
scf(3);
clf(3);
subplot(311);
plot2d(W(1:_W_/2),spectr(1:_spectr_/2));
xtitle("Спектр сигнала до фильтрации","f", "A")
subplot(312)
plot2d(W(1:_W_/2),abs(SpectrButt(1:_spectr_/2)))
xtitle("Спектр сигнала после фильтра Баттерворта","f", "A")
subplot(313)
plot2d(W(1:_W_/2),abs(SpectrCheb2(1:_W_/2)))
xtitle("Спектр сигнала после фильтра Чебышёва II рода","f", "A")
//Обратноe прeобразованиe обработанного спeктра
reverseButt=real(ifft(SpectrButt)); //Восстановление из спектра
reverseCheb2=real(ifft(SpectrCheb2));
reverseButt=reverseButt/max(reverseButt); //Нормировка по максимальному значению
reverseCheb2=reverseCheb2/max(reverseCheb2);
clf(4);
scf(4);
subplot(311);
plot2d(x, (y/max(abs(y))));
xtitle("Сигнал до фильтрации","t", "А");
subplot(312);
plot2d(x, reverseButt)
xtitle("Сигнал после фильтра Баттерворта","t", "А")
subplot(313)
plot2d(x, reverseCheb2)
xtitle("Сигнал после фильтра Чебышёва II рода","t", "А")
//Переходим к цифровым фильтрам
D_hsButt=iir(4,'lp','butt',[0.3 0], [0 0]);
//                                      2            3            4
//   0.1671793 + 0.6687171z + 1.0030756z + 0.6687171z + 0.1671793z
//   --------------------------------------------------------------
//                                          2            3   4
//       0.0301189 + 0.1826757z + 0.6799785z + 0.7820952z + z
D_hsCheb2=iir(4,'lp','cheb2',[0.3 0], [0.1 0.1]);
//                                      2            3            4
//   0.2727433 + 0.6738200z + 0.8994179z + 0.6738200z + 0.2727433z
//   --------------------------------------------------------------
//                                          2            3   4
//       0.0780474 + 0.2898983z + 0.8045473z + 0.6200515z + z
[F_Butt,D_butt]=repfreq(D_hsButt,0,0.5);
clf(5);scf(5);
subplot(211);
plot2d(F_Butt*100,abs(D_butt));
xtitle("АЧХ дискретного фильтра Баттерворта","Частота (Гц)", "Амплитуда");
//https://help.scilab.org/doc/5.5.2/en_US/dbphi.html
[db,phi]=dbphi(D_butt);
subplot(212);
plot2d(F_Butt*100, phi);
xtitle("ФЧХ дискретного фильтра Баттерворта","Частота (Гц)", "Фаза (Град)");
[F_Cheb2,D_Cheb2]=repfreq(D_hsCheb2,0,0.5);
clf(6);scf(6);
subplot(211);
plot2d(F_Cheb2*100,abs(D_Cheb2));
xtitle("АЧХ дискретного фильтра Чебышёва II рода","Частота (Гц)", "Амплитуда");
[db,phi]=dbphi(D_Cheb2);
subplot(212);
plot2d(F_Cheb2*100,phi)
xtitle("ФЧХ дискретного фильтра Чебышёва II рода","Частота (Гц)", "Фаза (Град)");


//каскадная схема из прямых форм для баттерфорда
//вторая строка коэффициентов с обратным знаком
//в соответсвии со схемой - учитываем знаки и берем обратные
c=factors(D_hsButt(3))
c1=coeff(c(1));
c2=coeff(c(2));
dButtsig=y
dButtsig2=y
G1=1
G2=1
for i=3:N
    dButtsig(i)=G1*y(i)+2*G1*y(i-1)+G1*y(i-2)-c1(2)*dButtsig(i-1)-c1(1)*dButtsig(i-2)
    dButtsig2(i)=G2*dButtsig(i)+2*G2*dButtsig(i-1)+G2*dButtsig(i-2)-c2(2)*dButtsig2(i-1)-c2(1)*dButtsig2(i-2)
end

//прямая схема (прямая форма II) каноническая для фильтра Чебышева II рода
//определение числителя и знаменателя
a = coeff(D_hsCheb2(3));
b = coeff(D_hsCheb2(2));
//
yc_d = zeros(1,length(y)+5);
v = zeros(1,length(y)+5);
//реализуем схему
for i=6:N
    v(i) = y(i) - (a(4)*v(i-1) + a(3)*v(i-2) + a(2)*v(i-3) + a(1)*v(i-4));
    yc_d(i) = b(5)*v(i-1) + b(4)*v(i-2) + b(3)*v(i-3) +..
              b(2)*v(i-4) + b(1)*v(i-5);
end

dButtsig2=dButtsig2/max(dButtsig2);
Yc_d=yc_d/max(yc_d);
//Спектр сигнала

//Yc_d = fftshift(abs(fft(yc_d)));
SdButtsig2=fft(dButtsig2)
SYc_d = fft(Yc_d);
clf(7)
scf(7)
subplot(311)
plot2d(W(1:_W_/2),spectr(1:_spectr_/2))
xtitle("Спектр сигнала до фильтрации","f", "А")
subplot(312)
plot2d(W(1:_W_/2),abs(SdButtsig2(1:_W_/2)))
xtitle("Спектр сигнала после фильтра Баттерворта каскадной формы","f", "А");
subplot(313)
plot2d(W(1:_W_/2),abs(SYc_d(1:_W_/2)))
xtitle("Спектр сигнала после фильтра Чебышева II рода","f", "А");

//Сигнал после фильтрации
clf(8)
scf(8)
subplot(311)
plot2d(x,(y/max(abs(y))))
xtitle("Исходный сигнал","t", "А")
subplot(312)
plot2d(x(1:length(dButtsig2)),dButtsig2)
xtitle("Сигнал после обработки цифровым фильтром Баттерворта каскадной формы","t", "А")
subplot(313)
//plot2d(x(1:length(Yc_d)),Yc_d)
//plot2d(x,(6:length(SYc_d)),Yc_d)
plot2d(x,(6:length(Yc_d)),Yc_d)
xtitle("Сигнал после обработки цифровым фильтром Чебышева","t", "А")